[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Getfem-commits] (no subject)
|
From: |
Tetsuo Koyama |
|
Subject: |
[Getfem-commits] (no subject) |
|
Date: |
Wed, 13 May 2020 08:49:15 -0400 (EDT) |
branch: devel-tetsuo-xml02
commit 17ac426352b171ee5cf2840d59b6b9601afac0b4
Author: Tetsuo Koyama <address@hidden>
AuthorDate: Wed May 13 21:46:11 2020 +0900
add class vtu_export class
---
doc/sphinx/source/replaces.txt | 3 +
doc/sphinx/source/userdoc/export.rst | 54 ++--
src/getfem/getfem_export.h | 172 ++++++++++++-
src/getfem_export.cc | 334 ++++++++++++++++++-------
tests/thermo_elasticity_electrical_coupling.cc | 47 +++-
5 files changed, 481 insertions(+), 129 deletions(-)
diff --git a/doc/sphinx/source/replaces.txt b/doc/sphinx/source/replaces.txt
index 5710048..00d83e8 100644
--- a/doc/sphinx/source/replaces.txt
+++ b/doc/sphinx/source/replaces.txt
@@ -2,6 +2,7 @@
.. |gnu| replace:: *GNU*
.. |c++| replace:: *C++*
.. |vtk| replace:: *VTK*
+.. |vtu| replace:: *VTU*
.. |opendx| replace:: *OpenDX*
.. |gmsh| replace:: *Gmsh*
.. |emc2| replace:: *emc2*
@@ -48,6 +49,7 @@
.. |gf_sl_a| replace:: ``getfem::slicer_action``
.. |gf_sl_ddb| replace:: ``getfem::mesh_slice_cv_dof_data_base``
.. |gf_vtk_export| replace:: ``getfem::vtk_export``
+.. |gf_vtu_export| replace:: ``getfem::vtu_export``
.. |gf_dx_export| replace:: ``getfem::dx_export``
.. |gf_pos_export| replace:: ``getfem::pos_export``
.. |gf_gasm| replace:: ``getfem::generic_assembly``
@@ -85,6 +87,7 @@
.. _user documentation: http://getfem.org/userdoc/index.html
.. _vocabulary: http://getfem.org/getfem_reference/index.html
.. _VTK: http://www.vtk.org
+.. _VTU: https://vtk.org/Wiki/VTK_XML_Formats
.. _MayaVi: http://mayavi.sourceforge.net
.. _OpenDX: http://www.opendx.org
.. _Gmsh: http://www.geuz.org/gmsh
diff --git a/doc/sphinx/source/userdoc/export.rst
b/doc/sphinx/source/userdoc/export.rst
index 945b656..4f16252 100644
--- a/doc/sphinx/source/userdoc/export.rst
+++ b/doc/sphinx/source/userdoc/export.rst
@@ -12,7 +12,7 @@ Export and view a solution
There are essentially four ways to view the result of getfem computations:
* Scilab, Octave or Matlab, with the interface.
-* The open-source Paraview or Mayavi or any other VTK files viewer.
+* The open-source Paraview or Mayavi or any other VTK (or VTU) files viewer.
* The open-source OpenDX program.
* The open-source Gmsh program.
@@ -43,7 +43,7 @@ and then, under Scilab, Octave or Matlab:
See the getfem-matlab interface documentation for more details.
-Two other file formats are supported for export: the `VTK`_ file format, the
+Two other file formats are supported for export: the `VTK`_ (and `VTU`_ )file
format, the
`OpenDX`_ file format and the `Gmsh`_ post-processing file format. Both can
export
either a |gf_m| or |gf_mf| , but also the more versatile |gf_smsl|.
@@ -178,20 +178,20 @@ In order to build a |gf_smsl| object during the slicing
operation, the ``stored_
getfem::slicer_half_space(base_node(0,0), base_node(1, 0), -1),
nrefine);
-The simplest way to use these slices is to export them to |vtk|, |opendx|, or
-|gmsh|. The file :file:`getfem/getfem_export.h` contains three classes:
-|gf_vtk_export|, |gf_dx_export| and |gf_pos_export|.
+The simplest way to use these slices is to export them to |vtk|, |vtu|,
|opendx|, or
+|gmsh|. The file :file:`getfem/getfem_export.h` contains four classes:
+|gf_vtk_export|, |gf_vtu_export|, |gf_dx_export| and |gf_pos_export|.
-Exporting |m|, |mf| or slices to VTK
-------------------------------------
+Exporting |m|, |mf| or slices to VTK (or VTU)
+----------------------------------------------
-First, it is important to know the limitation of VTK data files: each file can
+First, it is important to know the limitation of VTK (or VTU) data files: each
file can
contain only one mesh, with at most one scalar field and one vector field and
one
-tensor field on this mesh (in that order). VTK files can handle data on
segment,
+tensor field on this mesh (in that order). VTK (or VTU) files can handle data
on segment,
triangles, quadrangles, tetrahedrons and hexahedrons. Although quadratic
triangles, segments etc are said to be supported, it is just equivalent to
using
-``nrefine=2`` when building a slice. VTK data file do support meshes with more
+``nrefine=2`` when building a slice. VTK (or VTU) data file do support meshes
with more
than one type of element (i.e. meshes with triangles and quadrangles, for
example).
@@ -199,26 +199,38 @@ For example, supposing that a |smsl| ``sl`` has already
been built::
// an optional the 2nd argument can be set to true to produce
// a text file instead of a binary file
- vtk_export exp("output.vtk");
- exp.exporting(sl); // will save the geometrical structure of the slice
- exp.write_point_data(mfp, P, "pressure"); // write a scalar field
- exp.write_point_data(mfu, U, "displacement"); // write a vector field
+
+ vtk_export vtk_exp("output.vtk");
+ vtk_exp.exporting(sl); // will save the geometrical structure of the slice
+ vtk_exp.write_point_data(mfp, P, "pressure"); // write a scalar field
+ vtk_exp.write_point_data(mfu, U, "displacement"); // write a vector field
+
+ vtu_export vtu_exp("output.vtu", true);
+ vtu_exp.exporting(sl); // will save the geometrical structure of the slice
+ vtu_exp.write_point_data(mfp, P, "pressure"); // write a scalar field
+ vtu_exp.write_point_data(mfu, U, "displacement"); // write a vector field
In this example, the fields ``P`` and ``U`` are interpolated on the slice
nodes,
-and then written into the VTK field. The vector fields should always be written
+and then written into the VTK (or VTU) field. The vector fields should always
be written
after the scalar fields (and the tensor fields should be written last).
It is also possible to export a |mf| without having to build a slice::
// an optional the 2nd argument can be set to true to produce
// a text file instead of a binary file
- vtk_export exp("output.vtk");
- exp.exporting(mfu);
- exp.write_point_data(mfp, P, "pressure"); // write a scalar field
- exp.write_point_data(mfu, U, "displacement"); // write a vector field
-Note however that with this approach, the ``vtk_export`` will map each
convex/fem
-of ``mfu`` to a VTK element type. As VTK does not handle elements of degree
+ vtk_export vtk_exp("output.vtk");
+ vtk_exp.exporting(mfu);
+ vtk_exp.write_point_data(mfp, P, "pressure"); // write a scalar field
+ vtk_exp.write_point_data(mfu, U, "displacement"); // write a vector field
+
+ vtu_export vtu_exp("output.vtu", true);
+ vtu_exp.exporting(mfu);
+ vtu_exp.write_point_data(mfp, P, "pressure"); // write a scalar field
+ vtu_exp.write_point_data(mfu, U, "displacement"); // write a vector field
+
+Note however that with this approach, the ``vtk_export`` (or ``vtu_export``)
will map each convex/fem
+of ``mfu`` to a VTK (or VTU) element type. As VTK (or VTU) does not handle
elements of degree
greater than 2, there will be a loss of precision for higher degree FEMs.
Exporting |m|, |mf| or slices to OpenDX
diff --git a/src/getfem/getfem_export.h b/src/getfem/getfem_export.h
index de1a654..42b0328 100644
--- a/src/getfem/getfem_export.h
+++ b/src/getfem/getfem_export.h
@@ -57,6 +57,18 @@ namespace getfem {
return s2;
}
+ template<class VECT> VECT remove_dof_unused(VECT V, dal::bit_vector
dof_used, size_type Q) {
+ size_type cnt = 0;
+ for (dal::bv_visitor d(dof_used); !d.finished(); ++d, ++cnt ) {
+ if (cnt !=d)
+ for (size_type q=0; q < Q; ++q) {
+ V[cnt*Q + q] = V[d*Q + q];
+ }
+ }
+ V.resize(Q*dof_used.card());
+ return V;
+ }
+
/** @brief VTK export.
export class to VTK ( http://www.kitware.com/vtk.html ) file format
@@ -219,15 +231,9 @@ namespace getfem {
if (&mf != &(*pmf)) {
interpolation(mf, *pmf, U, V);
} else gmm::copy(U,V);
- size_type cnt = 0;
- for (dal::bv_visitor d(pmf_dof_used); !d.finished(); ++d, ++cnt) {
- if (cnt != d)
- for (size_type q=0; q < Q; ++q) {
- V[cnt*Q + q] = V[d*Q + q];
- }
- }
- V.resize(Q*pmf_dof_used.card());
- write_dataset_(V, name, qdim);
+ std::vector<scalar_type> W(Q*pmf_dof_used.card());
+ gmm::copy(remove_dof_unused(V, pmf_dof_used, Q), W);
+ write_dataset_(W, name, qdim);
}
}
@@ -286,6 +292,154 @@ namespace getfem {
write_separ();
}
+ /** @brief VTU export.
+
+ export class to VTU file format
+ (Serial vtkUnstructuredGrid (unstructured))
+ */
+ class vtu_export {
+ protected:
+ std::ostream &os;
+ char header[256]; // hard limit in vtu
+ bool ascii;
+ const stored_mesh_slice *psl;
+ std::unique_ptr<mesh_fem> pmf;
+ dal::bit_vector pmf_dof_used;
+ std::vector<unsigned> pmf_mapping_type;
+ std::ofstream real_os;
+ dim_type dim_;
+ enum { EMPTY, HEADER_WRITTEN, STRUCTURE_WRITTEN,
+ IN_POINT_DATA, FOOTER_WRITTEN } state;
+
+ template<class T> void write_val(T v);
+ template<class V> void write_vec(V p, size_type qdim);
+ template<class IT> void write_3x3tensor(IT p);
+ void write_separ();
+ public:
+ vtu_export(const std::string& fname, bool ascii_= false);
+ vtu_export(std::ostream &os_, bool ascii_ = false);
+ ~vtu_export(); /* the file is not complete until the destructor
+ has been executed */
+ void exporting(const mesh& m);
+ void exporting(const mesh_fem& mf);
+ void write_mesh(bool only_mesh = true);
+
+ /** append a new scalar or vector field defined on mf to the .vtu file. If
+ you are exporting a slice, or if mf != get_exported_mesh_fem(), U will
+ be interpolated on the slice, or on get_exported_mesh_fem().
+
+ Note that vectors should be written AFTER scalars, and tensors
+ after vectors
+
+ NO SPACE ALLOWED in 'name' */
+ template<class VECT> void write_point_data(const getfem::mesh_fem &mf,
+ const VECT& U0,
+ const std::string& name);
+
+ private:
+ void init();
+ void check_header();
+ void check_footer();
+ void write_mesh_structure_from_mesh_fem();
+ void switch_to_cell_data();
+ void switch_to_point_data();
+ template<class VECT> void write_dataset_(const VECT& U,
+ const std::string& name,
+ size_type qdim,
+ bool cell_data=false);
+ };
+
+ template<class T> void vtu_export::write_val(T v) {
+ os << " " << v;
+ }
+
+ template<class IT> void vtu_export::write_vec(IT p, size_type qdim) {
+ float v[3];
+ for (size_type i=0; i < qdim; ++i) {
+ v[i] = float(p[i]);
+ }
+ for (size_type i=qdim; i < 3; ++i) v[i] = 0.0f;
+ write_val(v[0]);write_val(v[1]);write_val(v[2]);
+ }
+
+ template<class IT> void vtu_export::write_3x3tensor(IT p) {
+ float v[3][3];
+ memset(v, 0, sizeof v);
+ for (size_type i=0; i < dim_; ++i) {
+ for (size_type j=0; j < dim_; ++j)
+ v[i][j] = float(p[i + j*dim_]);
+ }
+ for (size_type i=0; i < 3; ++i) {
+ for (size_type j=0; j < 3; ++j) {
+ write_val(v[i][j]);
+ }
+ if (ascii) os << "\n";
+ }
+ }
+
+ template<class VECT>
+ void vtu_export::write_point_data(const getfem::mesh_fem &mf, const VECT& U,
+ const std::string& name) {
+ size_type Q = (gmm::vect_size(U) / mf.nb_dof()) * mf.get_qdim();
+ size_type qdim = mf.get_qdim();
+ if (psl) {
+ std::vector<scalar_type> Uslice(Q*psl->nb_points());
+ psl->interpolate(mf, U, Uslice);
+ write_dataset_(Uslice, name, qdim);
+ } else {
+ std::vector<scalar_type> V(pmf->nb_dof() * Q);
+ if (&mf != &(*pmf)) {
+ interpolation(mf, *pmf, U, V);
+ } else gmm::copy(U,V);
+ std::vector<scalar_type> W(Q*pmf_dof_used.card());
+ gmm::copy(remove_dof_unused(V, pmf_dof_used, Q), W);
+ write_dataset_(W, name, qdim);
+ }
+ }
+
+ template<class VECT>
+ void vtu_export::write_dataset_(const VECT& U, const std::string& name,
+ size_type qdim, bool cell_data) {
+ write_mesh(false);
+ size_type nb_val = 0;
+ if (cell_data) {
+ switch_to_cell_data();
+ nb_val = psl ? psl->linked_mesh().convex_index().card()
+ : pmf->linked_mesh().convex_index().card();
+ } else {
+ switch_to_point_data();
+ nb_val = psl ? psl->nb_points() : pmf_dof_used.card();
+ }
+ size_type Q = qdim;
+ if (Q == 1) Q = gmm::vect_size(U) / nb_val;
+ GMM_ASSERT1(gmm::vect_size(U) == nb_val*Q,
+ "inconsistency in the size of the dataset: "
+ << gmm::vect_size(U) << " != " << nb_val << "*" << Q);
+ if (Q == 1) {
+ os << "<DataArray type=\"Float32\" Name=\"" << remove_spaces(name);
+ os << "\" format=\"ascii\">\n";
+ for (size_type i=0; i < nb_val; ++i) {
+ write_val(float(U[i]));
+ }
+ } else if (Q <= 3) {
+ os << "<DataArray type=\"Float32\" Name=\"" << remove_spaces(name);
+ os << "\" NumberOfComponents=\"3\" format=\"ascii\">\n";
+ for (size_type i=0; i < nb_val; ++i) {
+ write_vec(U.begin() + i*Q, Q);
+ }
+ } else if (Q == gmm::sqr(dim_)) {
+ /* tensors : coef are supposed to be stored in FORTRAN order
+ in the VTU file, they are written with C (row major) order
+ */
+ os << "<DataArray type=\"Float32\" Name=\"" << remove_spaces(name);
+ os << "\" NumberOfComponents=\"9\" format=\"ascii\">";
+ for (size_type i=0; i < nb_val; ++i) {
+ write_3x3tensor(U.begin() + i*Q);
+ }
+ } else GMM_ASSERT1(false, "vtu does not accept vectors of dimension > 3");
+ write_separ();
+ os << "</DataArray>\n";
+ }
/** @brief A (quite large) class for exportation of data to IBM OpenDX.
diff --git a/src/getfem_export.cc b/src/getfem_export.cc
index f3de086..4bbd1cc 100644
--- a/src/getfem_export.cc
+++ b/src/getfem_export.cc
@@ -124,6 +124,93 @@ namespace getfem
return vtktypes[t];
}
+ pfem select_finite_element_for_vtk(const mesh_fem& mf, size_type cv) {
+ bgeot::pgeometric_trans pgt = mf.linked_mesh().trans_of_convex(cv);
+ pfem pf = mf.fem_of_element(cv);
+
+ if (pf != fem_descriptor("FEM_Q2_INCOMPLETE(2)") &&
+ pf != fem_descriptor("FEM_Q2_INCOMPLETE(3)") &&
+ pf != fem_descriptor("FEM_PYRAMID_Q2_INCOMPLETE") &&
+ pf != fem_descriptor("FEM_PYRAMID_Q2_INCOMPLETE_DISCONTINUOUS") &&
+ pf != fem_descriptor("FEM_PRISM_INCOMPLETE_P2") &&
+ pf != fem_descriptor("FEM_PRISM_INCOMPLETE_P2_DISCONTINUOUS")) {
+ bool discontinuous = false;
+ for (unsigned i=0; i < pf->nb_dof(cv); ++i) {
+ /* could be a better test for discontinuity .. */
+ if (!dof_linkable(pf->dof_types()[i])) { discontinuous = true; break; }
+ }
+
+ pfem classical_pf1 = discontinuous ? classical_discontinuous_fem(pgt, 1)
+ : classical_fem(pgt, 1);
+
+ short_type degree = 1;
+ if ((pf != classical_pf1 && pf->estimated_degree() > 1) ||
+ pgt->structure() != pgt->basic_structure())
+ degree = 2;
+
+ pf = discontinuous ? classical_discontinuous_fem(pgt, degree, 0, true) :
+ classical_fem(pgt, degree, true);
+ }
+
+ return pf;
+ }
+
+ /* try to check if a quad or hexahedric cell is "rectangular" and oriented
+ along the axes */
+ template<typename C> static bool check_voxel(const C& c) {
+ scalar_type h[3];
+ unsigned N = c[0].size();
+ if (c.size() != (1U << N)) return false;
+ const base_node P0 = c[0];
+ h[0] = c[1][0] - P0[0];
+ h[1] = c[2][0] - P0[0];
+ if (c.size() != 4) h[2] = c[4][0] - P0[0];
+ for (unsigned i=1; i < c.size(); ++i) {
+ const base_node d = c[i] - P0;
+ for (unsigned j=0; j < N; ++j)
+ if (gmm::abs(d[j]) > 1e-7*h[j] && gmm::abs(d[j] - h[j]) > 1e-7*h[j])
+ return false;
+ }
+ return true;
+ }
+
+
+ vtk_mapping_type select_vtk_mapping_type(mesh m, size_type cv, size_type
dim, size_type nbd) {
+ vtk_mapping_type t = NO_VTK_MAPPING;
+ switch (dim) {
+ case 0: t = N1_TO_VTK_VERTEX; break;
+ case 1:
+ if (nbd == 2) t = N2_TO_VTK_LINE;
+ else if (nbd == 3) t = N3_TO_VTK_QUADRATIC_EDGE;
+ break;
+ case 2:
+ if (nbd == 3) t = N3_TO_VTK_TRIANGLE;
+ else if (nbd == 4)
+ t = check_voxel(m.points_of_convex(cv)) ? N4_TO_VTK_PIXEL
+ : N4_TO_VTK_QUAD;
+ else if (nbd == 6) t = N6_TO_VTK_QUADRATIC_TRIANGLE;
+ else if (nbd == 8) t = N8_TO_VTK_QUADRATIC_QUAD;
+ else if (nbd == 9) t = N9_TO_VTK_BIQUADRATIC_QUAD;
+ break;
+ case 3:
+ if (nbd == 4) t = N4_TO_VTK_TETRA;
+ else if (nbd == 10) t = N10_TO_VTK_QUADRATIC_TETRA;
+ else if (nbd == 8)
+ t = check_voxel(m.points_of_convex(cv)) ? N8_TO_VTK_VOXEL
+ : N8_TO_VTK_HEXAHEDRON;
+ else if (nbd == 20) t = N20_TO_VTK_QUADRATIC_HEXAHEDRON;
+ else if (nbd == 27) t = N27_TO_VTK_TRIQUADRATIC_HEXAHEDRON;
+ else if (nbd == 5) t = N5_TO_VTK_PYRAMID;
+ else if (nbd == 13) t = N13_TO_VTK_QUADRATIC_PYRAMID;
+ else if (nbd == 14) t = N14_TO_VTK_QUADRATIC_PYRAMID;
+ else if (nbd == 6) t = N6_TO_VTK_WEDGE;
+ else if (nbd == 15) t = N15_TO_VTK_QUADRATIC_WEDGE;
+ else if (nbd == 18) t = N18_TO_VTK_BIQUADRATIC_QUADRATIC_WEDGE;
+ break;
+ }
+ return t;
+ }
+
vtk_export::vtk_export(std::ostream &os_, bool ascii_)
: os(os_), ascii(ascii_) { init(); }
@@ -173,26 +260,6 @@ namespace getfem
}
- /* try to check if a quad or hexahedric cell is "rectangular" and oriented
- along the axes */
- template<typename C> static bool check_voxel(const C& c) {
- scalar_type h[3];
- unsigned N = c[0].size();
- if (c.size() != (1U << N)) return false;
- const base_node P0 = c[0];
- h[0] = c[1][0] - P0[0];
- h[1] = c[2][0] - P0[0];
- if (c.size() != 4) h[2] = c[4][0] - P0[0];
- for (unsigned i=1; i < c.size(); ++i) {
- const base_node d = c[i] - P0;
- for (unsigned j=0; j < N; ++j)
- if (gmm::abs(d[j]) > 1e-7*h[j] && gmm::abs(d[j] - h[j]) > 1e-7*h[j])
- return false;
- }
- return true;
- }
-
-
void vtk_export::exporting(const stored_mesh_slice& sl) {
psl = &sl; dim_ = dim_type(sl.dim());
GMM_ASSERT1(psl->dim() <= 3, "attempt to export a " << int(dim_)
@@ -202,7 +269,7 @@ namespace getfem
void vtk_export::exporting(const mesh& m) {
dim_ = m.dim();
GMM_ASSERT1(dim_ <= 3, "attempt to export a " << int(dim_)
- << "D slice (not supported)");
+ << "D mesh (not supported)");
pmf = std::make_unique<mesh_fem>(const_cast<mesh&>(m), dim_type(1));
for (dal::bv_visitor cv(m.convex_index()); !cv.finished(); ++cv) {
bgeot::pgeometric_trans pgt = m.trans_of_convex(cv);
@@ -215,41 +282,14 @@ namespace getfem
void vtk_export::exporting(const mesh_fem& mf) {
dim_ = mf.linked_mesh().dim();
GMM_ASSERT1(dim_ <= 3, "attempt to export a " << int(dim_)
- << "D slice (not supported)");
+ << "D mesh_fem (not supported)");
if (&mf != pmf.get())
pmf = std::make_unique<mesh_fem>(mf.linked_mesh());
/* initialize pmf with finite elements suitable for VTK (which only knows
isoparametric FEMs of order 1 and 2) */
for (dal::bv_visitor cv(mf.convex_index()); !cv.finished(); ++cv) {
- bgeot::pgeometric_trans pgt = mf.linked_mesh().trans_of_convex(cv);
- pfem pf = mf.fem_of_element(cv);
-
- if (pf == fem_descriptor("FEM_Q2_INCOMPLETE(2)") ||
- pf == fem_descriptor("FEM_Q2_INCOMPLETE(3)") ||
- pf == fem_descriptor("FEM_PYRAMID_Q2_INCOMPLETE") ||
- pf == fem_descriptor("FEM_PYRAMID_Q2_INCOMPLETE_DISCONTINUOUS") ||
- pf == fem_descriptor("FEM_PRISM_INCOMPLETE_P2") ||
- pf == fem_descriptor("FEM_PRISM_INCOMPLETE_P2_DISCONTINUOUS"))
- pmf->set_finite_element(cv, pf);
- else {
- bool discontinuous = false;
- for (unsigned i=0; i < pf->nb_dof(cv); ++i) {
- /* could be a better test for discontinuity .. */
- if (!dof_linkable(pf->dof_types()[i])) { discontinuous = true;
break; }
- }
-
- pfem classical_pf1 = discontinuous ? classical_discontinuous_fem(pgt,
1)
- : classical_fem(pgt, 1);
-
- short_type degree = 1;
- if ((pf != classical_pf1 && pf->estimated_degree() > 1) ||
- pgt->structure() != pgt->basic_structure())
- degree = 2;
-
- pmf->set_finite_element(cv, discontinuous ?
- classical_discontinuous_fem(pgt, degree, 0,
true) :
- classical_fem(pgt, degree, true));
- }
+ pfem pf = select_finite_element_for_vtk(mf, cv);
+ pmf->set_finite_element(cv, pf);
}
/* find out which dof will be exported to VTK */
@@ -257,53 +297,20 @@ namespace getfem
pmf_mapping_type.resize(pmf->convex_index().last_true() + 1, unsigned(-1));
pmf_dof_used.sup(0, pmf->nb_basic_dof());
for (dal::bv_visitor cv(pmf->convex_index()); !cv.finished(); ++cv) {
- vtk_mapping_type t = NO_VTK_MAPPING;
+ size_type dim = pmf->fem_of_element(cv)->dim();
size_type nbd = pmf->fem_of_element(cv)->nb_dof(cv);
- switch (pmf->fem_of_element(cv)->dim()) {
- case 0: t = N1_TO_VTK_VERTEX; break;
- case 1:
- if (nbd == 2) t = N2_TO_VTK_LINE;
- else if (nbd == 3) t = N3_TO_VTK_QUADRATIC_EDGE;
- break;
- case 2:
- if (nbd == 3) t = N3_TO_VTK_TRIANGLE;
- else if (nbd == 4)
- t = check_voxel(m.points_of_convex(cv)) ? N4_TO_VTK_PIXEL
- : N4_TO_VTK_QUAD;
- else if (nbd == 6) t = N6_TO_VTK_QUADRATIC_TRIANGLE;
- else if (nbd == 8) t = N8_TO_VTK_QUADRATIC_QUAD;
- else if (nbd == 9) t = N9_TO_VTK_BIQUADRATIC_QUAD;
- break;
- case 3:
- if (nbd == 4) t = N4_TO_VTK_TETRA;
- else if (nbd == 10) t = N10_TO_VTK_QUADRATIC_TETRA;
- else if (nbd == 8)
- t = check_voxel(m.points_of_convex(cv)) ? N8_TO_VTK_VOXEL
- : N8_TO_VTK_HEXAHEDRON;
- else if (nbd == 20) t = N20_TO_VTK_QUADRATIC_HEXAHEDRON;
- else if (nbd == 27) t = N27_TO_VTK_TRIQUADRATIC_HEXAHEDRON;
- else if (nbd == 5) t = N5_TO_VTK_PYRAMID;
- else if (nbd == 13) t = N13_TO_VTK_QUADRATIC_PYRAMID;
- else if (nbd == 14) t = N14_TO_VTK_QUADRATIC_PYRAMID;
- else if (nbd == 6) t = N6_TO_VTK_WEDGE;
- else if (nbd == 15) t = N15_TO_VTK_QUADRATIC_WEDGE;
- else if (nbd == 18) t = N18_TO_VTK_BIQUADRATIC_QUADRATIC_WEDGE;
- break;
- }
+ vtk_mapping_type t = select_vtk_mapping_type(m, cv, dim, nbd);
GMM_ASSERT1(t != -1, "semi internal error. Could not map " <<
name_of_fem(pmf->fem_of_element(cv))
<< " to a VTK cell type");
pmf_mapping_type[cv] = t;
const std::vector<unsigned> &dmap = select_vtk_dof_mapping(t);
- //cout << "nbd = " << nbd << ", t = " << t << ", dmap = "<<dmap << "\n";
GMM_ASSERT1(dmap.size() <= pmf->nb_basic_dof_of_element(cv),
"inconsistency in vtk_dof_mapping");
for (unsigned i=0; i < dmap.size(); ++i)
pmf_dof_used.add(pmf->ind_basic_dof_of_element(cv)[dmap[i]]);
}
- // cout << "mf.nb_dof = " << mf.nb_dof() << ", pmf->nb_dof="
- // << pmf->nb_dof() << ", dof_used = " << pmf_dof_used.card() << "\n";
}
@@ -449,6 +456,161 @@ namespace getfem
}
}
+ vtu_export::vtu_export(std::ostream &os_, bool ascii_)
+ : os(os_), ascii(ascii_) { init(); }
+
+ vtu_export::vtu_export(const std::string& fname, bool ascii_)
+ : os(real_os), ascii(ascii_),
+ real_os(fname.c_str(), !ascii ? std::ios_base::binary | std::ios_base::out
:
+ std::ios_base::out) {
+ GMM_ASSERT1(real_os, "impossible to write to vtu file '" << fname << "'");
+ init();
+ }
+
+ vtu_export::~vtu_export(){
+ check_footer();
+ }
+
+ void vtu_export::init() {
+ /* TODO: add ascii format */
+ GMM_ASSERT1(ascii, "vtu support only ascii format.");
+ strcpy(header, "Exported by getfem++");
+ state = EMPTY;
+ psl = 0; dim_ = dim_type(-1);
+ check_header();
+ }
+
+ void vtu_export::switch_to_point_data() {
+ if (state != IN_POINT_DATA) {
+ os << "<PointData>\n";
+ state = IN_POINT_DATA;
+ };
+ }
+
+ void vtu_export::switch_to_cell_data() {
+ }
+
+ void vtu_export::exporting(const mesh& m) {
+ dim_ = m.dim();
+ GMM_ASSERT1(dim_ <= 3, "attempt to export a " << int(dim_)
+ << "D mesh (not supported)");
+ pmf = std::make_unique<mesh_fem>(const_cast<mesh&>(m), dim_type(1));
+ for (dal::bv_visitor cv(m.convex_index()); !cv.finished(); ++cv) {
+ bgeot::pgeometric_trans pgt = m.trans_of_convex(cv);
+ pfem pf = getfem::classical_fem(pgt, pgt->complexity() > 1 ? 2 : 1);
+ pmf->set_finite_element(cv, pf);
+ }
+ exporting(*pmf);
+ }
+
+ void vtu_export::exporting(const mesh_fem& mf) {
+ dim_ = mf.linked_mesh().dim();
+ GMM_ASSERT1(dim_ <= 3, "attempt to export a " << int(dim_)
+ << "D mesh_fem (not supported)");
+ if (&mf != pmf.get())
+ pmf = std::make_unique<mesh_fem>(mf.linked_mesh());
+ /* initialize pmf with finite elements suitable for VTK (which only knows
+ isoparametric FEMs of order 1 and 2) */
+ for (dal::bv_visitor cv(mf.convex_index()); !cv.finished(); ++cv) {
+ pfem pf = select_finite_element_for_vtk(mf, cv);
+ pmf->set_finite_element(cv, pf);
+ }
+ /* find out which dof will be exported to VTU */
+
+ const mesh &m = pmf->linked_mesh();
+ pmf_mapping_type.resize(pmf->convex_index().last_true() + 1, unsigned(-1));
+ pmf_dof_used.sup(0, pmf->nb_basic_dof());
+ for (dal::bv_visitor cv(pmf->convex_index()); !cv.finished(); ++cv) {
+ size_type dim = pmf->fem_of_element(cv)->dim();
+ size_type nbd = pmf->fem_of_element(cv)->nb_dof(cv);
+ vtk_mapping_type t = select_vtk_mapping_type(m, cv, dim, nbd);
+ GMM_ASSERT1(t != -1, "semi internal error. Could not map " <<
+ name_of_fem(pmf->fem_of_element(cv))
+ << " to a VTK cell type");
+ pmf_mapping_type[cv] = t;
+
+ const std::vector<unsigned> &dmap = select_vtk_dof_mapping(t);
+ GMM_ASSERT1(dmap.size() <= pmf->nb_basic_dof_of_element(cv),
+ "inconsistency in vtk_dof_mapping");
+ for (unsigned i=0; i < dmap.size(); ++i)
+ pmf_dof_used.add(pmf->ind_basic_dof_of_element(cv)[dmap[i]]);
+ }
+ os << "<Piece NumberOfPoints=\"" << pmf_dof_used.card();
+ os << "\" NumberOfCells=\"" << pmf->convex_index().card() << "\">\n";
+ }
+
+ void vtu_export::check_header() {
+ if (state >= HEADER_WRITTEN) return;
+ os << "<?xml version=\"1.0\"?>\n";
+ os << "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\">\n";
+ os << "<!--" << header << "-->\n";
+ os << "<UnstructuredGrid>\n";
+ state = HEADER_WRITTEN;
+ }
+
+ void vtu_export::check_footer() {
+ if (state >= FOOTER_WRITTEN) return;
+ if (state == IN_POINT_DATA) os << "</PointData>\n";
+ os << "</Piece>\n";
+ os << "</UnstructuredGrid>\n";
+ os << "</VTKFile>\n";
+ state = FOOTER_WRITTEN;
+ }
+
+ void vtu_export::write_separ()
+ { if (ascii) os << "\n"; }
+
+ void vtu_export::write_mesh(bool only_mesh) {
+ write_mesh_structure_from_mesh_fem();
+ if (only_mesh == true) check_footer();
+ }
+
+ void vtu_export::write_mesh_structure_from_mesh_fem() {
+ if (state >= STRUCTURE_WRITTEN) return;
+ os << "<Points>\n";
+ os << "<DataArray type=\"Float32\" Name=\"Points\" ";
+ os << "NumberOfComponents=\"3\" format=\"ascii\">\n";
+ std::vector<int> dofmap(pmf->nb_dof());
+ int cnt = 0;
+ for (dal::bv_visitor d(pmf_dof_used); !d.finished(); ++d) {
+ dofmap[d] = cnt++;
+ base_node P = pmf->point_of_basic_dof(d);
+ write_vec(P.const_begin(),P.size());
+ write_separ();
+ }
+ size_type nb_cell_values = 0;
+ for (dal::bv_visitor cv(pmf->convex_index()); !cv.finished(); ++cv)
+ nb_cell_values += (1 +
select_vtk_dof_mapping(pmf_mapping_type[cv]).size());
+ os << "</DataArray>\n";
+ os << "</Points>\n";
+ os << "<Cells>\n";
+ os << "<DataArray type=\"Int64\" Name=\"connectivity\"
format=\"ascii\">\n";
+ for (dal::bv_visitor cv(pmf->convex_index()); !cv.finished(); ++cv) {
+ const std::vector<unsigned> &dmap =
select_vtk_dof_mapping(pmf_mapping_type[cv]);
+ for (size_type i=0; i < dmap.size(); ++i)
+ write_val(int(dofmap[pmf->ind_basic_dof_of_element(cv)[dmap[i]]]));
+ write_separ();
+ }
+ os << "</DataArray>\n";
+ os << "<DataArray type=\"Int64\" Name=\"offsets\" format=\"ascii\">\n";
+ cnt = 0;
+ for (dal::bv_visitor cv(pmf->convex_index()); !cv.finished(); ++cv) {
+ const std::vector<unsigned> &dmap =
select_vtk_dof_mapping(pmf_mapping_type[cv]);
+ cnt += int(dmap.size());
+ write_val(cnt);
+ write_separ();
+ }
+ os << "</DataArray>\n";
+ os << "<DataArray type=\"Int64\" Name=\"types\" format=\"ascii\">\n";
+ for (dal::bv_visitor cv(pmf->convex_index()); !cv.finished(); ++cv) {
+ write_val(select_vtk_type(pmf_mapping_type[cv]));
+ write_separ();
+ }
+ os << "</DataArray>\n";
+ os << "</Cells>\n";
+
+ state = STRUCTURE_WRITTEN;
+ }
/* -------------------------------------------------------------
* OPENDX export
diff --git a/tests/thermo_elasticity_electrical_coupling.cc
b/tests/thermo_elasticity_electrical_coupling.cc
index 397e8f0..d545937 100644
--- a/tests/thermo_elasticity_electrical_coupling.cc
+++ b/tests/thermo_elasticity_electrical_coupling.cc
@@ -188,9 +188,14 @@ int main(int argc, char *argv[]) {
}
if (export_mesh) {
- getfem::vtk_export exp("mesh.vtk", false);
- exp.exporting(mesh);
- exp.write_mesh();
+ getfem::vtk_export vtk_exp("mesh.vtk", true);
+ vtk_exp.exporting(mesh);
+ vtk_exp.write_mesh();
+
+ getfem::vtu_export vtu_exp("mesh.vtu", true);
+ vtu_exp.exporting(mesh);
+ vtu_exp.write_mesh();
+
cout << "\nYou can view the mesh for instance with\n";
cout << "mayavi2 -d mesh.vtk -f ExtractEdges -m Surface\n" << endl;
}
@@ -304,21 +309,37 @@ int main(int argc, char *argv[]) {
plain_vector CO(mfvm.nb_dof() * 2);
getfem::ga_interpolation_Lagrange_fem(md, "-"+sigmaeps+"*Grad_V", mfvm, CO);
- getfem::vtk_export exp("displacement_with_von_mises.vtk", false);
- exp.exporting(mfu);
- exp.write_point_data(mfu, U, "elastostatic displacement");
- exp.write_point_data(mfvm, VM, "Von Mises stress");
+ getfem::vtk_export vtk_exp("displacement_with_von_mises.vtk", false);
+ vtk_exp.exporting(mfu);
+ vtk_exp.write_point_data(mfu, U, "elastostatic displacement");
+ vtk_exp.write_point_data(mfvm, VM, "Von Mises stress");
+
+ getfem::vtu_export vtu_exp("displacement_with_von_mises.vtu", true);
+ vtu_exp.exporting(mfu);
+ vtu_exp.write_point_data(mfu, U, "elastostatic displacement");
+ vtu_exp.write_point_data(mfvm, VM, "Von Mises stress");
+
cout << "\nYou can view solutions with for instance:\n\nmayavi2 "
"-d displacement_with_von_mises.vtk -f WarpVector -m Surface\n" << endl;
- getfem::vtk_export exp2("temperature.vtk", false);
- exp2.exporting(mft);
- exp2.write_point_data(mft, THETA, "Temperature");
+ getfem::vtk_export vtk_exp2("temperature.vtk", false);
+ vtk_exp2.exporting(mft);
+ vtk_exp2.write_point_data(mft, THETA, "Temperature");
+
+ getfem::vtu_export vtu_exp2("temperature.vtu", true);
+ vtu_exp2.exporting(mft);
+ vtu_exp2.write_point_data(mft, THETA, "Temperature");
+
cout << "mayavi2 -d temperature.vtk -f WarpScalar -m Surface\n" << endl;
- getfem::vtk_export exp3("electric_potential.vtk", false);
- exp3.exporting(mft);
- exp3.write_point_data(mft, V, "Electric potential");
+ getfem::vtk_export vtk_exp3("electric_potential.vtk", false);
+ vtk_exp3.exporting(mft);
+ vtk_exp3.write_point_data(mft, V, "Electric potential");
+
+ getfem::vtu_export vtu_exp3("electric_potential.vtu", true);
+ vtu_exp3.exporting(mft);
+ vtu_exp3.write_point_data(mft, V, "Electric potential");
+
cout << "mayavi2 -d electric_potential.vtk -f WarpScalar -m Surface\n"
<< endl;