Rheolef  7.2
an efficient C++ finite element environment
field_seq_visu_vtk_paraview.cc
Go to the documentation of this file.
1 //
22 // paraview vtk visualization
23 //
24 // author: Pierre.Saramito@imag.fr
25 //
26 // date: 12 may 1997 update: 23 oct 2011
27 //
28 #include "rheolef/field_expr.h"
29 #include "rheolef/piola_util.h"
30 #include "rheolef/rheostream.h"
31 #include "rheolef/iorheo.h"
32 #include "rheolef/iofem.h"
33 #include "rheolef/interpolate.h"
34 #include "rheolef/render_option.h"
35 
36 using namespace std;
37 namespace rheolef {
38 
39 // ----------------------------------------------------------------------------
40 // field puts
41 // ----------------------------------------------------------------------------
42 // extern:
43 template <class T> odiststream& field_put_vtk (odiststream&, const field_basic<T,sequential>&);
44 
45 template <class T>
46 odiststream&
48 {
49  //
50  // 0) prerequises
51  //
52  using namespace std;
53  typedef typename field_basic<T,sequential>::float_type float_type;
55  typedef point_basic<size_type> ilat;
56  ostream& os = ops.os();
57  //
58  // 1) set render options
59  //
60  render_option popt;
61  popt.valued = uh.get_space().valued();
62  popt.fill = iorheo::getfill(os); // isocontours or color fill
63  popt.elevation = iorheo::getelevation(os);
64  popt.color = iorheo::getcolor(os);
65  popt.gray = iorheo::getgray(os);
66  popt.black_and_white = iorheo::getblack_and_white(os);
67  popt.showlabel = iorheo::getshowlabel(os);
68  popt.stereo = iorheo::getstereo(os);
69  popt.volume = iorheo::getvolume(os);
70  popt.iso = true;
71  popt.cut = iorheo::getcut(os);
72  popt.grid = iorheo::getgrid(os);
73  popt.format = iorheo::getimage_format(os);
74  popt.mark = iorheo::getmark(os);
75  bool is_scalar = (popt.valued == "scalar");
76  if (popt.mark == "") popt.mark = popt.valued;
77  popt.style = (popt.valued == "vector") ? (iorheo::getvelocity(os) ? "velocity" : "deformation") : "none";
78  popt.scale = iorheo::getvectorscale(os);
79  popt.origin = iofem::getorigin(os);
80  popt.normal = iofem::getnormal(os);
81  popt.resolution = iofem::getresolution(os);
82  popt.n_isovalue = iorheo::getn_isovalue(os);
83  popt.n_isovalue_negative = iorheo::getn_isovalue_negative(os);
84  popt.isovalue = iorheo::getisovalue(os); // isovalue is always a Float
85  popt.label = iorheo::getlabel(os);
86  const geo_basic<float_type,sequential>& omega = uh.get_geo();
87  size_type dim = omega.dimension();
88  size_type map_dim = omega.map_dimension();
89  popt.view_2d = (dim == 2);
90  popt.view_map = (dim > map_dim);
91 #if (_RHEOLEF_PARAVIEW_VERSION_MAJOR >= 5) && (_RHEOLEF_PARAVIEW_VERSION_MINOR >= 5)
92  // paraview version >= 5.5 has high order elements
93  popt.high_order = (omega.order() > 1 || uh.get_space().degree() > 1);
94 #else
95  popt.high_order = false;
96 #endif
97 #if (_RHEOLEF_PARAVIEW_VERSION_MAJOR == 5) && (_RHEOLEF_PARAVIEW_VERSION_MINOR == 7)
98  // paraview version == 5.7 has opacity bug
99  popt.have_opacity_bug = true;
100 #else
101  popt.have_opacity_bug = false;
102 #endif
103  popt.xmin = uh.get_geo().xmin();
104  popt.xmax = uh.get_geo().xmax();
105  field_basic<T,sequential> uh_scalar;
106  if (popt.valued == "scalar") {
107  uh_scalar = uh;
108  } else {
109  size_type k = uh.get_space().degree();
110  space_basic<T,sequential> T0h (uh.get_geo(), "P"+std::to_string(k));
111  uh_scalar = interpolate(T0h, norm(uh));
112  }
113  // TODO: u_min=min(uh.u.min,uh.b.min) instead of interpolate norm ?
114  popt.f_min = uh_scalar.min();
115  popt.f_max = uh_scalar.max();
116  //
117  // 2) output data
118  //
119  bool verbose = iorheo::getverbose(os);
120  bool clean = iorheo::getclean(os);
121  bool execute = iorheo::getexecute(os);
122  string basename = iorheo::getbasename(os);
123  string outfile_fmt = "";
124  string tmp = get_tmpdir() + "/";
125  if (!clean) tmp = "";
126  string filelist;
127  string filename = tmp+basename + ".vtk";
128  filelist = filelist + " " + filename;
129  ofstream vtk_os (filename.c_str());
130  odiststream vtk (vtk_os);
131  if (verbose) clog << "! file \"" << filename << "\" created.\n";
132  field_put_vtk (vtk, uh);
133  vtk.close();
134  //
135  // 3) create python data file
136  //
137  std::string py_name = filename = tmp+basename + ".py";
138  filelist = filelist + " " + filename;
139  ofstream py (filename.c_str());
140  if (verbose) clog << "! file \"" << filename << "\" created.\n";
141  point xmin = uh.get_geo().xmin(),
142  xmax = uh.get_geo().xmax();
143  py << popt
144  << endl
145  << "paraview_field_" << popt.valued << "(paraview, \"" << tmp+basename << "\", opt)" << endl
146  << endl
147  ;
148  py.close();
149  //
150  // 4) run pyton
151  //
152  int status = 0;
153  string command;
154  if (execute) {
155  string prog = (popt.format == "") ? "paraview --script=" : "pvbatch --use-offscreen-rendering ";
156  command = "LANG=C PYTHONPATH=" + string(_RHEOLEF_PKGDATADIR) + " " + prog + py_name;
157  if (popt.format != "") command = "DISPLAY=:0.0 " + command;
158  if (popt.stereo && popt.format == "") command = command + " --stereo";
159  if (verbose) clog << "! " << command << endl;
160  status = system (command.c_str());
161  }
162  //
163  // 4) clear vtk data
164  //
165  if (clean) {
166  command = "/bin/rm -f " + filelist;
167  if (verbose) clog << "! " << command << endl;
168  status = system (command.c_str());
169  }
170  return ops;
171 }
172 // --------------------------------------------------------------------
173 // plane cut : via vtk-paraview script
174 // --------------------------------------------------------------------
175 template <class T> idiststream& geo_get_vtk (idiststream&, geo_basic<T,sequential>&);
176 
177 template <class T>
178 field_basic<T,sequential>
180  const field_basic<T,sequential>& uh,
181  const point_basic<T>& origin,
182  const point_basic<T>& normal)
183 {
184  //
185  // 1) prerequises
186  //
187  using namespace std;
188  typedef typename field_basic<T,sequential>::float_type float_type;
190  typedef typename geo_basic<float_type,sequential>::node_type node_type;
191  typedef point_basic<size_type> ilat;
192  ostream& os = std::cout;
193  bool verbose = iorheo::getverbose(os);
194  bool clean = iorheo::getclean(os);
195  bool execute = iorheo::getexecute(os);
196  string basename = iorheo::getbasename(os);
197  string tmp = get_tmpdir() + "/";
198  if (!clean) tmp = "";
199  //
200  // 2) output data
201  //
202  string filelist;
203  string vtk_name = tmp+basename + ".vtk";
204  filelist = filelist + " " + vtk_name;
205  ofstream vtk_os (vtk_name.c_str());
206  odiststream vtk (vtk_os);
207  if (verbose) clog << "! file \"" << vtk_name << "\" created.\n";
208  field_put_vtk (vtk, uh);
209  vtk.close();
210  //
211  // create python script
212  //
213  //
214  // 3) create python data file
215  //
216  string py_name = tmp+basename + "-cut.py";
217  string vtk_cut_name = tmp+basename + "-cut.vtk";
218  filelist = filelist + " " + py_name + " " + vtk_cut_name;
219  ofstream py (py_name.c_str());
220  if (verbose) clog << "! file \"" << py_name << "\" created.\n";
221  py << setprecision(numeric_limits<T>::digits10)
222  << "from paraview.simple import *" << endl
223  << "reader = LegacyVTKReader(FileNames=['" << vtk_name << "'])" << endl
224  << "slice = Slice(SliceType=\"Plane\")" << endl
225  << "slice.SliceOffsetValues = [0.0]" << endl
226  << "slice.SliceType.Origin = " << render_option::python(origin) << endl
227  << "slice.SliceType.Normal = " << render_option::python(normal) << endl
228  << "writer = paraview.simple.CreateWriter(\"" << vtk_cut_name << "\", slice)" << endl
229  << "writer.FileType = 'Ascii'" << endl
230  << "writer.UpdatePipeline()" << endl
231  ;
232  py.close();
233  //
234  // 3) run pyton
235  //
236  int status = 0;
237  string command;
238  command = "LANG=C PYTHONPATH=" + string(_RHEOLEF_PKGDATADIR) + " pvpython < " + py_name;
239  if (verbose) clog << "! " << command << endl;
240  status = system (command.c_str());
241  //
242  // 4) load vtk cutted mesh
243  //
244  ifstream vtk_polydata (vtk_cut_name.c_str());
245  check_macro (vtk_polydata, "field: vtk polydata file \"" << vtk_cut_name << "\" not found.");
246  if (verbose) clog << "! load `" << vtk_cut_name << "'." << endl;
247  idiststream ips_vtk_polydata (vtk_polydata);
248  geo_basic<T,sequential> gamma_cut;
249  geo_get_vtk (ips_vtk_polydata, gamma_cut);
250  gamma_cut.set_name (basename + "-cut");
251  check_macro (gamma_cut.n_node() > 0, "empty mesh & plane intersection: HINT check normal and origin");
252  //
253  // clean
254  //
255  if (clean) {
256  command = "/bin/rm -f " + filelist;
257  if (verbose) clog << "! " << command << endl;
258  status = system (command.c_str());
259  }
260  //
261  // project geometry from 2d/3d on plane in 1d/2d
262  //
263  // 1D) x := dot(OM, tangent)
264  //
265  // 2D) x1 := dot(OM, t1)
266  // x2 := dot(OM, t2)
267  //
268  bool use_projection = true;
269  if (use_projection) {
270  switch (uh.get_geo().dimension()) {
271  case 2: {
272  // project in 1D space
273  gamma_cut.set_dimension (1);
274  point_basic<T> tangent = point_basic<T>(normal[1], -normal[0]);
275  disarray<node_type,sequential> node = gamma_cut.get_nodes();
276  for (size_type i = 0, n = node.size(); i < n; i++) {
277  node[i] = point_basic<T> (dot(node[i], tangent), 0);
278  }
279  gamma_cut.set_nodes (node);
280  break;
281  }
282  case 3: {
283  gamma_cut.set_dimension (2);
284  point_basic<T> t1 = point_basic<T>(1, 0, 0);
285  if (norm(vect(t1, normal)) == Float(0)) t1 = point_basic<T>(0, 1, 0);
286  t1 = t1 - dot(t1,normal)*normal;
287  t1 = t1 / norm(t1);
288  point_basic<T> t2 = vect(normal,t1);
289  disarray<node_type,sequential> node = gamma_cut.get_nodes();
290  for (size_type i = 0, n = node.size(); i < n; i++) {
291  node[i] = point_basic<T> (dot(node[i] - origin, t1), dot(node[i] - origin, t2), 0);
292  }
293  gamma_cut.set_nodes (node);
294  break;
295  }
296  default: error_macro ("unexpected dimension="<< uh.get_geo().dimension());
297  }
298  }
299  // -------------------------------
300  // get field values
301  // polydata header:
302  // POINT_DATA <n_node>
303  // SCALARS scalars float
304  // LOOKUP_TABLE default
305  // -------------------------------
306  string data_type;
307  bool founded = false;
308  while (vtk_polydata.good()) {
309  vtk_polydata >> data_type;
310  if ((data_type == "POINT_DATA") ||
311  (data_type == "CELL_DATA" )) {
312  break;
313  }
314  }
315  check_macro (data_type == "POINT_DATA" || data_type == "CELL_DATA",
316  "paraview_plane_cut: unexpected vtk polydata file");
317  scatch (vtk_polydata, "default");
318  string approx = (data_type == "POINT_DATA") ? "P1" : "P0";
319  space_basic<T,sequential> V_cut (gamma_cut, approx);
320  field_basic<T,sequential> u_cut (V_cut);
321  u_cut.set_u().get_values (ips_vtk_polydata);
322  return u_cut;
323 }
324 // --------------------------------------------------------------------
325 // isovalue surface : via vtk-paraview script
326 // --------------------------------------------------------------------
327 template <class T> idiststream& geo_get_vtk (idiststream&, geo_basic<T,sequential>&);
328 
329 template <class T>
330 geo_basic<T,sequential>
332 {
333  //
334  // 1) prerequises
335  //
336  using namespace std;
337  typedef typename field_basic<T,sequential>::float_type float_type;
339  typedef typename geo_basic<float_type,sequential>::node_type node_type;
340  typedef point_basic<size_type> ilat;
341  ostream& os = std::cout;
342  bool verbose = iorheo::getverbose(os);
343  bool clean = iorheo::getclean(os);
344  bool execute = iorheo::getexecute(os);
345  T isovalue = iorheo::getisovalue(os);
346  string basename = iorheo::getbasename(os);
347  string tmp = get_tmpdir() + "/";
348  if (!clean) tmp = "";
349  //
350  // 2) output data
351  //
352  string filelist;
353  string vtk_name = tmp+basename + ".vtk";
354  filelist = filelist + " " + vtk_name;
355  ofstream vtk_os (vtk_name.c_str());
356  odiststream vtk (vtk_os);
357  if (verbose) clog << "! file \"" << vtk_name << "\" created.\n";
358  field_put_vtk (vtk, uh);
359  vtk.close();
360  //
361  // create python script
362  //
363  //
364  // 3) create python data file
365  //
366  string py_name = tmp+basename + "-iso.py";
367  string vtk_iso_name = tmp+basename + "-iso.vtk";
368  filelist = filelist + " " + py_name + " " + vtk_iso_name;
369  ofstream py (py_name.c_str());
370  if (verbose) clog << "! file \"" << py_name << "\" created.\n";
371  py << setprecision(numeric_limits<T>::digits10)
372  << "from paraview.simple import *" << endl
373  << "reader = LegacyVTKReader(FileNames=['" << vtk_name << "'])" << endl
374  << "iso_surface = Contour(PointMergeMethod=\"Uniform Binning\")" << endl
375  << "iso_surface.ContourBy = ['POINTS', 'scalar']" << endl
376  << "iso_surface.Isosurfaces = [" << isovalue << "]" << endl
377  << "writer = paraview.simple.CreateWriter(\"" << vtk_iso_name << "\", iso_surface)" << endl
378  << "writer.FileType = 'Ascii'" << endl
379  << "writer.UpdatePipeline()" << endl
380  ;
381  py.close();
382  //
383  // 3) run pyton
384  //
385  int status = 0;
386  string command;
387  command = "LANG=C PYTHONPATH=" + string(_RHEOLEF_PKGDATADIR) + " pvbatch " + py_name;
388  if (verbose) clog << "! " << command << endl;
389  status = system (command.c_str());
390  //
391  // 4) load vtk isosurface mesh
392  //
393  ifstream vtk_polydata (vtk_iso_name.c_str());
394  check_macro (vtk_polydata, "field: vtk polydata file \"" << vtk_iso_name << "\" not found.");
395  if (verbose) clog << "! load `" << vtk_iso_name << "'." << endl;
396  idiststream ips_vtk_polydata (vtk_polydata);
397  geo_basic<T,sequential> gamma_iso;
398  geo_get_vtk (ips_vtk_polydata, gamma_iso);
399  gamma_iso.set_name (basename + "-iso");
400  gamma_iso.set_dimension (uh.get_geo().dimension());
401  check_macro (gamma_iso.n_node() > 0, "empty mesh & plane intersection: HINT check normal and origin");
402  //
403  // clean
404  //
405  if (clean) {
406  command = "/bin/rm -f " + filelist;
407  if (verbose) clog << "! " << command << endl;
408  status = system (command.c_str());
409  }
410  return gamma_iso;
411 }
412 // ----------------------------------------------------------------------------
413 // instanciation in library
414 // ----------------------------------------------------------------------------
415 #define _RHEOLEF_instanciation(T) \
416 template odiststream& \
417 visu_vtk_paraview<Float> ( \
418  odiststream&, \
419  const field_basic<Float,sequential>&); \
420 template field_basic<Float,sequential> \
421 paraview_plane_cut ( \
422  const field_basic<Float,sequential>&, \
423  const point_basic<Float>&, \
424  const point_basic<Float>&); \
425 template geo_basic<Float,sequential> \
426 paraview_extract_isosurface ( \
427  const field_basic<Float,sequential>&);
428 
430 #undef _RHEOLEF_instanciation
431 
432 }// namespace rheolef
field::size_type size_type
Definition: branch.cc:430
see the Float page for the full documentation
see the point page for the full documentation
see the disarray page for the full documentation
Definition: disarray.h:497
const space_type & get_space() const
Definition: field.h:270
const geo_type & get_geo() const
Definition: field.h:271
T min() const
Definition: field.h:786
vec< T, M > & set_u()
Definition: field.h:284
T max() const
Definition: field.h:802
typename float_traits< T >::type float_type
Definition: field.h:228
size_type n_node() const
Definition: geo.h:1172
void set_name(std::string name)
void set_nodes(const disarray< node_type, sequential > &x)
const disarray< node_type, sequential > & get_nodes() const
Definition: geo.h:1177
generic mesh with rerefence counting
Definition: geo.h:1089
idiststream: see the diststream page for the full documentation
Definition: diststream.h:336
odiststream: see the diststream page for the full documentation
Definition: diststream.h:137
std::ostream & os()
Definition: diststream.h:247
#define _RHEOLEF_PKGDATADIR
Definition: config.h:234
#define error_macro(message)
Definition: dis_macros.h:49
Expr1::float_type T
Definition: field_expr.h:230
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
#define _RHEOLEF_instanciation(T)
verbose clean transpose logscale grid shrink ball stereo iso volume skipvtk deformation fastfieldload lattice reader_on_stdin color format format format format format format format format format format format format format format format format vtk
string command
Definition: mkgeo_ball.sh:136
rheolef::details::is_vec dot
This file is part of Rheolef.
geo_basic< T, sequential > paraview_extract_isosurface(const field_basic< T, sequential > &uh)
T norm(const vec< T, M > &x)
norm(x): see the expression page for the full documentation
Definition: vec.h:387
bool scatch(std::istream &in, const std::string &ch, bool full_match=true)
scatch: see the rheostream page for the full documentation
Definition: scatch.icc:44
field_basic< T, sequential > paraview_plane_cut(const field_basic< T, sequential > &uh, const point_basic< T > &origin, const point_basic< T > &normal)
field_basic< T, M > interpolate(const space_basic< T, M > &V2h, const field_basic< T, M > &u1h)
see the interpolate page for the full documentation
Definition: interpolate.cc:233
point_basic< T > vect(const point_basic< T > &v, const point_basic< T > &w)
Definition: point.h:264
odiststream & field_put_vtk(odiststream &, const field_basic< T, sequential > &, std::string, bool)
odiststream & visu_vtk_paraview(odiststream &, const field_basic< T, sequential > &)
std::string get_tmpdir()
get_tmpdir: see the rheostream page for the full documentation
Definition: rheostream.cc:54
idiststream & geo_get_vtk(idiststream &ips, geo_basic< T, sequential > &omega)
details::field_expr_v2_nonlinear_terminal_function< details::normal_pseudo_function< Float > > normal()
normal: see the expression page for the full documentation
std::string python(const point_basic< T > &x, size_t d=3)
point_basic< size_t > resolution
Definition: render_option.h:48