[Getfem-commits] (no subject)

[Tetsuo Koyama]

branch: devel-tetsuo-xml
commit 09bcdfedf16b555a1f3e133edb9eb897d75c22b5
Author: Tetsuo Koyama <address@hidden>
AuthorDate: Wed Apr 29 03:08:07 2020 +0000

    :heavy_plus_sign: VTU export
---
 src/getfem/getfem_export.h                     |  29 ++++
 src/getfem_export.cc                           | 176 ++++++++++++++++++++++++-
 tests/test_export.cc                           |  69 ----------
 tests/thermo_elasticity_electrical_coupling.cc |  11 +-
 4 files changed, 209 insertions(+), 76 deletions(-)

diff --git a/src/getfem/getfem_export.h b/src/getfem/getfem_export.h
index f7b601b..5de5257 100644
--- a/src/getfem/getfem_export.h
+++ b/src/getfem/getfem_export.h
@@ -294,15 +294,44 @@ namespace getfem {
   class vtu_export {
   protected:
     std::ostream &os;
+    char header[256]; // hard limit in vtu
     bool ascii;
+    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, FOOTER_WRITTEN } state;
+
+    template<class T> void write_val(T v);
+    template<class V> void write_vec(V p, size_type qdim);
+    void write_separ();
   public:
     vtu_export(const std::string& fname, bool ascii_= false);
     vtu_export(std::ostream &os_, bool ascii_ = false);
     void exporting(const mesh& m);
+    void exporting(const mesh_fem& mf);
     void write_mesh();
+  private:
+    void init();
+    void check_header();
+    void check_footer();
+    void write_mesh_structure_from_mesh_fem();
   };
 
+  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]);
+  }
+
   /** @brief A (quite large) class for exportation of data to IBM OpenDX.
 
                      http://www.opendx.org/
diff --git a/src/getfem_export.cc b/src/getfem_export.cc
index e3ae7ba..9082b30 100644
--- a/src/getfem_export.cc
+++ b/src/getfem_export.cc
@@ -202,7 +202,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,7 +215,7 @@ 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
@@ -450,20 +450,188 @@ namespace getfem
   }
 
   vtu_export::vtu_export(std::ostream &os_, bool ascii_)
-    : os(os_), ascii(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();
+  }
+
+  void vtu_export::init() {
+    strcpy(header, "Exported by getfem++");
+    state = EMPTY;
   }
 
   void vtu_export::exporting(const mesh& m) {
-     (void) 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) {
+      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));
+      }
+    }
+    /* find out which dof will be exported to VTK */
+
+    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) {
+      vtk_mapping_type t = NO_VTK_MAPPING;
+      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;
+      }
+      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";
+  }
+
+  void vtu_export::check_header() {
+    if (state >= HEADER_WRITTEN) return;
+    os << "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" 
byte_order=\"BigEndian\">\n";
+    state = HEADER_WRITTEN;
   }
 
+  void vtu_export::check_footer() {
+    if (state >= FOOTER_WRITTEN) return;
+    os << "</VTKFile>\n";
+    state = FOOTER_WRITTEN;
+  }
+
+  void vtu_export::write_separ()
+  { if (ascii) os << "\n"; }
+
   void vtu_export::write_mesh() {
+    write_mesh_structure_from_mesh_fem();
+  }
+
+  void vtu_export::write_mesh_structure_from_mesh_fem() {
+    if (state >= STRUCTURE_WRITTEN) return;
+    check_header();
+    os << "<UnstructuredGrid>\n";
+    os << "<Piece \"NumberOfPoints=\"" << pmf_dof_used.card() << 
"\"NumberOfCells=\"" << pmf->convex_index().card() << "\n";
+    os << "<Points>\n";
+    os << "<DataArray type=\"Float32\" NumberOfComponents=\"" << dim_<< "\" 
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=\"Int32\" 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=\"Int32\" 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";
+    os << "</Piece>\n";
+    os << "</UnstructuredGrid>\n";
+    check_footer();
+
+    state = STRUCTURE_WRITTEN;
   }
 
   /* -------------------------------------------------------------
diff --git a/tests/test_export.cc b/tests/test_export.cc
deleted file mode 100644
index 7cacfee..0000000
--- a/tests/test_export.cc
+++ /dev/null
@@ -1,69 +0,0 @@
-/*===========================================================================
-
- Copyright (C) 2020-2020 Tetsuo Koyama
-
- This file is a part of GetFEM
-
- GetFEM  is  free software;  you  can  redistribute  it  and/or modify it
- under  the  terms  of the  GNU  Lesser General Public License as published
- by  the  Free Software Foundation;  either version 3 of the License,  or
- (at your option) any later version along with the GCC Runtime Library
- Exception either version 3.1 or (at your option) any later version.
- This program  is  distributed  in  the  hope  that it will be useful,  but
- WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
- or  FITNESS  FOR  A PARTICULAR PURPOSE.  See the GNU Lesser General Public
- License and GCC Runtime Library Exception for more details.
- You  should  have received a copy of the GNU Lesser General Public License
- along  with  this program;  if not, write to the Free Software Foundation,
- Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301, USA.
-
-===========================================================================*/
-#include <iostream>
-#include <string>
-#include <boost/property_tree/ptree.hpp>
-#include <boost/property_tree/xml_parser.hpp>
-#include <boost/foreach.hpp>
-#include <boost/lexical_cast.hpp>
-#include "getfem/getfem_regular_meshes.h"
-#include "getfem/getfem_export.h"
-
-using bgeot::base_node;
-using namespace boost::property_tree;
-
-int main(void) {
-
-  getfem::mesh m0;
-
-  base_node org(1);
-  org[0] = 0.0;
-
-  bgeot::base_small_vector h(1);
-  h[0] = 1;
-  std::vector<bgeot::base_small_vector> vect(1);
-  vect[0] = bgeot::base_small_vector(h);
-
-  std::vector<int> ref(1);
-  ref[0] = 1;
-
-  getfem::parallelepiped_regular_simplex_mesh(m0, 1, org, vect.begin(), 
ref.begin());
-
-  getfem::vtk_export vtk_exp("m0.vtk", true);
-  vtk_exp.exporting(m0);
-  vtk_exp.write_mesh();
-
-  getfem::vtu_export vtu_exp("m0.vtu", true);
-  vtu_exp.exporting(m0);
-  vtu_exp.write_mesh();
-
-  ptree pt;
-  read_xml("m0.vtu", pt);
-
-  if (boost::optional<std::string> str = 
pt.get_optional<std::string>("VTKFile.<xmlattr>.type")) {
-    GMM_ASSERT1(str.get() == "UnstructuredGrid", "UnstructuredGrid attribute 
is not incorrect.");
-  } else {
-    GMM_ASSERT1(false, "UnstructuredGrid attribute is not exist.");
-  }
-
-  return 0;
-}
-
diff --git a/tests/thermo_elasticity_electrical_coupling.cc 
b/tests/thermo_elasticity_electrical_coupling.cc
index 397e8f0..18a8ebe 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;
   }