Rheolef  7.2
an efficient C++ finite element environment
geo_seq_put_gmsh.cc
Go to the documentation of this file.
1 //
22 // gmsh output
23 //
24 // author: Pierre.Saramito@imag.fr
25 //
26 // date: 26 mars 2012
27 //
28 #include "rheolef/geo.h"
29 #include "rheolef/space_numbering.h"
30 #include "rheolef/piola_util.h"
31 #include "rheolef/rheostream.h"
32 #include "rheolef/iorheo.h"
33 using namespace std;
34 namespace rheolef {
35 
36 // ----------------------------------------------------------------------------
37 // one element puts
38 // ----------------------------------------------------------------------------
39 template <class T>
40 static
41 void
42 put_edge (ostream& gmsh, const geo_element& K, const basis_basic<T>& my_numb, const geo_basic<T,sequential>& omega)
43 {
45  static size_type order2gmsh_type [11] = {0, 1, 8, 26, 27, 28, 62, 63, 64, 65, 66 };
46  size_type my_order = my_numb.degree();
47  check_macro (my_order <= 10, "gmsh output: element 'e' order > 10 not yet supported");
48  std::vector<size_type> inod;
49  space_numbering::dis_idof (my_numb, omega.sizes(), K, inod);
50  gmsh << K.dis_ie()+1 << " " << order2gmsh_type [my_order] << " 2 99 2"; // TODO: domains
51  for (size_type iloc = 0, nloc = inod.size(); iloc < nloc; iloc++) {
52  gmsh << " " << inod[iloc]+1;
53  }
54  gmsh << endl;
55 }
56 template <class T>
57 static
58 void
59 put_triangle (ostream& gmsh, const geo_element& K, const basis_basic<T>& my_numb, const geo_basic<T,sequential>& omega)
60 {
61  typedef typename geo_basic<T,sequential>::size_type size_type;
62  static size_type order2gmsh_type [11] = {0, 2, 9, 21, 23, 25, 42, 43, 44, 45, 46};
63  size_type my_order = my_numb.degree();
64  // TODO: permutations of internal nodes for order >= 4
65  check_macro (my_order <= 3, "gmsh output: element 't' order > 10 not yet supported");
66  std::vector<size_type> inod;
67  space_numbering::dis_idof (my_numb, omega.sizes(), K, inod);
68  gmsh << K.dis_ie()+1 << " " << order2gmsh_type [my_order] << " 2 99 2"; // TODO: domains
69  for (size_type iloc = 0, nloc = inod.size(); iloc < nloc; iloc++) {
70  gmsh << " " << inod[iloc]+1;
71  }
72  gmsh << endl;
73 }
74 #ifdef TODO
75 template <class T>
76 static
77 void
78 put_quadrangle (ostream& gmsh, const geo_element& K, const basis_basic<T>& my_numb, const geo_basic<T,sequential>& omega)
79 {
80  typedef typename geo_basic<T,sequential>::size_type size_type;
81  typedef point_basic<size_type> ilat;
82  std::vector<size_type> inod;
83  space_numbering::dis_idof (my_numb, omega.sizes(), K, inod);
84  size_type my_order = my_numb.degree();
85  for (size_type i = 0; i < my_order; i++) {
86  for (size_type j = 0; j < my_order; j++) {
87  size_type loc_inod00 = reference_element_q::ilat2loc_inod (my_order, ilat(i, j));
88  size_type loc_inod10 = reference_element_q::ilat2loc_inod (my_order, ilat(i+1, j));
89  size_type loc_inod11 = reference_element_q::ilat2loc_inod (my_order, ilat(i+1, j+1));
90  size_type loc_inod01 = reference_element_q::ilat2loc_inod (my_order, ilat(i, j+1));
91  gmsh << "4\t" << inod[loc_inod00] << " "
92  << inod[loc_inod10] << " "
93  << inod[loc_inod11] << " "
94  << inod[loc_inod01] << endl;
95  }
96  }
97 }
98 template <class T>
99 static
100 void
101 put_tetrahedron (ostream& gmsh, const geo_element& K, const basis_basic<T>& my_numb, const geo_basic<T,sequential>& omega)
102 {
103  typedef typename geo_basic<T,sequential>::size_type size_type;
104  typedef point_basic<size_type> ilat;
105  std::vector<size_type> inod;
106  space_numbering::dis_idof (my_numb, omega.sizes(), K, inod);
107  size_type my_order = my_numb.degree();
108  for (size_type i = 0; i < my_order; i++) {
109  for (size_type j = 0; i+j < my_order; j++) {
110  for (size_type k = 0; i+j+k < my_order; k++) {
111  size_type loc_inod000 = reference_element_T::ilat2loc_inod (my_order, ilat(i, j, k));
112  size_type loc_inod100 = reference_element_T::ilat2loc_inod (my_order, ilat(i+1, j, k));
113  size_type loc_inod010 = reference_element_T::ilat2loc_inod (my_order, ilat(i, j+1, k));
114  size_type loc_inod001 = reference_element_T::ilat2loc_inod (my_order, ilat(i, j, k+1));
115  gmsh << "4\t" << inod[loc_inod000] << " "
116  << inod[loc_inod100] << " "
117  << inod[loc_inod010] << " "
118  << inod[loc_inod001] << endl;
119  if (i+j+k+2 > my_order) continue;
120  // complete the ijk-th cube: 4 more tetras
121  size_type loc_inod110 = reference_element_T::ilat2loc_inod (my_order, ilat(i+1, j+1, k));
122  size_type loc_inod101 = reference_element_T::ilat2loc_inod (my_order, ilat(i+1, j, k+1));
123  size_type loc_inod011 = reference_element_T::ilat2loc_inod (my_order, ilat(i, j+1, k+1));
124  gmsh << "4\t" << inod[loc_inod100] << " " // face in x0 & x2 direction
125  << inod[loc_inod101] << " "
126  << inod[loc_inod010] << " "
127  << inod[loc_inod001] << endl
128  << "4\t" << inod[loc_inod010] << " " // face in x1 & x2 direction
129  << inod[loc_inod011] << " "
130  << inod[loc_inod001] << " "
131  << inod[loc_inod101] << endl
132  << "4\t" << inod[loc_inod100] << " "
133  << inod[loc_inod101] << " "
134  << inod[loc_inod110] << " "
135  << inod[loc_inod010] << endl
136  << "4\t" << inod[loc_inod010] << " "
137  << inod[loc_inod110] << " "
138  << inod[loc_inod011] << " "
139  << inod[loc_inod101] << endl;
140  // the last 6th sub-tetra that fully fills the ijk-th cube
141  if (i+j+k+3 > my_order) continue;
142  size_type loc_inod111 = reference_element_T::ilat2loc_inod (my_order, ilat(i+1, j+1, k+1));
143  gmsh << "4\t" << inod[loc_inod111] << " " // face in x0 & x2 direction
144  << inod[loc_inod101] << " "
145  << inod[loc_inod011] << " "
146  << inod[loc_inod110] << endl;
147  }
148  }
149  }
150 }
151 static
152 void
153 raw_put_prism (ostream& gmsh,
154  size_t i000, size_t i100, size_t i010,
155  size_t i001, size_t i101, size_t i011)
156 {
157  // gmsh prism has swaped x & y axis order: 00z 10z 01z replaced by 00z 01z 10z
158  gmsh << "6\t" << i000 << " "
159  << i010 << " "
160  << i100 << " "
161  << i001 << " "
162  << i011 << " "
163  << i101 << endl;
164 }
165 template <class T>
166 static
167 void
168 put_prism (ostream& gmsh, const geo_element& K, const basis_basic<T>& my_numb, const geo_basic<T,sequential>& omega, const disarray<point_basic<Float>,sequential>& my_node)
169 {
170  typedef typename geo_basic<T,sequential>::size_type size_type;
171  typedef point_basic<size_type> ilat;
172  std::vector<size_type> inod;
173  space_numbering::dis_idof (my_numb, omega.sizes(), K, inod);
174  size_type my_order = my_numb.degree();
175  for (size_type k = 0; k < my_order; k++) {
176  for (size_type j = 0; j < my_order; j++) {
177  for (size_type i = 0; i+j < my_order; i++) {
178  size_type loc_inod000 = reference_element_P::ilat2loc_inod (my_order, ilat(i, j, k));
179  size_type loc_inod100 = reference_element_P::ilat2loc_inod (my_order, ilat(i+1, j, k));
180  size_type loc_inod010 = reference_element_P::ilat2loc_inod (my_order, ilat(i, j+1, k));
181  size_type loc_inod001 = reference_element_P::ilat2loc_inod (my_order, ilat(i, j, k+1));
182  size_type loc_inod101 = reference_element_P::ilat2loc_inod (my_order, ilat(i+1, j, k+1));
183  size_type loc_inod011 = reference_element_P::ilat2loc_inod (my_order, ilat(i, j+1, k+1));
184  raw_put_prism (gmsh,
185  inod[loc_inod000],
186  inod[loc_inod100],
187  inod[loc_inod010],
188  inod[loc_inod001],
189  inod[loc_inod101],
190  inod[loc_inod011]);
191  if (i+j+1 >= my_order) continue;
192  size_type loc_inod110 = reference_element_P::ilat2loc_inod (my_order, ilat(i+1, j+1, k));
193  size_type loc_inod111 = reference_element_P::ilat2loc_inod (my_order, ilat(i+1, j+1, k+1));
194  raw_put_prism (gmsh,
195  inod[loc_inod100],
196  inod[loc_inod110],
197  inod[loc_inod010],
198  inod[loc_inod101],
199  inod[loc_inod111],
200  inod[loc_inod011]);
201  }
202  }
203  }
204 }
205 template <class T>
206 static
207 void
208 put_hexahedron (ostream& gmsh, const geo_element& K, const basis_basic<T>& my_numb, const geo_basic<T,sequential>& omega)
209 {
210  typedef typename geo_basic<T,sequential>::size_type size_type;
211  typedef point_basic<size_type> ilat;
212  std::vector<size_type> inod;
213  space_numbering::dis_idof (my_numb, omega.sizes(), K, inod);
214  size_type my_order = my_numb.degree();
215  for (size_type i = 0; i < my_order; i++) {
216  for (size_type j = 0; j < my_order; j++) {
217  for (size_type k = 0; k < my_order; k++) {
218  size_type loc_inod000 = reference_element_H::ilat2loc_inod (my_order, ilat(i, j, k));
219  size_type loc_inod100 = reference_element_H::ilat2loc_inod (my_order, ilat(i+1, j, k));
220  size_type loc_inod110 = reference_element_H::ilat2loc_inod (my_order, ilat(i+1, j+1, k));
221  size_type loc_inod010 = reference_element_H::ilat2loc_inod (my_order, ilat(i, j+1, k));
222  size_type loc_inod001 = reference_element_H::ilat2loc_inod (my_order, ilat(i, j, k+1));
223  size_type loc_inod101 = reference_element_H::ilat2loc_inod (my_order, ilat(i+1, j, k+1));
224  size_type loc_inod011 = reference_element_H::ilat2loc_inod (my_order, ilat(i, j+1, k+1));
225  size_type loc_inod111 = reference_element_H::ilat2loc_inod (my_order, ilat(i+1, j+1, k+1));
226  gmsh << "8\t" << inod[loc_inod000] << " "
227  << inod[loc_inod100] << " "
228  << inod[loc_inod110] << " "
229  << inod[loc_inod010] << " "
230  << inod[loc_inod001] << " "
231  << inod[loc_inod101] << " "
232  << inod[loc_inod111] << " "
233  << inod[loc_inod011] << endl;
234  }
235  }
236  }
237 }
238 #endif // TODO
239 template <class T>
240 static
241 void
242 put (ostream& gmsh, const geo_element& K, const basis_basic<T>& my_numb, const geo_basic<T,sequential>& omega, const disarray<point_basic<Float>,sequential>& my_node)
243 {
244  switch (K.variant()) {
245 #ifdef TODO
246  case reference_element::p: gmsh << "1\t" << K[0] << endl; break;
247 #endif // TODO
248  case reference_element::e: put_edge (gmsh, K, my_numb, omega); break;
249  case reference_element::t: put_triangle (gmsh, K, my_numb, omega); break;
250 #ifdef TODO
251  case reference_element::q: put_quadrangle (gmsh, K, my_numb, omega); break;
252  case reference_element::T: put_tetrahedron (gmsh, K, my_numb, omega); break;
253  case reference_element::P: put_prism (gmsh, K, my_numb, omega, my_node); break;
254  case reference_element::H: put_hexahedron (gmsh, K, my_numb, omega); break;
255 #endif // TODO
256  default: error_macro ("unsupported element variant `" << K.name() <<"'");
257  }
258 }
259 // ----------------------------------------------------------------------------
260 // geo puts
261 // ----------------------------------------------------------------------------
262 
263 template <class T>
264 odiststream&
265 geo_put_gmsh (odiststream& ops, const geo_basic<T,sequential>& omega, const basis_basic<T>& my_numb, const disarray<point_basic<T>,sequential>& my_node)
266 {
267  //
268  // 0) pre-requises
269  //
270  typedef typename geo_basic<T,sequential>::size_type size_type;
271  size_type my_order = my_numb.degree();
272  ostream& gmsh = ops.os();
273  //
274  // 1) put header
275  //
276  gmsh << setprecision (std::numeric_limits<T>::digits10)
277  << "$MeshFormat" << endl
278  << "2.2 0 8" << endl
279  << "$EndMeshFormat" << endl;
280  // TODO: add domains: scan by domain and add for earch element to a domain list
281  //
282  // 2) put nodes
283  //
284  gmsh << "$Nodes" << endl
285  << my_node.size() << endl;
286  for (size_type inod = 0, nnod = my_node.size(); inod < nnod; inod++) {
287  gmsh << inod+1 << " " << my_node[inod] << endl;
288  }
289  gmsh << "$EndNodes" << endl;
290  //
291  // 3) put elements
292  //
293  size_type map_dim = omega.map_dimension();
294  gmsh << "$Elements" << endl
295  << omega.size() << endl;
296  for (size_type ie = 0, ne = omega.size(); ie < ne; ie++) {
297  const geo_element& K = omega.get_geo_element (map_dim, ie);
298  put (gmsh, K, my_numb, omega, my_node);
299  }
300  gmsh << "$EndElements" << endl;
301  return ops;
302 }
303 template <class T>
304 odiststream&
305 geo_put_gmsh (odiststream& ops, const geo_basic<T,sequential>& omega)
306 {
307  basis_basic<T> my_numb ("P" + std::to_string(omega.order()));
308  return geo_put_gmsh (ops, omega, my_numb, omega.get_nodes());
309 }
310 // ----------------------------------------------------------------------------
311 // instanciation in library
312 // ----------------------------------------------------------------------------
313 template odiststream& geo_put_gmsh<Float> (odiststream&, const geo_basic<Float,sequential>&, const basis_basic<Float>&, const disarray<point_basic<Float>,sequential>&);
314 template odiststream& geo_put_gmsh<Float> (odiststream&, const geo_basic<Float,sequential>&);
315 
316 }// namespace rheolef
field::size_type size_type
Definition: branch.cc:430
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
verbose clean transpose logscale grid shrink ball stereo iso volume skipvtk deformation fastfieldload lattice reader_on_stdin color format format format format format format gmsh
This file is part of Rheolef.