rheolef  6.6
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
const basis_basic< T > & get_basis() const
Definition: numbering.h:76
point_basic< T > vect(const point_basic< T > &v, const point_basic< T > &w)
Definition: point.h:151
field_nonlinear_expr< field_expr_terminal_function< normal_pseudo_function< Float > > > normal()
Definition: operators2.h:22
idiststream & geo_get_vtk(idiststream &, geo_basic< T, sequential > &)
#define error_macro(message)
Definition: compiler.h:100
define_sequential_odiststream_macro(char) define_sequential_odiststream_macro(int) define_sequential_odiststream_macro(unsigned int) define_sequential_odiststream_macro(long int) define_sequential_odiststream_macro(long unsigned int) define_sequential_odiststream_macro(float) define_sequential_odiststream_macro(double) define_sequential_odiststream_macro(long double) define_sequential_odiststream_macro(char *const ) define_sequential_odiststream_macro(st idiststream)()
Definition: diststream.h:228
field - piecewise polynomial finite element field
Definition: field.h:225
const geo_type & get_geo() const
Definition: field.h:272
std::string get_tmpdir()
Definition: rheostream.cc:33
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:55
point - vertex of a mesh
Definition: point.h:22
STL namespace.
vec< T, M > & set_u()
Definition: field.h:281
const array< node_type, sequential > & get_nodes() const
Definition: geo.h:941
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
bool scatch(istream &in, const string &ch, bool full_match)
Definition: rheostream.cc:220
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_t d
field_nonlinear_expr< field_nonlinear_expr_uf< details::sqrt_,field_expr_terminal_field< T, M > > > sqrt(const field_basic< T, M > &x)
template odiststream & visu_vtk_mayavi< Float >(odiststream &, const field_basic< Float, sequential > &)
#define _RHEOLEF_PKGDATADIR
Definition: config.h:402
size_type n_node() const
Definition: geo.h:935
odiststream & visu_vtk_mayavi(odiststream &, const field_basic< T, sequential > &)
odiststream & field_put_vtk(odiststream &, const field_basic< T, sequential > &, std::string, bool)
void set_nodes(const array< node_type, sequential > &x)
Float epsilon
T dot(const vec< T, M > &x, const int &y)
const space_type & get_space() const
Definition: field.h:271
reference_element_H::size_type size_type
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