GetFEM  5.4.2
getfem_export.h
Go to the documentation of this file.
1 /* -*- c++ -*- (enables emacs c++ mode) */
2 /*===========================================================================
3 
4  Copyright (C) 2001-2020 Yves Renard, Julien Pommier
5 
6  This file is a part of GetFEM
7 
8  GetFEM is free software; you can redistribute it and/or modify it
9  under the terms of the GNU Lesser General Public License as published
10  by the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version along with the GCC Runtime Library
12  Exception either version 3.1 or (at your option) any later version.
13  This program is distributed in the hope that it will be useful, but
14  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16  License and GCC Runtime Library Exception for more details.
17  You should have received a copy of the GNU Lesser General Public License
18  along with this program; if not, write to the Free Software Foundation,
19  Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
20 
21  As a special exception, you may use this file as it is a part of a free
22  software library without restriction. Specifically, if other files
23  instantiate templates or use macros or inline functions from this file,
24  or you compile this file and link it with other files to produce an
25  executable, this file does not by itself cause the resulting executable
26  to be covered by the GNU Lesser General Public License. This exception
27  does not however invalidate any other reasons why the executable file
28  might be covered by the GNU Lesser General Public License.
29 
30 ===========================================================================*/
31 
32 /**@file getfem_export.h
33  @author Yves Renard <[email protected]>,
34  @author Julien Pommier <[email protected]>
35  @date October 15, 2001.
36  @brief Export solutions to various formats.
37 */
38 #ifndef GETFEM_EXPORT_H__
39 #define GETFEM_EXPORT_H__
40 
41 #include "getfem_interpolation.h"
42 #include "getfem_mesh_slice.h"
43 #include <list>
44 
45 namespace getfem {
46 
47  /* ********************************************************************* */
48  /* */
49  /* Save a solution in a file with a Pk interpolation. */
50  /* */
51  /* ********************************************************************* */
52 
53  inline std::string remove_spaces(const std::string &s) {
54  std::string s2(s);
55  for (unsigned i=0; i < s.size(); ++i)
56  if (s2[i] <= ' ') s2[i] = '_';
57  return s2;
58  }
59 
60  /** @brief VTK/VTU export.
61 
62  export class to VTK/VTU file format
63  ( http://www.kitware.com/vtk.html,
64  legacy and serial vtkUnstructuredGrid)
65 
66  A vtk_export can store multiple scalar/vector fields.
67  */
68  class vtk_export {
69  protected:
70  std::ostream &os;
71  char header[256]; // hard limit in vtk/vtu
72  bool ascii;
73  bool vtk; // true for vtk, false for vtu
74  const stored_mesh_slice *psl;
75  std::unique_ptr<mesh_fem> pmf;
76  dal::bit_vector pmf_dof_used;
77  std::vector<unsigned> pmf_mapping_type;
78  std::ofstream real_os;
79  dim_type dim_;
80  bool reverse_endian;
81  std::vector<unsigned char> vals;
82  enum { EMPTY, HEADER_WRITTEN, STRUCTURE_WRITTEN, IN_CELL_DATA,
83  IN_POINT_DATA } state;
84 
85  template<class T> void write_val(T v);
86  template<class V> void write_vec(V p, size_type qdim);
87  template<class IT> void write_3x3tensor(IT p);
88  void write_separ();
89  void clear_vals();
90  void write_vals();
91 
92  public:
93  vtk_export(const std::string& fname, bool ascii_ = false, bool vtk_= true);
94  vtk_export(std::ostream &os_, bool ascii_ = false, bool vtk_ = true);
95  ~vtk_export(); // a vtu file is completed upon calling the destructor
96  /** should be called before write_*_data */
97  void exporting(const mesh& m);
98  void exporting(const mesh_fem& mf);
99  void exporting(const stored_mesh_slice& sl);
100 
101  /** the header is the second line of text in the exported file,
102  you can put whatever you want -- call this before any write_dataset
103  or write_mesh */
104  void set_header(const std::string& s);
105  void write_mesh();
106 
107  /** append a new scalar or vector field defined on mf to the .vtk/.vtu file.
108  If you are exporting a slice, or if mf != get_exported_mesh_fem(), U
109  will be interpolated on the slice, or on get_exported_mesh_fem().
110 
111  Note that vectors should be written AFTER scalars, and tensors
112  after vectors
113 
114  NO SPACE ALLOWED in 'name' */
115  template<class VECT> void write_point_data(const getfem::mesh_fem &mf,
116  const VECT& U0,
117  const std::string& name);
118 
119  /** append a new scalar or vector field to .vtk file. The Uslice vector is
120  the field interpolated on the exported mesh_slice This function should
121  not be used if you are not exporting a slice! NO SPACE ALLOWED in
122  'name' */
123  template<class VECT> void write_sliced_point_data(const VECT& Uslice,
124  const std::string& name,
125  size_type qdim=1);
126  /** export data which is constant over each element. You should not use
127  this function if you are exporting a slice. U should have
128  convex_index().card() elements. */
129 
130  template<class VECT> void write_cell_data(const VECT& U,
131  const std::string& name,
132  size_type qdim = 1);
133  /** export a data_set correspounding to measures of quality for each convex
134  of the supplied mesh (which should have the same number of convex than
135  the one used in the vtk_export)
136 
137  If a slice is being exported, the convex quality is written as
138  point_data (TO IMPROVE ONEDAY), if a mesh/mesh_fem is being exported,
139  it is written as cell_data
140  */
141  void write_mesh_quality(const mesh &m);
142  void write_normals();
143  const stored_mesh_slice& get_exported_slice() const;
144  const mesh_fem& get_exported_mesh_fem() const;
145  private:
146  void init();
147  void check_header();
148  void write_mesh_structure_from_slice();
149  void write_mesh_structure_from_mesh_fem();
150  void switch_to_cell_data();
151  void switch_to_point_data();
152  template<class VECT> void write_dataset_(const VECT& U,
153  const std::string& name,
154  size_type qdim,
155  bool cell_data=false);
156  };
157 
158  template<class T> void vtk_export::write_val(T v) {
159  if (ascii) os << " " << v;
160  else {
161  char *p = (char*)&v;
162  if (vtk) {
163  if (reverse_endian)
164  for (size_type i=0; i < sizeof(v)/2; ++i)
165  std::swap(p[i], p[sizeof(v)-i-1]);
166  os.write(p, sizeof(T));
167  } else {
168  union { T value; unsigned char bytes[sizeof(T)]; } UNION;
169  UNION.value = v;
170  for (size_type i=0; i < sizeof(T); i++)
171  vals.push_back(UNION.bytes[i]);
172  }
173  }
174  }
175 
176  template<class IT> void vtk_export::write_vec(IT p, size_type qdim) {
177  float v[3];
178  for (size_type i=0; i < qdim; ++i) {
179  v[i] = float(p[i]);
180  }
181  for (size_type i=qdim; i < 3; ++i) v[i] = 0.0f;
182  write_val(v[0]);write_val(v[1]);write_val(v[2]);
183  }
184 
185  template<class IT> void vtk_export::write_3x3tensor(IT p) {
186  float v[3][3];
187  memset(v, 0, sizeof v);
188  for (size_type i=0; i < dim_; ++i) {
189  for (size_type j=0; j < dim_; ++j)
190  v[i][j] = float(p[i + j*dim_]);
191  }
192  for (size_type i=0; i < 3; ++i) {
193  for (size_type j=0; j < 3; ++j) {
194  write_val(v[i][j]);
195  }
196  if (ascii) os << "\n";
197  }
198  }
199 
200  template<class VECT>
201  void vtk_export::write_point_data(const getfem::mesh_fem &mf, const VECT& U,
202  const std::string& name) {
203  size_type Q = (gmm::vect_size(U) / mf.nb_dof()) * mf.get_qdim();
204  size_type qdim = mf.get_qdim();
205  if (psl) {
206  std::vector<scalar_type> Uslice(Q*psl->nb_points());
207  psl->interpolate(mf, U, Uslice);
208  write_dataset_(Uslice, name, qdim);
209  } else {
210  std::vector<scalar_type> V(pmf->nb_dof() * Q);
211  if (&mf != &(*pmf)) {
212  interpolation(mf, *pmf, U, V);
213  } else gmm::copy(U,V);
214  size_type cnt = 0;
215  for (dal::bv_visitor d(pmf_dof_used); !d.finished(); ++d, ++cnt) {
216  if (cnt != d)
217  for (size_type q=0; q < Q; ++q) {
218  V[cnt*Q + q] = V[d*Q + q];
219  }
220  }
221  V.resize(Q*pmf_dof_used.card());
222  write_dataset_(V, name, qdim);
223  }
224  }
225 
226  template<class VECT>
227  void vtk_export::write_cell_data(const VECT& U, const std::string& name,
228  size_type qdim) {
229  write_dataset_(U, name, qdim, true);
230  }
231 
232  template<class VECT>
234  const std::string& name,
235  size_type qdim) {
236  write_dataset_(U, name, qdim, false);
237  }
238 
239  template<class VECT>
240  void vtk_export::write_dataset_(const VECT& U, const std::string& name,
241  size_type qdim, bool cell_data) {
242  write_mesh();
243  size_type nb_val = 0;
244  if (cell_data) {
245  switch_to_cell_data();
246  nb_val = psl ? psl->linked_mesh().convex_index().card()
247  : pmf->linked_mesh().convex_index().card();
248  } else {
249  switch_to_point_data();
250  nb_val = psl ? psl->nb_points() : pmf_dof_used.card();
251  }
252  size_type Q = qdim;
253  if (Q == 1) Q = gmm::vect_size(U) / nb_val;
254  GMM_ASSERT1(gmm::vect_size(U) == nb_val*Q,
255  "inconsistency in the size of the dataset: "
256  << gmm::vect_size(U) << " != " << nb_val << "*" << Q);
257  if (vtk) write_separ();
258  if (!vtk && !ascii) write_val(float(gmm::vect_size(U)));
259  if (Q == 1) {
260  if (vtk)
261  os << "SCALARS " << remove_spaces(name) << " float 1\n"
262  << "LOOKUP_TABLE default\n";
263  else
264  os << "<DataArray type=\"Float32\" Name=\"" << remove_spaces(name) << "\" "
265  << (ascii ? "format=\"ascii\">\n" : "format=\"binary\">\n");
266  for (size_type i=0; i < nb_val; ++i)
267  write_val(float(U[i]));
268  } else if (Q <= 3) {
269  if (vtk)
270  os << "VECTORS " << remove_spaces(name) << " float\n";
271  else
272  os << "<DataArray type=\"Float32\" Name=\"" << remove_spaces(name) << "\" "
273  << "NumberOfComponents=\"3\" "
274  << (ascii ? "format=\"ascii\">\n" : "format=\"binary\">\n");
275  for (size_type i=0; i < nb_val; ++i)
276  write_vec(U.begin() + i*Q, Q);
277  } else if (Q == gmm::sqr(dim_)) {
278  /* tensors : coef are supposed to be stored in FORTRAN order
279  in the VTK/VTU file, they are written with C (row major) order
280  */
281  if (vtk)
282  os << "TENSORS " << remove_spaces(name) << " float\n";
283  else
284  os << "<DataArray type=\"Float32\" Name=\"" << remove_spaces(name)
285  << "\" NumberOfComponents=\"9\" "
286  << (ascii ? "format=\"ascii\">\n" : "format=\"binary\">\n");
287  for (size_type i=0; i < nb_val; ++i)
288  write_3x3tensor(U.begin() + i*Q);
289  } else
290  GMM_ASSERT1(false, std::string(vtk ? "vtk" : "vtu")
291  + " does not accept vectors of dimension > 3");
292  write_vals();
293  if (vtk) write_separ();
294  if (!vtk) os << "\n" << "</DataArray>\n";
295  }
296 
297 
298  class vtu_export : public vtk_export {
299  public:
300  vtu_export(const std::string& fname, bool ascii_ = false) : vtk_export(fname, ascii_, false) {}
301  vtu_export(std::ostream &os_, bool ascii_ = false) : vtk_export(os_, ascii_, false) {}
302  };
303 
304  /** @brief A (quite large) class for exportation of data to IBM OpenDX.
305 
306  http://www.opendx.org/
307 
308  This class is more capable than the VTK export, as it is
309  possible to export many different meshes/slices, with their
310  edges, datasets, and create series of dataset for animations
311  etc, in a single '.dx' file.
312 
313  Moreover, it is able to reopen a '.dx' file and append new data
314  into it. Hence it is possible, if many time-steps are to be
315  saved, to view intermediate results in OpenDX during the
316  computation.
317  */
318  class dx_export {
319  std::ostream &os;
320  char header[256];
321  bool ascii;
322  const stored_mesh_slice *psl;
323  bool psl_use_merged; /* flag enabled if we merge the points of
324  psl before export */
325  std::unique_ptr<mesh_fem> pmf;
326  dal::bit_vector pmf_dof_used;
327  std::vector<unsigned> pmf_cell_type;
328  std::fstream real_os;
329  dim_type dim_, connections_dim;
330  struct dxSeries {
331  std::string name;
332  std::list<std::string> members;
333  };
334  struct dxObject {
335  std::string name;
336  std::string mesh;
337  };
338  struct dxMesh {
339  unsigned flags;
340  typedef enum { NONE=0, WITH_EDGES=1, STRUCTURE_WRITTEN=2 } flags_t;
341  std::string name;
342  dxMesh() : flags(NONE) {}
343  };
344  std::list<dxObject> objects;
345  std::list<dxSeries> series;
346  std::list<dxMesh> meshes;
347  bool header_written;
348  public:
349  dx_export(const std::string& fname, bool ascii_ = false,
350  bool append_ = false);
351  dx_export(std::ostream &os_, bool ascii_ = false);
352  ~dx_export(); /* the file is not complete until the destructor
353  has been executed */
354  void exporting(const mesh& m, std::string name = std::string());
355  void exporting(const mesh_fem& mf, std::string name = std::string());
356  void exporting(const stored_mesh_slice& sl, bool merge_points = true,
357  std::string name = std::string());
358  /** append edges information (if you want to draw the mesh and are
359  using a refined slice. Should be called just after exporting(..) */
360  void exporting_mesh_edges(bool with_slice_edge = true);
361 
362  /** the header is the second line of text in the exported file,
363  you can put whatever you want -- call this before any write_dataset
364  or write_mesh */
365  void set_header(const std::string& s);
366  void write_mesh();
367  /** add an object (typically the name of a data field) to a
368  'series', i.e. an aggregate of consecutive objects. Using
369  'series' is useful for animations in opendx
370 
371  If 'field_name' corresponds to a data_set whose mesh edges have
372  been exported, a second series called serie_name + '_edges'
373  will be filled, which will allow you to view the mesh edges.
374  */
375  void serie_add_object(const std::string& serie_name,
376  const std::string& object_name);
377  void serie_add_object(const std::string& serie_name)
378  { serie_add_object(serie_name, current_data_name()); }
379  /** return the name of current mesh (use exporting(...) to change
380  the current mesh) */
381  std::string current_mesh_name() { return current_mesh().name; }
382  /** return the name of last written data_set */
383  std::string current_data_name() { return current_data().name; }
384  template<class VECT> void
385  write_point_data(const getfem::mesh_fem &mf,
386  const VECT& U0, std::string name = std::string());
387  template<class VECT> void
388  write_sliced_point_data(const VECT& Uslice,
389  std::string name = std::string());
390  /* TOBEDONE !!!!!!!!!!!
391  template<class VECT> void
392  write_cell_data(const VECT& U, std::string name = std::string());
393  void write_mesh_quality(const mesh &m);*/
394  void write_normals();
395  const stored_mesh_slice& get_exported_slice() const;
396  const mesh_fem& get_exported_mesh_fem() const;
397 
398  private:
399  void init();
400  void reread_metadata();
401  void update_metadata(std::ios::pos_type);
402  void write_series();
403  void serie_add_object_(const std::string &serie_name,
404  const std::string &object_name);
405  void write_separ();
406  std::string default_name(std::string s, int count,
407  const char *default_prefix) {
408  if (s.size() == 0) {
409  std::stringstream ss; ss << default_prefix << count; return ss.str();
410  } else return s;
411  }
412  template<class T> void write_val(T v) {
413  if (ascii) os << " " << v;
414  else os.write((char*)&v, sizeof(T));
415  }
416  static const char* endianness() {
417  static int i=0x12345678;
418  char *p = (char*)&i;
419  if (*p == 0x12) return "msb";
420  else if (*p == 0x78) return "lsb";
421  else return "this is very strange..";
422  }
423  bool new_mesh(std::string &name);
424  std::list<dxMesh>::iterator get_mesh(const std::string& name,
425  bool raise_error = true);
426  std::list<dxObject>::iterator get_object(const std::string& name,
427  bool raise_error = true);
428  dxMesh &current_mesh() {
429  if (meshes.size()) return meshes.back();
430  else GMM_ASSERT1(false, "no mesh!");
431  }
432  dxObject &current_data() {
433  if (objects.size()) return objects.back();
434  else GMM_ASSERT1(false, "no data!");
435  }
436 
437  std::string name_of_pts_array(const std::string &meshname)
438  { return meshname + std::string("_pts"); }
439  std::string name_of_conn_array(const std::string &meshname)
440  { return meshname + std::string("_conn"); }
441  std::string name_of_edges_array(const std::string &meshname)
442  { return meshname + std::string("_edges"); }
443  void check_header();
444  const char *dxname_of_convex_structure(bgeot::pconvex_structure cvs);
445  void write_convex_attributes(bgeot::pconvex_structure cvs);
446  void write_mesh_structure_from_slice();
447  void write_mesh_structure_from_mesh_fem();
448  void write_mesh_edges_from_slice(bool with_slice_edge);
449  void write_mesh_edges_from_mesh();
450  template <class VECT>
451  void smooth_field(const VECT& U, base_vector &sU);
452  template<class VECT>
453  void write_dataset_(const VECT& U, std::string name, bool cell_data=false);
454  };
455 
456  template <class VECT>
457  void dx_export::smooth_field(const VECT& U, base_vector &sU) {
458  size_type Q = gmm::vect_size(U) / psl->nb_points();
459  sU.clear(); sU.resize(Q*psl->nb_merged_nodes());
460  for (size_type i=0; i < psl->nb_merged_nodes(); ++i) {
461  for (size_type j=0; j < psl->merged_point_cnt(i); ++j)
462  for (size_type q=0; q < Q; ++q)
463  sU[i*Q+q] += U[psl->merged_point_nodes(i)[j].pos*Q+q];
464  for (size_type q=0; q < Q; ++q)
465  sU[i*Q+q] /= double(psl->merged_point_cnt(i));
466  }
467  }
468 
469  template<class VECT>
470  void dx_export::write_point_data(const getfem::mesh_fem &mf, const VECT& U,
471  std::string name) {
472  size_type Q = (gmm::vect_size(U) / mf.nb_dof())*mf.get_qdim();
473  if (psl) {
474  std::vector<scalar_type> Uslice(Q*psl->nb_points());
475  psl->interpolate(mf, U, Uslice);
476  write_sliced_point_data(Uslice,name);
477  } else {
478  std::vector<scalar_type> V(pmf->nb_dof() * Q);
479  if (&mf != &(*pmf)) {
480  interpolation(mf, *pmf, U, V);
481  } else gmm::copy(U,V);
482  size_type cnt = 0;
483  for (dal::bv_visitor d(pmf_dof_used); !d.finished(); ++d, ++cnt) {
484  if (cnt != d)
485  for (size_type q=0; q < Q; ++q) {
486  V[cnt*Q + q] = V[d*Q + q];
487  }
488  }
489  V.resize(Q*pmf_dof_used.card());
490  write_dataset_(V, name);
491  }
492  }
493 
494  template<class VECT> void
495  dx_export::write_sliced_point_data(const VECT& Uslice, std::string name) {
496  if (!psl_use_merged)
497  write_dataset_(Uslice, name, false);
498  else {
499  base_vector Umerged; smooth_field(Uslice,Umerged);
500  write_dataset_(Umerged, name, false);
501  }
502  }
503 
504  template<class VECT> void
505  dx_export::write_dataset_(const VECT& U, std::string name, bool cell_data) {
506  write_mesh();
507  objects.push_back(dxObject());
508  name = default_name(name, int(objects.size()), "gf_field");
509  objects.back().name = name;
510  objects.back().mesh = current_mesh_name();
511  size_type nb_val = 0;
512  if (cell_data) {
513  nb_val = psl ? psl->linked_mesh().convex_index().card()
514  : pmf->linked_mesh().convex_index().card();
515  } else {
516  nb_val = psl ? (psl_use_merged ? psl->nb_merged_nodes() : psl->nb_points())
517  : pmf_dof_used.card();
518  }
519  size_type Q = gmm::vect_size(U) / nb_val;
520  GMM_ASSERT1(gmm::vect_size(U) == nb_val*Q,
521  "inconsistency in the size of the dataset: "
522  << gmm::vect_size(U) << " != " << nb_val << "*" << Q);
523 
524  os << "\nobject \"" << name << "_data\" class array type float rank ";
525  if (Q == 1) { os << "0"; } /* scalar data */
526  else if (Q == 4) { os << "2 shape 2 2"; } /* or 2x2 tensor data */
527  else if (Q == 9) { os << "2 shape 3 3"; } /* or 2x2 tensor data */
528  else { os << "1 shape " << Q; } /* fallback: vector data */
529  os << " items " << nb_val;
530  if (!ascii) os << " " << endianness() << " binary";
531  os << " data follows" << endl;
532  for (size_type i=0; i < nb_val*Q; ++i) {
533  write_val(float(U[i]));
534  if (((i+1) % (Q > 1 ? Q : 10)) == 0) write_separ();
535  }
536  write_separ();
537 
538  if (!cell_data)
539  os << "\n attribute \"dep\" string \"positions\"\n";
540  else os << "\n attribute \"dep\" string \"connections\"\n";
541  os << "\n";
542 
543  if (current_mesh().flags & dxMesh::WITH_EDGES) {
544  os << "\nobject \"" << name << "_edges\" class field\n"
545  << " component \"positions\" value \""
546  << name_of_pts_array(current_mesh_name()) << "\"\n"
547  << " component \"connections\" value \""
548  << name_of_conn_array(name_of_edges_array(current_mesh_name()))
549  << "\"\n"
550  << " component \"data\" value \"" << name << "_data\"\n";
551  }
552 
553  /* write footer */
554  os << "\nobject \"" << name << "\" class field\n"
555  << " component \"positions\" value \""
556  << name_of_pts_array(current_mesh_name()) << "\"\n"
557  << " component \"connections\" value \""
558  << name_of_conn_array(current_mesh_name()) << "\"\n"
559  << " component \"data\" value \"" << name << "_data\"\n";
560  }
561 
562  /** @brief POS export.
563 
564  export class to Gmsh post-processing file format.
565 
566  ( http://geuz.org/gmsh )
567 
568  A pos_export can store multiple scalar/vector/tensor fields.
569  */
570 
571  class pos_export {
572  protected:
573  std::ostream& os;
574  char header[256];
575 
576  std::vector<std::vector<float> > pos_pts;
577  std::vector<unsigned> pos_cell_type;
578  std::vector<std::vector<unsigned> > pos_cell_dof;
579 
580  std::unique_ptr<mesh_fem> pmf;
581  const stored_mesh_slice *psl;
582 
583  size_type view;
584  dim_type dim;
585  enum { EMPTY, HEADER_WRITTEN, STRUCTURE_WRITTEN, IN_CELL_DATA } state;
586  std::ofstream real_os;
587 
588  public:
589  typedef enum {
590  POS_PT = 0, // point
591  POS_LN = 1, // line
592  POS_TR = 2, // triangles
593  POS_QU = 3, // quadrangles
594  POS_SI = 4, // tetrahedra
595  POS_HE = 5, // hexahedra
596  POS_PR = 6, // prisms
597  POS_PY = 7 // pyramids
598  } pos_cell_types;
599 
600  pos_export(const std::string& fname);
601  pos_export(std::ostream& osname);
602 
603  void set_header(const std::string& s);
604 
605  void exporting(const mesh& m);
606  void exporting(const mesh_fem& mf);
607  void exporting(const stored_mesh_slice& sl);
608 
609  void write(const mesh& m, const std::string& name="");
610  void write(const mesh_fem& mf, const std::string& name="");
611  void write(const stored_mesh_slice& sl, const std::string& name="");
612 
613  template <class VECT>
614  void write(const mesh_fem& mf,const VECT& U, const std::string& name);
615  template <class VECT>
616  void write(const stored_mesh_slice& sl,const VECT& U, const std::string& name);
617 
618  private:
619  void init();
620  void check_header();
621 
622  template <class VECT>
623  void write(const VECT& V, const size_type qdim_v);
624 
625  template <class VECT>
626  void write_cell(const int& t, const std::vector<unsigned>& dof,
627  const VECT& val);
628  };
629 
630  template <class VECT>
631  void pos_export::write(const mesh_fem& mf,const VECT& U,
632  const std::string& name){
633  check_header();
634  exporting(mf);
635 
636  os << "View \"" << name.c_str() <<"\" {\n";
637 
638  size_type nb_points = mf.nb_dof()/mf.get_qdim();
639  size_type qdim_u = gmm::vect_size(U)/nb_points;
640  if (psl){
641  std::vector<scalar_type> Uslice(psl->nb_points()*qdim_u);
642  psl->interpolate(mf, U, Uslice);
643  qdim_u = gmm::vect_size(Uslice)/psl->nb_points();
644  write(Uslice, qdim_u);
645  }else {
646  std::vector<scalar_type> V(pmf->nb_dof()*qdim_u);
647  if (&mf != &(*pmf)) {
648  interpolation(mf, *pmf, U, V);
649  } else gmm::copy(U,V);
650  /*for (dal::bv_visitor d(pmf_dof_used); !d.finished(); ++d, ++cnt) {
651  if (cnt != d)
652  for (size_type q=0; q < Q; ++q) {
653  V[cnt*Q + q] = V[d*Q + q];
654  }
655  }
656  V.resize(Q*pmf_dof_used.card());*/
657  nb_points = pmf->nb_dof()/pmf->get_qdim();
658  qdim_u = gmm::vect_size(V)/nb_points;
659  write(V, qdim_u);
660  }
661 
662  os << "};\n";
663  os << "View[" << view << "].ShowScale = 1;\n";
664  os << "View[" << view << "].ShowElement = 0;\n";
665  os << "View[" << view << "].DrawScalars = 1;\n";
666  os << "View[" << view << "].DrawVectors = 1;\n";
667  os << "View[" << view++ << "].DrawTensors = 1;\n";
668  }
669 
670  template <class VECT>
671  void pos_export::write(const stored_mesh_slice& sl,const VECT& V,
672  const std::string& name){
673  check_header();
674  exporting(sl);
675 
676  os << "View \"" << name.c_str() <<"\" {\n";
677 
678  size_type qdim_v = gmm::vect_size(V)/psl->nb_points();
679  write(V, qdim_v);
680 
681  os << "};\n";
682  os << "View[" << view << "].ShowScale = 1;\n";
683  os << "View[" << view << "].ShowElement = 0;\n";
684  os << "View[" << view << "].DrawScalars = 1;\n";
685  os << "View[" << view << "].DrawVectors = 1;\n";
686  os << "View[" << view++ << "].DrawTensors = 1;\n";
687  }
688 
689  template <class VECT>
690  void pos_export::write(const VECT& V, const size_type qdim_v) {
691  int t;
692  std::vector<unsigned> cell_dof;
693  std::vector<scalar_type> cell_dof_val;
694  for (size_type cell = 0; cell < pos_cell_type.size(); ++cell) {
695  t = pos_cell_type[cell];
696  cell_dof = pos_cell_dof[cell];
697  cell_dof_val.resize(cell_dof.size()*qdim_v, scalar_type(0));
698  for (size_type i=0; i< cell_dof.size(); ++i)
699  for (size_type j=0; j< qdim_v; ++j)
700  cell_dof_val[i*qdim_v+j] = scalar_type(V[cell_dof[i]*qdim_v+j]);
701  write_cell(t,cell_dof,cell_dof_val);
702  }
703  }
704 
705  template <class VECT>
706  void pos_export::write_cell(const int& t, const std::vector<unsigned>& dof,
707  const VECT& val) {
708  size_type qdim_cell = val.size()/dof.size();
709  size_type dim3D = size_type(-1);
710  if (1==qdim_cell){
711  dim3D = size_type(1);
712  os << "S";
713  } else if (2==qdim_cell || 3==qdim_cell){
714  dim3D = size_type(3);
715  os << "V";
716  } else if (4<=qdim_cell && qdim_cell<=9){
717  dim3D = size_type(9);
718  os << "T";
719  }
720  switch (t){
721  case POS_PT: os << "P("; break; // point
722  case POS_LN: os << "L("; break; // line
723  case POS_TR: os << "T("; break; // triangle
724  case POS_QU: os << "Q("; break; // quadrangle
725  case POS_SI: os << "S("; break; // tetrahedra (simplex)
726  case POS_HE: os << "H("; break; // hexahedra
727  case POS_PR: os << "I("; break; // prism
728  case POS_PY: os << "Y("; break; // pyramid
729  }
730  for (size_type i=0; i<dof.size(); ++i){
731  for(size_type j=0; j<dim; ++j){
732  if(0!=i || 0!=j) os << ",";
733  os << pos_pts[dof[i]][j];
734  }
735  for (size_type j=dim; j<3; ++j){
736  os << ",0.00";
737  }
738  }
739 
740  os << "){";
741  for (size_type i=0; i<dof.size(); ++i){
742  for(size_type j=0; j<qdim_cell; ++j){
743  if(0!=i || 0!=j) os << ",";
744  os << val[i*qdim_cell+j];
745  }
746  for (size_type j=qdim_cell; j< dim3D; ++j){
747  os << ",0.00";
748  }
749  }
750  os << "};\n";
751  }
752 } /* end of namespace getfem. */
753 
754 #endif /* GETFEM_EXPORT_H__ */
getfem::vtk_export::exporting
void exporting(const mesh &m)
should be called before write_*_data
Definition: getfem_export.cc:269
getfem::mesh_fem::get_qdim
virtual dim_type get_qdim() const
Return the Q dimension.
Definition: getfem_mesh_fem.h:312
getfem::vtk_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:387
getfem::dx_export::current_data_name
std::string current_data_name()
return the name of last written data_set
Definition: getfem_export.h:383
getfem::stored_mesh_slice
The output of a getfem::mesh_slicer which has been recorded.
Definition: getfem_mesh_slice.h:47
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
bgeot::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:49
getfem::vtk_export::write_sliced_point_data
void write_sliced_point_data(const VECT &Uslice, const std::string &name, size_type qdim=1)
append a new scalar or vector field to .vtk file.
Definition: getfem_export.h:233
getfem::mesh_fem
Describe a finite element method linked to a mesh.
Definition: getfem_mesh_fem.h:148
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
getfem::interpolation
void interpolation(const mesh_fem &mf_source, const mesh_fem &mf_target, const VECTU &U, VECTV &V, int extrapolation=0, double EPS=1E-10, mesh_region rg_source=mesh_region::all_convexes(), mesh_region rg_target=mesh_region::all_convexes())
interpolation/extrapolation of (mf_source, U) on mf_target.
Definition: getfem_interpolation.h:693
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::vtk_export::write_point_data
void write_point_data(const getfem::mesh_fem &mf, const VECT &U0, const std::string &name)
append a new scalar or vector field defined on mf to the .vtk/.vtu file.
Definition: getfem_export.h:201
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::stored_mesh_slice::interpolate
void interpolate(const getfem::mesh_fem &mf, const V1 &UU, V2 &V) const
Interpolation of a mesh_fem on a slice.
Definition: getfem_mesh_slice.h:245
getfem_mesh_slice.h
Define the class getfem::stored_mesh_slice.
getfem::pos_export
POS export.
Definition: getfem_export.h:571
getfem::vtk_export::write_cell_data
void write_cell_data(const VECT &U, const std::string &name, size_type qdim=1)
export data which is constant over each element.
Definition: getfem_export.h:227
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_interpolation.h
Interpolation of fields from a mesh_fem onto another.
getfem::vtk_export::write_mesh_quality
void write_mesh_quality(const mesh &m)
export a data_set correspounding to measures of quality for each convex of the supplied mesh (which s...
Definition: getfem_export.cc:636
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::merged_point_cnt
size_type merged_point_cnt(size_type i_merged) const
Return the number of nodes that were merged to form the merged one.
Definition: getfem_mesh_slice.h:168
getfem::stored_mesh_slice::linked_mesh
const mesh & linked_mesh() const
return a pointer to the original mesh
Definition: getfem_mesh_slice.h:110
getfem::dx_export
A (quite large) class for exportation of data to IBM OpenDX.
Definition: getfem_export.h:318
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::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
getfem::vtk_export
VTK/VTU export.
Definition: getfem_export.h:68