GetFEM  5.4.2
getfem_export.cc
1 /*===========================================================================
2 
3  Copyright (C) 2001-2020 Yves Renard
4 
5  This file is a part of GetFEM
6 
7  GetFEM is free software; you can redistribute it and/or modify it
8  under the terms of the GNU Lesser General Public License as published
9  by the Free Software Foundation; either version 3 of the License, or
10  (at your option) any later version along with the GCC Runtime Library
11  Exception either version 3.1 or (at your option) any later version.
12  This program is distributed in the hope that it will be useful, but
13  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
14  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
15  License and GCC Runtime Library Exception for more details.
16  You should have received a copy of the GNU Lesser General Public License
17  along with this program; if not, write to the Free Software Foundation,
18  Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
19 
20 ===========================================================================*/
21 
22 #include <iomanip>
23 #include "getfem/dal_singleton.h"
24 #include "getfem/bgeot_comma_init.h"
25 #include "getfem/getfem_export.h"
26 
27 namespace getfem
28 {
29 
30  /* try to check if a quad or hexahedric cell is "rectangular" and oriented
31  along the axes */
32  template<typename C> static bool check_voxel(const C& c) {
33  scalar_type h[3];
34  unsigned N = c[0].size();
35  if (c.size() != (1U << N)) return false;
36  const base_node P0 = c[0];
37  h[0] = c[1][0] - P0[0];
38  h[1] = c[2][0] - P0[0];
39  if (c.size() != 4) h[2] = c[4][0] - P0[0];
40  for (unsigned i=1; i < c.size(); ++i) {
41  const base_node d = c[i] - P0;
42  for (unsigned j=0; j < N; ++j)
43  if (gmm::abs(d[j]) > 1e-7*h[j] && gmm::abs(d[j] - h[j]) > 1e-7*h[j])
44  return false;
45  }
46  return true;
47  }
48 
49  /* Base64 encoding (https://en.wikipedia.org/wiki/Base64) */
50  std::string base64_encode(const std::vector<unsigned char>& src)
51  {
52  const std::string table("ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789+/");
53  std::string dst;
54  for (size_type i = 0; i < src.size(); ++i) {
55  switch (i % 3) {
56  case 0:
57  dst.push_back(table[(src[i] & 0xFC) >> 2]);
58  if (i + 1 == src.size()) {
59  dst.push_back(table[(src[i] & 0x03) << 4]);
60  dst.push_back('=');
61  dst.push_back('=');
62  }
63  break;
64  case 1:
65  dst.push_back(table[((src[i - 1] & 0x03) << 4) | ((src[i + 0] & 0xF0) >> 4)]);
66  if (i + 1 == src.size()) {
67  dst.push_back(table[(src[i] & 0x0F) << 2]);
68  dst.push_back('=');
69  }
70  break;
71  case 2:
72  dst.push_back(table[((src[i - 1] & 0x0F) << 2) | ((src[i + 0] & 0xC0) >> 6)]);
73  dst.push_back(table[src[i] & 0x3F]);
74  break;
75  }
76  }
77  return dst;
78  }
79 
80  /* -------------------------------------------------------------
81  * VTK export
82  * ------------------------------------------------------------- */
83 
84  struct gf2vtk_dof_mapping : public std::vector<std::vector<unsigned> > {};
85  struct gf2vtk_vtk_type : public std::vector<int> {};
86 
87  typedef enum { VTK_VERTEX = 1,
88  VTK_LINE = 3,
89  VTK_TRIANGLE = 5,
90  VTK_PIXEL = 8,
91  VTK_QUAD = 9,
92  VTK_TETRA = 10,
93  VTK_VOXEL = 11,
94  VTK_HEXAHEDRON = 12,
95  VTK_WEDGE = 13,
96  VTK_PYRAMID = 14,
97  VTK_QUADRATIC_EDGE = 21,
98  VTK_QUADRATIC_TRIANGLE = 22,
99  VTK_QUADRATIC_QUAD = 23,
100  VTK_QUADRATIC_TETRA = 24,
101  VTK_QUADRATIC_HEXAHEDRON = 25,
102  VTK_QUADRATIC_WEDGE = 26,
103  VTK_QUADRATIC_PYRAMID = 27,
104  VTK_BIQUADRATIC_QUAD = 28,
105  VTK_TRIQUADRATIC_HEXAHEDRON = 29,
106  VTK_BIQUADRATIC_QUADRATIC_WEDGE = 32 } vtk_cell_type;
107 
108  typedef enum { NO_VTK_MAPPING,
109  N1_TO_VTK_VERTEX,
110  N2_TO_VTK_LINE,
111  N3_TO_VTK_TRIANGLE,
112  N4_TO_VTK_PIXEL,
113  N4_TO_VTK_QUAD,
114  N4_TO_VTK_TETRA,
115  N8_TO_VTK_VOXEL,
116  N8_TO_VTK_HEXAHEDRON,
117  N6_TO_VTK_WEDGE,
118  N5_TO_VTK_PYRAMID,
119  N3_TO_VTK_QUADRATIC_EDGE,
120  N6_TO_VTK_QUADRATIC_TRIANGLE,
121  N8_TO_VTK_QUADRATIC_QUAD,
122  N10_TO_VTK_QUADRATIC_TETRA,
123  N20_TO_VTK_QUADRATIC_HEXAHEDRON,
124  N15_TO_VTK_QUADRATIC_WEDGE,
125  N13_TO_VTK_QUADRATIC_PYRAMID,
126  N14_TO_VTK_QUADRATIC_PYRAMID,
127  N9_TO_VTK_BIQUADRATIC_QUAD,
128  N27_TO_VTK_TRIQUADRATIC_HEXAHEDRON,
129  N18_TO_VTK_BIQUADRATIC_QUADRATIC_WEDGE } vtk_mapping_type;
130 
131 
132  void init_gf2vtk() {
133  // see https://www.kitware.com/products/books/VTKUsersGuide.pdf
134  // and https://www.kitware.com/products/books/VTKTextbook.pdf
135  // (there are some conflicts between the two)
136  gf2vtk_dof_mapping &vtkmaps = dal::singleton<gf2vtk_dof_mapping>::instance();
137  gf2vtk_vtk_type &vtktypes = dal::singleton<gf2vtk_vtk_type>::instance();
138  vtkmaps.resize(25);
139  vtktypes.resize(25);
140 
141  vtktypes[N1_TO_VTK_VERTEX] = VTK_VERTEX;
142  vtkmaps [N1_TO_VTK_VERTEX] = {0};
143  vtktypes[N2_TO_VTK_LINE] = VTK_LINE;
144  vtkmaps [N2_TO_VTK_LINE] = {0, 1};
145  vtktypes[N3_TO_VTK_TRIANGLE] = VTK_TRIANGLE;
146  vtkmaps [N3_TO_VTK_TRIANGLE] = {0, 1, 2};
147  vtktypes[N4_TO_VTK_PIXEL] = VTK_PIXEL;
148  vtkmaps [N4_TO_VTK_PIXEL] = {0, 1, 2, 3};
149  vtktypes[N4_TO_VTK_QUAD] = VTK_QUAD;
150  vtkmaps [N4_TO_VTK_QUAD] = {0, 1, 3, 2};
151  vtktypes[N4_TO_VTK_TETRA] = VTK_TETRA;
152  vtkmaps [N4_TO_VTK_TETRA] = {0, 1, 2, 3};
153  vtktypes[N8_TO_VTK_VOXEL] = VTK_VOXEL;
154  vtkmaps [N8_TO_VTK_VOXEL] = {0, 1, 2, 3, 4, 5, 6, 7};
155  vtktypes[N8_TO_VTK_HEXAHEDRON] = VTK_HEXAHEDRON;
156  vtkmaps [N8_TO_VTK_HEXAHEDRON] = {0, 1, 3, 2, 4, 5, 7, 6};
157  vtktypes[N6_TO_VTK_WEDGE] = VTK_WEDGE;
158  vtkmaps [N6_TO_VTK_WEDGE] = {0, 1, 2, 3, 4, 5};
159  vtktypes[N5_TO_VTK_PYRAMID] = VTK_PYRAMID;
160  vtkmaps [N5_TO_VTK_PYRAMID] = {0, 1, 3, 2, 4};
161  vtktypes[N3_TO_VTK_QUADRATIC_EDGE] = VTK_QUADRATIC_EDGE;
162  vtkmaps [N3_TO_VTK_QUADRATIC_EDGE] = {0, 2, 1};
163  vtktypes[N6_TO_VTK_QUADRATIC_TRIANGLE] = VTK_QUADRATIC_TRIANGLE;
164  vtkmaps [N6_TO_VTK_QUADRATIC_TRIANGLE] = {0, 2, 5, 1, 4, 3};
165  vtktypes[N8_TO_VTK_QUADRATIC_QUAD] = VTK_QUADRATIC_QUAD;
166  vtkmaps [N8_TO_VTK_QUADRATIC_QUAD] = {0, 2, 7, 5, 1, 4, 6, 3};
167  vtktypes[N10_TO_VTK_QUADRATIC_TETRA] = VTK_QUADRATIC_TETRA;
168  vtkmaps [N10_TO_VTK_QUADRATIC_TETRA] = {0, 2, 5, 9, 1, 4, 3, 6, 7, 8};
169  vtktypes[N20_TO_VTK_QUADRATIC_HEXAHEDRON] = VTK_QUADRATIC_HEXAHEDRON;
170  vtkmaps [N20_TO_VTK_QUADRATIC_HEXAHEDRON] = {0, 2, 7, 5, 12, 14, 19, 17, 1, 4, 6, 3, 13, 16, 18, 15, 8, 9, 11, 10};
171  vtktypes[N15_TO_VTK_QUADRATIC_WEDGE] = VTK_QUADRATIC_WEDGE;
172  vtkmaps [N15_TO_VTK_QUADRATIC_WEDGE] = {0, 2, 5, 9, 11, 14, 1, 4, 3, 10, 13, 12, 6, 7, 8};
173  // = {0, 5, 2, 9, 14, 11, 3, 4, 1, 12, 13, 10, 6, 8, 7};
174  vtktypes[N13_TO_VTK_QUADRATIC_PYRAMID] = VTK_QUADRATIC_PYRAMID;
175  vtkmaps [N13_TO_VTK_QUADRATIC_PYRAMID] = {0, 2, 7, 5, 12, 1, 4, 6, 3, 8, 9, 11, 10};
176  vtktypes[N14_TO_VTK_QUADRATIC_PYRAMID] = VTK_QUADRATIC_PYRAMID;
177  vtkmaps [N14_TO_VTK_QUADRATIC_PYRAMID] = {0, 2, 8, 6, 13, 1, 5, 7, 3, 9, 10, 12, 11};
178  vtktypes[N9_TO_VTK_BIQUADRATIC_QUAD] = VTK_BIQUADRATIC_QUAD;
179  vtkmaps [N9_TO_VTK_BIQUADRATIC_QUAD] = {0, 2, 8, 6, 1, 5, 7, 3, 4};
180  vtktypes[N27_TO_VTK_TRIQUADRATIC_HEXAHEDRON] = VTK_TRIQUADRATIC_HEXAHEDRON;
181  vtkmaps [N27_TO_VTK_TRIQUADRATIC_HEXAHEDRON] = {0, 2, 8, 6, 18, 20, 26, 24, 1, 5, 7, 3, 19, 23, 25, 21, 9, 11, 17, 15, 12, 14, 10, 16, 4, 22};
182  vtktypes[N18_TO_VTK_BIQUADRATIC_QUADRATIC_WEDGE] = VTK_BIQUADRATIC_QUADRATIC_WEDGE;
183  vtkmaps [N18_TO_VTK_BIQUADRATIC_QUADRATIC_WEDGE] = {0, 2, 5, 12, 14, 17, 1, 4, 3, 13, 16, 15, 6, 8, 11, 7, 10, 9};
184  }
185 
186  static const std::vector<unsigned> &
187  select_vtk_dof_mapping(unsigned t) {
188  gf2vtk_dof_mapping &vtkmaps = dal::singleton<gf2vtk_dof_mapping>::instance();
189  if (vtkmaps.size() == 0) init_gf2vtk();
190  return vtkmaps[t];
191  }
192 
193  int select_vtk_type(unsigned t) {
194  gf2vtk_vtk_type &vtktypes = dal::singleton<gf2vtk_vtk_type>::instance();
195  if (vtktypes.size() == 0) init_gf2vtk();
196  return vtktypes[t];
197  }
198 
199 
200  vtk_export::vtk_export(std::ostream &os_, bool ascii_, bool vtk_)
201  : os(os_), ascii(ascii_), vtk(vtk_) { init(); }
202 
203  vtk_export::vtk_export(const std::string& fname, bool ascii_, bool vtk_)
204  : os(real_os), ascii(ascii_), vtk(vtk_),
205  real_os(fname.c_str(), !ascii ? std::ios_base::binary | std::ios_base::out
206  : std::ios_base::out) {
207  GMM_ASSERT1(real_os, "impossible to write to file '" << fname << "'");
208  init();
209  }
210 
211  vtk_export::~vtk_export(){
212  if (!vtk) {
213  if (state == IN_CELL_DATA) os << "</CellData>\n";
214  if (state == IN_POINT_DATA) os << "</PointData>\n";
215  os << "</Piece>\n";
216  os << "</UnstructuredGrid>\n";
217  os << "</VTKFile>\n";
218  }
219  }
220 
221  void vtk_export::init() {
222  strcpy(header, "Exported by GetFEM");
223  psl = 0; dim_ = dim_type(-1);
224  static int test_endian = 0x01234567;
225  reverse_endian = (*((char*)&test_endian) == 0x67);
226  state = EMPTY;
227  clear_vals();
228  }
229 
230  void vtk_export::switch_to_cell_data() {
231  if (state != IN_CELL_DATA) {
232  if (vtk) {
233  size_type nb_cells=0;
234  if (psl) for (auto i : {0,1,2,3}) nb_cells += psl->nb_simplexes(i);
235  else nb_cells = pmf->convex_index().card();
236  write_separ();
237  os << "CELL_DATA " << nb_cells << "\n";
238  write_separ();
239  } else {
240  if (state == IN_POINT_DATA) os << "</PointData>\n";
241  os << "<CellData>\n";
242  }
243  state = IN_CELL_DATA;
244  }
245  }
246 
247  void vtk_export::switch_to_point_data() {
248  if (state != IN_POINT_DATA) {
249  if (vtk) {
250  write_separ();
251  os << "POINT_DATA " << (psl ? psl->nb_points()
252  : pmf_dof_used.card()) << "\n";
253  write_separ();
254  } else {
255  if (state == IN_CELL_DATA) os << "</CellData>\n";
256  os << "<PointData>\n";
257  }
258  state = IN_POINT_DATA;
259  }
260  }
261 
262 
263  void vtk_export::exporting(const stored_mesh_slice& sl) {
264  psl = &sl; dim_ = dim_type(sl.dim());
265  GMM_ASSERT1(psl->dim() <= 3, "attempt to export a " << int(dim_)
266  << "D slice (not supported)");
267  }
268 
269  void vtk_export::exporting(const mesh& m) {
270  dim_ = m.dim();
271  GMM_ASSERT1(dim_ <= 3, "attempt to export a " << int(dim_)
272  << "D mesh (not supported)");
273  pmf = std::make_unique<mesh_fem>(const_cast<mesh&>(m), dim_type(1));
274  for (dal::bv_visitor cv(m.convex_index()); !cv.finished(); ++cv) {
275  bgeot::pgeometric_trans pgt = m.trans_of_convex(cv);
276  pfem pf = getfem::classical_fem(pgt, pgt->complexity() > 1 ? 2 : 1);
277  pmf->set_finite_element(cv, pf);
278  }
279  exporting(*pmf);
280  }
281 
282  void vtk_export::exporting(const mesh_fem& mf) {
283  dim_ = mf.linked_mesh().dim();
284  GMM_ASSERT1(dim_ <= 3, "attempt to export a " << int(dim_)
285  << "D mesh_fem (not supported)");
286  if (&mf != pmf.get())
287  pmf = std::make_unique<mesh_fem>(mf.linked_mesh());
288  /* initialize pmf with finite elements suitable for VTK (which only knows
289  isoparametric FEMs of order 1 and 2) */
290  for (dal::bv_visitor cv(mf.convex_index()); !cv.finished(); ++cv) {
291  bgeot::pgeometric_trans pgt = mf.linked_mesh().trans_of_convex(cv);
292  pfem pf = mf.fem_of_element(cv);
293 
294  if (pf == fem_descriptor("FEM_Q2_INCOMPLETE(2)") ||
295  pf == fem_descriptor("FEM_Q2_INCOMPLETE(3)") ||
296  pf == fem_descriptor("FEM_PYRAMID_Q2_INCOMPLETE") ||
297  pf == fem_descriptor("FEM_PYRAMID_Q2_INCOMPLETE_DISCONTINUOUS") ||
298  pf == fem_descriptor("FEM_PRISM_INCOMPLETE_P2") ||
299  pf == fem_descriptor("FEM_PRISM_INCOMPLETE_P2_DISCONTINUOUS"))
300  pmf->set_finite_element(cv, pf);
301  else {
302  bool discontinuous = false;
303  for (unsigned i=0; i < pf->nb_dof(cv); ++i) {
304  /* could be a better test for discontinuity .. */
305  if (!dof_linkable(pf->dof_types()[i])) { discontinuous = true; break; }
306  }
307 
308  pfem classical_pf1 = discontinuous ? classical_discontinuous_fem(pgt, 1)
309  : classical_fem(pgt, 1);
310 
311  short_type degree = 1;
312  if ((pf != classical_pf1 && pf->estimated_degree() > 1) ||
313  pgt->structure() != pgt->basic_structure())
314  degree = 2;
315 
316  pmf->set_finite_element(cv, discontinuous ?
317  classical_discontinuous_fem(pgt, degree, 0, true) :
318  classical_fem(pgt, degree, true));
319  }
320  }
321  /* find out which dof will be exported to VTK/VTU */
322 
323  const mesh &m = pmf->linked_mesh();
324  pmf_mapping_type.resize(pmf->convex_index().last_true() + 1, unsigned(-1));
325  pmf_dof_used.sup(0, pmf->nb_basic_dof());
326  for (dal::bv_visitor cv(pmf->convex_index()); !cv.finished(); ++cv) {
327  vtk_mapping_type t = NO_VTK_MAPPING;
328  size_type nbd = pmf->fem_of_element(cv)->nb_dof(cv);
329  switch (pmf->fem_of_element(cv)->dim()) {
330  case 0: t = N1_TO_VTK_VERTEX; break;
331  case 1:
332  if (nbd == 2) t = N2_TO_VTK_LINE;
333  else if (nbd == 3) t = N3_TO_VTK_QUADRATIC_EDGE;
334  break;
335  case 2:
336  if (nbd == 3) t = N3_TO_VTK_TRIANGLE;
337  else if (nbd == 4)
338  t = check_voxel(m.points_of_convex(cv)) ? N4_TO_VTK_PIXEL
339  : N4_TO_VTK_QUAD;
340  else if (nbd == 6) t = N6_TO_VTK_QUADRATIC_TRIANGLE;
341  else if (nbd == 8) t = N8_TO_VTK_QUADRATIC_QUAD;
342  else if (nbd == 9) t = N9_TO_VTK_BIQUADRATIC_QUAD;
343  break;
344  case 3:
345  if (nbd == 4) t = N4_TO_VTK_TETRA;
346  else if (nbd == 10) t = N10_TO_VTK_QUADRATIC_TETRA;
347  else if (nbd == 8)
348  t = check_voxel(m.points_of_convex(cv)) ? N8_TO_VTK_VOXEL
349  : N8_TO_VTK_HEXAHEDRON;
350  else if (nbd == 20) t = N20_TO_VTK_QUADRATIC_HEXAHEDRON;
351  else if (nbd == 27) t = N27_TO_VTK_TRIQUADRATIC_HEXAHEDRON;
352  else if (nbd == 5) t = N5_TO_VTK_PYRAMID;
353  else if (nbd == 13) t = N13_TO_VTK_QUADRATIC_PYRAMID;
354  else if (nbd == 14) t = N14_TO_VTK_QUADRATIC_PYRAMID;
355  else if (nbd == 6) t = N6_TO_VTK_WEDGE;
356  else if (nbd == 15) t = N15_TO_VTK_QUADRATIC_WEDGE;
357  else if (nbd == 18) t = N18_TO_VTK_BIQUADRATIC_QUADRATIC_WEDGE;
358  break;
359  }
360  GMM_ASSERT1(t != -1, "semi internal error. Could not map " <<
361  name_of_fem(pmf->fem_of_element(cv))
362  << " to a VTK cell type");
363  pmf_mapping_type[cv] = t;
364 
365  const std::vector<unsigned> &dmap = select_vtk_dof_mapping(t);
366  //cout << "nbd = " << nbd << ", t = " << t << ", dmap = "<<dmap << "\n";
367  GMM_ASSERT1(dmap.size() <= pmf->nb_basic_dof_of_element(cv),
368  "inconsistency in vtk_dof_mapping");
369  for (unsigned i=0; i < dmap.size(); ++i)
370  pmf_dof_used.add(pmf->ind_basic_dof_of_element(cv)[dmap[i]]);
371  }
372  // cout << "mf.nb_dof = " << mf.nb_dof() << ", pmf->nb_dof="
373  // << pmf->nb_dof() << ", dof_used = " << pmf_dof_used.card() << "\n";
374  }
375 
376 
377  const stored_mesh_slice& vtk_export::get_exported_slice() const {
378  GMM_ASSERT1(psl, "no slice!");
379  return *psl;
380  }
381 
382  const mesh_fem& vtk_export::get_exported_mesh_fem() const {
383  GMM_ASSERT1(pmf.get(), "no mesh_fem!");
384  return *pmf;
385  }
386 
387  void vtk_export::set_header(const std::string& s)
388  {
389  strncpy(header, s.c_str(), 256);
390  header[255] = 0;
391  }
392 
393  void vtk_export::check_header() {
394  if (state >= HEADER_WRITTEN) return;
395  if (vtk) {
396  os << "# vtk DataFile Version 2.0\n";
397  os << header << "\n";
398  os << (ascii ? "ASCII\n" : "BINARY\n");
399  } else {
400  os << "<?xml version=\"1.0\"?>\n";
401  os << "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" ";
402  os << "byte_order=\"" << (reverse_endian ? "LittleEndian" : "BigEndian") << "\">\n";
403  os << "<!--" << header << "-->\n";
404  os << "<UnstructuredGrid>\n";
405  }
406  state = HEADER_WRITTEN;
407  }
408 
409  void vtk_export::write_separ()
410  { if (ascii) os << "\n"; }
411 
412  void vtk_export::clear_vals()
413  { if (!vtk && !ascii) vals.clear(); }
414 
415  void vtk_export::write_vals() {
416  if (!vtk && !ascii) {
417  os << base64_encode(vals);
418  clear_vals();
419  }
420  }
421 
422  void vtk_export::write_mesh() {
423  if (psl) write_mesh_structure_from_slice();
424  else write_mesh_structure_from_mesh_fem();
425  }
426 
427  /* export the slice data as an unstructured mesh composed of simplexes */
428  void vtk_export::write_mesh_structure_from_slice() {
429  /* element type code for (linear) simplexes of dimensions 0,1,2,3 in VTK */
430  static int vtk_simplex_code[4] = { VTK_VERTEX, VTK_LINE, VTK_TRIANGLE, VTK_TETRA };
431  if (state >= STRUCTURE_WRITTEN) return;
432  check_header();
433  /* count total number of simplexes, and total number of entries */
434  size_type cells_cnt = 0, splx_cnt = 0;
435  for (size_type ic=0; ic < psl->nb_convex(); ++ic) {
436  for (const slice_simplex &s : psl->simplexes(ic))
437  cells_cnt += s.dim() + 2;
438  splx_cnt += psl->simplexes(ic).size();
439  }
440  if (vtk) {
441  /* possible improvement: detect structured grids */
442  os << "DATASET UNSTRUCTURED_GRID\n";
443  os << "POINTS " << psl->nb_points() << " float\n";
444  } else {
445  os << "<Piece NumberOfPoints=\"" << psl->nb_points();
446  os << "\" NumberOfCells=\"" << splx_cnt << "\">\n";
447  os << "<Points>\n";
448  os << "<DataArray type=\"Float32\" Name=\"Points\" ";
449  os << "NumberOfComponents=\"3\" ";
450  os << (ascii ? "format=\"ascii\">\n" : "format=\"binary\">\n");
451  if (!ascii) write_val(int(sizeof(float)*psl->nb_points()*3));
452  }
453  /*
454  points are not merge, vtk is mostly fine with that (except for
455  transparency and normals at vertices of 3D elements: all simplex faces
456  are considered to be "boundary" faces by vtk)
457  */
458  for (size_type ic=0; ic < psl->nb_convex(); ++ic) {
459  for (size_type i=0; i < psl->nodes(ic).size(); ++i)
460  write_vec(psl->nodes(ic)[i].pt.begin(),psl->nodes(ic)[i].pt.size());
461  write_separ();
462  }
463  write_vals();
464  if (!vtk) {
465  os << (ascii ? "" : "\n") << "</DataArray>\n";
466  os << "</Points>\n";
467  }
468 
469  /* CELLS section */
470  size_type nodes_cnt = 0;
471  if (vtk) {
472  write_separ(); os << "CELLS " << splx_cnt << " " << cells_cnt << "\n";
473  } else {
474  os << "<Cells>\n";
475  os << "<DataArray type=\"Int32\" Name=\"connectivity\" ";
476  os << (ascii ? "format=\"ascii\">\n" : "format=\"binary\">\n");
477  if (!ascii) {
478  int size = 0;
479  for (size_type ic=0; ic < psl->nb_convex(); ++ic) {
480  for (const slice_simplex &s : psl->simplexes(ic))
481  size += int((s.dim()+1)*sizeof(int));
482  }
483  write_val(size);
484  }
485  }
486  for (size_type ic=0; ic < psl->nb_convex(); ++ic) {
487  for (const slice_simplex &s : psl->simplexes(ic)) {
488  if (vtk)
489  write_val(int(s.dim()+1));
490  for (size_type j=0; j < s.dim()+1; ++j)
491  write_val(int(s.inodes[j] + nodes_cnt));
492  write_separ();
493  }
494  nodes_cnt += psl->nodes(ic).size();
495  }
496  assert(nodes_cnt == psl->nb_points()); // sanity check
497  write_vals();
498  if (vtk) {
499  write_separ(); os << "CELL_TYPES " << splx_cnt << "\n";
500  } else {
501  os << (ascii ? "" : "\n") << "</DataArray>\n";
502  os << "<DataArray type=\"Int32\" Name=\"offsets\" ";
503  os << (ascii ? "format=\"ascii\">\n" : "format=\"binary\">\n");
504  if (!ascii) {
505  int size = 0;
506  for (size_type ic=0; ic < psl->nb_convex(); ++ic)
507  size += int(psl->simplexes(ic).size()*sizeof(int));
508  write_val(size);
509  }
510  }
511  int cnt = 0;
512  for (size_type ic=0; ic < psl->nb_convex(); ++ic) {
513  for (const slice_simplex &s : psl->simplexes(ic))
514  if (vtk) {
515  write_val(int(vtk_simplex_code[s.dim()]));
516  } else {
517  cnt += int(s.dim()+1);
518  write_val(cnt);
519  }
520  write_separ();
521  splx_cnt -= psl->simplexes(ic).size();
522  }
523  write_vals();
524  assert(splx_cnt == 0); // sanity check
525  if (!vtk) {
526  os << (ascii ? "" : "\n") << "</DataArray>\n";
527  os << "<DataArray type=\"Int32\" Name=\"types\" ";
528  os << (ascii ? "format=\"ascii\">\n" : "format=\"binary\">\n");
529  if (!ascii) {
530  int size = 0;
531  for (size_type ic=0; ic < psl->nb_convex(); ++ic)
532  size += int(psl->simplexes(ic).size()*sizeof(int));
533  write_val(size);
534  }
535  for (size_type ic=0; ic < psl->nb_convex(); ++ic)
536  for (const slice_simplex &s : psl->simplexes(ic))
537  write_val(int(vtk_simplex_code[s.dim()]));
538  write_vals();
539  os << "\n" << "</DataArray>\n";
540  os << "</Cells>\n";
541  }
542  state = STRUCTURE_WRITTEN;
543  }
544 
545 
546  void vtk_export::write_mesh_structure_from_mesh_fem() {
547  if (state >= STRUCTURE_WRITTEN) return;
548  check_header();
549  if (vtk) {
550  /* possible improvement: detect structured grids */
551  os << "DATASET UNSTRUCTURED_GRID\n";
552  os << "POINTS " << pmf_dof_used.card() << " float\n";
553  } else {
554  os << "<Piece NumberOfPoints=\"" << pmf_dof_used.card()
555  << "\" NumberOfCells=\"" << pmf->convex_index().card() << "\">\n";
556  os << "<Points>\n";
557  os << "<DataArray type=\"Float32\" Name=\"Points\" ";
558  os << "NumberOfComponents=\"3\" ";
559  os << (ascii ? "format=\"ascii\">\n" : "format=\"binary\">\n");
560  if (!ascii) write_val(int(sizeof(float)*pmf_dof_used.card()*3));
561  }
562  std::vector<int> dofmap(pmf->nb_dof());
563  int cnt = 0;
564  for (dal::bv_visitor d(pmf_dof_used); !d.finished(); ++d) {
565  dofmap[d] = cnt++;
566  base_node P = pmf->point_of_basic_dof(d);
567  write_vec(P.const_begin(),P.size());
568  write_separ();
569  }
570  write_vals();
571 
572  size_type nb_cell_values = 0;
573  for (dal::bv_visitor cv(pmf->convex_index()); !cv.finished(); ++cv)
574  nb_cell_values += (1 + select_vtk_dof_mapping(pmf_mapping_type[cv]).size());
575 
576  if (vtk) {
577  write_separ();
578  os << "CELLS " << pmf->convex_index().card() << " " << nb_cell_values << "\n";
579  } else {
580  os << (ascii ? "" : "\n") << "</DataArray>\n";
581  os << "</Points>\n";
582  os << "<Cells>\n";
583  os << "<DataArray type=\"Int32\" Name=\"connectivity\" ";
584  os << (ascii ? "format=\"ascii\">\n" : "format=\"binary\">\n");
585  if (!ascii) {
586  int size = 0;
587  for (dal::bv_visitor cv(pmf->convex_index()); !cv.finished(); ++cv) {
588  const std::vector<unsigned> &dmap = select_vtk_dof_mapping(pmf_mapping_type[cv]);
589  size += int(sizeof(int)*dmap.size());
590  }
591  write_val(size);
592  }
593  }
594 
595  for (dal::bv_visitor cv(pmf->convex_index()); !cv.finished(); ++cv) {
596  const std::vector<unsigned> &dmap = select_vtk_dof_mapping(pmf_mapping_type[cv]);
597  if (vtk) write_val(int(dmap.size()));
598  for (size_type i=0; i < dmap.size(); ++i)
599  write_val(int(dofmap[pmf->ind_basic_dof_of_element(cv)[dmap[i]]]));
600  write_separ();
601  }
602  write_vals();
603 
604  if (vtk) {
605  write_separ();
606  os << "CELL_TYPES " << pmf->convex_index().card() << "\n";
607  } else {
608  os << (ascii ? "" : "\n") << "</DataArray>\n";
609  os << "<DataArray type=\"Int32\" Name=\"offsets\" ";
610  os << (ascii ? "format=\"ascii\">\n" : "format=\"binary\">\n");
611  if (!ascii)
612  write_val(int(pmf->convex_index().card()*sizeof(int)));
613  cnt = 0;
614  for (dal::bv_visitor cv(pmf->convex_index()); !cv.finished(); ++cv) {
615  const std::vector<unsigned> &dmap = select_vtk_dof_mapping(pmf_mapping_type[cv]);
616  cnt += int(dmap.size());
617  write_val(cnt);
618  }
619  write_vals();
620  os << "\n" << "</DataArray>\n";
621  os << "<DataArray type=\"Int32\" Name=\"types\" ";
622  os << (ascii ? "format=\"ascii\">\n" : "format=\"binary\">\n");
623  if (!ascii)
624  write_val(int(pmf->convex_index().card()*sizeof(int)));
625  }
626  for (dal::bv_visitor cv(pmf->convex_index()); !cv.finished(); ++cv) {
627  write_val(int(select_vtk_type(pmf_mapping_type[cv])));
628  if (vtk) write_separ();
629  }
630  write_vals();
631  if (!vtk) os << "\n" << "</DataArray>\n" << "</Cells>\n";
632 
633  state = STRUCTURE_WRITTEN;
634  }
635 
636  void vtk_export::write_mesh_quality(const mesh &m) {
637  if (psl) {
638  mesh_fem mf(const_cast<mesh&>(m),1);
639  mf.set_classical_finite_element(0);
640  std::vector<scalar_type> q(mf.nb_dof());
641  for (size_type d=0; d < mf.nb_dof(); ++d) {
642  q[d] = m.convex_quality_estimate(mf.first_convex_of_basic_dof(d));
643  }
644  write_point_data(mf, q, "convex_quality");
645  } else {
646  std::vector<scalar_type> q(pmf->convex_index().card());
647  for (dal::bv_visitor cv(pmf->convex_index()); !cv.finished(); ++cv) {
648  q[cv] = m.convex_quality_estimate(cv);
649  }
650  write_cell_data(q, "convex_quality");
651  }
652  }
653 
654 
655  /* -------------------------------------------------------------
656  * OPENDX export
657  * ------------------------------------------------------------- */
658 
659  dx_export::dx_export(std::ostream &os_, bool ascii_)
660  : os(os_), ascii(ascii_) { init(); }
661 
662  dx_export::dx_export(const std::string& fname, bool ascii_, bool append_)
663  : os(real_os), ascii(ascii_) {
664  real_os.open(fname.c_str(),
665  std::ios_base::openmode(std::ios_base::in |
666  std::ios_base::out |
667  (append_ ? std::ios_base::ate : std::ios_base::trunc)));
668  GMM_ASSERT1(real_os.good(), "impossible to write to dx file '"
669  << fname << "'");
670  init();
671  if (append_) { reread_metadata(); header_written = true; }
672  }
673 
674  dx_export::~dx_export() {
675  std::ios::pos_type p = os.tellp();
676  write_series();
677  os << "\n# --end of getfem export\nend\n";
678  update_metadata(p);
679  }
680 
681  void dx_export::init() {
682  strcpy(header, "Exported by GetFEM");
683  psl = 0; dim_ = dim_type(-1); connections_dim = dim_type(-1);
684  psl_use_merged = false;
685  header_written = false;
686  }
687 
688  void dx_export::write_separ()
689  { if (ascii) os << "\n"; }
690 
691  template<typename T> static typename std::list<T>::iterator
692  get_from_name(std::list<T> &c,
693  const std::string& name, bool raise_error) {
694  for (typename std::list<T>::iterator it = c.begin();
695  it != c.end(); ++it) {
696  if (it->name == name) return it;
697  }
698  GMM_ASSERT1(!raise_error, "object not found in dx file: " << name);
699  return c.end();
700  }
701 
702  std::list<dx_export::dxMesh>::iterator
703  dx_export::get_mesh(const std::string& name, bool raise_error) {
704  return get_from_name(meshes,name,raise_error);
705  }
706  std::list<dx_export::dxObject>::iterator
707  dx_export::get_object(const std::string& name, bool raise_error) {
708  return get_from_name(objects,name,raise_error);
709  }
710 
711 
712  bool dx_export::new_mesh(std::string &name) {
713  name = default_name(name, int(meshes.size()), "mesh");
714  std::list<dxMesh>::iterator it = get_mesh(name, false);
715  if (it != meshes.end()) {
716  if (&(*it) != &current_mesh())
717  std::swap(current_mesh(),*it);
718  return false;
719  } else {
720  meshes.push_back(dxMesh()); meshes.back().name = name;
721  return true;
722  }
723  }
724 
725  void dx_export::exporting(const stored_mesh_slice& sl, bool merge_points,
726  std::string name) {
727  if (!new_mesh(name)) return;
728  psl_use_merged = merge_points;
729  if (merge_points) sl.merge_nodes();
730  psl = &sl; dim_ = dim_type(sl.dim());
731  GMM_ASSERT1(psl->dim() <= 3, "4D slices and more are not supported");
732  for (dim_type d = 0; d <= psl->dim(); ++d) {
733  if (psl->nb_simplexes(d)) {
734  if (connections_dim == dim_type(-1)) connections_dim = d;
735  else GMM_ASSERT1(false, "Cannot export a slice containing "
736  "simplexes of different dimensions");
737  }
738  }
739  GMM_ASSERT1(connections_dim != dim_type(-1), "empty slice!");
740  }
741 
742 
743  void dx_export::exporting(const mesh_fem& mf, std::string name) {
744  name = default_name(name, int(meshes.size()), "mesh");
745  if (!new_mesh(name)) return;
746  const mesh &m = mf.linked_mesh();
747  GMM_ASSERT1(mf.linked_mesh().convex_index().card() != 0,
748  "won't export an empty mesh");
749 
750  dim_ = m.dim();
751  GMM_ASSERT1(dim_ <= 3, "4D meshes and more are not supported");
752  if (&mf != pmf.get())
753  pmf = std::make_unique<mesh_fem>(const_cast<mesh&>(m), dim_type(1));
754  bgeot::pgeometric_trans pgt = m.trans_of_convex(m.convex_index().first_true());
755  GMM_ASSERT1(dxname_of_convex_structure
756  (basic_structure(pgt->structure())) != 0,
757  "DX Cannot handle " <<
758  bgeot::name_of_geometric_trans(pgt) << ", use slices");
759  /* initialize pmf with finite elements suitable for OpenDX */
760  for (dal::bv_visitor cv(mf.convex_index()); !cv.finished(); ++cv) {
761  bgeot::pgeometric_trans pgt2 = mf.linked_mesh().trans_of_convex(cv);
762  GMM_ASSERT1(basic_structure(pgt->structure()) ==
763  basic_structure(pgt2->structure()),
764  "Cannot export this mesh to opendx, it contains "
765  "different convex types. Slice it first.");
766  pfem pf = mf.fem_of_element(cv);
767  bool discontinuous = false;
768  for (unsigned i=0; i < pf->nb_dof(cv); ++i) {
769  /* could be a better test for discontinuity .. */
770  if (!dof_linkable(pf->dof_types()[i])) { discontinuous = true; break; }
771  }
772  pfem classical_pf1 = discontinuous ? classical_discontinuous_fem(pgt, 1)
773  : classical_fem(pgt, 1);
774  pmf->set_finite_element(cv, classical_pf1);
775  }
776  pmf_dof_used.add(0, pmf->nb_basic_dof());
777  connections_dim = dim_type(pmf->nb_basic_dof_of_element(m.convex_index().first_true()));
778  }
779 
780  void dx_export::exporting(const mesh& m, std::string name) {
781  dim_ = m.dim();
782  GMM_ASSERT1(dim_ <= 3, "4D meshes and more are not supported");
783  pmf = std::make_unique<mesh_fem>(const_cast<mesh&>(m), dim_type(1));
784  pmf->set_classical_finite_element(1);
785  exporting(*pmf, name);
786  }
787 
788  void dx_export::write_series() {
789  for (std::list<dxSeries>::const_iterator it = series.begin();
790  it != series.end(); ++it) {
791  if (it->members.size() == 0) continue;
792  size_type count = 0;
793  os << "\nobject \"" << it->name << "\" class series\n";
794  for (std::list<std::string>::const_iterator ito = it->members.begin();
795  ito != it->members.end(); ++ito, ++count) {
796  os << " member " << count << " \"" << (*ito) << "\"\n";
797  }
798  }
799  }
800 
801  void dx_export::serie_add_object_(const std::string &serie_name,
802  const std::string &object_name) {
803  std::list<dxSeries>::iterator it = series.begin();
804  while (it != series.end() && it->name != serie_name) ++it;
805  if (it == series.end()) {
806  series.push_back(dxSeries()); it = series.end(); --it;
807  it->name = serie_name;
808  }
809  it->members.push_back(object_name);
810  }
811 
812  void dx_export::serie_add_object(const std::string &serie_name,
813  const std::string &object_name) {
814  /* create a series for edge data if possible (the cost is null
815  and it may be useful) */
816  std::list<dxObject>::iterator ito = get_object(object_name, false);
817  if (ito != objects.end()) {
818  std::list<dxMesh>::iterator itm = get_mesh(ito->mesh);
819  if (itm != meshes.end() && (itm->flags & dxMesh::WITH_EDGES)) {
820  serie_add_object_(serie_name + "_edges",
821  object_name + "_edges");
822  }
823  }
824  /* fill the real serie */
825  serie_add_object_(serie_name, object_name);
826  }
827 
828  void dx_export::set_header(const std::string& s)
829  { strncpy(header, s.c_str(), 256); header[255] = 0; }
830 
831  void dx_export::check_header() {
832  if (header_written) return;
833  header_written = true;
834  os << "# data file for IBM OpenDX, generated by GetFem++ v "
835  << GETFEM_VERSION << "\n";
836  os << "# " << header << "\n";
837  }
838 
839  void dx_export::update_metadata(std::ios::pos_type pos_series) {
840  os.seekp(0,std::ios::end);
841  os << "# This file contains the following objects\n";
842  std::ios::pos_type pos_end = os.tellp();
843  for (std::list<dxSeries>::const_iterator it = series.begin();
844  it != series.end(); ++it) {
845  os << "#S \"" << it->name << "\" which contains:\n";
846  for (std::list<std::string>::const_iterator its = it->members.begin();
847  its != it->members.end(); ++its)
848  os << "#+ \"" << *its << "\"\n";
849  }
850  for (std::list<dxObject>::const_iterator it = objects.begin();
851  it != objects.end(); ++it) {
852  os << "#O \"" << it->name << "\" \"" << it->mesh << "\"\n";
853  }
854  for (std::list<dxMesh>::const_iterator it = meshes.begin();
855  it != meshes.end(); ++it) {
856  os << "#M \"" << it->name << "\" " << it->flags << "\n";
857  }
858  os << "#E \"THE_END\" " << std::setw(20) << pos_series << std::setw(20) << pos_end << "\n";
859  }
860 
861  void dx_export::reread_metadata() {
862  char line[512];
863  real_os.seekg(0, std::ios::end);
864  int count=0; char c;
865  unsigned long lu_end, lu_series;
866  do {
867  real_os.seekg(-1, std::ios::cur);
868  c = char(real_os.peek());
869  } while (++count < 512 && c != '#');
870  real_os.getline(line, sizeof line);
871  if (sscanf(line, "#E \"THE_END\" %lu %lu", &lu_series, &lu_end) != 2)
872  GMM_ASSERT1(false, "this file was not generated by getfem, "
873  "cannot append data to it!\n");
874  real_os.seekg(lu_end, std::ios::beg);
875  do {
876  char name[512]; unsigned n;
877  int pos;
878  real_os.getline(line, sizeof line);
879  if (sscanf(line, "#%c \"%512[^\"]\"%n", &c, name, &pos) < 1)
880  GMM_ASSERT1(false, "corrupted file! your .dx file is broken\n");
881  if (c == 'S') {
882  series.push_back(dxSeries()); series.back().name = name;
883  } else if (c == '+') {
884  series.back().members.push_back(name);
885  } else if (c == 'O') {
886  objects.push_back(dxObject()); objects.back().name = name;
887  sscanf(line+pos, " \"%512[^\"]\"", name); objects.back().mesh = name;
888  } else if (c == 'M') {
889  meshes.push_back(dxMesh()); meshes.back().name = name;
890  sscanf(line+pos, "%u", &n); meshes.back().flags = n;
891  } else if (c == 'E') {
892  break;
893  } else GMM_ASSERT1(false, "corrupted file! your .dx file is broken\n");
894  } while (1);
895  real_os.seekp(lu_series, std::ios::beg);
896  }
897 
898  void dx_export::write_convex_attributes(bgeot::pconvex_structure cvs) {
899  const char *s_elem_type = dxname_of_convex_structure(cvs);
900  if (!s_elem_type)
901  GMM_WARNING1("OpenDX won't handle this kind of convexes");
902  os << "\n attribute \"element type\" string \"" << s_elem_type << "\"\n"
903  << " attribute \"ref\" string \"positions\"\n\n";
904  }
905 
906  const char *dx_export::dxname_of_convex_structure(bgeot::pconvex_structure cvs) {
907  const char *s_elem_type = 0;
908  switch (cvs->dim()) {
909  /* TODO: do something for point data */
910  case 1: s_elem_type = "lines"; break;
911  case 2:
912  if (cvs->nb_points() == 3)
913  s_elem_type = "triangles";
914  else if (cvs->nb_points() == 4)
915  s_elem_type = "quads";
916  break;
917  case 3:
918  if (cvs->nb_points() == 4)
919  s_elem_type = "tetrahedra";
920  else if (cvs->nb_points() == 8)
921  s_elem_type = "cubes";
922  break;
923  }
924  return s_elem_type;
925  }
926 
927  void dx_export::write_mesh() {
928  check_header();
929  if (current_mesh().flags & dxMesh::STRUCTURE_WRITTEN) return;
930  if (psl) write_mesh_structure_from_slice();
931  else write_mesh_structure_from_mesh_fem();
932 
933  os << "\nobject \"" << current_mesh_name() << "\" class field\n"
934  << " component \"positions\" value \""
935  << name_of_pts_array(current_mesh_name()) << "\"\n"
936  << " component \"connections\" value \""
937  << name_of_conn_array(current_mesh_name()) << "\"\n";
938  current_mesh().flags |= dxMesh::STRUCTURE_WRITTEN;
939  }
940 
941  /* export the slice data as an unstructured mesh composed of simplexes */
942  void dx_export::write_mesh_structure_from_slice() {
943  os << "\nobject \"" << name_of_pts_array(current_mesh_name())
944  << "\" class array type float rank 1 shape "
945  << int(psl->dim())
946  << " items " << (psl_use_merged ? psl->nb_merged_nodes() : psl->nb_points());
947  if (!ascii) os << " " << endianness() << " binary";
948  os << " data follows\n";
949  if (psl_use_merged) {
950  for (size_type i=0; i < psl->nb_merged_nodes(); ++i) {
951  for (size_type k=0; k < psl->dim(); ++k)
952  write_val(float(psl->merged_point(i)[k]));
953  write_separ();
954  }
955  } else {
956  for (size_type ic=0; ic < psl->nb_convex(); ++ic) {
957  for (size_type i=0; i < psl->nodes(ic).size(); ++i)
958  for (size_type k=0; k < psl->dim(); ++k)
959  write_val(float(psl->nodes(ic)[i].pt[k]));
960  write_separ();
961  }
962  }
963 
964  os << "\nobject \"" << name_of_conn_array(current_mesh_name())
965  << "\" class array type int rank 1 shape "
966  << int(connections_dim+1)
967  << " items " << psl->nb_simplexes(connections_dim);
968  if (!ascii) os << " " << endianness() << " binary";
969  os << " data follows\n";
970 
971  size_type nodes_cnt = 0; /* <- a virer , global_index le remplace */
972  for (size_type ic=0; ic < psl->nb_convex(); ++ic) {
973  for (const slice_simplex &s : psl->simplexes(ic)) {
974  if (s.dim() == connections_dim) {
975  for (size_type j=0; j < s.dim()+1; ++j) {
976  size_type k;
977  if (psl_use_merged)
978  k = psl->merged_index(ic, s.inodes[j]);
979  else k = psl->global_index(ic, s.inodes[j]);
980  write_val(int(k));
981  }
982  write_separ();
983  }
984  }
985  nodes_cnt += psl->nodes(ic).size();
986  }
987 
988  write_convex_attributes(bgeot::simplex_structure(connections_dim));
989  assert(nodes_cnt == psl->nb_points()); // sanity check
990  }
991 
992  void dx_export::write_mesh_structure_from_mesh_fem() {
993  os << "\nobject \"" << name_of_pts_array(current_mesh_name())
994  << "\" class array type float rank 1 shape "
995  << int(pmf->linked_mesh().dim())
996  << " items " << pmf->nb_dof();
997  if (!ascii) os << " " << endianness() << " binary";
998  os << " data follows\n";
999 
1000  /* possible improvement: detect structured grids */
1001  for (size_type d = 0; d < pmf->nb_basic_dof(); ++d) {
1002  const base_node P = pmf->point_of_basic_dof(d);
1003  for (size_type k=0; k < dim_; ++k)
1004  write_val(float(P[k]));
1005  write_separ();
1006  }
1007 
1008  os << "\nobject \"" << name_of_conn_array(current_mesh_name())
1009  << "\" class array type int rank 1 shape "
1010  << int(connections_dim)
1011  << " items " << pmf->convex_index().card();
1012  if (!ascii) os << " " << endianness() << " binary";
1013  os << " data follows\n";
1014 
1015  for (dal::bv_visitor cv(pmf->convex_index()); !cv.finished(); ++cv) {
1016  for (size_type i=0; i < connections_dim; ++i)
1017  write_val(int(pmf->ind_basic_dof_of_element(cv)[i]));
1018  write_separ();
1019  }
1020  write_convex_attributes(basic_structure(pmf->linked_mesh().structure_of_convex(pmf->convex_index().first_true())));
1021  }
1022 
1023  void dx_export::exporting_mesh_edges(bool with_slice_edges) {
1024  write_mesh();
1025  if (current_mesh().flags & dxMesh::WITH_EDGES) return;
1026  if (psl) write_mesh_edges_from_slice(with_slice_edges);
1027  else write_mesh_edges_from_mesh();
1028  current_mesh().flags |= dxMesh::WITH_EDGES;
1029  os << "\nobject \"" << name_of_edges_array(current_mesh_name())
1030  << "\" class field\n"
1031  << " component \"positions\" value \""
1032  << name_of_pts_array(current_mesh_name()) << "\"\n"
1033  << " component \"connections\" value \""
1034  << name_of_conn_array(name_of_edges_array(current_mesh_name()))
1035  << "\"\n";
1036  }
1037 
1038  void dx_export::write_mesh_edges_from_slice(bool with_slice_edges) {
1039  std::vector<size_type> edges;
1040  dal::bit_vector slice_edges;
1041  psl->get_edges(edges, slice_edges, psl_use_merged);
1042  if (with_slice_edges) slice_edges.clear();
1043  os << "\nobject \""
1044  << name_of_conn_array(name_of_edges_array(current_mesh_name()))
1045  << "\" class array type int rank 1 shape 2"
1046  << " items " << edges.size()/2 - slice_edges.card();
1047  if (!ascii) os << " " << endianness() << " binary";
1048  os << " data follows\n";
1049  for (size_type i=0; i < edges.size()/2; ++i) {
1050  if (!slice_edges.is_in(i)) {
1051  write_val(int(edges[2*i]));
1052  write_val(int(edges[2*i+1]));
1053  }
1054  if ((i+1)%10 == 0) write_separ();
1055  }
1056  write_separ();
1057  write_convex_attributes(bgeot::simplex_structure(1));
1058  }
1059 
1060  void dx_export::write_mesh_edges_from_mesh() {
1061  bgeot::mesh_structure ms(pmf->linked_mesh()); ms.to_edges();
1062  os << "\nobject \""
1063  << name_of_conn_array(name_of_edges_array(current_mesh_name()))
1064  << "\" class array type int rank 1 shape 2"
1065  << " items " << ms.convex_index().card();
1066  if (!ascii) os << " " << endianness() << " binary";
1067  os << " data follows\n";
1068  for (dal::bv_visitor cv(ms.convex_index()); !cv.finished(); ++cv) {
1069  write_val(int(ms.ind_points_of_convex(cv)[0]));
1070  write_val(int(ms.ind_points_of_convex(cv)[1]));
1071  if ((cv+1)%20 == 0) write_separ();
1072  }
1073  write_separ();
1074  write_convex_attributes(bgeot::simplex_structure(1));
1075  }
1076 
1077 
1078  /* -------------------------------------------------------------
1079  * POS export (Gmsh post-processing format)
1080  * ------------------------------------------------------------- */
1081  struct gf2pos_dof_mapping : public std::vector<std::vector<unsigned> > {};
1082 
1083  static const std::vector<unsigned>& getfem_to_pos_dof_mapping(int t) {
1084  gf2pos_dof_mapping &dm = dal::singleton<gf2pos_dof_mapping>::instance();
1085  if (dm.size() == 0) {
1086  dm.resize(8);
1087  dm[pos_export::POS_PT] = {0};
1088  dm[pos_export::POS_LN] = {0, 1};
1089  dm[pos_export::POS_TR] = {0, 1, 2};
1090  dm[pos_export::POS_QU] = {0, 1, 3, 2};
1091  dm[pos_export::POS_SI] = {0, 1, 2, 3};
1092  dm[pos_export::POS_HE] = {0, 1, 3, 2, 4, 5, 7, 6};
1093  dm[pos_export::POS_PR] = {0, 1, 2, 3, 4, 5};
1094  dm[pos_export::POS_PY] = {0, 1, 3, 2, 4};
1095  }
1096  return dm[t];
1097  }
1098 
1099  pos_export::pos_export(std::ostream& osname): os(osname) {
1100  init();
1101  }
1102 
1103  pos_export::pos_export(const std::string& fname)
1104  : os(real_os), real_os(fname.c_str()) {
1105  GMM_ASSERT1(real_os, "impossible to write to pos file '" << fname << "'");
1106  init();
1107  }
1108 
1109  void pos_export::init() {
1110  strcpy(header, "Exported by GetFEM");
1111  state = EMPTY;
1112  view = 0;
1113  }
1114 
1115  void pos_export::set_header(const std::string& s){
1116  strncpy(header, s.c_str(), 256);
1117  header[255] = 0;
1118  }
1119 
1120  void pos_export::check_header() {
1121  if (state >= HEADER_WRITTEN) return;
1122  os << "/* " << header << " */\n";
1123  os << "General.FastRedraw = 0;\n";
1124  os << "General.ColorScheme = 1;\n";
1125  state = HEADER_WRITTEN;
1126  }
1127 
1128  void pos_export::exporting(const mesh& m) {
1129  if (state >= STRUCTURE_WRITTEN) return;
1130  dim = dim_type(m.dim());
1131  GMM_ASSERT1(int(dim) <= 3, "attempt to export a "
1132  << int(dim) << "D mesh (not supported)");
1133  pmf = std::make_unique<mesh_fem>(const_cast<mesh&>(m), dim_type(1));
1134  pmf->set_classical_finite_element(1);
1135  exporting(*pmf);
1136  state = STRUCTURE_WRITTEN;
1137  }
1138 
1139  void pos_export::exporting(const mesh_fem& mf) {
1140  if (state >= STRUCTURE_WRITTEN) return;
1141  dim = dim_type(mf.linked_mesh().dim());
1142  GMM_ASSERT1(int(dim) <= 3, "attempt to export a "
1143  << int(dim) << "D mesh_fem (not supported)");
1144  if (&mf != pmf.get())
1145  pmf = std::make_unique<mesh_fem>(mf.linked_mesh(), dim_type(1));
1146 
1147  /* initialize pmf with finite elements suitable for Gmsh */
1148  for (dal::bv_visitor cv(mf.convex_index()); !cv.finished(); ++cv) {
1149  bgeot::pgeometric_trans pgt = mf.linked_mesh().trans_of_convex(cv);
1150  pfem pf = mf.fem_of_element(cv);
1151 
1152  bool discontinuous = false;
1153  for (unsigned i=0; i < pf->nb_dof(cv); ++i) {
1154  // could be a better test for discontinuity ...
1155  if (!dof_linkable(pf->dof_types()[i])) { discontinuous = true; break; }
1156  }
1157  pfem classical_pf1 = discontinuous ?
1158  classical_discontinuous_fem(pgt, 1) : classical_fem(pgt, 1);
1159  pmf->set_finite_element(cv, classical_pf1);
1160  }
1161  psl = NULL;
1162 
1163  /* find out which dof will be exported to Gmsh */
1164  for (dal::bv_visitor cv(pmf->convex_index()); !cv.finished(); ++cv) {
1165  int t = -1;
1166  switch (pmf->fem_of_element(cv)->nb_dof(cv)){
1167  case 1: t = POS_PT; break;
1168  case 2: t = POS_LN; break;
1169  case 3: t = POS_TR; break;
1170  case 4:
1171  if ( 2 == pmf->fem_of_element(cv)->dim() ) t = POS_QU;
1172  else if (3 == pmf->fem_of_element(cv)->dim()) t = POS_SI;
1173  break;
1174  case 6: t = POS_PR; break;
1175  case 8: t = POS_HE; break;
1176  case 5: t = POS_PY; break;
1177  }
1178  GMM_ASSERT1(t != -1, "semi internal error. Could not map "
1179  << name_of_fem(pmf->fem_of_element(cv))
1180  << " to a POS primitive type");
1181  pos_cell_type.push_back(unsigned(t));
1182 
1183  const std::vector<unsigned>& dmap = getfem_to_pos_dof_mapping(t);
1184  GMM_ASSERT1(dmap.size() <= pmf->nb_basic_dof_of_element(cv),
1185  "inconsistency in pos_dof_mapping");
1186  std::vector<unsigned> cell_dof;
1187  cell_dof.resize(dmap.size(),unsigned(-1));
1188  for (size_type i=0; i < dmap.size(); ++i)
1189  cell_dof[i] = unsigned(pmf->ind_basic_dof_of_element(cv)[dmap[i]]);
1190  pos_cell_dof.push_back(cell_dof);
1191  }
1192  for (size_type i=0; i< pmf->nb_basic_dof(); ++i){
1193  std::vector<float> pt;
1194  pt.resize(dim,float(0));
1195  for (size_type j=0; j<dim; ++j)
1196  pt[j] = float(pmf->point_of_basic_dof(i)[j]);
1197  pos_pts.push_back(pt);
1198  }
1199  state = STRUCTURE_WRITTEN;
1200  }
1201 
1202  void pos_export::exporting(const stored_mesh_slice& sl) {
1203  if (state >= STRUCTURE_WRITTEN) return;
1204  psl = &sl;
1205  dim = dim_type(sl.dim());
1206  GMM_ASSERT1(int(dim) <= 3, "attempt to export a "
1207  << int(dim) << "D slice (not supported)");
1208 
1209  for (size_type ic=0, pcnt=0; ic < psl->nb_convex(); ++ic) {
1210  for (const slice_simplex &s : psl->simplexes(ic)) {
1211  int t = -1;
1212  switch (s.dim()){
1213  case 0: t = POS_PT; break;
1214  case 1: t = POS_LN; break;
1215  case 2: t = POS_TR; break;
1216  case 3: t = POS_SI; break;
1217  }
1218  GMM_ASSERT1(t != -1, "semi internal error.. could not map to a POS primitive type");
1219  pos_cell_type.push_back(unsigned(t));
1220 
1221  const std::vector<unsigned>& dmap = getfem_to_pos_dof_mapping(t);
1222  GMM_ASSERT1(dmap.size() <= size_type(t+1), "inconsistency in pos_dof_mapping");
1223 
1224  std::vector<unsigned> cell_dof;
1225  cell_dof.resize(dmap.size(),unsigned(-1));
1226  for (size_type i=0; i < dmap.size(); ++i)
1227  cell_dof[i] = unsigned(s.inodes[dmap[i]] + pcnt);
1228  pos_cell_dof.push_back(cell_dof);
1229  }
1230  for (const slice_node &n : psl->nodes(ic)) {
1231  std::vector<float> pt;
1232  pt.resize(dim,float(0));
1233  for (size_type i=0; i<dim; ++i)
1234  pt[i] = float(n.pt[i]);
1235  pos_pts.push_back(pt);
1236  }
1237  pcnt += psl->nodes(ic).size();
1238  }
1239  state = STRUCTURE_WRITTEN;
1240  }
1241 
1242  void pos_export::write(const mesh& m, const std::string &name){
1243  if (state >= IN_CELL_DATA) return;
1244  GMM_ASSERT1(int(m.dim()) <= 3, "attempt to export a "
1245  << int(m.dim()) << "D mesh (not supported)");
1246  pmf = std::make_unique<mesh_fem>(const_cast<mesh&>(m), dim_type(1));
1247  pmf->set_classical_finite_element(1);
1248  write(*pmf,name);
1249  state = IN_CELL_DATA;
1250  }
1251 
1252  void pos_export::write(const mesh_fem& mf, const std::string &name){
1253  if (state >= IN_CELL_DATA) return;
1254  check_header();
1255  exporting(mf);
1256 
1257  if (""==name) os << "View \"mesh " << view <<"\" {\n";
1258  else os << "View \"" << name <<"\" {\n";
1259 
1260  int t;
1261  std::vector<unsigned> cell_dof;
1262  std::vector<float> cell_dof_val;
1263  for (size_type cell = 0; cell < pos_cell_type.size(); ++cell) {
1264  t = pos_cell_type[cell];
1265  cell_dof = pos_cell_dof[cell];
1266  cell_dof_val.resize(cell_dof.size(),float(0));
1267  write_cell(t,cell_dof,cell_dof_val);
1268  }
1269 
1270  os << "};\n";
1271  os << "View[" << view << "].ShowScale = 0;\n";
1272  os << "View[" << view << "].ShowElement = 1;\n";
1273  os << "View[" << view << "].DrawScalars = 0;\n";
1274  os << "View[" << view << "].DrawVectors = 0;\n";
1275  os << "View[" << view++ << "].DrawTensors = 0;\n";
1276  state = IN_CELL_DATA;
1277  }
1278 
1279  void pos_export::write(const stored_mesh_slice& sl, const std::string &name){
1280  if (state >= IN_CELL_DATA) return;
1281  check_header();
1282  exporting(sl);
1283 
1284  if (""==name) os << "View \"mesh " << view <<"\" {\n";
1285  else os << "View \"" << name <<"\" {\n";
1286 
1287  int t;
1288  std::vector<unsigned> cell_dof;
1289  std::vector<float> cell_dof_val;
1290  for (size_type cell = 0; cell < pos_cell_type.size(); ++cell) {
1291  t = pos_cell_type[cell];
1292  cell_dof = pos_cell_dof[cell];
1293  cell_dof_val.resize(cell_dof.size(),float(0));
1294  write_cell(t,cell_dof,cell_dof_val);
1295  }
1296 
1297  os << "};\n";
1298  os << "View[" << view << "].ShowScale = 0;\n";
1299  os << "View[" << view << "].ShowElement = 1;\n";
1300  os << "View[" << view << "].DrawScalars = 0;\n";
1301  os << "View[" << view << "].DrawVectors = 0;\n";
1302  os << "View[" << view++ << "].DrawTensors = 0;\n";
1303  state = IN_CELL_DATA;
1304  }
1305 } /* end of namespace getfem. */
dal::singleton::instance
static T & instance()
Instance from the current thread.
Definition: dal_singleton.h:165
getfem_export.h
Export solutions to various formats.
bgeot::name_of_geometric_trans
std::string name_of_geometric_trans(pgeometric_trans p)
Get the string name of a geometric transformation.
Definition: bgeot_geometric_trans.cc:1173
getfem::stored_mesh_slice::nb_merged_nodes
size_type nb_merged_nodes() const
Return the number of merged nodes in slice.
Definition: getfem_mesh_slice.h:149
getfem::mesh_fem::first_convex_of_basic_dof
virtual size_type first_convex_of_basic_dof(size_type d) const
Shortcut for convex_to_dof(d)[0].
Definition: getfem_mesh_fem.cc:258
bgeot::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:49
dal_singleton.h
A simple singleton implementation.
getfem::mesh_fem::fem_of_element
virtual pfem fem_of_element(size_type cv) const
Return the basic fem associated with an element (if no fem is associated, the function will crash!...
Definition: getfem_mesh_fem.h:444
getfem::classical_fem
pfem classical_fem(bgeot::pgeometric_trans pgt, short_type k, bool complete=false)
Give a pointer on the structures describing the classical polynomial fem of degree k on a given conve...
Definition: getfem_fem.cc:4141
getfem::mesh::convex_quality_estimate
scalar_type convex_quality_estimate(size_type ic) const
Return an estimate of the convex quality (0 <= Q <= 1).
Definition: getfem_mesh.cc:450
bgeot::basic_structure
pconvex_structure basic_structure(pconvex_structure cv)
Original structure (if concerned)
Definition: bgeot_convex_structure.h:172
getfem::mesh_fem
Describe a finite element method linked to a mesh.
Definition: getfem_mesh_fem.h:148
getfem::stored_mesh_slice::get_edges
void get_edges(std::vector< size_type > &edges, dal::bit_vector &slice_edges, bool from_merged_nodes) const
Extract the list of mesh edges.
Definition: getfem_mesh_slice.cc:231
getfem::stored_mesh_slice::nb_points
size_type nb_points() const
Return the number of nodes in the slice.
Definition: getfem_mesh_slice.h:118
bgeot::short_type
gmm::uint16_type short_type
used as the common short type integer in the library
Definition: bgeot_config.h:72
bgeot::simplex_structure
pconvex_structure simplex_structure(dim_type nc)
Give a pointer on the structures of a simplex of dimension d.
Definition: bgeot_convex_structure.cc:148
getfem
GEneric Tool for Finite Element Methods.
Definition: getfem_accumulated_distro.h:46
bgeot::mesh_structure::convex_index
const dal::bit_vector & convex_index() const
Return the list of valid convex IDs.
Definition: bgeot_mesh_structure.h:91
getfem::dx_export::set_header
void set_header(const std::string &s)
the header is the second line of text in the exported file, you can put whatever you want – call this...
Definition: getfem_export.cc:828
getfem::mesh_fem::convex_index
const dal::bit_vector & convex_index() const
Get the set of convexes where a finite element has been assigned.
Definition: getfem_mesh_fem.h:191
getfem::stored_mesh_slice::nb_convex
size_type nb_convex() const
return the number of convexes of the original mesh referenced in the slice
Definition: getfem_mesh_slice.h:98
getfem::dx_export::serie_add_object
void serie_add_object(const std::string &serie_name, const std::string &object_name)
add an object (typically the name of a data field) to a 'series', i.e.
Definition: getfem_export.cc:812
getfem::dof_linkable
bool dof_linkable(pdof_description)
Says if the dof is linkable.
Definition: getfem_fem.cc:615
getfem::pfem
std::shared_ptr< const getfem::virtual_fem > pfem
type of pointer on a fem description
Definition: getfem_fem.h:244
getfem::mesh_fem::linked_mesh
const mesh & linked_mesh() const
Return a reference to the underlying mesh.
Definition: getfem_mesh_fem.h:279
getfem::stored_mesh_slice::dim
size_type dim() const
return the slice dimension
Definition: getfem_mesh_slice.h:108
getfem::stored_mesh_slice::nodes
const mesh_slicer::cs_nodes_ct & nodes(size_type ic) const
Return the list of nodes for the 'ic'th convex of the slice.
Definition: getfem_mesh_slice.h:120
bgeot::mesh_structure
Mesh structure definition.
Definition: bgeot_mesh_structure.h:73
bgeot::pconvex_structure
std::shared_ptr< const convex_structure > pconvex_structure
Pointer on a convex structure description.
Definition: bgeot_convex_structure.h:54
getfem::mesh
Describe a mesh (collection of convexes (elements) and points).
Definition: getfem_mesh.h:95
getfem::stored_mesh_slice::merged_point
const base_node merged_point(size_type i_merged) const
Return the physical position of the merged node.
Definition: getfem_mesh_slice.h:155
bgeot::pgeometric_trans
std::shared_ptr< const bgeot::geometric_trans > pgeometric_trans
pointer type for a geometric transformation
Definition: bgeot_geometric_trans.h:186
getfem::name_of_fem
std::string name_of_fem(pfem p)
get the string name of a fem descriptor.
Definition: getfem_fem.cc:4241
getfem::dx_export::current_mesh_name
std::string current_mesh_name()
return the name of current mesh (use exporting(...) to change the current mesh)
Definition: getfem_export.h:381
getfem::stored_mesh_slice::simplexes
const mesh_slicer::cs_simplexes_ct & simplexes(size_type ic) const
Return the list of simplexes for the 'ic'th convex of the slice.
Definition: getfem_mesh_slice.h:124
getfem::mesh_fem::set_classical_finite_element
void set_classical_finite_element(size_type cv, dim_type fem_degree, bool complete=false)
Set a classical (i.e.
Definition: getfem_mesh_fem.cc:174
bgeot_comma_init.h
convenient initialization of vectors via overload of "operator,".
getfem::classical_discontinuous_fem
pfem classical_discontinuous_fem(bgeot::pgeometric_trans pg, short_type k, scalar_type alpha=0, bool complete=false)
Give a pointer on the structures describing the classical polynomial discontinuous fem of degree k on...
Definition: getfem_fem.cc:4146
getfem::mesh_fem::nb_dof
virtual size_type nb_dof() const
Return the total number of degrees of freedom.
Definition: getfem_mesh_fem.h:562
getfem::fem_descriptor
pfem fem_descriptor(const std::string &name)
get a fem descriptor from its string name.
Definition: getfem_fem.cc:4232
getfem::stored_mesh_slice::nb_simplexes
void nb_simplexes(std::vector< size_type > &c) const
return the simplex count, in an array.
Definition: getfem_mesh_slice.h:114
getfem::dx_export::exporting_mesh_edges
void exporting_mesh_edges(bool with_slice_edge=true)
append edges information (if you want to draw the mesh and are using a refined slice.
Definition: getfem_export.cc:1023