Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
field_seq_visu_gnuplot.cc
Go to the documentation of this file.
1
21//
22// gnuplot geo visualisation
23//
24// author: Pierre.Saramito@imag.fr
25//
26// date: 16 sept 2011
27//
28# include "rheolef/field.h"
29# include "rheolef/field_wdof_indirect.h"
30# include "rheolef/field_evaluate.h"
31# include "rheolef/piola_util.h"
32# include "rheolef/rheostream.h"
34# include "rheolef/rounder.h"
35# include "rheolef/interpolate.h"
36
37namespace rheolef {
38
39template <class T>
40odiststream& visu_gnuplot (odiststream& ops, const geo_basic<T,sequential>& omega);
41
42// ----------------------------------------------------------------------------
43// puts for one element
44// ----------------------------------------------------------------------------
45template<class T>
46static
47void
48put_edge (
49 std::ostream& gdat,
50 const geo_basic<T,sequential>& omega,
51 const geo_element& K,
52 const field_basic<T,sequential>& uh,
53 const fem_on_pointset<T>& fops,
54 size_t my_order,
55 bound_type<T>& bbox)
56{
57 using namespace std;
58 typedef typename geo_basic<T,sequential>::size_type size_type;
59 typedef point_basic<size_type> ilat;
60 size_type dim = omega.dimension();
61 const Eigen::Matrix<piola<T>,Eigen::Dynamic,1>& piola = fops.get_piola_on_pointset().get_piola (omega, K);
62 Eigen::Matrix<T,Eigen::Dynamic,1> value;
63 field_evaluate (uh, fops, omega, K, value);
64 for (size_type loc_idof = 0, loc_ndof = piola.size(); loc_idof < loc_ndof; loc_idof++) {
65 bbox.update (piola[loc_idof].F, value[loc_idof]);
66 }
67 for (size_type i = 0; i <= my_order; i++) {
68 size_type loc_idof = reference_element_e::ilat2loc_inod (my_order, ilat(i));
69 piola[loc_idof].F.put (gdat, dim); gdat << " " << value[loc_idof] << endl;
70 }
71 gdat << endl;
72 gdat << endl;
73}
74template<class T>
75static
76void
77put_triangle (
78 std::ostream& gdat,
79 const geo_basic<T,sequential>& omega,
80 const geo_element& K,
81 const field_basic<T,sequential>& uh,
82 const fem_on_pointset<T>& fops,
83 size_t my_order,
84 bound_type<T>& bbox)
85{
86 using namespace std;
87 typedef typename geo_basic<T,sequential>::size_type size_type;
88 typedef point_basic<size_type> ilat;
89 size_type dim = omega.dimension();
90 const Eigen::Matrix<piola<T>,Eigen::Dynamic,1>& piola = fops.get_piola_on_pointset().get_piola (omega, K);
91 Eigen::Matrix<T,Eigen::Dynamic,1> value;
92 field_evaluate (uh, fops, omega, K, value);
93 for (size_type loc_idof = 0, loc_ndof = piola.size(); loc_idof < loc_ndof; loc_idof++) {
94 bbox.update (piola[loc_idof].F, value[loc_idof]);
95 }
96 for (size_type j = 0; j <= my_order; j++) {
97 for (size_type i1 = 0; i1 <= my_order; i1++) {
98 size_type i = std::min(i1, my_order-j);
99 size_type iloc = reference_element_t::ilat2loc_inod (my_order, ilat(i, j));
100 piola[iloc].F.put (gdat, dim); gdat << " " << value[iloc] << endl;
101 }
102 gdat << endl;
103 }
104 gdat << endl << endl;
105}
106template<class T>
107static
108void
109put_quadrangle (
110 std::ostream& gdat,
111 const geo_basic<T,sequential>& omega,
112 const geo_element& K,
113 const field_basic<T,sequential>& uh,
114 const fem_on_pointset<T>& fops,
115 size_t my_order,
116 bound_type<T>& bbox)
117{
118 using namespace std;
119 typedef typename geo_basic<T,sequential>::size_type size_type;
120 typedef point_basic<size_type> ilat;
121 size_type dim = omega.dimension();
122 const Eigen::Matrix<piola<T>,Eigen::Dynamic,1>& piola = fops.get_piola_on_pointset().get_piola (omega, K);
123 Eigen::Matrix<T,Eigen::Dynamic,1> value;
124 field_evaluate (uh, fops, omega, K, value);
125 for (size_type loc_idof = 0, loc_ndof = piola.size(); loc_idof < loc_ndof; loc_idof++) {
126 bbox.update (piola[loc_idof].F, value[loc_idof]);
127 }
128 for (size_type j = 0; j < my_order+1; j++) {
129 for (size_type i = 0; i < my_order+1; i++) {
130 size_type loc_idof00 = reference_element_q::ilat2loc_inod (my_order, ilat(i, j));
131 piola[loc_idof00].F.put (gdat, dim); gdat << " " << value[loc_idof00] << endl;
132 }
133 gdat << endl;
134 }
135 gdat << endl << endl;
136}
137template<class T>
138void
140 std::ostream& gdat,
141 const geo_basic<T,sequential>& omega,
142 const geo_element& K,
144 const fem_on_pointset<T>& fops,
145 size_t my_order,
146 bound_type<T>& bbox)
147{
148 switch (K.variant()) {
149 case reference_element::e: put_edge (gdat, omega, K, uh, fops, my_order, bbox); break;
150 case reference_element::t: put_triangle (gdat, omega, K, uh, fops, my_order, bbox); break;
151 case reference_element::q: put_quadrangle (gdat, omega, K, uh, fops, my_order, bbox); break;
152 default: error_macro ("unsupported element variant `"<<K.variant()<<"'");
153 }
154}
155// ----------------------------------------------------------------------------
156// scalar field puts
157// ----------------------------------------------------------------------------
158template <class T>
159odiststream&
161{
162 using namespace std;
163 typedef typename field_basic<T,sequential>::float_type float_type;
165 typedef point_basic<size_type> ilat;
166 ostream& os = ods.os();
167 bool verbose = iorheo::getverbose(os);
168 bool clean = iorheo::getclean(os);
169 bool execute = iorheo::getexecute(os);
170 bool fill = iorheo::getfill(os); // show grid or fill elements
171 bool elevation = iorheo::getelevation(os);
172 bool color = iorheo::getcolor(os);
173 bool gray = iorheo::getgray(os);
174 bool black_and_white = iorheo::getblack_and_white(os);
175 bool reader_on_stdin = iorheo::getreader_on_stdin(os);
176 string format = iorheo::getimage_format(os);
177 string basename = iorheo::getbasename(os);
178 size_type subdivide = iorheo::getsubdivide(os);
179 size_type n_isovalue = iorheo::getn_isovalue(os);
180 size_type n_isovalue_negative = iorheo::getn_isovalue_negative(os);
181 string outfile_fmt = "";
182 string tmp = get_tmpdir() + "/";
183 if (!clean) tmp = "";
184
185 const geo_basic<float_type,sequential>& omega = uh.get_geo();
186 size_type dim = omega.dimension();
187 size_type map_dim = omega.map_dimension();
188 size_type nv = omega.sizes().ownership_by_dimension[0].size();
189 size_type nedg = omega.sizes().ownership_by_dimension[1].size();
190 size_type nfac = omega.sizes().ownership_by_dimension[2].size();
191 size_type nvol = omega.sizes().ownership_by_dimension[3].size();
192 size_type ne = omega.sizes().ownership_by_dimension[map_dim].size();
193
194 const basis_basic<float_type>& b = uh.get_space().get_basis();
195 if (subdivide == 0) { // subdivide is unset: use default
196 subdivide = std::max(omega.order(), subdivide);
197 subdivide = std::max(b.degree (), subdivide);
198 }
199 basis_basic<T> subdivide_pointset ("P"+std::to_string(subdivide));
200 piola_on_pointset<T> pops; pops.initialize (omega.get_piola_basis(), subdivide_pointset, integrate_option());
201 fem_on_pointset<T> fops; fops.initialize (b, pops);
202
203 bound_type<T> bbox;
204 bbox.xmin = omega.xmin();
205 bbox.xmax = omega.xmax();
206 bbox.umin = uh.min();
207 bbox.umax = uh.max();
208
209 std::vector<T> values;
210 if (n_isovalue_negative != 0) {
211 for (size_t i = 0; i <= n_isovalue_negative; i++) {
212 values.push_back (bbox.umin*(n_isovalue_negative - i)/n_isovalue_negative);
213 }
214 for (size_t i = 1; i <= n_isovalue - n_isovalue_negative; i++) {
215 values.push_back (bbox.umax*i/(n_isovalue - n_isovalue_negative));
216 }
217 }
218 string filelist;
219 //
220 // output .gdat
221 //
222 string filename = tmp+basename + ".gdat";
223 string gdatname = filename;
224 filelist = filelist + " " + filename;
225 ofstream gdat (filename.c_str());
226 if (verbose) clog << "! file \"" << filename << "\" created.\n";
227 gdat << setprecision(numeric_limits<float_type>::digits10);
228 size_type used_dim = (fill ? map_dim : 1);
229 for (size_type ie = 0, ne = omega.size(used_dim); ie < ne; ie++) {
230 const geo_element& K = omega.get_geo_element(used_dim,ie);
231 put (gdat, omega, K, uh, fops, subdivide, bbox);
232 }
233 gdat.close();
234 //
235 // rounding bounds
236 //
237 T eps = 1e-7;
238 bbox.umin = floorer(eps) (bbox.umin);
239 bbox.umax = ceiler(eps) (bbox.umax);
240 if (fabs(bbox.umax - bbox.umin) < eps) { // empty range: constant field
241 bbox.umax = 1.1*bbox.umax;
242 bbox.umin = 0.9*bbox.umin;
243 }
244 //
245 // output .plot
246 //
247 filename = tmp+basename + ".plot";
248 filelist = filelist + " " + filename;
249 ofstream plot (filename.c_str());
250 if (verbose) clog << "! file \"" << filename << "\" created.\n";
251
252 plot << "#!gnuplot" << endl
253 << setprecision(numeric_limits<float_type>::digits10);
254 if (format != "") {
255 outfile_fmt = basename + "." + format;
256 string terminal = format;
257 if (terminal == "ps") {
258 terminal = "postscript eps";
259 if (color) terminal += " color";
260 }
261 if (terminal == "jpg") terminal = "jpeg";
262 if (terminal == "jpeg" || terminal == "png" || terminal == "gif") {
263 terminal += " crop";
264 }
265 plot << "set terminal " << terminal << endl
266 << "set output \"" << outfile_fmt << "\"" << endl;
267 }
268 plot << "umin = " << bbox.umin << endl
269 << "umax = " << bbox.umax << endl
270 << "n_isovalue = " << n_isovalue << endl
271 << "uincr = 1.0*(umax-umin)/n_isovalue" << endl;
272 if (dim == 2) {
273 if (bbox.xmin[0] >= bbox.xmax[0]) plot << "#";
274 plot << "set xrange [" << bbox.xmin[0] << ":" << bbox.xmax[0] << "]" << endl;
275 if (bbox.xmin[1] >= bbox.xmax[1]) plot << "#";
276 plot << "set yrange [" << bbox.xmin[1] << ":" << bbox.xmax[1] << "]" << endl;
277 if (bbox.umin >= bbox.umax) plot << "#";
278 plot << "set zrange [umin:umax]" << endl;
279 if (bbox.xmin[0] >= bbox.xmax[0] || bbox.xmin[1] >= bbox.xmax[1]) {
280 plot << "set size square" << endl;
281 } else {
282 plot << "set size ratio -1 # equal scales" << endl
283 << "set noxtics" << endl
284 << "set noytics" << endl;
285 }
286 if (elevation) {
287 plot << "set xyplane at umin-0.1*(umax-umin)" << endl;
288 } else {
289 plot << "set view map" << endl;
290 }
291 if (map_dim == 2) {
292 plot << "unset surface" << endl
293 << "set contour base" << endl;
294 if (values.size() != 0) {
295 plot << "set cntrparam levels discrete ";
296 for (size_t i = 0, n = values.size(); i < n; i++) {
297 plot << values[i];
298 if (i+1 != n) plot << ", ";
299 }
300 plot << endl;
301 } else {
302 plot << "eps = (umax-umin)*1e-4"<< endl;
303 plot << "uincr_eps = 1.0*(umax-umin-2*eps)/n_isovalue"<< endl;
304 plot << "set cntrparam levels incremental umin+eps,uincr_eps,umax-eps"<< endl;
305 }
306 }
307 if (black_and_white) {
308 plot << "set cbtics out scale 0.5" << endl;
309 }
310 plot << "set cbtics uincr" << endl
311 << "set cbrange [umin:umax]" << endl;
312 if (gray) {
313 plot << "set palette gray" << endl;
314 } else if (color) {
315 plot << "set palette rgbformulae 33,13,-4" << endl;
316 } else { // bw
317 plot << "set palette rgbformulae 0,0,0" << endl;
318 }
319 plot << "set palette maxcolors n_isovalue+1" << endl
320 << "set nokey" << endl;
321 } else if (dim == 3 && map_dim == 2) {
322 // field defined on a 3D surface
323 point_basic<T> dx = 0.1*(omega.xmax() - omega.xmin());
324 T dx_max = max(dx[0],max(dx[1],dx[2]));
325 if (dx_max == 0) dx_max = 0.1;
326 dx[0] = max(dx[0],dx_max);
327 if (omega.dimension() >= 2) dx[1] = max(dx[1],dx_max);
328 if (omega.dimension() == 3) dx[2] = max(dx[2],dx_max);
329 point_basic<T> xmin = omega.xmin() - dx;
330 point_basic<T> xmax = omega.xmax() + dx;
331 // TODO: visu contours discrets comme en 2d classique
332 plot << "set xrange [" << xmin[0] << ":" << xmax[0] << "]" << endl
333 << "set yrange [" << xmin[1] << ":" << xmax[1] << "]" << endl
334 << "set zrange [" << xmin[2] << ":" << xmax[2] << "]" << endl
335 << "set xyplane at " << xmin[2] << endl
336 << "set view equal xyz # equal scales" << endl
337 << "set view 70,120" << endl
338 << "n_isovalue = " << n_isovalue << endl
339 << "n_subdivide = 40 " << endl
340 << "uincr = 1.0*(umax-umin)/n_isovalue" << endl
341 << "set cbtics uincr" << endl
342 << "set pm3d interpolate n_subdivide,n_subdivide corners2color mean" << endl
343 << "set palette rgbformulae 33,13,-4 maxcolors n_isovalue" << endl;
344 if (format != "") {
345 plot << "set noxlabel" << endl
346 << "set noylabel" << endl
347 << "set nozlabel" << endl;
348 } else {
349 plot << "set xlabel \"x\"" << endl
350 << "set ylabel \"y\"" << endl
351 << "set zlabel \"z\"" << endl;
352 }
353 }
354 if (dim == 1) {
355 if (color) {
356 plot << "set colors classic" << endl;
357 }
358 plot << "plot \"" << gdatname << "\" notitle with lines lw 2";
359 if (gray || black_and_white) {
360 plot << " lc 0" << endl;
361 }
362 plot << endl;
363 } else {
364 if (!fill && dim == 2) {
365 plot << "plot";
366 } else {
367 plot << "splot";
368 }
369 plot << " \"" << gdatname << "\" notitle";
370 if (map_dim == 2) {
371 if (black_and_white) {
372 plot << " with lines palette lw 2" << endl;
373 } else {
374 plot << " with lines palette" << endl;
375 }
376 } else { // a 2d line
377 plot << " with lines palette lw 2" << endl;
378 }
379 }
380 //
381 // end of plot
382 //
383 if (format == "" && !reader_on_stdin) {
384 plot << "pause -1 \"<return>\"\n";
385 }
386 plot.close();
387 //
388 // run gnuplot
389 //
390 int status = 0;
391 string command;
392 if (execute) {
393 command = "gnuplot ";
394 if (reader_on_stdin) command += "-persist ";
395 command += tmp + basename + ".plot";
396 if (verbose) clog << "! " << command << endl;
397 cin.sync();
398 status = system (command.c_str());
399 if (format != "") {
400 check_macro (file_exists (outfile_fmt), "! file \"" << outfile_fmt << "\" creation failed");
401 if (verbose) clog << "! file \"" << outfile_fmt << "\" created" << endl;
402 }
403 }
404 //
405 // clear gnuplot data
406 //
407 if (clean) {
408 command = "/bin/rm -f " + filelist;
409 if (verbose) clog << "! " << command << endl;
410 status = system (command.c_str());
411 }
412 return ods;
413}
414// ----------------------------------------------------------------------------
415// vector field puts
416// ----------------------------------------------------------------------------
417template <class T>
418odiststream&
420{
422 Float vscale = iorheo::getvectorscale(ods.os());
423 const geo_basic<T,sequential>& omega = uh.get_geo();
424 // TODO: non-isoparam => interpolate
425 check_macro (omega.get_piola_basis().degree() == uh.get_space().get_basis().degree(),
426 "gnuplot vector: unsupported non-isoparametric approx " << uh.get_space().get_basis().name());
427 size_type n_elt = omega.size();
428 size_type d = omega.dimension();
429 T diam_omega = norm (omega.xmax() - omega.xmin()); // characteristic length in omega
430 T h_moy = diam_omega/pow(n_elt,1./d);
431 space_basic<T,sequential> Xh (uh.get_geo(), "P"+std::to_string(uh.get_space().degree()));
432 field_basic<T,sequential> norm_uh = interpolate(Xh, norm(uh));
433 T norm_max_uh = norm_uh.max_abs(); // get max vector length
434 if (norm_max_uh + 1 == 1) norm_max_uh = 1;
435#ifdef TODO
436 T scale = vscale*(h_moy/norm_max_uh);
437#endif // TODO
438 size_type n_comp = uh.get_space().get_basis().size();
439 disarray<point_basic<T>, sequential> x = omega.get_nodes();
440 for (size_type inod = 0, nnod = x.size(); inod < nnod; inod++) {
442 for (size_type i_comp = 0; i_comp < n_comp; i_comp++) {
443 size_type idof = n_comp*inod + i_comp;
444 vi[i_comp] = uh.dof(idof);
445 }
446 x[inod] += vscale*vi;
447 }
448 geo_basic<T,sequential> deformed_omega = omega;
449 deformed_omega.set_nodes(x);
450 space_basic<T,sequential> deformed_Vh (deformed_omega, norm_uh.get_space().get_basis().name());
451 field_basic<T,sequential> deformed_norm_uh (deformed_Vh);
452 std::copy (norm_uh.begin_dof(), norm_uh.end_dof(), deformed_norm_uh.begin_dof());
453 visu_gnuplot (ods, deformed_norm_uh);
454 return ods;
455}
456// ----------------------------------------------------------------------------
457// switch
458// ----------------------------------------------------------------------------
459template <class T>
460odiststream&
462{
463 switch (uh.get_space().valued_tag()) {
464 case space_constant::scalar: visu_gnuplot_scalar (ods, uh); break;
465 case space_constant::vector: visu_gnuplot_vector (ods, uh); break;
466 default: error_macro ("do not known how to print " << uh.valued() << "-valued field");
467 }
468 return ods;
469}
470// ----------------------------------------------------------------------------
471// instanciation in library
472// ----------------------------------------------------------------------------
473#define _RHEOLEF_instanciate(T) \
474template odiststream& visu_gnuplot<T> (odiststream&, const field_basic<T,sequential>&);
475
477#undef _RHEOLEF_instanciate
478} // rheolef namespace
field::size_type size_type
Definition branch.cc:430
see the Float page for the full documentation
see the disarray page for the full documentation
Definition disarray.h:497
void initialize(const basis_basic< T > &fem_basis, const piola_on_pointset< T > &pops)
T max_abs() const
Definition field.h:834
T & dof(size_type idof)
Definition field.h:738
iterator begin_dof()
Definition field.h:595
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
iterator end_dof()
Definition field.h:603
const std::string & valued() const
Definition field.h:274
typename float_traits< T >::type float_type
Definition field.h:228
generic mesh with rerefence counting
Definition geo.h:1089
see the geo_element page for the full documentation
variant_type variant() const
see the integrate_option page for the full documentation
odiststream: see the diststream page for the full documentation
Definition diststream.h:137
std::ostream & os()
Definition diststream.h:247
void initialize(const basis_basic< T > &piola_basis, const quadrature< T > &quad, const integrate_option &iopt)
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
static const variant_type q
static const variant_type e
static const variant_type t
the finite element space
Definition space.h:382
rheolef::std value
#define _RHEOLEF_instanciate(T)
Definition csr_seq.cc:491
#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 gray
This file is part of Rheolef.
void field_evaluate(const field_basic< T, M > &uh, const basis_on_pointset< T > &bops, reference_element hat_K, const std::vector< size_t > &dis_idof, Eigen::Matrix< T, Eigen::Dynamic, 1 > &value)
odiststream & visu_gnuplot_vector(odiststream &ods, const field_basic< T, sequential > &uh)
std::string get_tmpdir()
get_tmpdir: see the rheostream page for the full documentation
Definition rheostream.cc:54
space_mult_list< T, M > pow(const space_basic< T, M > &X, size_t n)
Definition space_mult.h:120
void put(std::ostream &out, std::string name, const tiny_matrix< T > &a)
Definition tiny_lu.h:155
odiststream & visu_gnuplot_scalar(odiststream &ods, const field_basic< T, sequential > &uh)
floorer_type< T > floorer(const T &prec)
Definition rounder.h:61
ceiler_type< T > ceiler(const T &prec)
Definition rounder.h:62
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
bool file_exists(const std::string &filename)
file_exists: see the rheostream page for the full documentation
Definition scatch.icc:34
odiststream & visu_gnuplot(odiststream &, const field_basic< T, sequential > &)
T norm(const vec< T, M > &x)
norm(x): see the expression page for the full documentation
Definition vec.h:387
STL namespace.