[Getfem-commits] (no subject)

[Tetsuo Koyama]

branch: devel-tetsuo-xml-binary
commit a79688f46f6ce867615dc824bdf18b1e0154400b
Author: Tetsuo Koyama <tkoyama010@gmail.com>
AuthorDate: Sun May 31 12:02:00 2020 +0000

    Add the scripting interface for vtu export
---
 interface/src/gf_mesh_fem_get.cc | 33 ++++++++++++++++++
 interface/src/gf_mesh_get.cc     | 25 ++++++++++++++
 interface/src/gf_slice_get.cc    | 72 ++++++++++++++++++++++++++++++++++++++++
 src/getfem/getfem_export.h       |  4 +--
 src/getfem_export.cc             | 16 ++++-----
 tests/helmholtz.cc               |  3 +-
 6 files changed, 141 insertions(+), 12 deletions(-)

diff --git a/interface/src/gf_mesh_fem_get.cc b/interface/src/gf_mesh_fem_get.cc
index 03abb0c..8a179b1 100644
--- a/interface/src/gf_mesh_fem_get.cc
+++ b/interface/src/gf_mesh_fem_get.cc
@@ -763,6 +763,39 @@ void gf_mesh_fem_get(getfemint::mexargs_in& m_in,
        }
        );
 
+    /*@GET ('export to vtu',@str filename, ... ['ascii'], U, 'name'...)
+    Export a @tmf and some fields to a vtu file.
+
+    The FEM and geometric transformations will be mapped to order 1
+    or 2 isoparametric Pk (or Qk) FEMs (as VTK(XML) does not handle higher
+    order elements). If you need to represent high-order FEMs or
+    high-order geometric transformations, you should consider
+    SLICE:GET('export to vtu').@*/
+    sub_command
+      ("export to vtu", 0, -1, 0, 0,
+       std::string fname = in.pop().to_string();
+       bool ascii = false;
+       while (in.remaining() && in.front().is_string()) {
+        std::string cmd2 = in.pop().to_string();
+        if (cmd_strmatch(cmd2, "ascii"))
+          ascii = true;
+        else THROW_BADARG("expecting 'ascii', got " << cmd2);
+       }
+       getfem::vtu_export exp(fname, ascii);
+       exp.exporting(*mf);
+       exp.write_mesh();
+       int count = 1;
+       while (in.remaining()) {
+        const getfem::mesh_fem *mf2 = mf;
+        if (in.remaining() >= 2 && is_meshfem_object(in.front()))
+          mf2 = to_meshfem_object(in.pop());
+        darray U = in.pop().to_darray();
+        in.last_popped().check_trailing_dimension(int(mf2->nb_dof()));
+        exp.write_point_data(*mf2, U, get_vtk_dataset_name(in, count));
+        count+=1;
+       }
+       );
+
 
     /*@GET ('export to dx',@str filename, ...['as', @str 
mesh_name][,'edges']['serie',@str serie_name][,'ascii'][,'append'], U, 
'name'...)
     Export a @tmf and some fields to an OpenDX file.
diff --git a/interface/src/gf_mesh_get.cc b/interface/src/gf_mesh_get.cc
index 6c9687f..9fcf011 100644
--- a/interface/src/gf_mesh_get.cc
+++ b/interface/src/gf_mesh_get.cc
@@ -1247,6 +1247,31 @@ void gf_mesh_get(getfemint::mexargs_in& m_in,
        );
 
 
+    /*@GET ('export to vtu', @str filename, ... [,'ascii'][,'quality'])
+    Exports a mesh to a VTK(XML) file .
+
+    If 'quality' is specified, an estimation of the quality of each
+    convex will be written to the file.
+
+    See also MESH_FEM:GET('export to vtu'), SLICE:GET('export to vtu').@*/
+    sub_command
+      ("export to vtu", 1, 3, 0, 1,
+       bool write_q = false; bool ascii = false;
+       std::string fname = in.pop().to_string();
+       while (in.remaining() && in.front().is_string()) {
+         std::string cmd2 = in.pop().to_string();
+         if (cmd_strmatch(cmd2, "ascii"))
+           ascii = true;
+         else if (cmd_strmatch(cmd2, "quality"))
+           write_q = true;
+         else THROW_BADARG("expecting 'ascii' or 'quality', got " << cmd2);
+       }
+       getfem::vtu_export exp(fname, ascii);
+       exp.exporting(*pmesh);
+       exp.write_mesh(); if (write_q) exp.write_mesh_quality(*pmesh);
+       );
+
+
     /*@GET ('export to dx', @str filename, ... 
[,'ascii'][,'append'][,'as',@str name,[,'serie',@str serie_name]][,'edges'])
     Exports a mesh to an OpenDX file.
 
diff --git a/interface/src/gf_slice_get.cc b/interface/src/gf_slice_get.cc
index 985072c..9a47f09 100644
--- a/interface/src/gf_slice_get.cc
+++ b/interface/src/gf_slice_get.cc
@@ -438,6 +438,78 @@ void gf_slice_get(getfemint::mexargs_in& m_in,
        );
 
 
+    /*@GET ('export to vtu', @str filename, ...)
+    Export a slice to VTK(XML).
+
+    Following the `filename`, you may use any of the following options:
+
+    - if 'ascii' is not used, the file will contain binary data
+      (non portable, but fast).
+    - if 'edges' is used, the edges of the original mesh will be
+      written instead of the slice content.
+
+    More than one dataset may be written, just list them. Each dataset
+    consists of either:
+
+    - a field interpolated on the slice (scalar, vector or tensor),
+      followed by an optional name.
+    - a mesh_fem and a field, followed by an optional name.
+
+    Examples:
+
+    - SLICE:GET('export to vtu', 'test.vtu', Usl, 'first_dataset', mf,
+      U2, 'second_dataset')
+    - SLICE:GET('export to vtu', 'test.vtu', 'ascii', mf,U2)
+    - SLICE:GET('export to vtu', 'test.vtu', 'edges', 'ascii', Uslice)@*/
+    sub_command
+      ("export to vtu",1, -1, 0, 0,
+       std::string fname = in.pop().to_string();
+       bool ascii = false;
+       bool edges = false;
+       while (in.remaining() && in.front().is_string()) {
+         std::string cmd2 = in.pop().to_string();
+         if (cmd_strmatch(cmd2, "ascii"))
+           ascii = true;
+         else if (cmd_strmatch(cmd2, "edges"))
+           edges = true;
+         else THROW_BADARG("expecting 'ascii' or 'edges', got " << cmd2);
+       }
+       getfem::vtu_export exp(fname, ascii);
+       getfem::stored_mesh_slice sl_edges;
+       const getfem::stored_mesh_slice *vtu_slice = sl;
+       getfem::mesh m_edges;
+       if (edges) {
+         vtu_slice = &sl_edges;
+         dal::bit_vector slice_edges;
+         getfem::mesh_slicer slicer(sl->linked_mesh());
+         getfem::slicer_build_edges_mesh action(m_edges,slice_edges);
+         slicer.push_back_action(action); slicer.exec(*sl);
+         sl_edges.build(m_edges, getfem::slicer_none());
+       }
+       exp.exporting(*vtu_slice);
+       exp.write_mesh();
+       int count = 1;
+       if (in.remaining()) {
+         do {
+           if (in.remaining() >= 2 && is_meshfem_object(in.front())) {
+             const getfem::mesh_fem &mf = *to_meshfem_object(in.pop());
+
+             darray U = in.pop().to_darray();
+             in.last_popped().check_trailing_dimension(int(mf.nb_dof()));
+
+             exp.write_point_data(mf,U,get_vtk_dataset_name(in, count));
+           } else if (in.remaining()) {
+             darray slU = in.pop().to_darray();
+             
in.last_popped().check_trailing_dimension(int(vtu_slice->nb_points()));
+
+             exp.write_sliced_point_data(slU,get_vtk_dataset_name(in, count));
+           } else THROW_BADARG("don't know what to do with this argument")
+               count+=1;
+         } while (in.remaining());
+       }
+       );
+
+
     /*@GET ('export to pov', @str filename)
       Export a the triangles of the slice to POV-RAY.@*/
     sub_command
diff --git a/src/getfem/getfem_export.h b/src/getfem/getfem_export.h
index b4436cf..d5888ff 100644
--- a/src/getfem/getfem_export.h
+++ b/src/getfem/getfem_export.h
@@ -286,8 +286,8 @@ namespace getfem {
 
   class vtu_export : public vtk_export {
   public:
-    vtu_export(const std::string& fname) : vtk_export(fname, true, false) {}
-    vtu_export(std::ostream &os_) : vtk_export(os_, true, false) {}
+    vtu_export(const std::string& fname, bool ascii_) : vtk_export(fname, 
ascii_, false) {}
+    vtu_export(std::ostream &os_, bool ascii_) : vtk_export(os_, ascii_, 
false) {}
   };
 
   /** @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 acfca88..63abc0c 100644
--- a/src/getfem_export.cc
+++ b/src/getfem_export.cc
@@ -189,7 +189,6 @@ namespace getfem
   }
 
   void vtk_export::init() {
-    GMM_ASSERT1(vtk || ascii, "vtu support only ascii format.");
     strcpy(header, "Exported by GetFEM");
     psl = 0; dim_ = dim_type(-1);
     static int test_endian = 0x01234567;
@@ -432,12 +431,12 @@ namespace getfem
     }
     for (size_type ic=0; ic < psl->nb_convex(); ++ic) {
       const getfem::mesh_slicer::cs_simplexes_ct& s = psl->simplexes(ic);
-      for (size_type i=0; i < s.size(); ++i) {
+      for (const auto &val : s) {
         if (vtk) {
-          write_val(int(s[i].dim()+1));
+          write_val(int(val.dim()+1));
         }
-        for (size_type j=0; j < s[i].dim()+1; ++j)
-          write_val(int(s[i].inodes[j] + nodes_cnt));
+        for (size_type j=0; j < val.dim()+1; ++j)
+          write_val(int(val.inodes[j] + nodes_cnt));
         write_separ();
       }
       nodes_cnt += psl->nodes(ic).size();
@@ -452,14 +451,13 @@ namespace getfem
     int cnt = 0;
     for (size_type ic=0; ic < psl->nb_convex(); ++ic) {
       const getfem::mesh_slicer::cs_simplexes_ct& s = psl->simplexes(ic);
-      for (size_type i=0; i < s.size(); ++i) {
+      for (const auto &val : s)
         if (vtk) {
-          write_val(int(vtk_simplex_code[s[i].dim()]));
+          write_val(int(vtk_simplex_code[val.dim()]));
         } else {
-          cnt += int(s[i].dim()+1);
+          cnt += int(val.dim()+1);
           write_val(cnt);
         }
-      }
       write_separ();
       splx_cnt -= s.size();
     }
diff --git a/tests/helmholtz.cc b/tests/helmholtz.cc
index 6b4a8d1..6cf28ac 100644
--- a/tests/helmholtz.cc
+++ b/tests/helmholtz.cc
@@ -235,7 +235,8 @@ int main(int argc, char *argv[]) {
     getfem::vtk_export vtk_exp(p.datafilename + ".vtk",
                           p.PARAM.int_value("VTK_EXPORT")==1, true);
     cout << "export to " << p.datafilename + ".vtu" << "..\n";
-    getfem::vtu_export vtu_exp(p.datafilename + ".vtu");
+    getfem::vtu_export vtu_exp(p.datafilename + ".vtu",
+                           p.PARAM.int_value("VTK_EXPORT")==0);
     getfem::stored_mesh_slice sl(p.mesh, p.mesh.nb_convex() < 2000 ? 8 : 6);
     vtk_exp.exporting(sl);
     vtk_exp.write_point_data(p.mf_u, gmm::real_part(U), "helmholtz_rfield");