Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
field_seq_visu_vtk_paraview.cc
Go to the documentation of this file.
1
21//
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
36using namespace std;
37namespace rheolef {
38
39// ----------------------------------------------------------------------------
40// field puts
41// ----------------------------------------------------------------------------
42// extern:
43template <class T> odiststream& field_put_vtk (odiststream&, const field_basic<T,sequential>&);
44
45template <class T>
46odiststream&
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();
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// --------------------------------------------------------------------
175template <class T> idiststream& geo_get_vtk (idiststream&, geo_basic<T,sequential>&);
176
177template <class T>
178field_basic<T,sequential>
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// --------------------------------------------------------------------
327template <class T> idiststream& geo_get_vtk (idiststream&, geo_basic<T,sequential>&);
328
329template <class T>
330geo_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) \
416template odiststream& \
417visu_vtk_paraview<Float> ( \
418 odiststream&, \
419 const field_basic<Float,sequential>&); \
420template field_basic<Float,sequential> \
421paraview_plane_cut ( \
422 const field_basic<Float,sequential>&, \
423 const point_basic<Float>&, \
424 const point_basic<Float>&); \
425template geo_basic<Float,sequential> \
426paraview_extract_isosurface ( \
427 const field_basic<Float,sequential>&);
428
430#undef _RHEOLEF_instanciation
431
432}// namespace rheolef
#define _RHEOLEF_instanciation(T, M, A)
Definition asr.cc:223
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
vec< T, M > & set_u()
Definition field.h:284
T min() const
Definition field.h:786
T max() const
Definition field.h:802
const geo_type & get_geo() const
Definition field.h:271
const space_type & get_space() const
Definition field.h:270
typename float_traits< T >::type float_type
Definition field.h:228
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
the finite element space
Definition space.h:382
#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)")
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
This file is part of Rheolef.
geo_basic< T, sequential > paraview_extract_isosurface(const field_basic< T, sequential > &uh)
odiststream & visu_vtk_paraview(odiststream &, const field_basic< T, sequential > &)
point_basic< T > vect(const point_basic< T > &v, const point_basic< T > &w)
Definition point.h:264
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
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)
field_basic< T, sequential > paraview_plane_cut(const field_basic< T, sequential > &uh, const point_basic< T > &origin, const point_basic< T > &normal)
details::field_expr_v2_nonlinear_terminal_function< details::normal_pseudo_function< Float > > normal()
normal: see the expression page for the full documentation
odiststream & field_put_vtk(odiststream &, const field_basic< T, sequential > &, std::string, bool)
rheolef::std enable_if ::type dot const Expr1 expr1, const Expr2 expr2 dot(const Expr1 &expr1, const Expr2 &expr2)
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
T norm(const vec< T, M > &x)
norm(x): see the expression page for the full documentation
Definition vec.h:387
STL namespace.
point_basic< size_t > resolution
static std::string python(const point &x, size_t d=3)