Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
field2gmsh_pos.cc
Go to the documentation of this file.
1
21//
22// Passage du format field au format bb de Bamg
23//
24// author: Pierre.Saramito@imag.fr
25//
26// date: 03/09/1997 ; 04/01/2018
27//
28/*Prog:field2bb
29NAME: @code{field2gmsh_pos} - convert from .field file format to gmsh .pos one
30@pindex field2gmsh_pos
31@fiindex @file{.pos} gmsh field
32@fiindex @file{.field} field
33@toindex @code{gmsh}
34SYNOPSIS:
35@example
36 field2gmsh_pos < @var{input}.field > @var{output}.pos
37@end example
38
39DESCRIPTION:
40 Convert a @file{.field} file into a gmsh @file{.pos} one.
41 The output goes to standart output.
42 This command is useful for mesh adaptation with @code{gmsh}.
43End:*/
44#include "rheolef/compiler.h" // after index_set to avoid boost
45#include "scatch.icc"
46#include <iostream>
47#include <iomanip>
48#include <string>
49#include <vector>
50#include <array>
51#include <limits>
52#include <unistd.h>
53using namespace std;
54using namespace rheolef;
55
56namespace rheolef {
57
58template<class T>
59struct point_basic: std::array<T,3> {
60 typedef std::array<T,3> base;
61 point_basic(T x0=T(), T x1=T(), T x2=T())
62 : base() {
63 base::operator[](0) = x0;
64 base::operator[](1) = x1;
65 base::operator[](2) = x2;
66 }
67};
69
70} // namespace rheolef
71
72#include "reference_element_aux.icc"
73
74namespace rheolef {
75
76// TODO: move somewhere in reference_element_xxx.cc
77const char
79 'p',
80 'e',
81 't',
82 'q',
83 'T',
84 'P',
85 'H'
86};
87
88struct geo_element: std::array<size_t,8> {
89 void setname (char name) { _variant = reference_element_variant(name); }
90 char name() const { return _reference_element_name [_variant]; }
91 size_t variant() const { return _variant; }
92 size_t dimension()const { return reference_element_dimension_by_variant(variant()); }
93 size_t n_vertex() const { return reference_element_n_vertex (variant()); }
94 size_t size() const { return n_vertex(); }
95 geo_element() : array<size_t,8>(), _variant(-1) {}
96protected:
97 size_t _variant;
98};
99struct my_geo: vector<geo_element> {
100// allocator:
101 my_geo();
102// data:
103 size_t dimension;
104 size_t map_dimension;
105 size_t order;
106 string sys_coord;
107 array<size_t,reference_element__max_variant> size_by_variant;
108 array<size_t,3> size_by_dimension;
109 vector<point> node;
110};
111inline
112my_geo::my_geo()
113: dimension(0),
114 map_dimension(0),
115 order(1),
116 sys_coord("cartesian"),
117 size_by_variant(),
118 size_by_dimension(),
119 node()
120{
121}
122
123} // namespace rheolef
124
125void
126get_geo (istream& in, my_geo& omega)
127{
128 static const char* label_variant [] = {
129 "nodes", "edges", "triangles", "quadrangles", "tetrahedra", "prisms", "hexahedra" };
130 check_macro (scatch(in,"\nmesh",true), "input stream does not contains a geo.");
131 size_t version;
132 in >> version;
133 check_macro (version == 4, "mesh format version 4 expected, but format version " << version << " founded");
134 // get header
135 check_macro (scatch(in,"header",true), "geo file format version 4: \"header\" keyword not found");
136 string label;
137 omega.size_by_dimension.fill (0);
138 while (in.good()) {
139 size_t value;
140 in >> label;
141 if (label == "end") break;
142 if (label == "dimension") { in >> omega.dimension; }
143 else if (label == "order") { in >> omega.order; }
144 else if (label == "coordinate_system") { in >> omega.sys_coord; }
145 else {
146 size_t variant = 0;
147 for (; variant < reference_element__max_variant; variant++) {
148 if (label == label_variant[variant]) break;
149 }
150 check_macro (variant < reference_element__max_variant, "unexpected header member: \""<<label<<"\"");
151 in >> omega.size_by_variant [variant];
152 size_t dim = reference_element_dimension_by_variant(variant);
153 omega.size_by_dimension [dim] += omega.size_by_variant [variant];
154 }
155 }
156 check_macro (omega.order == 1, "unsupported geo order = "<<omega.order<<" > 1");
157 for (omega.map_dimension = 3; omega.map_dimension + 1 >= 1; omega.map_dimension--) {
158 if (omega.size_by_dimension [omega.map_dimension] != 0) break;
159 }
160 in >> label;
161 check_macro (label == "header", "geo file format version 4: \"end header\" keyword not found");
162 // get nodes
163 omega.node.resize (omega.size_by_variant[reference_element__p]);
164 for (size_t inod = 0; inod < omega.node.size(); ++inod) {
165 for (size_t i = 0; i < omega.dimension; ++i) {
166 in >> omega.node[inod][i];
167 }
168 }
169 // get elements
170 omega.resize (omega.size_by_dimension[omega.map_dimension]);
171 for (size_t ie = 0; ie < omega.size(); ++ie) {
172 char name;
173 in >> name;
174 omega[ie].setname(name);
175 for (size_t j = 0; j < omega[ie].size(); ++j) {
176 in >> omega[ie][j];
177 }
178 }
179}
180int main() {
181 //-------------------------
182 // input field
183 //-------------------------
184 scatch (cin, "field", true);
185 size_t version;
186 cin >> version;
187 check_macro (version == 1, "version " << version << " field format not yet supported");
188 size_t nbpts;
189 cin >> nbpts;
190 string geo_name;
191 string approx;
192 cin >> geo_name
193 >> approx;
194 check_macro (approx == "P1", "approx " << approx << " field not yet supported");
195 vector<double> uh (nbpts);
196 for (size_t i=0; i<nbpts; i++) {
197 cin >> uh[i];
198 }
199 //----------------------------
200 // input geo
201 //----------------------------
202 my_geo omega;
203 if (file_exists (geo_name + ".geo")) {
204 ifstream in ((geo_name + ".geo").c_str());
205 get_geo (in, omega);
206 } else if (file_exists (geo_name + ".geo.gz")) {
207 // use fifo:
208 // https://stackoverflow.com/questions/40740914/using-istream-to-read-from-named-pipe
209 bool do_verbose = true;
210 string fifo_name = tmpnam(0);
211#ifdef TO_CLEAN
212 string command = "rm -f " + fifo_name;
213 if (do_verbose) cerr << "! " << command << endl;
214 check_macro (system (command.c_str()) == 0, "adapt: unix command failed");
215#endif // TO_CLEAN
216 unlink (fifo_name.c_str()); // remove
217 mkfifo (fifo_name.c_str(), S_IWUSR | S_IRUSR); // | S_IRGRP | S_IROTH);
218 string command = "zcat " + geo_name + ".geo.gz > " + fifo_name + " &";
219 if (do_verbose) cerr << "! " << command << endl;
220 check_macro (system (command.c_str()) == 0, "adapt: unix command failed");
221 ifstream in (fifo_name.c_str());
222 get_geo (in, omega);
223 unlink (fifo_name.c_str()); // remove
224#ifdef TO_CLEAN
225 command = "rm -f " + fifo_name;
226 if (do_verbose) cerr << "! " << command << endl;
227 check_macro (system (command.c_str()) == 0, "adapt: unix command failed");
228#endif // TO_CLEAN
229 } else {
230 error_macro ("file \""<<geo_name<< ".geo[.gz]\" not found");
231 }
232 //----------------------------
233 // output pos
234 //----------------------------
235 /*
236 type #list-of-coords #list-of-values
237 --------------------------------------------------------------------
238 Scalar point SP 3 1 * nb-time-steps
239 Vector point VP 3 3 * nb-time-steps
240 Tensor point TP 3 9 * nb-time-steps
241 Scalar line SL 6 2 * nb-time-steps
242 Vector line VL 6 6 * nb-time-steps
243 Tensor line TL 6 18 * nb-time-steps
244 Scalar triangle ST 9 3 * nb-time-steps
245 Vector triangle VT 9 9 * nb-time-steps
246 Tensor triangle TT 9 27 * nb-time-steps
247 Scalar quadrangle SQ 12 4 * nb-time-steps
248 Vector quadrangle VQ 12 12 * nb-time-steps
249 Tensor quadrangle TQ 12 36 * nb-time-steps
250 Scalar tetrahedron SS 12 4 * nb-time-steps
251 Vector tetrahedron VS 12 12 * nb-time-steps
252 Tensor tetrahedron TS 12 36 * nb-time-steps
253 Scalar hexahedron SH 24 8 * nb-time-steps
254 Vector hexahedron VH 24 24 * nb-time-steps
255 Tensor hexahedron TH 24 72 * nb-time-steps
256 Scalar prism SI 18 6 * nb-time-steps
257 Vector prism VI 18 18 * nb-time-steps
258 Tensor prism TI 18 54 * nb-time-steps
259 Scalar pyramid SY 15 5 * nb-time-steps
260 Vector pyramid VY 15 15 * nb-time-steps
261 Tensor pyramid TY 15 45 * nb-time-steps
262 */
263 cout << setprecision(numeric_limits<double>::digits10)
264 << "View \"scalar\" {" << endl;
265 static char gmsh_pos_type [reference_element__max_variant];
266 gmsh_pos_type [reference_element__p] = 'P';
267 gmsh_pos_type [reference_element__e] = 'L';
268 gmsh_pos_type [reference_element__t] = 'T';
269 gmsh_pos_type [reference_element__q] = 'Q';
270 gmsh_pos_type [reference_element__T] = 'S';
271 gmsh_pos_type [reference_element__H] = 'H';
272 gmsh_pos_type [reference_element__P] = 'I';
273#ifdef TODDO
274 switch (uh.get_space().valued_tag()) {
275 // ---------------------------
277 // ---------------------------
278#endif // TODO
279 for (size_t ie = 0, ne = omega.size(); ie < ne; ie++) {
280 const geo_element& K = omega[ie];
281 cout << "S" << gmsh_pos_type[K.variant()] << "(";
282 for (size_t iloc = 0, nloc = K.size(); iloc < nloc; iloc++) {
283 const point& x = omega.node[K[iloc]];
284 for (size_t i_comp = 0; i_comp < 3; i_comp++) {
285 cout << x[i_comp];
286 if (i_comp != 2) cout << ",";
287 }
288 if (iloc != nloc-1) cout << ",";
289 }
290 cout << "){";
291 for (size_t iloc = 0, nloc = K.size(); iloc < nloc; iloc++) {
292 cout << uh[K[iloc]];
293 if (iloc != nloc-1) cout << ",";
294 }
295 cout << "};" << endl;
296 }
297#ifdef TODO
298 }
299 // ---------------------------
301 // ---------------------------
302 break;
303 }
304 // ---------------------------
306 // ---------------------------
307 field_component_const<T,sequential> uh_comp [3][3];
308 size_t d = uh.get_geo().dimension();
309 for (size_t i_comp = 0; i_comp < d; i_comp++) {
310 for (size_t j_comp = 0; j_comp < d; j_comp++) {
311 uh_comp[i_comp][j_comp].proxy_assign(uh(i_comp,j_comp));
312 }}
313 tensor_basic<T> value;
314#define HAVE_ID_DEFAULT_TENSOR
315#ifdef HAVE_ID_DEFAULT_TENSOR
316 // default (3,3) is 1 instead of 0
317 value (0,0) = value (1,1) = value (2,2) = 1;
318#endif
319 for (size_t ie = 0, ne = omega.size(); ie < ne; ie++) {
320 const geo_element& K = omega[ie];
321 gmsh << "T" << gmsh_pos_type[K.variant()] << "(";
322 for (size_t iloc = 0, nloc = K.size(); iloc < nloc; iloc++) {
323 const point_basic<T>& x = omega.node(K[iloc]);
324 for (size_t i_comp = 0; i_comp < 3; i_comp++) {
325 gmsh << x[i_comp];
326 if (i_comp != 2) gmsh << ",";
327 }
328 if (iloc != nloc-1) gmsh << ",";
329 }
330 gmsh << "){";
331 for (size_t iloc = 0, nloc = K.size(); iloc < nloc; iloc++) {
332 size_t inod = K[iloc];
333 const point_basic<T>& x = omega.node(inod);
334 for (size_t i_comp = 0; i_comp < d; i_comp++) {
335 for (size_t j_comp = 0; j_comp < d; j_comp++) {
336 value(i_comp,j_comp) = uh_comp [i_comp][j_comp].dof (inod);
337 }}
338 for (size_t i_comp = 0; i_comp < 3; i_comp++) {
339 for (size_t j_comp = 0; j_comp < 3; j_comp++) {
340 gmsh << value(i_comp,j_comp);
341 if (i_comp != 2 || j_comp != 2 || iloc != nloc-1) gmsh << ",";
342 }}
343 }
344 gmsh << "};" << endl;
345 }
346 break;
347 }
348 default: error_macro ("put_gmsh: do not known how to print " << uh.valued() << "-valued field");
349 }
350#endif // TODO
351 cout << "};" << endl;
352}
see the point page for the full documentation
see the geo_element page for the full documentation
size_t variant() const
size_type size() const
size_t dimension() const
void setname(char name)
variant_type variant() const
size_t n_vertex() const
point_basic(T x0=T(), T x1=T(), T x2=T())
std::array< T, 3 > base
point_basic< Float > point
Definition point.h:163
constexpr size_t reference_element__max_variant
const char _reference_element_name[reference_element__max_variant]
#define error_macro(message)
Definition dis_macros.h:49
const size_t dimension
Definition edge.icc:64
int main()
void get_geo(istream &in, my_geo &omega)
Expr1::float_type T
Definition field_expr.h:230
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.
bool scatch(std::istream &in, const std::string &ch, bool full_match=true)
scatch: see the rheostream page for the full documentation
Definition scatch.icc:44
bool file_exists(const std::string &filename)
file_exists: see the rheostream page for the full documentation
Definition scatch.icc:34
STL namespace.