rheolef  6.5
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>
25 odiststream&
27 {
28  using namespace std;
29  typedef typename field_basic<T,sequential>::float_type float_type;
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 = iorheo::getiso(os);
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>
157 field_basic<T,sequential>
159  const field_basic<T,sequential>& uh,
160  const point_basic<T>& origin,
161  const point_basic<T>& normal)
162 {
163  using namespace std;
164  typedef typename field_basic<T,sequential>::float_type float_type;
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 #ifdef TO_CLEAN
176  point_basic<float_type> origin = iofem::getorigin(os);
177  point_basic<float_type> normal = iofem::getnormal(os);
178 #endif // TO_CLEAN
179  string filelist;
180  string vtk_name = tmp+basename + ".vtk";
181  filelist = filelist + " " + vtk_name;
182  ofstream vtk_os (vtk_name.c_str());
183  odiststream vtk (vtk_os);
184  if (verbose) clog << "! file \"" << vtk_name << "\" created.\n";
185  field_put_vtk (vtk, uh);
186  vtk.close();
187  string py_name = tmp+basename + "-cut.py";
188  string vtk_cut_name = tmp+basename + "-cut.vtk";
189  filelist = filelist + " " + py_name + " " + vtk_cut_name;
190  ofstream py (py_name.c_str());
191  if (verbose) clog << "! file \"" << py_name << "\" created.\n";
192  py << "#!/usr/bin/env mayavi2" << endl
193  << "# This is a mayavi script for the visualization of " << basename << ".vtk" << endl
194  << "# automatically generated by rheolef." << endl
195  << endl
196  ;
197  py << "from sys import version_info" << endl
198  << "# mayavi2.py file has changed from directory since python 2.7:" << endl
199  << "if version_info[0] >= 2 and version_info[1] >= 7 and version_info[2] >= 2:" << endl
200  << " from mayavi.scripts import mayavi2" << endl
201  << "else:" << endl
202  << " from enthought.mayavi.scripts import mayavi2" << endl
203  << "mayavi2.standalone(globals()) # start mayvi2 graphic interface" << endl
204  << "from mayavi2_rheolef import * # load rheolef specific functions" << endl
205  << endl
206  << setprecision(numeric_limits<Float>::digits10)
207  << "input = \"" << vtk_name << "\"" << endl
208  << "output = \"" << vtk_cut_name << "\"" << endl
209  << "origin = " << python(origin) << endl
210  << "normal = " << python(normal) << endl
211  << "tol = " << sqrt(numeric_limits<Float>::epsilon()) << endl
212  << "plane_cutter_on_file(input, output, origin, normal, tol)" << endl
213  ;
214  py.close();
215  int status = 0;
216  string command;
217  command = "LANG=C PYTHONPATH=" + string(_RHEOLEF_PKGDATADIR) + " python " + py_name;
218  if (verbose) clog << "! " << command << endl;
219  status = system (command.c_str());
220  ifstream vtk_polydata (vtk_cut_name.c_str());
221  check_macro (vtk_polydata, "field: vtk polydata file \"" << vtk_cut_name << "\" not found.");
222  if (verbose) clog << "! load `" << vtk_cut_name << "'." << endl;
223  idiststream ips_vtk_polydata (vtk_polydata);
224  geo_basic<T,sequential> gamma_cut;
225  geo_get_vtk (ips_vtk_polydata, gamma_cut);
226  gamma_cut.set_name (basename + "-cut");
227  check_macro (gamma_cut.n_node() > 0, "empty mesh & plane intersection: HINT check normal and origin");
228  if (clean) {
229  command = "/bin/rm -f " + filelist;
230  if (verbose) clog << "! " << command << endl;
231  status = system (command.c_str());
232  }
233  // project geometry from 2d/3d on plane in 1d/2d
234  bool use_projection = true;
235  if (use_projection) {
236  switch (uh.get_geo().dimension()) {
237  case 2: {
238  gamma_cut.set_dimension (1);
239  point_basic<T> tangent = point_basic<T>(normal[1], -normal[0]);
240  array<node_type,sequential> node = gamma_cut.get_nodes();
241  for (size_type i = 0, n = node.size(); i < n; i++) {
242  node[i] = point_basic<T> (dot(node[i], tangent), 0);
243  }
244  gamma_cut.set_nodes (node);
245  break;
246  }
247  case 3: {
248  gamma_cut.set_dimension (2);
249  point_basic<T> t1 = point_basic<T>(1, 0, 0);
250  if (norm(vect(t1, normal)) == Float(0)) t1 = point_basic<T>(0, 1, 0);
251  t1 = t1 - dot(t1,normal)*normal;
252  t1 = t1 / norm(t1);
253  point_basic<T> t2 = vect(normal,t1);
254  array<node_type,sequential> node = gamma_cut.get_nodes();
255  for (size_type i = 0, n = node.size(); i < n; i++) {
256  node[i] = point_basic<T> (dot(node[i] - origin, t1), dot(node[i] - origin, t2), 0);
257  }
258  gamma_cut.set_nodes (node);
259  break;
260  }
261  default: error_macro ("unexpected dimension="<< uh.get_geo().dimension());
262  }
263  }
264  // POINT_DATA <n_node>
265  scatch (vtk_polydata, "POINT_DATA");
266  scatch (vtk_polydata, "SCALARS");
267  scatch (vtk_polydata, "default");
268  space_basic<T,sequential> V_cut (gamma_cut, "P1");
269  field_basic<T,sequential> u_cut (V_cut);
270  u_cut.set_u().get_values (ips_vtk_polydata);
271  return u_cut;
272 }
273 template odiststream& visu_vtk_mayavi<Float> (odiststream&, const field_basic<Float,sequential>&);
274 template field_basic<Float,sequential> plane_cut (const field_basic<Float,sequential>&, const point_basic<Float>&, const point_basic<Float>&);
275 
276 }// namespace rheolef
277