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