rheolef  6.6
geo_seq_put_vtk.cc
Go to the documentation of this file.
1 #include "rheolef/geo.h"
2 #include "rheolef/piola.h"
3 #include "rheolef/rheostream.h"
4 #include "rheolef/iorheo.h"
5 using namespace std;
6 namespace rheolef {
7 
8 template <class T>
9 static
10 void
11 put_edge (ostream& vtk, const geo_element& K, const numbering<T,sequential>& my_numb, const geo_basic<T,sequential>& omega)
12 {
14  typedef point_basic<size_type> ilat;
15  std::vector<size_type> inod;
16  my_numb.dis_idof (omega.sizes(), K, inod);
17  size_type my_order = my_numb.degree();
18  for (size_type i = 0; i < my_order; i++) {
19  size_type loc_inod0 = reference_element_e::ilat2loc_inod (my_order, ilat(i));
20  size_type loc_inod1 = reference_element_e::ilat2loc_inod (my_order, ilat(i+1));
21  vtk << "2\t" << inod[loc_inod0] << " " << inod[loc_inod1] << endl;
22  }
23 }
24 template <class T>
25 static
26 void
27 put_triangle (ostream& vtk, const geo_element& K, const numbering<T,sequential>& my_numb, const geo_basic<T,sequential>& omega)
28 {
30  typedef point_basic<size_type> ilat;
31  std::vector<size_type> inod;
32  my_numb.dis_idof (omega.sizes(), K, inod);
33  size_type my_order = my_numb.degree();
34  for (size_type i = 0; i < my_order; i++) {
35  for (size_type j = 0; i+j < my_order; j++) {
36  size_type loc_inod00 = reference_element_t::ilat2loc_inod (my_order, ilat(i, j));
37  size_type loc_inod10 = reference_element_t::ilat2loc_inod (my_order, ilat(i+1, j));
38  size_type loc_inod01 = reference_element_t::ilat2loc_inod (my_order, ilat(i, j+1));
39  vtk << "3\t" << inod[loc_inod00] << " "
40  << inod[loc_inod10] << " "
41  << inod[loc_inod01] << endl;
42  if (i+j+1 >= my_order) continue;
43  size_type loc_inod11 = reference_element_t::ilat2loc_inod (my_order, ilat(i+1, j+1));
44  vtk << "3\t" << inod[loc_inod10] << " "
45  << inod[loc_inod11] << " "
46  << inod[loc_inod01] << endl;
47  }
48  }
49 }
50 template <class T>
51 static
52 void
53 put_quadrangle (ostream& vtk, const geo_element& K, const numbering<T,sequential>& my_numb, const geo_basic<T,sequential>& omega)
54 {
56  typedef point_basic<size_type> ilat;
57  std::vector<size_type> inod;
58  my_numb.dis_idof (omega.sizes(), K, inod);
59  size_type my_order = my_numb.degree();
60  for (size_type i = 0; i < my_order; i++) {
61  for (size_type j = 0; j < my_order; j++) {
62  size_type loc_inod00 = reference_element_q::ilat2loc_inod (my_order, ilat(i, j));
63  size_type loc_inod10 = reference_element_q::ilat2loc_inod (my_order, ilat(i+1, j));
64  size_type loc_inod11 = reference_element_q::ilat2loc_inod (my_order, ilat(i+1, j+1));
65  size_type loc_inod01 = reference_element_q::ilat2loc_inod (my_order, ilat(i, j+1));
66  vtk << "4\t" << inod[loc_inod00] << " "
67  << inod[loc_inod10] << " "
68  << inod[loc_inod11] << " "
69  << inod[loc_inod01] << endl;
70  }
71  }
72 }
73 template <class T>
74 static
75 void
76 put_tetrahedron (ostream& vtk, const geo_element& K, const numbering<T,sequential>& my_numb, const geo_basic<T,sequential>& omega)
77 {
79  typedef point_basic<size_type> ilat;
80  std::vector<size_type> inod;
81  my_numb.dis_idof (omega.sizes(), K, inod);
82  size_type my_order = my_numb.degree();
83  for (size_type i = 0; i < my_order; i++) {
84  for (size_type j = 0; i+j < my_order; j++) {
85  for (size_type k = 0; i+j+k < my_order; k++) {
86  size_type loc_inod000 = reference_element_T::ilat2loc_inod (my_order, ilat(i, j, k));
87  size_type loc_inod100 = reference_element_T::ilat2loc_inod (my_order, ilat(i+1, j, k));
88  size_type loc_inod010 = reference_element_T::ilat2loc_inod (my_order, ilat(i, j+1, k));
89  size_type loc_inod001 = reference_element_T::ilat2loc_inod (my_order, ilat(i, j, k+1));
90  vtk << "4\t" << inod[loc_inod000] << " "
91  << inod[loc_inod100] << " "
92  << inod[loc_inod010] << " "
93  << inod[loc_inod001] << endl;
94  if (i+j+k+2 > my_order) continue;
95  size_type loc_inod110 = reference_element_T::ilat2loc_inod (my_order, ilat(i+1, j+1, k));
96  size_type loc_inod101 = reference_element_T::ilat2loc_inod (my_order, ilat(i+1, j, k+1));
97  size_type loc_inod011 = reference_element_T::ilat2loc_inod (my_order, ilat(i, j+1, k+1));
98  vtk << "4\t" << inod[loc_inod100] << " " // face in x0 & x2 direction
99  << inod[loc_inod101] << " "
100  << inod[loc_inod010] << " "
101  << inod[loc_inod001] << endl
102  << "4\t" << inod[loc_inod010] << " " // face in x1 & x2 direction
103  << inod[loc_inod011] << " "
104  << inod[loc_inod001] << " "
105  << inod[loc_inod101] << endl
106  << "4\t" << inod[loc_inod100] << " "
107  << inod[loc_inod101] << " "
108  << inod[loc_inod110] << " "
109  << inod[loc_inod010] << endl
110  << "4\t" << inod[loc_inod010] << " "
111  << inod[loc_inod110] << " "
112  << inod[loc_inod011] << " "
113  << inod[loc_inod101] << endl;
114  if (i+j+k+3 > my_order) continue;
115  size_type loc_inod111 = reference_element_T::ilat2loc_inod (my_order, ilat(i+1, j+1, k+1));
116  vtk << "4\t" << inod[loc_inod111] << " " // face in x0 & x2 direction
117  << inod[loc_inod101] << " "
118  << inod[loc_inod011] << " "
119  << inod[loc_inod110] << endl;
120  }
121  }
122  }
123 }
124 static
125 void
126 raw_put_prism (ostream& vtk,
127  size_t i000, size_t i100, size_t i010,
128  size_t i001, size_t i101, size_t i011)
129 {
130  vtk << "6\t" << i000 << " "
131  << i010 << " "
132  << i100 << " "
133  << i001 << " "
134  << i011 << " "
135  << i101 << endl;
136 }
137 template <class T>
138 static
139 void
140 put_prism (ostream& vtk, const geo_element& K, const numbering<T,sequential>& my_numb, const geo_basic<T,sequential>& omega, const array<point_basic<Float>,sequential>& my_node)
141 {
143  typedef point_basic<size_type> ilat;
144  std::vector<size_type> inod;
145  my_numb.dis_idof (omega.sizes(), K, inod);
146  size_type my_order = my_numb.degree();
147  for (size_type k = 0; k < my_order; k++) {
148  for (size_type j = 0; j < my_order; j++) {
149  for (size_type i = 0; i+j < my_order; i++) {
150  size_type loc_inod000 = reference_element_P::ilat2loc_inod (my_order, ilat(i, j, k));
151  size_type loc_inod100 = reference_element_P::ilat2loc_inod (my_order, ilat(i+1, j, k));
152  size_type loc_inod010 = reference_element_P::ilat2loc_inod (my_order, ilat(i, j+1, k));
153  size_type loc_inod001 = reference_element_P::ilat2loc_inod (my_order, ilat(i, j, k+1));
154  size_type loc_inod101 = reference_element_P::ilat2loc_inod (my_order, ilat(i+1, j, k+1));
155  size_type loc_inod011 = reference_element_P::ilat2loc_inod (my_order, ilat(i, j+1, k+1));
156  raw_put_prism (vtk,
157  inod[loc_inod000],
158  inod[loc_inod100],
159  inod[loc_inod010],
160  inod[loc_inod001],
161  inod[loc_inod101],
162  inod[loc_inod011]);
163  if (i+j+1 >= my_order) continue;
164  size_type loc_inod110 = reference_element_P::ilat2loc_inod (my_order, ilat(i+1, j+1, k));
165  size_type loc_inod111 = reference_element_P::ilat2loc_inod (my_order, ilat(i+1, j+1, k+1));
166  raw_put_prism (vtk,
167  inod[loc_inod100],
168  inod[loc_inod110],
169  inod[loc_inod010],
170  inod[loc_inod101],
171  inod[loc_inod111],
172  inod[loc_inod011]);
173  }
174  }
175  }
176 }
177 template <class T>
178 static
179 void
180 put_hexahedron (ostream& vtk, const geo_element& K, const numbering<T,sequential>& my_numb, const geo_basic<T,sequential>& omega)
181 {
183  typedef point_basic<size_type> ilat;
184  std::vector<size_type> inod;
185  my_numb.dis_idof (omega.sizes(), K, inod);
186  size_type my_order = my_numb.degree();
187  for (size_type i = 0; i < my_order; i++) {
188  for (size_type j = 0; j < my_order; j++) {
189  for (size_type k = 0; k < my_order; k++) {
190  size_type loc_inod000 = reference_element_H::ilat2loc_inod (my_order, ilat(i, j, k));
191  size_type loc_inod100 = reference_element_H::ilat2loc_inod (my_order, ilat(i+1, j, k));
192  size_type loc_inod110 = reference_element_H::ilat2loc_inod (my_order, ilat(i+1, j+1, k));
193  size_type loc_inod010 = reference_element_H::ilat2loc_inod (my_order, ilat(i, j+1, k));
194  size_type loc_inod001 = reference_element_H::ilat2loc_inod (my_order, ilat(i, j, k+1));
195  size_type loc_inod101 = reference_element_H::ilat2loc_inod (my_order, ilat(i+1, j, k+1));
196  size_type loc_inod011 = reference_element_H::ilat2loc_inod (my_order, ilat(i, j+1, k+1));
197  size_type loc_inod111 = reference_element_H::ilat2loc_inod (my_order, ilat(i+1, j+1, k+1));
198  vtk << "8\t" << inod[loc_inod000] << " "
199  << inod[loc_inod100] << " "
200  << inod[loc_inod110] << " "
201  << inod[loc_inod010] << " "
202  << inod[loc_inod001] << " "
203  << inod[loc_inod101] << " "
204  << inod[loc_inod111] << " "
205  << inod[loc_inod011] << endl;
206  }
207  }
208  }
209 }
210 template <class T>
211 static
212 void
213 put (ostream& vtk, const geo_element& K, const numbering<T,sequential>& my_numb, const geo_basic<T,sequential>& omega, const array<point_basic<Float>,sequential>& my_node)
214 {
215  switch (K.variant()) {
216  case reference_element::p: vtk << "1\t" << K[0] << endl; break;
217  case reference_element::e: put_edge (vtk, K, my_numb, omega); break;
218  case reference_element::t: put_triangle (vtk, K, my_numb, omega); break;
219  case reference_element::q: put_quadrangle (vtk, K, my_numb, omega); break;
220  case reference_element::T: put_tetrahedron (vtk, K, my_numb, omega); break;
221  case reference_element::P: put_prism (vtk, K, my_numb, omega, my_node); break;
222  case reference_element::H: put_hexahedron (vtk, K, my_numb, omega); break;
223  default: error_macro ("unsupported element variant `" << K.name() <<"'");
224  }
225 }
226 
227 #ifndef VTK_EMPTY_CELL
228 // from vtk-3.1.2/common/vtkCellType.h:
229 #define VTK_EMPTY_CELL 0
230 #define VTK_VERTEX 1
231 #define VTK_POLY_VERTEX 2
232 #define VTK_LINE 3
233 #define VTK_POLY_LINE 4
234 #define VTK_TRIANGLE 5
235 #define VTK_TRIANGLE_STRIP 6
236 #define VTK_POLYGON 7
237 #define VTK_PIXEL 8
238 #define VTK_QUAD 9
239 #define VTK_TETRA 10
240 #define VTK_VOXEL 11
241 #define VTK_HEXAHEDRON 12
242 #define VTK_WEDGE 13
243 #define VTK_PYRAMID 14
244 
245 #define VTK_PARAMETRIC_CURVE 51
246 #define VTK_PARAMETRIC_SURFACE 52
247 #endif // VTK
248 
249 // TODO: P3 geo with P2 field: draw on the P3 lattice and re-interpolate the field !
250 template <class T>
251 odiststream&
252 geo_put_vtk (odiststream& ops, const geo_basic<T,sequential>& omega, const numbering<T,sequential>& my_numb, const array<point_basic<T>,sequential>& my_node, bool append_data)
253 {
255  size_type my_order = my_numb.degree();
256  ostream& vtk = ops.os();
257  check_macro (my_order >= omega.order(), "order="<<omega.order()<<" > field degree="<<my_order);
258  vtk << "# vtk DataFile Version 1.0" << endl
259  << "Unstructured Grid" << endl
260  << "ASCII" << endl
261  << "DATASET UNSTRUCTURED_GRID" << endl;
262  vtk << "POINTS " << my_node.size() << " float" << endl;
263  for (size_type inod = 0, nnod = my_node.size(); inod < nnod; inod++) {
264  vtk << my_node[inod] << endl;
265  }
266  size_type map_dim = omega.map_dimension();
267  boost::array<size_type,reference_element::max_variant> size_by_variant;
268  size_by_variant.assign (0);
269  for (size_type ie = 0, ne = omega.size(); ie < ne; ie++) {
270  const geo_element& K = omega.get_geo_element (map_dim, ie);
271  size_by_variant [K.variant()]++;
272  }
273  size_type ncell = 0;
274  size_type ndata = 0;
275  boost::array<size_type,reference_element::max_variant> loc_ncell;
276  boost::array<size_type,reference_element::max_variant> loc_ndata;
277  for (size_type variant = reference_element::first_variant_by_dimension(map_dim);
278  variant < reference_element:: last_variant_by_dimension(map_dim); variant++) {
279  size_type d = reference_element::dimension(variant);
280  size_type n = reference_element::n_vertex (variant);
281  loc_ncell [variant] = pow(my_order,d);
282  loc_ndata [variant] = (n+1)*loc_ncell [variant];
283  ncell += loc_ncell[variant]*size_by_variant [variant];
284  ndata += loc_ndata[variant]*size_by_variant [variant];
285  }
286  vtk << "CELLS " << ncell << " " << ndata << endl;
287  for (size_type ie = 0, ne = omega.size(); ie < ne; ie++) {
288  const geo_element& K = omega.get_geo_element (map_dim, ie);
289  put (vtk, K, my_numb, omega, my_node);
290  }
291  boost::array<size_type,reference_element::max_variant> cell_type;
292  cell_type [reference_element::p] = VTK_VERTEX;
293  cell_type [reference_element::e] = VTK_LINE;
294  cell_type [reference_element::t] = VTK_TRIANGLE;
295  cell_type [reference_element::q] = VTK_QUAD;
296  cell_type [reference_element::T] = VTK_TETRA;
297  cell_type [reference_element::P] = VTK_WEDGE;
298  cell_type [reference_element::H] = VTK_HEXAHEDRON;
299  vtk << "CELL_TYPES " << ncell << endl;
300  for (size_type ie = 0, ne = omega.size(); ie < ne; ie++) {
301  const geo_element& K = omega.get_geo_element (map_dim, ie);
302  for (size_type k = 0; k < loc_ncell[K.variant()]; k++) {
303  vtk << cell_type [K.variant()] << endl;
304  }
305  }
306  if (! append_data) return ops;
307  std::string data_name = "mesh";
308  vtk << "POINT_DATA " << my_node.size() << endl
309  << "SCALARS " << data_name << " float" << endl
310  << "LOOKUP_TABLE default" << endl;
311  for (size_type inod = 0, nnod = my_node.size(); inod < nnod; inod++) {
312  vtk << "0" << endl;
313  }
314  vtk << endl;
315 
316  return ops;
317 }
318 template <class T>
319 odiststream&
320 geo_put_vtk (odiststream& ops, const geo_basic<T,sequential>& omega, size_t my_order, const array<point_basic<T>,sequential>& my_node, bool append_data)
321 {
322  numbering<T,sequential> my_numb ("P" + itos(my_order));
323  return geo_put_vtk (ops, omega, my_numb, my_node, append_data);
324 }
325 template <class T>
326 odiststream&
327 geo_put_vtk (odiststream& ops, const geo_basic<T,sequential>& omega, size_t my_order, const array<point_basic<T>,sequential>& my_node)
328 {
329  return geo_put_vtk (ops, omega, my_order, my_node, true);
330 }
331 template <class T>
332 odiststream&
333 geo_put_vtk (odiststream& ops, const geo_basic<T,sequential>& omega, bool append_data)
334 {
335  return geo_put_vtk (ops, omega, omega.order(), omega.get_nodes(), append_data);
336 }
337 template <class T>
338 odiststream&
340 {
341  return geo_put_vtk (ops, omega, omega.order(), omega.get_nodes(), true);
342 }
343 template odiststream& geo_put_vtk<Float> (odiststream&, const geo_basic<Float,sequential>&, const numbering<Float,sequential>&, const array<point_basic<Float>,sequential>&, bool);
344 template odiststream& geo_put_vtk<Float> (odiststream&, const geo_basic<Float,sequential>&, size_t, const array<point_basic<Float>,sequential>&, bool);
345 template odiststream& geo_put_vtk<Float> (odiststream&, const geo_basic<Float,sequential>&, size_t, const array<point_basic<Float>,sequential>&);
346 template odiststream& geo_put_vtk<Float> (odiststream&, const geo_basic<Float,sequential>&);
347 
348 }// namespace rheolef
#define error_macro(message)
Definition: compiler.h:100
#define VTK_WEDGE
#define VTK_TETRA
size_type size(size_type dim) const
Definition: geo.h:1025
idiststream, odiststream - distributed interface for large data streams
Definition: diststream.h:68
static void raw_put_prism(ostream &vtk, size_t i000, size_t i100, size_t i010, size_t i001, size_t i101, size_t i011)
#define VTK_LINE
point - vertex of a mesh
Definition: point.h:22
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)
const geo_size & sizes() const
Definition: geo.h:924
#define VTK_TRIANGLE
STL namespace.
void dis_idof(const geo_size &gs, const geo_element &K, std::vector< size_type > &dis_idof) const
Definition: numbering.h:171
const array< node_type, sequential > & get_nodes() const
Definition: geo.h:941
irheostream, orheostream - large data streams
Definition: compiler.h:7
string itos(string::size_type i)
Definition: rheostream.cc:38
variant_type variant() const
#define VTK_HEXAHEDRON
#define check_macro(ok_condition, message)
Definition: compiler.h:104
const_reference get_geo_element(size_type dim, size_type ige) const
Definition: geo.h:926
size_t d
template odiststream & geo_put_vtk< Float >(odiststream &, const geo_basic< Float, sequential > &)
size_type map_dimension() const
Definition: geo.h:912
static void put_tetrahedron(ostream &vtk, const geo_element &K, const numbering< T, sequential > &my_numb, const geo_basic< T, sequential > &omega)
#define VTK_QUAD
static void put_triangle(ostream &vtk, const geo_element &K, const numbering< T, sequential > &my_numb, const geo_basic< T, sequential > &omega)
static void put_hexahedron(ostream &vtk, const geo_element &K, const numbering< T, sequential > &my_numb, const geo_basic< T, sequential > &omega)
void put_prism(ostream &out, size_t my_order)
#define VTK_VERTEX
static void put_edge(ostream &vtk, const geo_element &K, const numbering< T, sequential > &my_numb, const geo_basic< T, sequential > &omega)
static void put_quadrangle(ostream &vtk, const geo_element &K, const numbering< T, sequential > &my_numb, const geo_basic< T, sequential > &omega)
geo_element - element of a mesh
size_type order() const
Definition: geo.h:918
static void put(ostream &vtk, const geo_element &K, const numbering< T, sequential > &my_numb, const geo_basic< T, sequential > &omega, const array< point_basic< Float >, sequential > &my_node)
odiststream & geo_put_vtk(odiststream &ops, const geo_basic< T, sequential > &omega)
reference_element_H::size_type size_type
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
size_type degree() const
Definition: numbering.h:118