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"
42put_edge (ostream&
gmsh,
const geo_element& K,
const basis_basic<T>& my_numb,
const geo_basic<T,sequential>& omega)
45 static size_type order2gmsh_type [11] = {0, 1, 8, 26, 27, 28, 62, 63, 64, 65, 66 };
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;
59put_triangle (ostream& gmsh, const geo_element& K, const basis_basic<T>& my_numb, const geo_basic<T,sequential>& omega)
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;
78put_quadrangle (ostream& gmsh, const geo_element& K, const basis_basic<T>& my_numb, const geo_basic<T,sequential>& omega)
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;
101put_tetrahedron (ostream& gmsh, const geo_element& K, const basis_basic<T>& my_numb, const geo_basic<T,sequential>& omega)
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;
153raw_put_prism (ostream& gmsh,
154 size_t i000, size_t i100, size_t i010,
155 size_t i001, size_t i101, size_t i011)
157 // gmsh prism has swaped x & y axis order: 00z 10z 01z replaced by 00z 01z 10z
158 gmsh << "6\
t" << i000 << " "
168put_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)
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));
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));
208put_hexahedron (ostream& gmsh, const geo_element& K, const basis_basic<T>& my_numb, const geo_basic<T,sequential>& omega)
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;
242put (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)
244 switch (K.variant()) {
246 case reference_element::p: gmsh << "1\
t" << K[0] << endl; break;
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;
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;
256 default: error_macro ("unsupported element
variant `
" << K.name() <<"'");
259// ----------------------------------------------------------------------------
261// ----------------------------------------------------------------------------
265geo_put_gmsh (odiststream& ops, const geo_basic<T,sequential>& omega, const basis_basic<T>& my_numb, const disarray<point_basic<T>,sequential>& my_node)
270 typedef typename geo_basic<T,sequential>::size_type size_type;
271 size_type my_order = my_numb.degree();
272 ostream& gmsh = ops.os();
276 gmsh << setprecision (std::numeric_limits<T>::digits10)
277 << "$MeshFormat" << endl
279 << "$EndMeshFormat" << endl;
280 // TODO: add domains: scan by domain and add for earch element to a domain list
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;
289 gmsh << "$EndNodes" << endl;
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);
300 gmsh << "$EndElements" << endl;
305geo_put_gmsh (odiststream& ops, const geo_basic<T,sequential>& omega)
307 basis_basic<T> my_numb ("P" + std::to_string(omega.order()));
308 return geo_put_gmsh (ops, omega, my_numb, omega.get_nodes());
310// ----------------------------------------------------------------------------
311// instanciation in library
312// ----------------------------------------------------------------------------
313template odiststream& geo_put_gmsh<Float> (odiststream&, const geo_basic<Float,sequential>&, const basis_basic<Float>&, const disarray<point_basic<Float>,sequential>&);
314template odiststream& geo_put_gmsh<Float> (odiststream&, const geo_basic<Float,sequential>&);
field::size_type size_type
rheolef::space_base_rep< T, M > t
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.