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