rheolef  6.5
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
349