[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Getfem-commits] (no subject)
|
From: |
Konstantinos Poulios |
|
Subject: |
[Getfem-commits] (no subject) |
|
Date: |
Sat, 23 May 2020 15:07:12 -0400 (EDT) |
branch: devel-logari81-xml
commit 9f956bb759bda182a25ae071f8e691c4a345bf12
Author: Konstantinos Poulios <address@hidden>
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;
}