Rheolef  7.1
an efficient C++ finite element environment
field2gmsh_pos.cc
Go to the documentation of this file.
1 //
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
29 NAME: @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}
34 SYNOPSIS:
35 @example
36  field2gmsh_pos < @var{input}.field > @var{output}.pos
37 @end example
38 
39 DESCRIPTION:
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}.
43 End:*/
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>
53 using namespace std;
54 using namespace rheolef;
55 
56 namespace rheolef {
57 
58 template<class T>
59 struct 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 
74 namespace rheolef {
75 
76 // TODO: move somewhere in reference_element_xxx.cc
77 const char
78 _reference_element_name [reference_element__max_variant] = {
79  'p',
80  'e',
81  't',
82  'q',
83  'T',
84  'P',
85  'H'
86 };
87 
88 struct 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) {}
96 protected:
97  size_t _variant;
98 };
99 struct 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 };
111 inline
112 my_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 
125 void
126 get_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 }
180 int 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  // ---------------------------
276  case space_constant::scalar: {
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  // ---------------------------
300  case space_constant::vector: {
301  // ---------------------------
302  break;
303  }
304  // ---------------------------
305  case space_constant::tensor: {
306  // ---------------------------
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 }
rheolef::geo_element::size
size_t size() const
Definition: field2gmsh_pos.cc:94
rheolef::geo_element::_variant
size_t _variant
Definition: field2gmsh_pos.cc:97
rheolef::point_basic
Definition: point.h:87
check_macro
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
rheolef::file_exists
bool file_exists(const std::string &filename)
file_exists: see the rheostream page for the full documentation
Definition: scatch.icc:42
rheolef::geo_element::dimension
size_t dimension() const
Definition: field2gmsh_pos.cc:92
rheolef::field_component_const
Definition: field.h:207
rheolef::tensor_basic
Definition: tensor.h:90
mkgeo_ball.order
order
Definition: mkgeo_ball.sh:343
rheolef::geo_element::size
size_type size() const
Definition: geo_element.h:168
rheolef::geo_element::geo_element
geo_element()
Definition: field2gmsh_pos.cc:95
rheolef-config.version
version
Definition: rheolef-config.in:126
rheolef::geo_element
see the geo_element page for the full documentation
Definition: geo_element.h:102
rheolef::geo_element::variant
variant_type variant() const
Definition: geo_element.h:161
rheolef::geo_element::variant
size_t variant() const
Definition: field2gmsh_pos.cc:91
rheolef::point_basic::point_basic
point_basic(T x0=T(), T x1=T(), T x2=T())
Definition: field2gmsh_pos.cc:61
rheolef::geo_element::n_vertex
size_t n_vertex() const
Definition: field2gmsh_pos.cc:93
rheolef::space_constant::tensor
Definition: space_constant.h:138
mkgeo_ball.variant
variant
Definition: mkgeo_ball.sh:149
rheolef::space_constant::scalar
Definition: space_constant.h:136
rheolef::geo_element::setname
void setname(char name)
Definition: field2gmsh_pos.cc:89
gmsh
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
Definition: iorheo-members.h:83
rheolef::scatch
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:52
scatch.icc
main
int main()
Definition: field2gmsh_pos.cc:180
dimension
const size_t dimension
Definition: edge.icc:64
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
rheolef::field_component_const::proxy_assign
field_component_const< T, M > & proxy_assign(const field_component_const< T, M > &uh_comp)
Definition: field_component.h:280
error_macro
#define error_macro(message)
Definition: dis_macros.h:49
rheolef::field_component_const::dof
const T & dof(size_type idof) const
Definition: field_component.h:169
point
see the point page for the full documentation
rheolef::point_basic::base
std::array< T, 3 > base
Definition: field2gmsh_pos.cc:60
mkgeo_grid.sys_coord
sys_coord
Definition: mkgeo_grid.sh:171
mkgeo_ball.d
d
Definition: mkgeo_ball.sh:154
rheolef::space_constant::vector
Definition: space_constant.h:137
rheolef::geo_element::name
char name() const
Definition: field2gmsh_pos.cc:90
get_geo
void get_geo(istream &in, my_geo &omega)
Definition: field2gmsh_pos.cc:126
mkgeo_ball.dim
dim
Definition: mkgeo_ball.sh:307
n_vertex
const size_t n_vertex
Definition: edge.icc:66
rheolef::std
Definition: vec_expr_v2.h:391
mkgeo_ball.command
command
Definition: mkgeo_ball.sh:136
mkgeo_contraction.name
name
Definition: mkgeo_contraction.sh:133
T
Expr1::float_type T
Definition: field_expr.h:218