rheolef  6.5
geo_seq_put_gmsh.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& gmsh, const geo_element& K, const numbering<T,sequential>& my_numb, const geo_basic<T,sequential>& omega)
12 {
14  static size_type order2gmsh_type [11] = {0, 1, 8, 26, 27, 28, 62, 63, 64, 65, 66 };
15  size_type my_order = my_numb.degree();
16  check_macro (my_order <= 10, "gmsh output: element 'e' order > 10 not yet supported");
17  std::vector<size_type> inod;
18  my_numb.dis_idof (omega.sizes(), K, inod);
19  gmsh << K.dis_ie()+1 << " " << order2gmsh_type [my_order] << " 2 99 2"; // TODO: domains
20  for (size_type iloc = 0, nloc = inod.size(); iloc < nloc; iloc++) {
21  gmsh << " " << inod[iloc]+1;
22  }
23  gmsh << endl;
24 }
25 template <class T>
26 static
27 void
28 put_triangle (ostream& gmsh, const geo_element& K, const numbering<T,sequential>& my_numb, const geo_basic<T,sequential>& omega)
29 {
30  typedef typename geo_basic<T,sequential>::size_type size_type;
31  static size_type order2gmsh_type [11] = {0, 2, 9, 21, 23, 25, 42, 43, 44, 45, 46};
32  size_type my_order = my_numb.degree();
33  // TODO: permutations of internal nodes for order >= 4
34  check_macro (my_order <= 3, "gmsh output: element 't' order > 10 not yet supported");
35  std::vector<size_type> inod;
36  my_numb.dis_idof (omega.sizes(), K, inod);
37  gmsh << K.dis_ie()+1 << " " << order2gmsh_type [my_order] << " 2 99 2"; // TODO: domains
38  for (size_type iloc = 0, nloc = inod.size(); iloc < nloc; iloc++) {
39  gmsh << " " << inod[iloc]+1;
40  }
41  gmsh << endl;
42 }
43 #ifdef TODO
44 template <class T>
45 static
46 void
47 put_quadrangle (ostream& gmsh, const geo_element& K, const numbering<T,sequential>& my_numb, const geo_basic<T,sequential>& omega)
48 {
49  typedef typename geo_basic<T,sequential>::size_type size_type;
50  typedef point_basic<size_type> ilat;
51  std::vector<size_type> inod;
52  my_numb.dis_idof (omega.sizes(), K, inod);
53  size_type my_order = my_numb.degree();
54  for (size_type i = 0; i < my_order; i++) {
55  for (size_type j = 0; j < my_order; j++) {
56  size_type loc_inod00 = reference_element_q::ilat2loc_inod (my_order, ilat(i, j));
57  size_type loc_inod10 = reference_element_q::ilat2loc_inod (my_order, ilat(i+1, j));
58  size_type loc_inod11 = reference_element_q::ilat2loc_inod (my_order, ilat(i+1, j+1));
59  size_type loc_inod01 = reference_element_q::ilat2loc_inod (my_order, ilat(i, j+1));
60  gmsh << "4\t" << inod[loc_inod00] << " "
61  << inod[loc_inod10] << " "
62  << inod[loc_inod11] << " "
63  << inod[loc_inod01] << endl;
64  }
65  }
66 }
67 template <class T>
68 static
69 void
70 put_tetrahedron (ostream& gmsh, const geo_element& K, const numbering<T,sequential>& my_numb, const geo_basic<T,sequential>& omega)
71 {
72  typedef typename geo_basic<T,sequential>::size_type size_type;
73  typedef point_basic<size_type> ilat;
74  std::vector<size_type> inod;
75  my_numb.dis_idof (omega.sizes(), K, inod);
76  size_type my_order = my_numb.degree();
77  for (size_type i = 0; i < my_order; i++) {
78  for (size_type j = 0; i+j < my_order; j++) {
79  for (size_type k = 0; i+j+k < my_order; k++) {
80  size_type loc_inod000 = reference_element_T::ilat2loc_inod (my_order, ilat(i, j, k));
81  size_type loc_inod100 = reference_element_T::ilat2loc_inod (my_order, ilat(i+1, j, k));
82  size_type loc_inod010 = reference_element_T::ilat2loc_inod (my_order, ilat(i, j+1, k));
83  size_type loc_inod001 = reference_element_T::ilat2loc_inod (my_order, ilat(i, j, k+1));
84  gmsh << "4\t" << inod[loc_inod000] << " "
85  << inod[loc_inod100] << " "
86  << inod[loc_inod010] << " "
87  << inod[loc_inod001] << endl;
88  if (i+j+k+2 > my_order) continue;
89  size_type loc_inod110 = reference_element_T::ilat2loc_inod (my_order, ilat(i+1, j+1, k));
90  size_type loc_inod101 = reference_element_T::ilat2loc_inod (my_order, ilat(i+1, j, k+1));
91  size_type loc_inod011 = reference_element_T::ilat2loc_inod (my_order, ilat(i, j+1, k+1));
92  gmsh << "4\t" << inod[loc_inod100] << " " // face in x0 & x2 direction
93  << inod[loc_inod101] << " "
94  << inod[loc_inod010] << " "
95  << inod[loc_inod001] << endl
96  << "4\t" << inod[loc_inod010] << " " // face in x1 & x2 direction
97  << inod[loc_inod011] << " "
98  << inod[loc_inod001] << " "
99  << inod[loc_inod101] << endl
100  << "4\t" << inod[loc_inod100] << " "
101  << inod[loc_inod101] << " "
102  << inod[loc_inod110] << " "
103  << inod[loc_inod010] << endl
104  << "4\t" << inod[loc_inod010] << " "
105  << inod[loc_inod110] << " "
106  << inod[loc_inod011] << " "
107  << inod[loc_inod101] << endl;
108  if (i+j+k+3 > my_order) continue;
109  size_type loc_inod111 = reference_element_T::ilat2loc_inod (my_order, ilat(i+1, j+1, k+1));
110  gmsh << "4\t" << inod[loc_inod111] << " " // face in x0 & x2 direction
111  << inod[loc_inod101] << " "
112  << inod[loc_inod011] << " "
113  << inod[loc_inod110] << endl;
114  }
115  }
116  }
117 }
118 static
119 void
120 raw_put_prism (ostream& gmsh,
121  size_t i000, size_t i100, size_t i010,
122  size_t i001, size_t i101, size_t i011)
123 {
124  gmsh << "6\t" << i000 << " "
125  << i010 << " "
126  << i100 << " "
127  << i001 << " "
128  << i011 << " "
129  << i101 << endl;
130 }
131 template <class T>
132 static
133 void
134 put_prism (ostream& gmsh, const geo_element& K, const numbering<T,sequential>& my_numb, const geo_basic<T,sequential>& omega, const array<point_basic<Float>,sequential>& my_node)
135 {
136  typedef typename geo_basic<T,sequential>::size_type size_type;
137  typedef point_basic<size_type> ilat;
138  std::vector<size_type> inod;
139  my_numb.dis_idof (omega.sizes(), K, inod);
140  size_type my_order = my_numb.degree();
141  for (size_type k = 0; k < my_order; k++) {
142  for (size_type j = 0; j < my_order; j++) {
143  for (size_type i = 0; i+j < my_order; i++) {
144  size_type loc_inod000 = reference_element_P::ilat2loc_inod (my_order, ilat(i, j, k));
145  size_type loc_inod100 = reference_element_P::ilat2loc_inod (my_order, ilat(i+1, j, k));
146  size_type loc_inod010 = reference_element_P::ilat2loc_inod (my_order, ilat(i, j+1, k));
147  size_type loc_inod001 = reference_element_P::ilat2loc_inod (my_order, ilat(i, j, k+1));
148  size_type loc_inod101 = reference_element_P::ilat2loc_inod (my_order, ilat(i+1, j, k+1));
149  size_type loc_inod011 = reference_element_P::ilat2loc_inod (my_order, ilat(i, j+1, k+1));
150  raw_put_prism (gmsh,
151  inod[loc_inod000],
152  inod[loc_inod100],
153  inod[loc_inod010],
154  inod[loc_inod001],
155  inod[loc_inod101],
156  inod[loc_inod011]);
157  if (i+j+1 >= my_order) continue;
158  size_type loc_inod110 = reference_element_P::ilat2loc_inod (my_order, ilat(i+1, j+1, k));
159  size_type loc_inod111 = reference_element_P::ilat2loc_inod (my_order, ilat(i+1, j+1, k+1));
160  raw_put_prism (gmsh,
161  inod[loc_inod100],
162  inod[loc_inod110],
163  inod[loc_inod010],
164  inod[loc_inod101],
165  inod[loc_inod111],
166  inod[loc_inod011]);
167  }
168  }
169  }
170 }
171 template <class T>
172 static
173 void
174 put_hexahedron (ostream& gmsh, const geo_element& K, const numbering<T,sequential>& my_numb, const geo_basic<T,sequential>& omega)
175 {
176  typedef typename geo_basic<T,sequential>::size_type size_type;
177  typedef point_basic<size_type> ilat;
178  std::vector<size_type> inod;
179  my_numb.dis_idof (omega.sizes(), K, inod);
180  size_type my_order = my_numb.degree();
181  for (size_type i = 0; i < my_order; i++) {
182  for (size_type j = 0; j < my_order; j++) {
183  for (size_type k = 0; k < my_order; k++) {
184  size_type loc_inod000 = reference_element_H::ilat2loc_inod (my_order, ilat(i, j, k));
185  size_type loc_inod100 = reference_element_H::ilat2loc_inod (my_order, ilat(i+1, j, k));
186  size_type loc_inod110 = reference_element_H::ilat2loc_inod (my_order, ilat(i+1, j+1, k));
187  size_type loc_inod010 = reference_element_H::ilat2loc_inod (my_order, ilat(i, j+1, k));
188  size_type loc_inod001 = reference_element_H::ilat2loc_inod (my_order, ilat(i, j, k+1));
189  size_type loc_inod101 = reference_element_H::ilat2loc_inod (my_order, ilat(i+1, j, k+1));
190  size_type loc_inod011 = reference_element_H::ilat2loc_inod (my_order, ilat(i, j+1, k+1));
191  size_type loc_inod111 = reference_element_H::ilat2loc_inod (my_order, ilat(i+1, j+1, k+1));
192  gmsh << "8\t" << inod[loc_inod000] << " "
193  << inod[loc_inod100] << " "
194  << inod[loc_inod110] << " "
195  << inod[loc_inod010] << " "
196  << inod[loc_inod001] << " "
197  << inod[loc_inod101] << " "
198  << inod[loc_inod111] << " "
199  << inod[loc_inod011] << endl;
200  }
201  }
202  }
203 }
204 #endif // TODO
205 template <class T>
206 static
207 void
208 put (ostream& gmsh, const geo_element& K, const numbering<T,sequential>& my_numb, const geo_basic<T,sequential>& omega, const array<point_basic<Float>,sequential>& my_node)
209 {
210  switch (K.variant()) {
211 #ifdef TODO
212  case reference_element::p: gmsh << "1\t" << K[0] << endl; break;
213 #endif // TODO
214  case reference_element::e: put_edge (gmsh, K, my_numb, omega); break;
215  case reference_element::t: put_triangle (gmsh, K, my_numb, omega); break;
216 #ifdef TODO
217  case reference_element::q: put_quadrangle (gmsh, K, my_numb, omega); break;
218  case reference_element::T: put_tetrahedron (gmsh, K, my_numb, omega); break;
219  case reference_element::P: put_prism (gmsh, K, my_numb, omega, my_node); break;
220  case reference_element::H: put_hexahedron (gmsh, K, my_numb, omega); break;
221 #endif // TODO
222  default: error_macro ("unsupported element variant `" << K.name() <<"'");
223  }
224 }
225 
226 template <class T>
227 odiststream&
228 geo_put_gmsh (odiststream& ops, const geo_basic<T,sequential>& omega, const numbering<T,sequential>& my_numb, const array<point_basic<T>,sequential>& my_node)
229 {
230  typedef typename geo_basic<T,sequential>::size_type size_type;
231  size_type my_order = my_numb.degree();
232  ostream& gmsh = ops.os();
233  gmsh << setprecision (std::numeric_limits<T>::digits10)
234  << "$MeshFormat" << endl
235  << "2.2 0 8" << endl
236  << "$EndMeshFormat" << endl;
237  gmsh << "$Nodes" << endl
238  << my_node.size() << endl;
239  for (size_type inod = 0, nnod = my_node.size(); inod < nnod; inod++) {
240  gmsh << inod+1 << " " << my_node[inod] << endl;
241  }
242  gmsh << "$EndNodes" << endl;
243  size_type map_dim = omega.map_dimension();
244  gmsh << "$Elements" << endl
245  << omega.size() << endl;
246  for (size_type ie = 0, ne = omega.size(); ie < ne; ie++) {
247  const geo_element& K = omega.get_geo_element (map_dim, ie);
248  put (gmsh, K, my_numb, omega, my_node);
249  }
250  gmsh << "$EndElements" << endl;
251  return ops;
252 }
253 template <class T>
254 odiststream&
255 geo_put_gmsh (odiststream& ops, const geo_basic<T,sequential>& omega)
256 {
257  numbering<T,sequential> my_numb ("P" + itos(omega.order()));
258  return geo_put_gmsh (ops, omega, my_numb, omega.get_nodes());
259 }
260 template odiststream& geo_put_gmsh<Float> (odiststream&, const geo_basic<Float,sequential>&, const numbering<Float,sequential>&, const array<point_basic<Float>,sequential>&);
261 template odiststream& geo_put_gmsh<Float> (odiststream&, const geo_basic<Float,sequential>&);
262 
263 }// namespace rheolef
264