rheolef  7.0
field_seq_visu_vtk_mayavi.cc
Go to the documentation of this file.
1 #include "rheolef/field.h"
2 #include "rheolef/piola.h"
3 #include "rheolef/rheostream.h"
4 #include "rheolef/iorheo.h"
5 #include "rheolef/iofem.h"
6 using namespace std;
7 namespace rheolef {
8 
9 template<class T>
10 static
11 std::string
12 python (const point_basic<T>& x, size_t d=3)
13 {
14  std::ostringstream os;
15  os << "(" << x[0];
16  for (size_t i = 1; i < d; i++)
17  os << ", " << x[i];
18  os << ")" << std::flush;
19  string buffer = os.str();
20  return buffer;
21 }
22 template <class T> odiststream& field_put_vtk (odiststream&, const field_basic<T,sequential>&);
23 
24 template <class T>
27 {
28  using namespace std;
31  typedef point_basic<size_type> ilat;
32  ostream& os = ops.os();
33  string valued = uh.get_space().valued();
34  bool is_scalar = (valued == "scalar");
35  bool verbose = iorheo::getverbose(os);
36  bool clean = iorheo::getclean(os);
37  bool execute = iorheo::getexecute(os);
38  bool fill = iorheo::getfill(os); // isocontours or color fill
39  bool elevation = iorheo::getelevation(os);
40  bool color = iorheo::getcolor(os);
41  bool gray = iorheo::getgray(os);
42  bool black_and_white = iorheo::getblack_and_white(os);
43  bool stereo = iorheo::getstereo(os);
44  bool volume = iorheo::getvolume(os);
45  bool iso = true;
46  bool cut = iorheo::getcut(os);
47  bool grid = iorheo::getgrid(os);
48  string format = iorheo::getimage_format(os);
49  string basename = iorheo::getbasename(os);
50  string mark = iorheo::getmark(os);
51  if (mark != "" && !is_scalar) mark = "|"+mark +"|";
52  bool velocity = iorheo::getvelocity(os);
53  bool deformation = iorheo::getdeformation(os);
54  Float vscale = iorheo::getvectorscale(os);
55  size_type n_isovalue = iorheo::getn_isovalue(os);
56  size_type n_isovalue_negative = iorheo::getn_isovalue_negative(os);
57  point_basic<float_type> origin = iofem::getorigin(os);
58  point_basic<float_type> normal = iofem::getnormal(os);
59  if (black_and_white) {
60  fill = false;
61  }
62  string outfile_fmt = "";
63  string tmp = get_tmpdir() + "/";
64  if (!clean) tmp = "";
65 
66  const geo_basic<float_type,sequential>& omega = uh.get_geo();
67  size_type dim = omega.dimension();
68  size_type map_dim = omega.map_dimension();
69  size_type nv = omega.sizes().ownership_by_dimension[0].size();
70  size_type nedg = omega.sizes().ownership_by_dimension[1].size();
71  size_type nfac = omega.sizes().ownership_by_dimension[2].size();
72  size_type nvol = omega.sizes().ownership_by_dimension[3].size();
73  size_type ne = omega.sizes().ownership_by_dimension[map_dim].size();
74 
75  const numbering<float_type,sequential>& fem = uh.get_space().get_numbering();
76  basis_on_pointset<float_type> b (fem.get_basis(), omega.get_piola_basis());
77  string filelist;
78  string filename = tmp+basename + ".vtk";
79  filelist = filelist + " " + filename;
80  ofstream vtk_os (filename.c_str());
81  odiststream vtk (vtk_os);
82  if (verbose) clog << "! file \"" << filename << "\" created.\n";
83  field_put_vtk (vtk, uh);
84  vtk.close();
85  std::string py_name = filename = tmp+basename + ".py";
86  filelist = filelist + " " + filename;
87  ofstream py (filename.c_str());
88  if (verbose) clog << "! file \"" << filename << "\" created.\n";
89  string style = (valued == "tensor" ? "tensor" : (velocity ? "velocity" : "deformation"));
90  py << "#!/usr/bin/env mayavi2" << endl
91  << "# This is a mayavi script for the visualization of " << basename << ".vtk" << endl
92  << "# automatically generated by rheolef." << endl
93  << endl
94  ;
95  py << "from sys import version_info" << endl
96  << "# mayavi2.py file has changed from directory since python 2.7:" << endl
97  << "if version_info[0] >= 2 and version_info[1] >= 7 and version_info[2] >= 2:" << endl
98  << " from mayavi.scripts import mayavi2" << endl
99  << "else:" << endl
100  << " from enthought.mayavi.scripts import mayavi2" << endl
101  << "mayavi2.standalone(globals()) # start mayvi2 graphic interface" << endl
102  << "from mayavi2_rheolef import * # load rheolef specific functions" << endl
103  << endl
104  ;
105  py << "opt = { \\" << endl
106  << " 'label' : '" << mark << "', \\" << endl
107  << " 'axis' : 1," << endl
108  << " 'color' : '" << (color ? "color" : (gray ? "gray" : "black_and_white")) << "'," << endl
109  << " 'stereo' : " << stereo << ", \\" << endl
110  ;
111  if (valued == "vector") {
112  py << " 'style' : '" << style << "', \\" << endl
113  ;
114  }
115  py << " 'scale' : " << vscale << ", \\" << endl
116  << " 'elevation' : " << elevation << ", \\" << endl
117  ;
118  if (valued == "scalar") {
119  py << " 'volume' : " << volume << ", \\" << endl
120  << " 'iso' : " << iso << ", \\" << endl
121  << " 'n_isovalue' : " << n_isovalue << ", \\" << endl
122  << " 'n_isovalue_negative' : " << n_isovalue_negative << ", \\" << endl
123  << " 'cut' : " << cut << ", \\" << endl
124  << " 'origin' : " << python(origin) << ", \\" << endl
125  << " 'normal' : " << python(normal) << ", \\" << endl
126  ;
127  }
128  py << " 'grid' : " << grid << ", \\" << endl
129  << " 'fill' : " << fill << ", \\" << endl
130  << " 'view_1d' : 0, \\" << endl
131  << " 'view_2d' : " << (dim == 2) << ", \\" << endl
132  << " 'view_map' : " << (dim > map_dim) << " \\" << endl
133  << " }" << endl
134  << endl
135  ;
136  py << "mayavi2_field_" << valued << "(mayavi, \"" << tmp+basename << "\", opt)" << endl
137  << endl
138  ;
139  py.close();
140  int status = 0;
141  string command;
142  if (execute) {
143  command = "LANG=C PYTHONPATH=" + string(_RHEOLEF_PKGDATADIR) + " mayavi2 -x " + py_name;
144  if (verbose) clog << "! " << command << endl;
145  status = system (command.c_str());
146  }
147  if (clean) {
148  command = "/bin/rm -f " + filelist;
149  if (verbose) clog << "! " << command << endl;
150  status = system (command.c_str());
151  }
152  return ops;
153 }
154 template <class T> idiststream& geo_get_vtk (idiststream&, geo_basic<T,sequential>&);
155 
156 template <class T>
159  const field_basic<T,sequential>& uh,
160  const point_basic<T>& origin,
161  const point_basic<T>& normal)
162 {
163  using namespace std;
166  typedef typename geo_basic<float_type,sequential>::node_type node_type;
167  typedef point_basic<size_type> ilat;
168  ostream& os = std::cout;
169  bool verbose = iorheo::getverbose(os);
170  bool clean = iorheo::getclean(os);
171  bool execute = iorheo::getexecute(os);
172  string basename = iorheo::getbasename(os);
173  string tmp = get_tmpdir() + "/";
174  if (!clean) tmp = "";
175  string filelist;
176  string vtk_name = tmp+basename + ".vtk";
177  filelist = filelist + " " + vtk_name;
178  ofstream vtk_os (vtk_name.c_str());
179  odiststream vtk (vtk_os);
180  if (verbose) clog << "! file \"" << vtk_name << "\" created.\n";
181  field_put_vtk (vtk, uh);
182  vtk.close();
183  string py_name = tmp+basename + "-cut.py";
184  string vtk_cut_name = tmp+basename + "-cut.vtk";
185  filelist = filelist + " " + py_name + " " + vtk_cut_name;
186  ofstream py (py_name.c_str());
187  if (verbose) clog << "! file \"" << py_name << "\" created.\n";
188  py << "#!/usr/bin/env mayavi2" << endl
189  << "# This is a mayavi script for the visualization of " << basename << ".vtk" << endl
190  << "# automatically generated by rheolef." << endl
191  << endl
192  ;
193  py << "from sys import version_info" << endl
194  << "# mayavi2.py file has changed from directory since python 2.7:" << endl
195  << "if version_info[0] >= 2 and version_info[1] >= 7 and version_info[2] >= 2:" << endl
196  << " from mayavi.scripts import mayavi2" << endl
197  << "else:" << endl
198  << " from enthought.mayavi.scripts import mayavi2" << endl
199  << "mayavi2.standalone(globals()) # start mayvi2 graphic interface" << endl
200  << "from mayavi2_rheolef import * # load rheolef specific functions" << endl
201  << endl
202  << setprecision(numeric_limits<Float>::digits10)
203  << "input = \"" << vtk_name << "\"" << endl
204  << "output = \"" << vtk_cut_name << "\"" << endl
205  << "origin = " << python(origin) << endl
206  << "normal = " << python(normal) << endl
207  << "tol = " << sqrt(numeric_limits<Float>::epsilon()) << endl
208  << "plane_cutter_on_file(input, output, origin, normal, tol)" << endl
209  ;
210  py.close();
211  int status = 0;
212  string command;
213  command = "LANG=C PYTHONPATH=" + string(_RHEOLEF_PKGDATADIR) + " python " + py_name;
214  if (verbose) clog << "! " << command << endl;
215  status = system (command.c_str());
216  ifstream vtk_polydata (vtk_cut_name.c_str());
217  check_macro (vtk_polydata, "field: vtk polydata file \"" << vtk_cut_name << "\" not found.");
218  if (verbose) clog << "! load `" << vtk_cut_name << "'." << endl;
219  idiststream ips_vtk_polydata (vtk_polydata);
220  geo_basic<T,sequential> gamma_cut;
221  geo_get_vtk (ips_vtk_polydata, gamma_cut);
222  gamma_cut.set_name (basename + "-cut");
223  check_macro (gamma_cut.n_node() > 0, "empty mesh & plane intersection: HINT check normal and origin");
224  if (clean) {
225  command = "/bin/rm -f " + filelist;
226  if (verbose) clog << "! " << command << endl;
227  status = system (command.c_str());
228  }
229  // project geometry from 2d/3d on plane in 1d/2d
230  bool use_projection = true;
231  if (use_projection) {
232  switch (uh.get_geo().dimension()) {
233  case 2: {
234  gamma_cut.set_dimension (1);
235  point_basic<T> tangent = point_basic<T>(normal[1], -normal[0]);
236  disarray<node_type,sequential> node = gamma_cut.get_nodes();
237  for (size_type i = 0, n = node.size(); i < n; i++) {
238  node[i] = point_basic<T> (dot(node[i], tangent), 0);
239  }
240  gamma_cut.set_nodes (node);
241  break;
242  }
243  case 3: {
244  gamma_cut.set_dimension (2);
245  point_basic<T> t1 = point_basic<T>(1, 0, 0);
246  if (norm(vect(t1, normal)) == Float(0)) t1 = point_basic<T>(0, 1, 0);
247  t1 = t1 - dot(t1,normal)*normal;
248  t1 = t1 / norm(t1);
249  point_basic<T> t2 = vect(normal,t1);
250  disarray<node_type,sequential> node = gamma_cut.get_nodes();
251  for (size_type i = 0, n = node.size(); i < n; i++) {
252  node[i] = point_basic<T> (dot(node[i] - origin, t1), dot(node[i] - origin, t2), 0);
253  }
254  gamma_cut.set_nodes (node);
255  break;
256  }
257  default: error_macro ("unexpected dimension="<< uh.get_geo().dimension());
258  }
259  }
260  // POINT_DATA <n_node>
261  string data_type;
262  while (true) {
263  vtk_polydata >> data_type;
264  if ((data_type == "POINT_DATA") ||
265  (data_type == "CELL_DATA" )) {
266  break;
267  }
268  }
269  scatch (vtk_polydata, "default");
270  string approx = (data_type == "POINT_DATA") ? "P1" : "P0";
271  space_basic<T,sequential> V_cut (gamma_cut, approx);
272  field_basic<T,sequential> u_cut (V_cut);
273  u_cut.set_u().get_values (ips_vtk_polydata);
274  return u_cut;
275 }
278 
279 }// namespace rheolef
point_basic< T > vect(const point_basic< T > &v, const point_basic< T > &w)
Definition: point.h:138
idiststream & geo_get_vtk(idiststream &, geo_basic< T, sequential > &)
#define error_macro(message)
Definition: compiler.h:100
details::field_expr_v2_nonlinear_terminal_function< details::normal_pseudo_function< Float > > normal()
const basis_basic< T > & get_basis() const
Definition: numbering.h:77
field - piecewise polynomial finite element field
std::string get_tmpdir()
Definition: rheostream.cc:35
void set_name(std::string name)
idiststream, odiststream - distributed interface for large data streams
Definition: diststream.h:68
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:56
point - vertex of a mesh
Definition: point.h:22
float_traits< value_type >::type float_type
STL namespace.
rheolef::details::is_vec dot
vec< T, M > & set_u()
Definition: field.h:293
template field_basic< Float, sequential > plane_cut(const field_basic< Float, sequential > &, const point_basic< Float > &, const point_basic< Float > &)
irheostream, orheostream - large data streams
Definition: compiler.h:7
verbose clean transpose logscale grid shrink ball stereo iso volume velocity elevation bezieradapt label color color
const geo_type & get_geo() const
Definition: field.h:284
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
template odiststream & visu_vtk_mayavi< Float >(odiststream &, const field_basic< Float, sequential > &)
#define _RHEOLEF_PKGDATADIR
Definition: config.h:392
bool scatch(std::istream &in, const std::string &ch, bool full_match=true)
odiststream & visu_vtk_mayavi(odiststream &, const field_basic< T, sequential > &)
odiststream & field_put_vtk(odiststream &, const field_basic< T, sequential > &, std::string, bool)
size_type n_node() const
Definition: geo.h:972
const disarray< node_type, sequential > & get_nodes() const
Definition: geo.h:978
const space_type & get_space() const
Definition: field.h:283
Float epsilon
void set_nodes(const disarray< node_type, sequential > &x)
static std::string python(const point_basic< T > &x, size_t d=3)
std::ostream & os()
Definition: diststream.h:166
verbose clean transpose logscale grid shrink ball stereo iso volume velocity elevation bezieradapt label color color format format format format format format format format format format format format format format format format vtk