branch: devel-logari81-xml commit 9f956bb759bda182a25ae071f8e691c4a345bf12 Author: Konstantinos Poulios <logar...@gmail.com> AuthorDate: Sat May 23 21:06:56 2020 +0200
Export to VTU file format (only ascii), based on the original implementation by Tetsuo Koyama --- doc/sphinx/source/replaces.txt | 5 +- doc/sphinx/source/userdoc/export.rst | 53 ++++++++-------- src/getfem/getfem_export.h | 57 ++++++++++++------ src/getfem_export.cc | 113 ++++++++++++++++++++++++++--------- 4 files changed, 159 insertions(+), 69 deletions(-) diff --git a/doc/sphinx/source/replaces.txt b/doc/sphinx/source/replaces.txt index 5710048..38ed62a 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`` @@ -84,7 +86,8 @@ .. _specific guides: http://getfem.org/index.html .. _user documentation: http://getfem.org/userdoc/index.html .. _vocabulary: http://getfem.org/getfem_reference/index.html -.. _VTK: http://www.vtk.org +.. _VTK: https://vtk.org/Wiki/VTK_XML_Formats +.. _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..1433153 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/VTU file viewer. * The open-source OpenDX program. * The open-source Gmsh program. @@ -43,9 +43,13 @@ 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 -`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|. +Four file formats are supported for export: the `VTK`_ and `VTU`_ file +formats, the`OpenDX`_ file format and the `Gmsh`_ post-processing file +format. All four can be used for exporting either a |gf_m| or |gf_mf|, and +all except `VTU`_ can be used for exporting the more versatile |gf_smsl|. +The corresponding four classes: |gf_vtk_export|, |gf_vtu_export|, +|gf_dx_export| and |gf_pos_export| are contained in the file +:file:`getfem/getfem_export.h`. Examples of use can be found in the examples of the tests directory. @@ -178,22 +182,15 @@ 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|, +|opendx|, or |gmsh|. -Exporting |m|, |mf| or slices to VTK ------------------------------------- +Exporting |m|, |mf| or slices to VTK/VTU +----------------------------------------- -First, it is important to know the limitation of VTK 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, -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 -than one type of element (i.e. meshes with triangles and quadrangles, for -example). +VTK/VTU files can handle data on segment, triangles, quadrangles, +tetrahedrons and hexahedrons of first or second degree. For example, supposing that a |smsl| ``sl`` has already been built:: @@ -204,11 +201,10 @@ For example, supposing that a |smsl| ``sl`` has already been built:: exp.write_point_data(mfp, P, "pressure"); // write a scalar field 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 -after the scalar fields (and the tensor fields should be written last). +In this example, the fields ``P`` and ``U`` are interpolated on the slice +nodes and then written into the VTK field. -It is also possible to export a |mf| without having to build a slice:: +It is also possible to export a |mf| ``mfu`` 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 @@ -217,9 +213,18 @@ It is also possible to export a |mf| without having to build a slice:: 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 -greater than 2, there will be a loss of precision for higher degree FEMs. +An |mf| ``mfu`` can also be exported in the VTU format with:: + + // VTU export is limitted to ascii output and cannot be used for slices + vtu_export exp("output.vtu); + exp.exporting(mfu); // will save the geometrical structure of the mesh_fem + 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 when exporing a |mf| with ``vtk_export`` or ``vtu_export`` +each convex/fem of ``mfu`` will be mapped to a VTK/VTU element type. As +VTK/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 eb24d4b..b4436cf 100644 --- a/src/getfem/getfem_export.h +++ b/src/getfem/getfem_export.h @@ -57,18 +57,20 @@ namespace getfem { return s2; } - /** @brief VTK export. + /** @brief VTK/VTU export. - export class to VTK ( http://www.kitware.com/vtk.html ) file format - (not the XML format, but the old format) + export class to VTK/VTU file format + ( http://www.kitware.com/vtk.html, + legacy and serial vtkUnstructuredGrid) A vtk_export can store multiple scalar/vector fields. */ class vtk_export { protected: std::ostream &os; - char header[256]; // hard limit in vtk + char header[256]; // hard limit in vtk/vtu bool ascii; + bool vtk; // true for vtk, false for vtu const stored_mesh_slice *psl; std::unique_ptr<mesh_fem> pmf; dal::bit_vector pmf_dof_used; @@ -85,9 +87,9 @@ namespace getfem { void write_separ(); public: - vtk_export(const std::string& fname, bool ascii_ = false); - vtk_export(std::ostream &os_, bool ascii_ = false); - + vtk_export(const std::string& fname, bool ascii_= false, bool vtk_=true); + vtk_export(std::ostream &os_, bool ascii_ = false, bool vtk_=true); + ~vtk_export(); // a vtu file is completed upon calling the destructor /** should be called before write_*_data */ void exporting(const mesh& m); void exporting(const mesh_fem& mf); @@ -99,9 +101,9 @@ namespace getfem { void set_header(const std::string& s); void write_mesh(); - /** append a new scalar or vector field defined on mf to the .vtk 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(). + /** append a new scalar or vector field defined on mf to the .vtk/.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 @@ -242,31 +244,52 @@ namespace getfem { GMM_ASSERT1(gmm::vect_size(U) == nb_val*Q, "inconsistency in the size of the dataset: " << gmm::vect_size(U) << " != " << nb_val << "*" << Q); - write_separ(); + if (vtk) write_separ(); if (Q == 1) { - os << "SCALARS " << remove_spaces(name) << " float 1\n"; - os << "LOOKUP_TABLE default\n"; + if (vtk) + os << "SCALARS " << remove_spaces(name) << " float 1\n" + << "LOOKUP_TABLE default\n"; + else + os << "<DataArray type=\"Float32\" Name=\"" << remove_spaces(name) + << "\" format=\"ascii\">\n"; for (size_type i=0; i < nb_val; ++i) { write_val(float(U[i])); } } else if (Q <= 3) { - os << "VECTORS " << remove_spaces(name) << " float\n"; + if (vtk) + os << "VECTORS " << remove_spaces(name) << " float\n"; + else + os << "<DataArray type=\"Float32\" Name=\"" << remove_spaces(name) + << "\" 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 VTK file, they are written with C (row major) order + in the VTK/VTU file, they are written with C (row major) order */ - os << "TENSORS " << remove_spaces(name) << " float\n"; + if (vtk) + os << "TENSORS " << remove_spaces(name) << " float\n"; + else + os << "<DataArray type=\"Float32\" Name=\"" << remove_spaces(name) + << "\" NumberOfComponents=\"9\" format=\"ascii\">"; for (size_type i=0; i < nb_val; ++i) { write_3x3tensor(U.begin() + i*Q); } - } else GMM_ASSERT1(false, "vtk does not accept vectors of dimension > 3"); + } else + GMM_ASSERT1(false, std::string(vtk ? "vtk" : "vtu") + + " does not accept vectors of dimension > 3"); write_separ(); + if (!vtk) os << "</DataArray>\n"; } + 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) {} + }; + /** @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 3b12008..d884c95 100644 --- a/src/getfem_export.cc +++ b/src/getfem_export.cc @@ -167,18 +167,29 @@ namespace getfem } - vtk_export::vtk_export(std::ostream &os_, bool ascii_) - : os(os_), ascii(ascii_) { init(); } + vtk_export::vtk_export(std::ostream &os_, bool ascii_, bool vtk_) + : os(os_), ascii(ascii_), vtk(vtk_) { init(); } - vtk_export::vtk_export(const std::string& fname, bool ascii_) - : os(real_os), ascii(ascii_), + vtk_export::vtk_export(const std::string& fname, bool ascii_, bool vtk_) + : os(real_os), ascii(ascii_), vtk(vtk_), 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 file '" << fname << "'"); init(); } + vtk_export::~vtk_export(){ + if (!vtk) { + if (state == IN_CELL_DATA) os << "</CellData>\n"; + if (state == IN_POINT_DATA) os << "</PointData>\n"; + os << "</Piece>\n"; + os << "</UnstructuredGrid>\n"; + os << "</VTKFile>\n"; + } + } + 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; @@ -187,34 +198,41 @@ namespace getfem } void vtk_export::switch_to_cell_data() { + GMM_ASSERT1(vtk, "Export of cell data to vtu not supported yet."); if (state != IN_CELL_DATA) { - state = IN_CELL_DATA; - write_separ(); - if (psl) { - os << "CELL_DATA " << psl->nb_simplexes(0) + psl->nb_simplexes(1) - + psl->nb_simplexes(2) + psl->nb_simplexes(3) << "\n"; + if (vtk) { + size_type nb_cells=0; + if (psl) for (auto i : {0,1,2,3}) nb_cells += psl->nb_simplexes(i); + else nb_cells = pmf->convex_index().card(); + write_separ(); + os << "CELL_DATA " << nb_cells << "\n"; + write_separ(); } else { - os << "CELL_DATA " << pmf->convex_index().card() << "\n"; + if (state == IN_POINT_DATA) os << "</PointData>\n"; + os << "<CellData>\n"; } - write_separ(); + state = IN_CELL_DATA; } } void vtk_export::switch_to_point_data() { if (state != IN_POINT_DATA) { - state = IN_POINT_DATA; - write_separ(); - if (psl) { - write_separ(); os << "POINT_DATA " << psl->nb_points() << "\n"; + if (vtk) { + write_separ(); + os << "POINT_DATA " << (psl ? psl->nb_points() + : pmf_dof_used.card()) << "\n"; + write_separ(); } else { - os << "POINT_DATA " << pmf_dof_used.card() << "\n"; + if (state == IN_CELL_DATA) os << "</CellData>\n"; + os << "<PointData>\n"; } - write_separ(); + state = IN_POINT_DATA; } } void vtk_export::exporting(const stored_mesh_slice& sl) { + GMM_ASSERT1(vtk, "Export of mesh slice to vtu not supported yet."); psl = &sl; dim_ = dim_type(sl.dim()); GMM_ASSERT1(psl->dim() <= 3, "attempt to export a " << int(dim_) << "D slice (not supported)"); @@ -272,7 +290,7 @@ namespace getfem classical_fem(pgt, degree, true)); } } - /* find out which dof will be exported to VTK */ + /* find out which dof will be exported to VTK/VTU */ const mesh &m = pmf->linked_mesh(); pmf_mapping_type.resize(pmf->convex_index().last_true() + 1, unsigned(-1)); @@ -346,9 +364,16 @@ namespace getfem void vtk_export::check_header() { if (state >= HEADER_WRITTEN) return; - os << "# vtk DataFile Version 2.0\n"; - os << header << "\n"; - if (ascii) os << "ASCII\n"; else os << "BINARY\n"; + if (vtk) { + os << "# vtk DataFile Version 2.0\n"; + os << header << "\n"; + os << (ascii ? "ASCII\n" : "BINARY\n"); + } else { + os << "<?xml version=\"1.0\"?>\n"; + os << "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\">\n"; + os << "<!--" << header << "-->\n"; + os << "<UnstructuredGrid>\n"; + } state = HEADER_WRITTEN; } @@ -417,9 +442,17 @@ namespace getfem void vtk_export::write_mesh_structure_from_mesh_fem() { if (state >= STRUCTURE_WRITTEN) return; check_header(); - /* possible improvement: detect structured grids */ - os << "DATASET UNSTRUCTURED_GRID\n"; - os << "POINTS " << pmf_dof_used.card() << " float\n"; + if (vtk) { + /* possible improvement: detect structured grids */ + os << "DATASET UNSTRUCTURED_GRID\n"; + os << "POINTS " << pmf_dof_used.card() << " float\n"; + } else { + os << "<Piece NumberOfPoints=\"" << pmf_dof_used.card() + << "\" NumberOfCells=\"" << pmf->convex_index().card() << "\">\n"; + 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) { @@ -433,21 +466,47 @@ namespace getfem for (dal::bv_visitor cv(pmf->convex_index()); !cv.finished(); ++cv) nb_cell_values += (1 + select_vtk_dof_mapping(pmf_mapping_type[cv]).size()); - write_separ(); os << "CELLS " << pmf->convex_index().card() << " " << nb_cell_values << "\n"; + if (vtk) { + write_separ(); + os << "CELLS " << pmf->convex_index().card() << " " << nb_cell_values << "\n"; + } else { + 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]); - write_val(int(dmap.size())); + if (vtk) write_val(int(dmap.size())); for (size_type i=0; i < dmap.size(); ++i) write_val(int(dofmap[pmf->ind_basic_dof_of_element(cv)[dmap[i]]])); write_separ(); } - write_separ(); os << "CELL_TYPES " << pmf->convex_index().card() << "\n"; + if (vtk) { + write_separ(); + os << "CELL_TYPES " << pmf->convex_index().card() << "\n"; + } else { + 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(); } + if (!vtk) + os << "</DataArray>\n" + << "</Cells>\n"; state = STRUCTURE_WRITTEN; }