rheolef  6.5
field_seq_visu_gnuplot.cc
Go to the documentation of this file.
1 # include "rheolef/field.h"
2 # include "rheolef/field_indirect.h"
3 # include "rheolef/field_expr.h"
4 # include "rheolef/field_expr_ops.h"
5 # include "rheolef/field_nonlinear_expr_ops.h"
6 # include "rheolef/interpolate.h"
7 # include "rheolef/field_evaluate.h"
8 # include "rheolef/piola.h"
9 # include "rheolef/rounder.h"
10 # include "rheolef/rheostream.h"
11 
12 namespace rheolef {
13 
14 template <class T>
15 odiststream& visu_gnuplot (odiststream& ops, const geo_basic<T,sequential>& omega);
16 
17 template <class T>
18 struct bound_type {
20  T umin, umax;
21  void update (const point_basic<T>& x, const T& u) {
22  for (size_t i = 0; i < 3; i++) {
23  xmin[i] = std::min(xmin[i], x[i]);
24  xmax[i] = std::max(xmax[i], x[i]);
25  }
26  umin = std::min(umin, u);
27  umax = std::max(umax, u);
28  }
29 };
30 template<class T>
31 static
32 void
33 put_edge (std::ostream& gdat, const geo_element& E, const geo_basic<T,sequential>& omega,
34  const field_basic<T,sequential>& uh,
35  const basis_on_pointset<T>& geo_on_pointset,
36  const basis_on_pointset<T>& field_on_pointset,
37  size_t my_order,
38  bound_type<T>& bbox)
39 {
40  using namespace std;
42  typedef point_basic<size_type> ilat;
43  size_type dim = omega.dimension();
44  std::vector<size_type> dis_idof1;
45  uh.get_space().dis_idof (E, dis_idof1);
46  std::vector<size_type> dis_inod;
47  omega.dis_inod (E, dis_inod);
48  std::vector<point_basic<T> > hat_xi;
49  uh.get_space().get_numbering().get_basis().hat_node (E.variant(), hat_xi);
50  for (size_type i = 0; i <= my_order; i++) {
51  size_type iloc = reference_element_e::ilat2loc_inod (my_order, ilat(i));
52  point_basic<T> xi = piola_transformation (omega, geo_on_pointset, E.variant(), dis_inod, iloc);
53  T ui = field_evaluate (uh, field_on_pointset, E.variant(), dis_idof1, iloc);
54  xi.put (gdat, dim); gdat << " " << ui << endl;
55  bbox.update (xi, ui);
56  }
57  gdat << endl << endl;
58 }
59 template<class T>
60 static
61 void
62 put_triangle (std::ostream& gdat, const geo_element& F, const geo_basic<T,sequential>& omega,
63  const field_basic<T,sequential>& uh,
64  const basis_on_pointset<T>& geo_on_pointset,
65  const basis_on_pointset<T>& field_on_pointset,
66  size_t my_order,
67  bound_type<T>& bbox)
68 {
69  using namespace std;
71  typedef point_basic<size_type> ilat;
72  size_type dim = omega.dimension();
73  std::vector<size_type> dis_idof1;
74  std::vector<size_type> dis_inod;
75  omega.dis_inod (F, dis_inod);
76  uh.get_space().dis_idof (F, dis_idof1);
77  size_type nloc = geo_on_pointset.size(F);
78  std::vector<point_basic<T> > xdof (nloc);
79  std::vector<T> udof (nloc);
80  for (size_type iloc = 0; iloc < nloc; iloc++) {
81  xdof[iloc] = piola_transformation (omega, geo_on_pointset, F.variant(), dis_inod, iloc);
82  udof[iloc] = field_evaluate (uh, field_on_pointset, F.variant(), dis_idof1, iloc);
83  bbox.update (xdof[iloc], udof[iloc]);
84  }
85  for (size_type j = 0; j <= my_order; j++) {
86  for (size_type i1 = 0; i1 <= my_order; i1++) {
87  size_type i = std::min(i1, my_order-j);
88  size_type iloc = reference_element_t::ilat2loc_inod (my_order, ilat(i, j));
89  xdof [iloc].put (gdat, dim); gdat << " " << udof[iloc] << endl;
90  }
91  gdat << endl;
92  }
93  gdat << endl << endl;
94 }
95 template<class T>
96 static
97 void
98 put_quadrangle (std::ostream& gdat, const geo_element& F, const geo_basic<T,sequential>& omega,
100 {
101  using namespace std;
103  typedef point_basic<size_type> ilat;
104  size_type dim = omega.dimension();
105  size_type degree = uh.get_space().get_numbering().get_basis().degree();
106  std::vector<size_type> dis_idof1;
107  uh.get_space().dis_idof (F, dis_idof1);
108  std::vector<size_type> dis_inod;
109  omega.dis_inod (F, dis_inod);
110  std::vector<point_basic<T> > xdof (dis_idof1.size());
111  for (size_type loc_idof = 0, loc_ndof = xdof.size(); loc_idof < loc_ndof; loc_idof++) {
112  xdof[loc_idof] = piola_transformation (omega, b, F.variant(), dis_inod, loc_idof);
113  }
114  for (size_type j = 0; j < degree+1; j++) {
115  for (size_type i = 0; i < degree+1; i++) {
116  size_type loc_idof00 = reference_element_q::ilat2loc_inod (degree, ilat(i, j));
117  xdof [loc_idof00].put (gdat, dim); gdat << " " << uh.dof(dis_idof1[loc_idof00]) << endl;
118  }
119  gdat << endl;
120  }
121  gdat << endl << endl;
122 }
123 template<class T>
124 static
125 void
126 put (std::ostream& gdat, const geo_element& K, const geo_basic<T,sequential>& omega,
127  const field_basic<T,sequential>& uh,
128  const basis_on_pointset<T>& geo_on_pointset,
129  const basis_on_pointset<T>& field_on_pointset,
130  size_t my_order,
131  bound_type<T>& bbox)
132 {
133  switch (K.variant()) {
134  case reference_element::e: put_edge (gdat, K, omega, uh, geo_on_pointset, field_on_pointset, my_order, bbox); break;
135  case reference_element::t: put_triangle (gdat, K, omega, uh, geo_on_pointset, field_on_pointset, my_order, bbox); break;
136  case reference_element::q: put_quadrangle (gdat, K, omega, uh, geo_on_pointset); break;
137  default: error_macro ("unsupported element variant `"<<K.variant()<<"'");
138  }
139 }
140 template <class T>
141 odiststream&
143 {
144  using namespace std;
145  typedef typename field_basic<T,sequential>::float_type float_type;
147  typedef point_basic<size_type> ilat;
148  ostream& os = ods.os();
149  bool verbose = iorheo::getverbose(os);
150  bool clean = iorheo::getclean(os);
151  bool execute = iorheo::getexecute(os);
152  bool fill = iorheo::getfill(os); // show grid or fill elements
153  bool elevation = iorheo::getelevation(os);
154  bool color = iorheo::getcolor(os);
155  bool gray = iorheo::getgray(os);
156  bool black_and_white = iorheo::getblack_and_white(os);
157  bool reader_on_stdin = iorheo::getreader_on_stdin(os);
158  string format = iorheo::getimage_format(os);
159  string basename = iorheo::getbasename(os);
160  size_type subdivide = iorheo::getsubdivide(os);
161  size_type n_isovalue = iorheo::getn_isovalue(os);
162  size_type n_isovalue_negative = iorheo::getn_isovalue_negative(os);
163  string outfile_fmt = "";
164  string tmp = get_tmpdir() + "/";
165  if (!clean) tmp = "";
166 
167  const geo_basic<float_type,sequential>& omega = uh.get_geo();
168  size_type dim = omega.dimension();
169  size_type map_dim = omega.map_dimension();
170  size_type nv = omega.sizes().ownership_by_dimension[0].size();
171  size_type nedg = omega.sizes().ownership_by_dimension[1].size();
172  size_type nfac = omega.sizes().ownership_by_dimension[2].size();
173  size_type nvol = omega.sizes().ownership_by_dimension[3].size();
174  size_type ne = omega.sizes().ownership_by_dimension[map_dim].size();
175 
176  const numbering<float_type,sequential>& fem = uh.get_space().get_numbering();
177  if (subdivide == 0) { // subdivide is unset: use default
178  subdivide = std::max(omega.order(), subdivide);
179  subdivide = std::max(fem.degree (), subdivide);
180  }
181  basis_basic<T> subdivide_pointset ("P"+itos(subdivide));
182  basis_on_pointset<T> geo_on_pointset (subdivide_pointset, omega.get_piola_basis());
183  basis_on_pointset<T> field_on_pointset (subdivide_pointset, fem.get_basis());
184 
185  bound_type<T> bbox;
186  bbox.xmin = omega.xmin();
187  bbox.xmax = omega.xmax();
188  bbox.umin = uh.min();
189  bbox.umax = uh.max();
190 
191  std::vector<T> values;
192  if (n_isovalue_negative != 0) {
193  for (size_t i = 0; i <= n_isovalue_negative; i++) {
194  values.push_back (bbox.umin*(n_isovalue_negative - i)/n_isovalue_negative);
195  }
196  for (size_t i = 1; i <= n_isovalue - n_isovalue_negative; i++) {
197  values.push_back (bbox.umax*i/(n_isovalue - n_isovalue_negative));
198  }
199  }
200  string filelist;
201  string filename = tmp+basename + ".gdat";
202  string gdatname = filename;
203  filelist = filelist + " " + filename;
204  ofstream gdat (filename.c_str());
205  if (verbose) clog << "! file \"" << filename << "\" created.\n";
206  gdat << setprecision(numeric_limits<float_type>::digits10);
207  size_type used_dim = (fill ? map_dim : 1);
208  for (size_type ie = 0, ne = omega.size(used_dim); ie < ne; ie++) {
209  const geo_element& K = omega.get_geo_element(used_dim,ie);
210  put (gdat, K, omega, uh, geo_on_pointset, field_on_pointset, subdivide, bbox);
211  }
212  gdat.close();
213  Float eps = 1e-7;
214  bbox.umin = floorer(eps) (bbox.umin);
215  bbox.umax = ceiler(eps) (bbox.umax);
216  filename = tmp+basename + ".plot";
217  filelist = filelist + " " + filename;
218  ofstream plot (filename.c_str());
219  if (verbose) clog << "! file \"" << filename << "\" created.\n";
220 
221  plot << "#!gnuplot" << endl
222  << setprecision(numeric_limits<float_type>::digits10);
223  if (format != "") {
224  outfile_fmt = basename + "." + format;
225  string terminal = format;
226  if (terminal == "ps") {
227  terminal = "postscript eps";
228  if (color) terminal += " color";
229  }
230  if (terminal == "jpg") terminal = "jpeg";
231  if (terminal == "jpeg" || terminal == "png" || terminal == "gif") {
232  terminal += " crop";
233  }
234  plot << "set terminal " << terminal << endl
235  << "set output \"" << outfile_fmt << "\"" << endl;
236  }
237  if (!black_and_white && map_dim == dim) {
238  plot << "set noborder" << endl;
239  }
240  if (dim == 2) {
241  plot << "umin = " << bbox.umin << endl
242  << "umax = " << bbox.umax << endl;
243  if (bbox.xmin[0] >= bbox.xmax[0]) plot << "#";
244  plot << "set xrange [" << bbox.xmin[0] << ":" << bbox.xmax[0] << "]" << endl;
245  if (bbox.xmin[1] >= bbox.xmax[1]) plot << "#";
246  plot << "set yrange [" << bbox.xmin[1] << ":" << bbox.xmax[1] << "]" << endl;
247  if (bbox.umin >= bbox.umax) plot << "#";
248  plot << "set zrange [umin:umax]" << endl;
249  if (bbox.xmin[0] >= bbox.xmax[0] || bbox.xmin[1] >= bbox.xmax[1]) {
250  plot << "set size square" << endl;
251  } else {
252  plot << "set size ratio -1 # equal scales" << endl
253  << "set noxtics" << endl
254  << "set noytics" << endl;
255  }
256  if (elevation) {
257  plot << "set xyplane at umin-0.1*(umax-umin)" << endl;
258  } else {
259  plot << "set view map" << endl;
260  }
261  if (color || gray) {
262  plot << "set pm3d interpolate 10,10 corners2color mean" << endl;
263  if (gray) {
264  plot << "set palette gray" << endl;
265  } else {
266  plot << "set palette rgbformulae 33,13,-4" << endl;
267  }
268  } else { // black-and-white
269  if (values.size() != 0) {
270  plot << "set cntrparam levels discrete ";
271  for (size_t i = 0, n = values.size(); i < n; i++) {
272  plot << values[i];
273  if (i != n-1) plot << ", ";
274  }
275  plot << endl;
276  } else {
277  plot << "set cntrparam levels incremental umin,(umax-umin)/10.0,umax"<< endl;
278  }
279  plot << "unset surface" << endl
280  << "set contour base" << endl
281  << "set nokey" << endl
282  << "set palette rgbformulae 0,0,0" << endl
283  << "unset colorbox" << endl;
284  }
285  } else if (dim == 3 && map_dim == 2) {
286  point_basic<T> dx = 0.1*(omega.xmax() - omega.xmin());
287  T dx_max = max(dx[0],max(dx[1],dx[2]));
288  if (dx_max == 0) dx_max = 0.1;
289  dx[0] = max(dx[0],dx_max);
290  if (omega.dimension() >= 2) dx[1] = max(dx[1],dx_max);
291  if (omega.dimension() == 3) dx[2] = max(dx[2],dx_max);
292  point_basic<T> xmin = omega.xmin() - dx;
293  point_basic<T> xmax = omega.xmax() + dx;
294  plot << "set xrange [" << xmin[0] << ":" << xmax[0] << "]" << endl
295  << "set yrange [" << xmin[1] << ":" << xmax[1] << "]" << endl
296  << "set zrange [" << xmin[2] << ":" << xmax[2] << "]" << endl
297  << "set xyplane at " << xmin[2] << endl
298  << "set view equal xyz # equal scales" << endl
299  << "set view 70,120" << endl
300  << "set pm3d interpolate 5,5 corners2color mean" << endl
301  << "set palette rgbformulae 33,13,-4" << endl;
302  if (format != "") {
303  plot << "set noxlabel" << endl
304  << "set noylabel" << endl
305  << "set nozlabel" << endl;
306  } else {
307  plot << "set xlabel \"x\"" << endl
308  << "set ylabel \"y\"" << endl
309  << "set zlabel \"z\"" << endl;
310  }
311  }
312  if (dim == 1) {
313  plot << "plot \"" << gdatname << "\" notitle with lines";
314  if (black_and_white) {
315  plot << " lc 0" << endl;
316  }
317  plot << endl;
318  } else {
319  if (!fill && dim == 2) {
320  plot << "plot";
321  } else {
322  plot << "splot";
323  }
324  plot << " \"" << gdatname << "\" notitle";
325  if (map_dim == 2) {
326  if (!black_and_white && fill) {
327  plot << " with pm3d" << endl;
328  } else {
329  plot << " with lines palette" << endl;
330  }
331  } else { // a 2d line
332  plot << " with lines palette lw 2" << endl;
333  }
334  }
335  if (format == "" && !reader_on_stdin) {
336  plot << "pause -1 \"<return>\"\n";
337  }
338  plot.close();
339  int status = 0;
340  string command;
341  if (execute) {
342  command = "gnuplot ";
343  if (reader_on_stdin) command += "-persist ";
344  command += tmp + basename + ".plot";
345  if (verbose) clog << "! " << command << endl;
346  cin.sync();
347  status = system (command.c_str());
348  if (format != "") {
349  check_macro (file_exists (outfile_fmt), "! file \"" << outfile_fmt << "\" creation failed");
350  if (verbose) clog << "! file \"" << outfile_fmt << "\" created" << endl;
351  }
352  }
353  if (clean) {
354  command = "/bin/rm -f " + filelist;
355  if (verbose) clog << "! " << command << endl;
356  status = system (command.c_str());
357  }
358  return ods;
359 }
360 template <class T>
361 odiststream&
363 {
365  const geo_basic<T,sequential>& omega = uh.get_geo();
366  check_macro (omega.get_piola_basis().name() == uh.get_space().get_numbering().name(),
367  "gnuplot vector: unsupported non-isoparametric approx " << uh.get_space().get_numbering().name());
368  size_type n_elt = omega.size();
369  size_type d = omega.dimension();
370  Float diam_omega = norm (omega.xmax() - omega.xmin()); // characteristic length in omega
371  Float h_moy = diam_omega/pow(n_elt,1./d);
373  field_basic<T,sequential> norm_uh = interpolate(Xh, norm(uh));
374  Float norm_max_uh = norm_uh.max_abs(); // get max vector length
375  if (norm_max_uh + 1 == 1) norm_max_uh = 1;
376  Float scale = h_moy/norm_max_uh;
377  size_type n_comp = uh.size();
378  std::vector<field_component_const<T,sequential> > uh_comp (n_comp);
379  for (size_type i_comp = 0; i_comp < n_comp; i_comp++) {
380  uh_comp[i_comp] = uh[i_comp];
381  }
382  array<point_basic<T>, sequential> x = omega.get_nodes();
383  for (size_type inod = 0, nnod = x.size(); inod < nnod; inod++) {
384  point_basic<T> vi;
385  for (size_type i_comp = 0; i_comp < n_comp; i_comp++) {
386  vi[i_comp] = uh[i_comp].dof(inod); // TODO: non-isoparam => interpolate
387  }
388  x[inod] += vi;
389  }
390  geo_basic<T,sequential> deformed_omega = omega;
391  deformed_omega.set_nodes(x);
392  space_basic<T,sequential> deformed_Vh (deformed_omega, norm_uh.get_space().get_numbering().name());
393  field_basic<T,sequential> deformed_norm_uh (deformed_Vh);
394  std::copy (norm_uh.begin_dof(), norm_uh.end_dof(), deformed_norm_uh.begin_dof());
395  visu_gnuplot (ods, deformed_norm_uh);
396  return ods;
397 }
398 template <class T>
399 odiststream&
401 {
402  switch (uh.get_space().valued_tag()) {
403  case space_constant::scalar: visu_gnuplot_scalar (ods, uh); break;
404  case space_constant::vector: visu_gnuplot_vector (ods, uh); break;
405  default: error_macro ("do not known how to print " << uh.valued() << "-valued field");
406  }
407  return ods;
408 }
409 template odiststream& visu_gnuplot<Float> (odiststream&, const field_basic<Float,sequential>&);
410 
411 } // rheolef namespace
412