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