From ad10273d4481bc9cd597b86eb24ed47e380e9f16 Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Thu, 23 May 2024 20:10:09 +0200 Subject: [PATCH 01/58] lib/mito/fem: quadrature fields should always be named --- lib/mito/fem/QuadratureField.h | 2 +- lib/mito/quadrature/Integrator.h | 2 +- tests/mito.lib/fem/quadrature_field.cc | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/lib/mito/fem/QuadratureField.h b/lib/mito/fem/QuadratureField.h index e0fe6619..7f1c0f2e 100644 --- a/lib/mito/fem/QuadratureField.h +++ b/lib/mito/fem/QuadratureField.h @@ -32,7 +32,7 @@ namespace mito::fem { {} private: - inline QuadratureField(int /*nElements*/, const pack_t & packing, std::string name = "") : + inline QuadratureField(int /*nElements*/, const pack_t & packing, std::string name) : _grid{ packing, packing.cells() }, _name(name) {} diff --git a/lib/mito/quadrature/Integrator.h b/lib/mito/quadrature/Integrator.h index 2d7f5477..fd629084 100644 --- a/lib/mito/quadrature/Integrator.h +++ b/lib/mito/quadrature/Integrator.h @@ -53,7 +53,7 @@ namespace mito::quadrature { public: Integrator(const manifold_type & manifold) : _manifold(manifold), - _coordinates(manifold.nElements()) + _coordinates(manifold.nElements(), "coordinates") { _computeQuadPointCoordinates(); } diff --git a/tests/mito.lib/fem/quadrature_field.cc b/tests/mito.lib/fem/quadrature_field.cc index d9cc5e3d..9902202c 100644 --- a/tests/mito.lib/fem/quadrature_field.cc +++ b/tests/mito.lib/fem/quadrature_field.cc @@ -29,7 +29,7 @@ TEST(Fem, QuadratureFields) constexpr int Q = 10; // a field of {Q} mito vectors - auto field = mito::fem::quadrature_field_t>(nElements); + auto field = mito::fem::quadrature_field_t>(nElements, "field"); // get a reference the (0, 1) vector auto & vector = field(0, 1); From 961b550d78c18da4bfafc7487ffe21b55c39ece0 Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Thu, 23 May 2024 20:11:09 +0200 Subject: [PATCH 02/58] lib/mito/fem: remove unused parameter in constructor of class {QuadratureField} --- lib/mito/fem/QuadratureField.h | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/lib/mito/fem/QuadratureField.h b/lib/mito/fem/QuadratureField.h index 7f1c0f2e..55883daf 100644 --- a/lib/mito/fem/QuadratureField.h +++ b/lib/mito/fem/QuadratureField.h @@ -27,12 +27,12 @@ namespace mito::fem { * constructor * @param[in] elements number of elements for which data are stored */ - inline QuadratureField(int nElements, std::string name = "") : - QuadratureField(nElements, pack_t{ { nElements, Q } }, name) + inline QuadratureField(int nElements, std::string name) : + QuadratureField(pack_t{ { nElements, Q } }, name) {} private: - inline QuadratureField(int /*nElements*/, const pack_t & packing, std::string name) : + inline QuadratureField(const pack_t & packing, std::string name) : _grid{ packing, packing.cells() }, _name(name) {} From 5d5bf6be6a54c6adf68e3a240e1caf2bc5f2df92 Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Thu, 23 May 2024 20:11:49 +0200 Subject: [PATCH 03/58] lib/mito/fem: constructor of class {QuadratureField} now takes packing as a temporary --- lib/mito/fem/QuadratureField.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/mito/fem/QuadratureField.h b/lib/mito/fem/QuadratureField.h index 55883daf..8fd417cc 100644 --- a/lib/mito/fem/QuadratureField.h +++ b/lib/mito/fem/QuadratureField.h @@ -32,7 +32,7 @@ namespace mito::fem { {} private: - inline QuadratureField(const pack_t & packing, std::string name) : + inline QuadratureField(const pack_t && packing, std::string name) : _grid{ packing, packing.cells() }, _name(name) {} From 63a775aeb48ba74f10a94087c763f0635a1e09c5 Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Thu, 23 May 2024 20:13:46 +0200 Subject: [PATCH 04/58] lib/mito/fem: consistency fixes in class {QuadratureField} --- lib/mito/fem/QuadratureField.h | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/lib/mito/fem/QuadratureField.h b/lib/mito/fem/QuadratureField.h index 8fd417cc..5b5078c0 100644 --- a/lib/mito/fem/QuadratureField.h +++ b/lib/mito/fem/QuadratureField.h @@ -14,13 +14,13 @@ namespace mito::fem { private: // conventionally packed grid for {e, q} - using pack_t = pyre::grid::canonical_t<2>; + using pack_type = pyre::grid::canonical_t<2>; // of Y on the heap - using storage_t = pyre::memory::heap_t; + using storage_type = pyre::memory::heap_t; // putting it all together - using grid_t = pyre::grid::grid_t; + using grid_type = pyre::grid::grid_t; // index - using index_t = pack_t::index_type; + using index_type = pack_type::index_type; public: /** @@ -28,11 +28,11 @@ namespace mito::fem { * @param[in] elements number of elements for which data are stored */ inline QuadratureField(int nElements, std::string name) : - QuadratureField(pack_t{ { nElements, Q } }, name) + QuadratureField(pack_type{ { nElements, Q } }, name) {} private: - inline QuadratureField(const pack_t && packing, std::string name) : + inline QuadratureField(const pack_type && packing, std::string name) : _grid{ packing, packing.cells() }, _name(name) {} @@ -66,13 +66,13 @@ namespace mito::fem { return operator[]({ e, q }); } - inline Y & operator[](const index_t & index) + inline Y & operator[](const index_type & index) { // all done return _grid[index]; } - inline const Y & operator[](const index_t & index) const + inline const Y & operator[](const index_type & index) const { // all done return _grid[index]; @@ -112,7 +112,7 @@ namespace mito::fem { private: // instantiate the grid - grid_t _grid; + grid_type _grid; // name of the field std::string _name; From 583110eaf2cffd1d3b7ca507540a4c73a752cfe6 Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Thu, 23 May 2024 20:15:38 +0200 Subject: [PATCH 05/58] lib/mito/fem: consistency fixes in class {QuadratureField} --- lib/mito/fem/QuadratureField.h | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/lib/mito/fem/QuadratureField.h b/lib/mito/fem/QuadratureField.h index 5b5078c0..83eb5f28 100644 --- a/lib/mito/fem/QuadratureField.h +++ b/lib/mito/fem/QuadratureField.h @@ -48,7 +48,7 @@ namespace mito::fem { * @param[in] q local index of the quadrature point in the element * @return the data */ - inline Y & operator()(int e, int q) + inline auto operator()(int e, int q) -> Y & { // all done return operator[]({ e, q }); @@ -60,19 +60,19 @@ namespace mito::fem { * @param[in] q local index of the quadrature point in the element * @return the data */ - inline const Y & operator()(int e, int q) const + inline auto operator()(int e, int q) const -> const Y & { // all done return operator[]({ e, q }); } - inline Y & operator[](const index_type & index) + inline auto operator[](const index_type & index) -> Y & { // all done return _grid[index]; } - inline const Y & operator[](const index_type & index) const + inline auto operator[](const index_type & index) const -> const Y & { // all done return _grid[index]; @@ -82,18 +82,18 @@ namespace mito::fem { * accessor for the number of elements * @return the number of elements */ - inline int n_elements() const { return _grid.layout().shape()[0]; } + inline auto n_elements() const -> int { return _grid.layout().shape()[0]; } /* * accessor for the number of quadrature points per element * @return the number of quadrature point per element */ - inline constexpr int n_quad_points() const noexcept { return Q; } + inline constexpr auto n_quad_points() const noexcept -> int { return Q; } /** * const accessor for name */ - inline std::string name() const noexcept { return _name; } + inline auto name() const noexcept -> std::string { return _name; } /** * setter method for name From 34b455ef774f551f16d7aa67aaa3e7db753a5b2d Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Thu, 23 May 2024 20:19:40 +0200 Subject: [PATCH 06/58] lib/mito/fem: use default constructor for class {QuadratureField} --- lib/mito/fem/QuadratureField.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/mito/fem/QuadratureField.h b/lib/mito/fem/QuadratureField.h index 83eb5f28..1904d7c4 100644 --- a/lib/mito/fem/QuadratureField.h +++ b/lib/mito/fem/QuadratureField.h @@ -39,7 +39,7 @@ namespace mito::fem { public: // destructor - ~QuadratureField() {} + ~QuadratureField() = default; public: /** From 679e7a6a796887f241f7a3e914008248f2d9e751 Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Thu, 23 May 2024 20:23:29 +0200 Subject: [PATCH 07/58] lib/mito/fem: in class {QuadratureField} remove redundant {inline} of {constexpr} function --- lib/mito/fem/QuadratureField.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/mito/fem/QuadratureField.h b/lib/mito/fem/QuadratureField.h index 1904d7c4..7a7b688e 100644 --- a/lib/mito/fem/QuadratureField.h +++ b/lib/mito/fem/QuadratureField.h @@ -88,7 +88,7 @@ namespace mito::fem { * accessor for the number of quadrature points per element * @return the number of quadrature point per element */ - inline constexpr auto n_quad_points() const noexcept -> int { return Q; } + constexpr auto n_quad_points() const noexcept -> int { return Q; } /** * const accessor for name From f5fca26cfd31c65eb4c22ee540824ed169eb087e Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Thu, 23 May 2024 20:26:48 +0200 Subject: [PATCH 08/58] lib/mito/fem: in class {QuadratureField} remove setter for {name} attribute Fields should have immutable names. --- lib/mito/fem/QuadratureField.h | 9 --------- 1 file changed, 9 deletions(-) diff --git a/lib/mito/fem/QuadratureField.h b/lib/mito/fem/QuadratureField.h index 7a7b688e..06d96317 100644 --- a/lib/mito/fem/QuadratureField.h +++ b/lib/mito/fem/QuadratureField.h @@ -95,15 +95,6 @@ namespace mito::fem { */ inline auto name() const noexcept -> std::string { return _name; } - /** - * setter method for name - */ - inline void name(std::string name) noexcept - { - _name = name; - return; - } - // support for ranged for loops (wrapping grid) inline auto begin() const { return std::cbegin(_grid); } inline auto end() const { return std::cend(_grid); } From fe1df66589925f3c5b4af4ab60044bd21ae4a608 Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Thu, 23 May 2024 20:27:05 +0200 Subject: [PATCH 09/58] lib/mito/fem: narrative improvements in class {QuadratureField} --- lib/mito/fem/QuadratureField.h | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/lib/mito/fem/QuadratureField.h b/lib/mito/fem/QuadratureField.h index 06d96317..812d6303 100644 --- a/lib/mito/fem/QuadratureField.h +++ b/lib/mito/fem/QuadratureField.h @@ -91,7 +91,7 @@ namespace mito::fem { constexpr auto n_quad_points() const noexcept -> int { return Q; } /** - * const accessor for name + * accessor for name */ inline auto name() const noexcept -> std::string { return _name; } @@ -102,10 +102,10 @@ namespace mito::fem { inline auto end() { return std::end(_grid); } private: - // instantiate the grid + // the underlying grid grid_type _grid; - // name of the field + // the name of the field std::string _name; }; From 552cc1c4959063fb809e212065b406745bc6306f Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Thu, 23 May 2024 20:35:48 +0200 Subject: [PATCH 10/58] lib/mito/fem: redesign of class {NodalField} based on {pyre::grid} --- lib/mito/fem/NodalField.h | 111 ++++++++++++++------------------------ lib/mito/fem/factories.h | 4 +- lib/mito/fem/forward.h | 6 +-- 3 files changed, 46 insertions(+), 75 deletions(-) diff --git a/lib/mito/fem/NodalField.h b/lib/mito/fem/NodalField.h index 196a8d10..f6de87ae 100644 --- a/lib/mito/fem/NodalField.h +++ b/lib/mito/fem/NodalField.h @@ -9,104 +9,75 @@ namespace mito::fem { - // TOFIX: order of template parameters - template + template class NodalField { + + private: + // conventionally packed grid for {n} + using pack_type = pyre::grid::canonical_t<1>; + // of Y on the heap + using storage_type = pyre::memory::heap_t; + // putting it all together + using grid_type = pyre::grid::grid_t; + // index + using index_type = pack_type::index_type; + public: - NodalField(int nodes, std::string name = "") : - _nodes(nodes), - _name(name), - _nodalField(_nodes * D, 0.0) - {} + // constructor + NodalField(int nodes, std::string name) : NodalField(pack_type{ { nodes } }, name) {} - ~NodalField() {} + // destructor + ~NodalField() = default; + + private: + inline NodalField(const pack_type && packing, std::string name) : + _grid{ packing, packing.cells() }, + _name(name) + {} + public: /** * Operator() */ - const T & operator()(const int a, const int i) const - { - assert(i < D); - return _nodalField[D * a + i]; - } + inline auto operator()(const int a) const -> const Y & { return _grid[a]; } - T & operator()(const int a, const int i) - { - assert(i < D); - return _nodalField[D * a + i]; - } + inline auto operator()(const int a) -> Y & { return _grid[a]; } /** - * Iterators + * accessor for the number of nodes */ - typename std::vector::iterator begin() { return std::begin(_nodalField); } - - typename std::vector::iterator end() { return std::end(_nodalField); } - - typename std::vector::const_iterator begin() const { return std::cbegin(_nodalField); } - - typename std::vector::const_iterator end() const { return std::cend(_nodalField); } + inline auto n_nodes() const -> int { return _grid.layout().shape()[0]; } /** - * Accessors + * accessor for name */ inline const std::string & name() const noexcept { return _name; } - inline int dim() const noexcept { return D; } - - inline int nodes() const noexcept { return _nodes; } - - inline int size() const noexcept { return _nodes * D; } - - /** - * Set the field to zero. - */ - inline void init() { std::ranges::fill(_nodalField, 0.0); } + // support for ranged for loops (wrapping grid) + inline auto begin() const { return std::cbegin(_grid); } + inline auto end() const { return std::cend(_grid); } + inline auto begin() { return std::begin(_grid); } + inline auto end() { return std::end(_grid); } private: - /** - * number of nodes - */ - int _nodes; + // the underlying grid + grid_type _grid; - /** - * name of the nodal field - */ + // the name of the field std::string _name; - - /** - * nodal field - */ - std::vector _nodalField; }; - template - std::ostream & operator<<(std::ostream & os, const NodalField & nodalField) + template + std::ostream & operator<<(std::ostream & os, const NodalField & nodalField) { - os << "Nodal field \"" << nodalField.name() << "\" : "; + os << "Nodal field \"" << nodalField.name() << "\" : " << std::endl; - if (std::size(nodalField) == 0) { - os << "[]"; - return os; - } - os << "[(" << nodalField(0, 0); - for (int d = 1; d < D; ++d) { - os << ", " << nodalField(0, d); - } - os << ")"; - - for (auto i = 1; i < nodalField.nodes(); ++i) { - os << ", (" << nodalField(i, 0); - for (int d = 1; d < D; ++d) { - os << ", " << nodalField(i, d); - } - os << ")"; + for (auto i = 0; i < nodalField.n_nodes(); ++i) { + os << "Node " << i << ": " << nodalField(i) << std::endl; } - os << "]"; - return os; } diff --git a/lib/mito/fem/factories.h b/lib/mito/fem/factories.h index dc5e0e14..32d29e67 100644 --- a/lib/mito/fem/factories.h +++ b/lib/mito/fem/factories.h @@ -17,10 +17,10 @@ namespace mito::fem { } // nodal field factory - template + template constexpr auto nodal_field(int nodes, std::string name) { - return nodal_field_t(nodes, name); + return nodal_field_t(nodes, name); } } diff --git a/lib/mito/fem/forward.h b/lib/mito/fem/forward.h index 6664dce2..7b94e071 100644 --- a/lib/mito/fem/forward.h +++ b/lib/mito/fem/forward.h @@ -10,12 +10,12 @@ namespace mito::fem { // class nodal field - template + template class NodalField; // nodal field alias - template - using nodal_field_t = NodalField; + template + using nodal_field_t = NodalField; // class quadrature field template From 8be72ff115ecc1f4856061edad9fbc9e1466c1d7 Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Fri, 24 May 2024 18:30:02 +0200 Subject: [PATCH 11/58] tests/mito.lib/fem: add test of nodal fields --- .cmake/mito_tests_mito_lib.cmake | 1 + tests/mito.lib/fem/nodal_field.cc | 34 +++++++++++++++++++++++++++++++ 2 files changed, 35 insertions(+) create mode 100644 tests/mito.lib/fem/nodal_field.cc diff --git a/.cmake/mito_tests_mito_lib.cmake b/.cmake/mito_tests_mito_lib.cmake index 272566f9..6969148c 100644 --- a/.cmake/mito_tests_mito_lib.cmake +++ b/.cmake/mito_tests_mito_lib.cmake @@ -37,6 +37,7 @@ mito_test_driver(tests/mito.lib/geometry/spherical_metric_space.cc) # fem mito_test_driver(tests/mito.lib/fem/quadrature_field.cc) +mito_test_driver(tests/mito.lib/fem/nodal_field.cc) # io mito_test_driver(tests/mito.lib/io/read_mesh_summit_2D.cc) diff --git a/tests/mito.lib/fem/nodal_field.cc b/tests/mito.lib/fem/nodal_field.cc new file mode 100644 index 00000000..5e71d407 --- /dev/null +++ b/tests/mito.lib/fem/nodal_field.cc @@ -0,0 +1,34 @@ +// -*- c++ -*- +// +// Copyright (c) 2020-2024, the MiTo Authors, all rights reserved +// + +#include +#include + + +TEST(Fem, NodalField) +{ + // say we have 20 nodes + int n_nodes = 20; + + // a field of mito vectors + auto field = mito::fem::nodal_field>(n_nodes, "field"); + + // expect that the number of nodes is 20 + EXPECT_TRUE(field.n_nodes() == n_nodes); + + // populate the field + int i = 0; + for (auto & vector : field) { + if (i == 0) { + vector = mito::vector_t<2>{ 0.0, 0.0 }; + } else { + vector = field(i - 1) + mito::vector_t<2>{ 1, 1 }; + } + ++i; + } + + // expect that editing {vector} has edited the {field} + EXPECT_TRUE((field(10) == mito::vector_t<2>{ 10.0, 10.0 })); +} \ No newline at end of file From bf6a896ce62796d19a7c983c447edd2e3afa5718 Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Fri, 24 May 2024 20:03:16 +0200 Subject: [PATCH 12/58] lib/mito/fem: fix compiler error Define {HAVE_TENSOR} and {HAVE_COMPACT_PACKINGS} before including {pyre}. --- lib/mito/fem/externals.h | 2 ++ 1 file changed, 2 insertions(+) diff --git a/lib/mito/fem/externals.h b/lib/mito/fem/externals.h index ef88065c..23813515 100644 --- a/lib/mito/fem/externals.h +++ b/lib/mito/fem/externals.h @@ -8,6 +8,8 @@ // externals +#define HAVE_TENSOR +#define HAVE_COMPACT_PACKINGS #include #include #include From 56e323faf1cd6b66b3e966df5f0cdbf7070077e8 Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Fri, 24 May 2024 20:03:49 +0200 Subject: [PATCH 13/58] lib/mito/fem: remove unnecessary include --- lib/mito/fem/externals.h | 1 - 1 file changed, 1 deletion(-) diff --git a/lib/mito/fem/externals.h b/lib/mito/fem/externals.h index 23813515..6aa6c979 100644 --- a/lib/mito/fem/externals.h +++ b/lib/mito/fem/externals.h @@ -12,7 +12,6 @@ #define HAVE_COMPACT_PACKINGS #include #include -#include // support #include "../journal.h" From a4c00bba4f61c61a2e21afa43598a59dfb39db96 Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Fri, 24 May 2024 20:05:09 +0200 Subject: [PATCH 14/58] lib/mito/io: draft capability to write a nodal field to vtu file --- .cmake/mito_tests_mito_lib.cmake | 1 + lib/mito/io/externals.h | 1 + lib/mito/io/vtk/externals.h | 1 + lib/mito/io/vtk/writer.h | 31 ++++++++++++++++ tests/mito.lib/fem/nodal_field_sphere.cc | 46 ++++++++++++++++++++++++ 5 files changed, 80 insertions(+) create mode 100644 tests/mito.lib/fem/nodal_field_sphere.cc diff --git a/.cmake/mito_tests_mito_lib.cmake b/.cmake/mito_tests_mito_lib.cmake index 6969148c..da5783e0 100644 --- a/.cmake/mito_tests_mito_lib.cmake +++ b/.cmake/mito_tests_mito_lib.cmake @@ -38,6 +38,7 @@ mito_test_driver(tests/mito.lib/geometry/spherical_metric_space.cc) # fem mito_test_driver(tests/mito.lib/fem/quadrature_field.cc) mito_test_driver(tests/mito.lib/fem/nodal_field.cc) +mito_test_driver(tests/mito.lib/fem/nodal_field_sphere.cc) # io mito_test_driver(tests/mito.lib/io/read_mesh_summit_2D.cc) diff --git a/lib/mito/io/externals.h b/lib/mito/io/externals.h index 663760eb..dcf10bf2 100644 --- a/lib/mito/io/externals.h +++ b/lib/mito/io/externals.h @@ -15,6 +15,7 @@ #include "../journal.h" #include "../geometry.h" #include "../mesh.h" +#include "../fem.h" // end of file diff --git a/lib/mito/io/vtk/externals.h b/lib/mito/io/vtk/externals.h index c4625e10..362a525c 100644 --- a/lib/mito/io/vtk/externals.h +++ b/lib/mito/io/vtk/externals.h @@ -11,6 +11,7 @@ #include #include #include +#include #include #include #include diff --git a/lib/mito/io/vtk/writer.h b/lib/mito/io/vtk/writer.h index 54f6364b..29d5bdd5 100644 --- a/lib/mito/io/vtk/writer.h +++ b/lib/mito/io/vtk/writer.h @@ -171,6 +171,37 @@ namespace mito::io::vtk { write(fileName, gridVtk); } + template + requires(cellT::dim == coordT::dim) + auto writer( + std::string fileName, const mito::mesh::mesh_t & mesh, + const geometry::coordinate_system_t & coordinate_system, + const fem::nodal_field_t & nodal_field + /*const auto fields::field_c & field*/) -> void + { + // create vtk unstructured grid + const auto gridVtk = createVtkUnstructuredGrid(mesh, coordinate_system); + + // initialize a vtk array + auto vtkArray = vtkSmartPointer::New(); + vtkArray->SetName(nodal_field.name().data()); + vtkArray->SetNumberOfComponents(Y::size); + vtkArray->SetNumberOfTuples(nodal_field.n_nodes()); + + // std::cout << gridVtk->GetPoints()->GetPoint(0)[0] << std::endl; + // std::cout << gridVtk->GetNumberOfPoints() << std::endl; + + for (int i = 0; i < gridVtk->GetNumberOfPoints(); ++i) { + vtkArray->SetTuple(i, nodal_field(i).begin()); + } + + // insert array into output mesh + gridVtk->GetPointData()->AddArray(vtkArray); + + // write grid to file + write(fileName, gridVtk); + } + } // namespace mito::io::vtk diff --git a/tests/mito.lib/fem/nodal_field_sphere.cc b/tests/mito.lib/fem/nodal_field_sphere.cc new file mode 100644 index 00000000..1df32061 --- /dev/null +++ b/tests/mito.lib/fem/nodal_field_sphere.cc @@ -0,0 +1,46 @@ +// -*- c++ -*- +// +// Copyright (c) 2020-2024, the MiTo Authors, all rights reserved +// + +#include +#include +#include +#include + + +// cartesian coordinates in 3D +using coordinates_t = mito::geometry::coordinates_t<3, mito::geometry::CARTESIAN>; + + +TEST(Fem, NodalFieldSphere) +{ + // the coordinate system + auto coord_system = mito::geometry::coordinate_system(); + + // read the mesh of a sphere + std::ifstream fileStream("sphere.summit"); + auto mesh = mito::io::summit::reader>(fileStream, coord_system); + + // say we have 319 nodes + int n_nodes = 319; + + // a field of mito vectors + auto field = mito::fem::nodal_field>(n_nodes, "field"); + + // populate the field + int i = 0; + for (auto & vector : field) { + if (i == 0) { + vector = mito::vector_t<3>{ 0.0, 0.0, 0.0 }; + } else { + vector = field(i - 1) + mito::vector_t<3>{ 1.0, 1.0, 1.0 }; + } + ++i; + } + +#ifdef WITH_VTK + // write mesh to vtk file + mito::io::vtk::writer("sphere_culo", mesh, coord_system, field); +#endif +} \ No newline at end of file From b77835119f8acee1c8acfd073fa89a976f5ae653 Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Sat, 25 May 2024 08:46:40 +0200 Subject: [PATCH 15/58] lib/mito/io: second draft implementation of writing a field to vtu file Here, we pass down to the writer a {fields::field} which is then evaluated on the nodes. --- lib/mito/io/vtk/writer.h | 22 +++++++++++---------- tests/mito.lib/fem/nodal_field_sphere.cc | 25 ++++++++---------------- 2 files changed, 20 insertions(+), 27 deletions(-) diff --git a/lib/mito/io/vtk/writer.h b/lib/mito/io/vtk/writer.h index 29d5bdd5..1bda32ae 100644 --- a/lib/mito/io/vtk/writer.h +++ b/lib/mito/io/vtk/writer.h @@ -171,28 +171,30 @@ namespace mito::io::vtk { write(fileName, gridVtk); } - template + template requires(cellT::dim == coordT::dim) auto writer( std::string fileName, const mito::mesh::mesh_t & mesh, const geometry::coordinate_system_t & coordinate_system, - const fem::nodal_field_t & nodal_field - /*const auto fields::field_c & field*/) -> void + const fields::field_c auto & field) -> void { // create vtk unstructured grid const auto gridVtk = createVtkUnstructuredGrid(mesh, coordinate_system); + auto n_nodes = gridVtk->GetNumberOfPoints(); + // initialize a vtk array auto vtkArray = vtkSmartPointer::New(); - vtkArray->SetName(nodal_field.name().data()); - vtkArray->SetNumberOfComponents(Y::size); - vtkArray->SetNumberOfTuples(nodal_field.n_nodes()); - - // std::cout << gridVtk->GetPoints()->GetPoint(0)[0] << std::endl; - // std::cout << gridVtk->GetNumberOfPoints() << std::endl; + vtkArray->SetName("name"); + vtkArray->SetNumberOfComponents(3); + vtkArray->SetNumberOfTuples(n_nodes); for (int i = 0; i < gridVtk->GetNumberOfPoints(); ++i) { - vtkArray->SetTuple(i, nodal_field(i).begin()); + auto point = gridVtk->GetPoints()->GetPoint(i); + auto coord = geometry::cartesian_coordinates_t<3>({ point[0], point[1], point[2] }); + + // std::cout << gridVtk->GetNumberOfPoints() << std::endl; + vtkArray->SetTuple(i, field(coord).begin()); } // insert array into output mesh diff --git a/tests/mito.lib/fem/nodal_field_sphere.cc b/tests/mito.lib/fem/nodal_field_sphere.cc index 1df32061..724119c3 100644 --- a/tests/mito.lib/fem/nodal_field_sphere.cc +++ b/tests/mito.lib/fem/nodal_field_sphere.cc @@ -22,25 +22,16 @@ TEST(Fem, NodalFieldSphere) std::ifstream fileStream("sphere.summit"); auto mesh = mito::io::summit::reader>(fileStream, coord_system); - // say we have 319 nodes - int n_nodes = 319; - - // a field of mito vectors - auto field = mito::fem::nodal_field>(n_nodes, "field"); - - // populate the field - int i = 0; - for (auto & vector : field) { - if (i == 0) { - vector = mito::vector_t<3>{ 0.0, 0.0, 0.0 }; - } else { - vector = field(i - 1) + mito::vector_t<3>{ 1.0, 1.0, 1.0 }; - } - ++i; - } + // the normal field to the submanifold + constexpr auto normal_field = mito::fields::field([](const coordinates_t & x) -> auto { + mito::scalar_t phi = std::atan2(x[1], x[0]); + mito::scalar_t theta = std::atan2(std::hypot(x[1], x[0]), x[2]); + return std::sin(theta) * std::cos(phi) * mito::e_0<3> + + std::sin(theta) * std::sin(phi) * mito::e_1<3> + std::cos(theta) * mito::e_2<3>; + }); #ifdef WITH_VTK // write mesh to vtk file - mito::io::vtk::writer("sphere_culo", mesh, coord_system, field); + mito::io::vtk::writer("sphere_field", mesh, coord_system, normal_field); #endif } \ No newline at end of file From 56f47ec5dc563e0e231909ffb68452defda29c83 Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Sat, 25 May 2024 09:04:57 +0200 Subject: [PATCH 16/58] lib/mito/io/vtk: remove unnecessary template specialization --- lib/mito/io/vtk/writer.h | 19 ++----------------- 1 file changed, 2 insertions(+), 17 deletions(-) diff --git a/lib/mito/io/vtk/writer.h b/lib/mito/io/vtk/writer.h index 1bda32ae..c48875bb 100644 --- a/lib/mito/io/vtk/writer.h +++ b/lib/mito/io/vtk/writer.h @@ -10,24 +10,9 @@ namespace mito::io::vtk { template - auto vtkCellPointer() -> vtkSmartPointer::type>; - - template <> - auto vtkCellPointer() -> vtkSmartPointer - { - return vtkSmartPointer::New(); - } - - template <> - auto vtkCellPointer() -> vtkSmartPointer - { - return vtkSmartPointer::New(); - } - - template <> - auto vtkCellPointer() -> vtkSmartPointer + auto vtkCellPointer() -> vtkSmartPointer::type> { - return vtkSmartPointer::New(); + return vtkSmartPointer::type>::New(); } template From 0684a60a3f1a2ff42c6e3856d7a0ef9373092c78 Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Sat, 25 May 2024 09:05:28 +0200 Subject: [PATCH 17/58] lib/mito/io/vtk: move {vtkCellPointer} to vtk cell header --- lib/mito/io/vtk/vtk_cell.h | 5 +++++ lib/mito/io/vtk/writer.h | 6 ------ 2 files changed, 5 insertions(+), 6 deletions(-) diff --git a/lib/mito/io/vtk/vtk_cell.h b/lib/mito/io/vtk/vtk_cell.h index a7f5bdbc..482b8dd1 100644 --- a/lib/mito/io/vtk/vtk_cell.h +++ b/lib/mito/io/vtk/vtk_cell.h @@ -28,6 +28,11 @@ namespace mito::io::vtk { using type = vtkTetra; }; + template + auto vtkCellPointer() -> vtkSmartPointer::type> + { + return vtkSmartPointer::type>::New(); + } } diff --git a/lib/mito/io/vtk/writer.h b/lib/mito/io/vtk/writer.h index c48875bb..fb6feb42 100644 --- a/lib/mito/io/vtk/writer.h +++ b/lib/mito/io/vtk/writer.h @@ -9,12 +9,6 @@ namespace mito::io::vtk { - template - auto vtkCellPointer() -> vtkSmartPointer::type> - { - return vtkSmartPointer::type>::New(); - } - template auto insertVtkPoint(const coordT &, vtkSmartPointer &) -> void; From f35cd45fe809a72045f52581bdf51f9d493406a6 Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Sat, 25 May 2024 10:01:25 +0200 Subject: [PATCH 18/58] lib/mito/fem: redesign class {NodalField} as a mapping between mesh nodes and field values --- .cmake/mito_tests_mito_lib.cmake | 1 - lib/mito/fem/NodalField.h | 67 ++++++++++++++---------- lib/mito/fem/api.h | 5 +- lib/mito/fem/externals.h | 1 + lib/mito/fem/factories.h | 6 +-- lib/mito/fem/forward.h | 6 +-- lib/mito/geometry/GeometricSimplex.h | 8 +++ lib/mito/mesh/forward.h | 8 +++ tests/mito.lib/fem/nodal_field.cc | 34 ------------ tests/mito.lib/fem/nodal_field_sphere.cc | 19 +++++-- 10 files changed, 80 insertions(+), 75 deletions(-) delete mode 100644 tests/mito.lib/fem/nodal_field.cc diff --git a/.cmake/mito_tests_mito_lib.cmake b/.cmake/mito_tests_mito_lib.cmake index da5783e0..576363c2 100644 --- a/.cmake/mito_tests_mito_lib.cmake +++ b/.cmake/mito_tests_mito_lib.cmake @@ -37,7 +37,6 @@ mito_test_driver(tests/mito.lib/geometry/spherical_metric_space.cc) # fem mito_test_driver(tests/mito.lib/fem/quadrature_field.cc) -mito_test_driver(tests/mito.lib/fem/nodal_field.cc) mito_test_driver(tests/mito.lib/fem/nodal_field_sphere.cc) # io diff --git a/lib/mito/fem/NodalField.h b/lib/mito/fem/NodalField.h index f6de87ae..229e8807 100644 --- a/lib/mito/fem/NodalField.h +++ b/lib/mito/fem/NodalField.h @@ -9,44 +9,57 @@ namespace mito::fem { - template + template class NodalField { private: - // conventionally packed grid for {n} - using pack_type = pyre::grid::canonical_t<1>; - // of Y on the heap - using storage_type = pyre::memory::heap_t; - // putting it all together - using grid_type = pyre::grid::grid_t; - // index - using index_type = pack_type::index_type; + // the node type + using node_type = geometry::node_t; + + // hash function for nodes + struct hash_function { + auto operator()(const node_type & node) const -> uintptr_t + { + // reinterpret the address of the node as a {uintptr_t} and return it + return uintptr_t(&node); + } + }; + + // a map from {key_type} to {Y} values + using map_type = std::unordered_map; public: // constructor - NodalField(int nodes, std::string name) : NodalField(pack_type{ { nodes } }, name) {} + NodalField(const mesh::mesh_c auto & mesh, std::string name) : + _map_nodes_to_values(), + _name(name) + { + // populate the map with the mesh nodes (map all nodes to default item {Y}) + for (const auto & cell : mesh.cells()) { + for (const auto & node : cell.nodes()) { + _map_nodes_to_values[node] = Y(); + } + } + } // destructor ~NodalField() = default; - private: - inline NodalField(const pack_type && packing, std::string name) : - _grid{ packing, packing.cells() }, - _name(name) - {} - public: /** * Operator() */ - inline auto operator()(const int a) const -> const Y & { return _grid[a]; } + inline auto operator()(const node_type & node) const -> const Y & + { + return _map_nodes_to_values[node]; + } - inline auto operator()(const int a) -> Y & { return _grid[a]; } + inline auto operator()(const node_type & node) -> Y & { return _map_nodes_to_values[node]; } /** * accessor for the number of nodes */ - inline auto n_nodes() const -> int { return _grid.layout().shape()[0]; } + inline auto n_nodes() const -> int { return _map_nodes_to_values.size(); } /** * accessor for name @@ -54,21 +67,21 @@ namespace mito::fem { inline const std::string & name() const noexcept { return _name; } // support for ranged for loops (wrapping grid) - inline auto begin() const { return std::cbegin(_grid); } - inline auto end() const { return std::cend(_grid); } - inline auto begin() { return std::begin(_grid); } - inline auto end() { return std::end(_grid); } + inline auto begin() const { return std::cbegin(_map_nodes_to_values); } + inline auto end() const { return std::cend(_map_nodes_to_values); } + inline auto begin() { return std::begin(_map_nodes_to_values); } + inline auto end() { return std::end(_map_nodes_to_values); } private: - // the underlying grid - grid_type _grid; + // the underlying mapping of nodes to nodal values + map_type _map_nodes_to_values; // the name of the field std::string _name; }; - template - std::ostream & operator<<(std::ostream & os, const NodalField & nodalField) + template + std::ostream & operator<<(std::ostream & os, const NodalField & nodalField) { os << "Nodal field \"" << nodalField.name() << "\" : " << std::endl; diff --git a/lib/mito/fem/api.h b/lib/mito/fem/api.h index 3951942c..5c43f432 100644 --- a/lib/mito/fem/api.h +++ b/lib/mito/fem/api.h @@ -14,9 +14,8 @@ namespace mito::fem { constexpr auto quadrature_field(int nElements, std::string name = ""); // nodal field factory - template - constexpr auto nodal_field(int nodes, std::string name = ""); - + template + constexpr auto nodal_field(const mesh::mesh_c auto & mesh, std::string name); } diff --git a/lib/mito/fem/externals.h b/lib/mito/fem/externals.h index 6aa6c979..33a95cf5 100644 --- a/lib/mito/fem/externals.h +++ b/lib/mito/fem/externals.h @@ -16,6 +16,7 @@ // support #include "../journal.h" #include "../tensor.h" +#include "../mesh.h" // end of file diff --git a/lib/mito/fem/factories.h b/lib/mito/fem/factories.h index 32d29e67..c67ee923 100644 --- a/lib/mito/fem/factories.h +++ b/lib/mito/fem/factories.h @@ -17,10 +17,10 @@ namespace mito::fem { } // nodal field factory - template - constexpr auto nodal_field(int nodes, std::string name) + template + constexpr auto nodal_field(const mesh::mesh_c auto & mesh, std::string name) { - return nodal_field_t(nodes, name); + return nodal_field_t(mesh, name); } } diff --git a/lib/mito/fem/forward.h b/lib/mito/fem/forward.h index 7b94e071..af28aacf 100644 --- a/lib/mito/fem/forward.h +++ b/lib/mito/fem/forward.h @@ -10,12 +10,12 @@ namespace mito::fem { // class nodal field - template + template class NodalField; // nodal field alias - template - using nodal_field_t = NodalField; + template + using nodal_field_t = NodalField; // class quadrature field template diff --git a/lib/mito/geometry/GeometricSimplex.h b/lib/mito/geometry/GeometricSimplex.h index 6fc0a2fd..e4e7e152 100644 --- a/lib/mito/geometry/GeometricSimplex.h +++ b/lib/mito/geometry/GeometricSimplex.h @@ -202,6 +202,14 @@ namespace mito::geometry { point_type _point; }; + // operator== for nodes + template + constexpr auto operator==( + const GeometricSimplex<0, D> & node_a, const GeometricSimplex<0, D> & node_b) -> bool + { + // two nodes are the same if they live at the same address + return &node_a == &node_b; + } } diff --git a/lib/mito/mesh/forward.h b/lib/mito/mesh/forward.h index 9591aa78..d152e30d 100644 --- a/lib/mito/mesh/forward.h +++ b/lib/mito/mesh/forward.h @@ -17,6 +17,14 @@ namespace mito::mesh { template using mesh_t = Mesh; + // concept of a mesh + template + concept mesh_c = requires(F c) { + // require that F only binds to {Mesh} specializations + [](const mesh_t &) { + }(c); + }; + // class boundary template class Boundary; diff --git a/tests/mito.lib/fem/nodal_field.cc b/tests/mito.lib/fem/nodal_field.cc deleted file mode 100644 index 5e71d407..00000000 --- a/tests/mito.lib/fem/nodal_field.cc +++ /dev/null @@ -1,34 +0,0 @@ -// -*- c++ -*- -// -// Copyright (c) 2020-2024, the MiTo Authors, all rights reserved -// - -#include -#include - - -TEST(Fem, NodalField) -{ - // say we have 20 nodes - int n_nodes = 20; - - // a field of mito vectors - auto field = mito::fem::nodal_field>(n_nodes, "field"); - - // expect that the number of nodes is 20 - EXPECT_TRUE(field.n_nodes() == n_nodes); - - // populate the field - int i = 0; - for (auto & vector : field) { - if (i == 0) { - vector = mito::vector_t<2>{ 0.0, 0.0 }; - } else { - vector = field(i - 1) + mito::vector_t<2>{ 1, 1 }; - } - ++i; - } - - // expect that editing {vector} has edited the {field} - EXPECT_TRUE((field(10) == mito::vector_t<2>{ 10.0, 10.0 })); -} \ No newline at end of file diff --git a/tests/mito.lib/fem/nodal_field_sphere.cc b/tests/mito.lib/fem/nodal_field_sphere.cc index 724119c3..d9eb8ac9 100644 --- a/tests/mito.lib/fem/nodal_field_sphere.cc +++ b/tests/mito.lib/fem/nodal_field_sphere.cc @@ -22,6 +22,9 @@ TEST(Fem, NodalFieldSphere) std::ifstream fileStream("sphere.summit"); auto mesh = mito::io::summit::reader>(fileStream, coord_system); + // a nodal field on the mesh + auto nodal_field = mito::fem::nodal_field, 3>(mesh, "normal"); + // the normal field to the submanifold constexpr auto normal_field = mito::fields::field([](const coordinates_t & x) -> auto { mito::scalar_t phi = std::atan2(x[1], x[0]); @@ -30,8 +33,16 @@ TEST(Fem, NodalFieldSphere) + std::sin(theta) * std::sin(phi) * mito::e_1<3> + std::cos(theta) * mito::e_2<3>; }); -#ifdef WITH_VTK - // write mesh to vtk file - mito::io::vtk::writer("sphere_field", mesh, coord_system, normal_field); -#endif + // fill information in nodal field + for (auto & [node, value] : nodal_field) { + // get the coordinates of the node + auto & coordinates = coord_system.coordinates(node.point()); + // compute the value of the normal field at those coordinates + value = normal_field(coordinates); + } + + // #ifdef WITH_VTK + // // write mesh to vtk file + // mito::io::vtk::writer("sphere_field", mesh, coord_system, normal_field); + // #endif } \ No newline at end of file From 839873364909006f6e613ad935449ed122acbdf0 Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Sat, 25 May 2024 10:09:05 +0200 Subject: [PATCH 19/58] lib/mito/io/vtk: add possibility to write 1D points --- lib/mito/io/vtk/writer.h | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/lib/mito/io/vtk/writer.h b/lib/mito/io/vtk/writer.h index fb6feb42..572a6b69 100644 --- a/lib/mito/io/vtk/writer.h +++ b/lib/mito/io/vtk/writer.h @@ -30,6 +30,15 @@ namespace mito::io::vtk { pointsVtk->InsertNextPoint(coord[0], coord[1], 0.); } + template <> + auto insertVtkPoint( + const geometry::coordinates_t<1, geometry::CARTESIAN> & coord, + vtkSmartPointer & pointsVtk) -> void + { + // add the point as new vtk point + pointsVtk->InsertNextPoint(coord[0], 0., 0.); + } + template auto createVtkUnstructuredGrid( const mito::mesh::mesh_t & mesh, From 8a652bb8eada0fcf08ca07079ce34d854fcc82fc Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Sun, 26 May 2024 20:25:26 +0200 Subject: [PATCH 20/58] lib/mito/fem: bugfix Fixed bug for which nodes in {NodalField} were duplicated across adjacent elements. When nodes are inserted in a mesh with mesh.insert({node0, node1, node2}), nodes are copied. Therefore, a hash function and {operator==} based on the nodes addresses will create copies for each instance of a node in a mesh (that is for the same node appearing in adjacent elements). For now, the problem is fixed by setting nodes equal if their underlying vertex and points are equal and by editing the hash function to trivially combine the hashes of vertex and point (this needs to be improved to reduce the chance of clash). However, these choices will be revisited when we implement discontinuous Galerkin. --- lib/mito/fem/NodalField.h | 8 ++++++-- lib/mito/geometry/GeometricSimplex.h | 4 ++-- 2 files changed, 8 insertions(+), 4 deletions(-) diff --git a/lib/mito/fem/NodalField.h b/lib/mito/fem/NodalField.h index 229e8807..d81143b5 100644 --- a/lib/mito/fem/NodalField.h +++ b/lib/mito/fem/NodalField.h @@ -20,8 +20,12 @@ namespace mito::fem { struct hash_function { auto operator()(const node_type & node) const -> uintptr_t { - // reinterpret the address of the node as a {uintptr_t} and return it - return uintptr_t(&node); + // TOFIX: this hash function is not safe against collisions: either figure out a + // proper hash function or turn nodes into something shareable, so that the hash + // can simply be a check on the address of the underlying resource. We'll figure + // out which one of two routes to pursue when we do discontinuous Galerkin + // build a hash based on the id of the node and that of the point + return node.vertex().id() + node.point().id(); } }; diff --git a/lib/mito/geometry/GeometricSimplex.h b/lib/mito/geometry/GeometricSimplex.h index e4e7e152..601c616c 100644 --- a/lib/mito/geometry/GeometricSimplex.h +++ b/lib/mito/geometry/GeometricSimplex.h @@ -207,8 +207,8 @@ namespace mito::geometry { constexpr auto operator==( const GeometricSimplex<0, D> & node_a, const GeometricSimplex<0, D> & node_b) -> bool { - // two nodes are the same if they live at the same address - return &node_a == &node_b; + // two nodes are the same if their vertex and point are the same + return node_a.vertex() == node_b.vertex() && node_a.point() == node_b.point(); } } From 204b39bb364dcc59cda6760d9a479f198213dae1 Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Sun, 26 May 2024 20:25:57 +0200 Subject: [PATCH 21/58] lib/mito/geometry: publish dimension of physical space in class {PointCloud} --- lib/mito/geometry/PointCloud.h | 3 +++ 1 file changed, 3 insertions(+) diff --git a/lib/mito/geometry/PointCloud.h b/lib/mito/geometry/PointCloud.h index 5ecf1fc5..36f0021a 100644 --- a/lib/mito/geometry/PointCloud.h +++ b/lib/mito/geometry/PointCloud.h @@ -11,6 +11,9 @@ namespace mito::geometry { template class PointCloud { + public: + static constexpr int dim = D; + private: // a point using point_type = point_t; From 616b178bd3544a5f22783607aef95431e93504c0 Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Sun, 26 May 2024 20:27:03 +0200 Subject: [PATCH 22/58] lib/mito/geometry: add concepts for a coordinate system and a point cloud --- lib/mito/geometry/forward.h | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/lib/mito/geometry/forward.h b/lib/mito/geometry/forward.h index 7da28690..aed99cc6 100644 --- a/lib/mito/geometry/forward.h +++ b/lib/mito/geometry/forward.h @@ -17,6 +17,14 @@ namespace mito::geometry { template class CoordinateSystem; + // concept of a coordinate system + template + concept coordinate_system_c = requires(F c) { + // require that F only binds to {CoordinateSystem} specializations + [](const CoordinateSystem &) { + }(c); + }; + // class point template class Point; @@ -29,6 +37,14 @@ namespace mito::geometry { // class point cloud template class PointCloud; + + // concept of a point cloud + template + concept point_cloud_c = requires(F c) { + // require that F only binds to {PointCloud} specializations + [](const PointCloud &) { + }(c); + }; } From 1f299a6aa138918068ae07c7d594f2d1623b379a Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Sun, 26 May 2024 21:00:45 +0200 Subject: [PATCH 23/58] lib/mito/io/vtk: split file {writer.h} into {MeshWriter} and {PointCloudWriter} classes This is necessary because writing nodal fields requires to store the vtu grid as state. --- lib/mito/io/vtk/MeshWriter.h | 163 ++++++++++++++ lib/mito/io/vtk/PointCloudWriter.h | 145 +++++++++++++ lib/mito/io/vtk/api.h | 36 ++++ lib/mito/io/vtk/factories.h | 32 +++ lib/mito/io/vtk/forward.h | 24 +++ lib/mito/io/vtk/public.h | 17 +- lib/mito/io/vtk/vtk_point.h | 45 ++++ lib/mito/io/vtk/writer.h | 198 ------------------ tests/mito.lib/fem/nodal_field_sphere.cc | 12 +- .../integration/write_tetra_mesh_to_vtk.cc | 2 +- tests/mito.lib/io/summit_to_summit_mesh_2D.cc | 2 +- tests/mito.lib/io/summit_to_vtk_mesh_2D.cc | 2 +- tests/mito.lib/io/summit_to_vtk_mesh_3D.cc | 2 +- tests/mito.lib/io/write_point_cloud_vtk.cc | 4 +- tests/mito.lib/mesh/ball.cc | 2 +- tests/mito.lib/mesh/filter_ball.cc | 6 +- tests/mito.lib/mesh/sphere.cc | 2 +- 17 files changed, 479 insertions(+), 215 deletions(-) create mode 100644 lib/mito/io/vtk/MeshWriter.h create mode 100644 lib/mito/io/vtk/PointCloudWriter.h create mode 100644 lib/mito/io/vtk/api.h create mode 100644 lib/mito/io/vtk/factories.h create mode 100644 lib/mito/io/vtk/forward.h create mode 100644 lib/mito/io/vtk/vtk_point.h delete mode 100644 lib/mito/io/vtk/writer.h diff --git a/lib/mito/io/vtk/MeshWriter.h b/lib/mito/io/vtk/MeshWriter.h new file mode 100644 index 00000000..8603c2ff --- /dev/null +++ b/lib/mito/io/vtk/MeshWriter.h @@ -0,0 +1,163 @@ +// -*- c++ -*- +// +// Copyright (c) 2020-2024, the MiTo Authors, all rights reserved +// + +// code guard +#pragma once + + +namespace mito::io::vtk { + + template + requires(meshT::cell_type::dim == coordSystemT::coordinates_type::dim) + class MeshWriter { + + private: + // the mesh type + using mesh_type = meshT; + // the coordinate system type + using coord_system_type = coordSystemT; + // the dimension of the physical space + static constexpr int D = mesh_type::cell_type::dim; + // the type of grid + using grid_type = vtkSmartPointer; + // map mesh points to the index of the vtk points. Points that are shared among + // multiple elements have the same index. + using point_to_index_type = std::unordered_map< + geometry::point_t, int, utilities::hash_function>>; + + private: + auto _create_vtk_grid(const mesh_type & mesh, const coord_system_type & coordinate_system) + -> grid_type + { + // vtk unstructured grid + auto gridVtk = vtkSmartPointer::New(); + // vtk points and cells + auto pointsVtk = vtkSmartPointer::New(); + + // global index assigned to each vtk point + auto indexPointVtk = 0; + + // loop over the cells + for (const auto & cell : mesh.cells()) { + + // create vtk cell + auto cellVtk = vtkCellPointer(); + + // local index for the points of the cell + auto indexLocalPointVtk = 0; + + // loop over the nodes of the cell + for (const auto & node : cell.nodes()) { + // retrieve the corresponding point + const auto & point = node.point(); + // if the point is not present in the map + if (!_point_to_index.contains(point)) { + // insert the new vtk point + insert_vtk_point(coordinate_system.coordinates(point), pointsVtk); + // add the point to the map with its global index + _point_to_index[point] = indexPointVtk; + // update global index for the vtk point + ++indexPointVtk; + } + // set the id of the point + cellVtk->GetPointIds()->SetId(indexLocalPointVtk, _point_to_index[point]); + // update local index for the points in the cell + ++indexLocalPointVtk; + } + + // insert the new cell + gridVtk->InsertNextCell(cellVtk->GetCellType(), cellVtk->GetPointIds()); + } + + // set the grid points + gridVtk->SetPoints(pointsVtk); + + // all done + return gridVtk; + } + + public: + MeshWriter( + std::string filename, const mesh_type & mesh, const coord_system_type & coord_system) : + _filename(filename), + _point_to_index(), + _grid(_create_vtk_grid(mesh, coord_system)) + {} + + // TOFIX: use concepts to say that Y is a tensor + template + auto record(const fem::nodal_field_t & field, std::string fieldname = "") -> void + { + // if no name was provided + if (fieldname == "") { + // use the name of the field + fieldname = field.name(); + } + + // get the number of nodes + auto n_nodes = field.n_nodes(); + + // check the number of nodes in the field equals the number of points in the grid + std::cout << n_nodes << "\t" << _grid->GetNumberOfPoints() << std::endl; + // assert(n_nodes == _grid->GetNumberOfPoints()); + // and equal to the number of indices that we have recorded so far + // assert(n_nodes == int(std::size(_point_to_index))); + std::cout << n_nodes << "\t" << int(std::size(_point_to_index)) << std::endl; + + // initialize a vtk array + auto vtkArray = vtkSmartPointer::New(); + vtkArray->SetName(fieldname.data()); + vtkArray->SetNumberOfComponents(Y::size); + vtkArray->SetNumberOfTuples(field.n_nodes()); + + // populate the array with the nodal values + for (auto & [node, value] : field) { + // get the index corresponding to the current point + auto i = _point_to_index.at(node.point()); + vtkArray->SetTuple(i, value.begin()); + } + + // insert array into output mesh + _grid->GetPointData()->AddArray(vtkArray); + + // all done + return; + } + + auto write() const -> void + { + // create a new writer + auto writer = vtkSmartPointer::New(); + // set the name of the output file + writer->SetFileName((_filename + ".vtu").data()); + + // sign the grid up for writing +#if VTK_MAJOR_VERSION <= 8 + writer->SetInput(_grid); +#else + writer->SetInputData(_grid); +#endif + // write the grid to file + writer->Write(); + + // all done + return; + } + + private: + // the name of the output file to be written + std::string _filename; + + // a map storing point -> index relation + point_to_index_type _point_to_index; + + // the grid + grid_type _grid; + }; + +} // namespace mito::io::vtk + + +// end of file diff --git a/lib/mito/io/vtk/PointCloudWriter.h b/lib/mito/io/vtk/PointCloudWriter.h new file mode 100644 index 00000000..b654a39a --- /dev/null +++ b/lib/mito/io/vtk/PointCloudWriter.h @@ -0,0 +1,145 @@ +// -*- c++ -*- +// +// Copyright (c) 2020-2024, the MiTo Authors, all rights reserved +// + +// code guard +#pragma once + + +namespace mito::io::vtk { + + template + requires(cloudT::dim == coordSystemT::coordinates_type::dim) + class PointCloudWriter { + + private: + // the point cloud type + using cloud_type = cloudT; + // the coordinate system type + using coord_system_type = coordSystemT; + // the dimension of the physical space + static constexpr int D = cloudT::dim; + // the type of grid + using grid_type = vtkSmartPointer; + // map mesh points to the index of the vtk points. Points that are shared among + // multiple elements have the same index. + using point_to_index_type = std::unordered_map< + geometry::point_t, int, utilities::hash_function>>; + + private: + auto _create_vtk_grid(const cloud_type & cloud, const coord_system_type & coordinate_system) + -> grid_type + { + // vtk unstructured grid + auto gridVtk = vtkSmartPointer::New(); + // vtk points and cells + auto pointsVtk = vtkSmartPointer::New(); + + // global index assigned to each vtk point + auto indexPointVtk = 0; + + // get point cloud compositions + const auto & points = cloud.points(); + + // iterate over the points + for (const auto & point : points) { + // if the point is not present in the map + if (!_point_to_index.contains(point)) { + insert_vtk_point(coordinate_system.coordinates(point), pointsVtk); + // add the point to the map with its global index + _point_to_index[point] = indexPointVtk; + // update global index for the vtk point + ++indexPointVtk; + } + } + + // set the grid points + gridVtk->SetPoints(pointsVtk); + + // all done + return gridVtk; + } + + public: + PointCloudWriter( + std::string filename, const cloud_type & cloud, + const coord_system_type & coord_system) : + _filename(filename), + _point_to_index(), + _grid(_create_vtk_grid(cloud, coord_system)) + {} + + // TOFIX: use concepts to say that Y is a tensor + template + auto record(const fem::nodal_field_t & field, std::string fieldname = "") -> void + { + // if no name was provided + if (fieldname == "") { + // use the name of the field + fieldname = field.name(); + } + + // get the number of nodes + auto n_nodes = field.n_nodes(); + + // check the number of nodes in the field equals the number of points in the grid + assert(n_nodes == _grid->GetNumberOfPoints()); + // and equal to the number of indices that we have recorded so far + assert(n_nodes == int(std::size(_point_to_index))); + + // initialize a vtk array + auto vtkArray = vtkSmartPointer::New(); + vtkArray->SetName(fieldname.data()); + vtkArray->SetNumberOfComponents(Y::size); + vtkArray->SetNumberOfTuples(field.n_nodes()); + + // populate the array with the nodal values + for (auto & [node, value] : field) { + // get the index corresponding to the current point + auto i = _point_to_index.at(node.point()); + vtkArray->SetTuple(i, value.begin()); + } + + // insert array into output mesh + _grid->GetPointData()->AddArray(vtkArray); + + // all done + return; + } + + auto write() const -> void + { + // create a new writer + auto writer = vtkSmartPointer::New(); + // set the name of the output file + writer->SetFileName((_filename + ".vtu").data()); + + // sign the grid up for writing +#if VTK_MAJOR_VERSION <= 8 + writer->SetInput(_grid); +#else + writer->SetInputData(_grid); +#endif + // write the grid to file + writer->Write(); + + // all done + return; + } + + private: + // the name of the output file to be written + std::string _filename; + + // a map storing point -> index relation + point_to_index_type _point_to_index; + + // the grid + grid_type _grid; + }; + +} // namespace mito::io::vtk + + +// end of file diff --git a/lib/mito/io/vtk/api.h b/lib/mito/io/vtk/api.h new file mode 100644 index 00000000..c68cd80d --- /dev/null +++ b/lib/mito/io/vtk/api.h @@ -0,0 +1,36 @@ +// -*- c++ -*- +// +// Copyright (c) 2020-2024, the MiTo Authors, all rights reserved +// + +// code guard +#pragma once + + +namespace mito::io::vtk { + + // mesh writer alias + template + requires(meshT::cell_type::dim == coordSystemT::coordinates_type::dim) + using mesh_writer_t = MeshWriter; + + // vtk mesh writer factory + template + requires(meshT::cell_type::dim == coordSystemT::coordinates_type::dim) + auto writer(std::string filename, const meshT & mesh, const coordSystemT & coord_system) + -> mesh_writer_t; + + // point cloud writer alias + template + requires(cloudT::dim == coordSystemT::coordinates_type::dim) + using cloud_writer_t = PointCloudWriter; + + // point cloud writer factory + template + requires(cloudT::dim == coordSystemT::coordinates_type::dim) + auto writer(std::string filename, const cloudT & cloud, const coordSystemT & coord_system) + -> cloud_writer_t; +} + + +// end of file diff --git a/lib/mito/io/vtk/factories.h b/lib/mito/io/vtk/factories.h new file mode 100644 index 00000000..21d48a20 --- /dev/null +++ b/lib/mito/io/vtk/factories.h @@ -0,0 +1,32 @@ +// -*- c++ -*- +// +// Copyright (c) 2020-2024, the MiTo Authors, all rights reserved +// + +// code guard +#pragma once + + +namespace mito::io::vtk { + + // vtk mesh writer factory + template + requires(meshT::cell_type::dim == coordSystemT::coordinates_type::dim) + auto writer(std::string filename, const meshT & mesh, const coordSystemT & coord_system) + -> mesh_writer_t + { + return MeshWriter(filename, mesh, coord_system); + } + + // vtk point cloud writer factory + template + requires(cloudT::dim == coordSystemT::coordinates_type::dim) + auto writer(std::string filename, const cloudT & cloud, const coordSystemT & coord_system) + -> cloud_writer_t + { + return PointCloudWriter(filename, cloud, coord_system); + } +} + + +// end of file diff --git a/lib/mito/io/vtk/forward.h b/lib/mito/io/vtk/forward.h new file mode 100644 index 00000000..8b478f15 --- /dev/null +++ b/lib/mito/io/vtk/forward.h @@ -0,0 +1,24 @@ +// -*- c++ -*- +// +// Copyright (c) 2020-2024, the MiTo Authors, all rights reserved +// + +// code guard +#pragma once + + +namespace mito::io::vtk { + + // class mesh writer + template + requires(meshT::cell_type::dim == coordSystemT::coordinates_type::dim) + class MeshWriter; + + // class point cloud writer + template + requires(cloudT::dim == coordSystemT::coordinates_type::dim) + class PointCloudWriter; +} + + +// end of file diff --git a/lib/mito/io/vtk/public.h b/lib/mito/io/vtk/public.h index c7956cd0..c1495e68 100644 --- a/lib/mito/io/vtk/public.h +++ b/lib/mito/io/vtk/public.h @@ -10,9 +10,22 @@ // external packages #include "externals.h" -// classes implementation +// get the forward declarations +#include "forward.h" + +// published type factories; this is the file you are looking for... +#include "api.h" + +// wrappers #include "vtk_cell.h" -#include "writer.h" +#include "vtk_point.h" + +// classes implementation +#include "MeshWriter.h" +#include "PointCloudWriter.h" + +// factories implementation +#include "factories.h" // end of file diff --git a/lib/mito/io/vtk/vtk_point.h b/lib/mito/io/vtk/vtk_point.h new file mode 100644 index 00000000..54bd4380 --- /dev/null +++ b/lib/mito/io/vtk/vtk_point.h @@ -0,0 +1,45 @@ +// -*- c++ -*- +// +// Copyright (c) 2020-2024, the MiTo Authors, all rights reserved +// + +// code guard +#pragma once + + +namespace mito::io::vtk { + + template + auto insert_vtk_point(const coordT &, vtkSmartPointer &) -> void; + + template <> + auto insert_vtk_point( + const geometry::coordinates_t<3, geometry::CARTESIAN> & coord, + vtkSmartPointer & pointsVtk) -> void + { + // add the point as new vtk point + pointsVtk->InsertNextPoint(coord[0], coord[1], coord[2]); + } + + template <> + auto insert_vtk_point( + const geometry::coordinates_t<2, geometry::CARTESIAN> & coord, + vtkSmartPointer & pointsVtk) -> void + { + // add the point as new vtk point + pointsVtk->InsertNextPoint(coord[0], coord[1], 0.); + } + + template <> + auto insert_vtk_point( + const geometry::coordinates_t<1, geometry::CARTESIAN> & coord, + vtkSmartPointer & pointsVtk) -> void + { + // add the point as new vtk point + pointsVtk->InsertNextPoint(coord[0], 0., 0.); + } + +} // namespace mito::io::vtk + + +// end of file diff --git a/lib/mito/io/vtk/writer.h b/lib/mito/io/vtk/writer.h deleted file mode 100644 index 572a6b69..00000000 --- a/lib/mito/io/vtk/writer.h +++ /dev/null @@ -1,198 +0,0 @@ -// -*- c++ -*- -// -// Copyright (c) 2020-2024, the MiTo Authors, all rights reserved -// - -// code guard -#pragma once - - -namespace mito::io::vtk { - - template - auto insertVtkPoint(const coordT &, vtkSmartPointer &) -> void; - - template <> - auto insertVtkPoint( - const geometry::coordinates_t<3, geometry::CARTESIAN> & coord, - vtkSmartPointer & pointsVtk) -> void - { - // add the point as new vtk point - pointsVtk->InsertNextPoint(coord[0], coord[1], coord[2]); - } - - template <> - auto insertVtkPoint( - const geometry::coordinates_t<2, geometry::CARTESIAN> & coord, - vtkSmartPointer & pointsVtk) -> void - { - // add the point as new vtk point - pointsVtk->InsertNextPoint(coord[0], coord[1], 0.); - } - - template <> - auto insertVtkPoint( - const geometry::coordinates_t<1, geometry::CARTESIAN> & coord, - vtkSmartPointer & pointsVtk) -> void - { - // add the point as new vtk point - pointsVtk->InsertNextPoint(coord[0], 0., 0.); - } - - template - auto createVtkUnstructuredGrid( - const mito::mesh::mesh_t & mesh, - const geometry::coordinate_system_t & coordinate_system) - -> vtkSmartPointer - requires(cellT::dim == coordT::dim) - { - // the dimension of physical space - constexpr int D = cellT::dim; - - // vtk unstructured grid - auto gridVtk = vtkSmartPointer::New(); - // vtk points and cells - auto pointsVtk = vtkSmartPointer::New(); - - // map mesh points to the index of the vtk points. Points that are shared among - // multiple elements have the same index. - std::unordered_map< - geometry::point_t, int, utilities::hash_function>> - mapPoints; - - // global index assigned to each vtk point - auto indexPointVtk = 0; - - // loop over the cells - for (const auto & cell : mesh.cells()) { - // create vtk cell - auto cellVtk = vtkCellPointer(); - - // local index for the points of the cell - auto indexLocalPointVtk = 0; - - // retrieve nodes of the cell - auto nodes = cell.nodes(); - - // loop over the nodes of the cell - for (const auto & node : nodes) { - // retrieve the corresponding point - const auto pPoint = node.point(); - // if the point is not present in the map - if (mapPoints.count(pPoint) == 0) { - // insert the new vtk point - insertVtkPoint(coordinate_system.coordinates(pPoint), pointsVtk); - // add the point to the map with its global index - mapPoints[pPoint] = indexPointVtk; - // update global index for the vtk point - ++indexPointVtk; - } - // set the id of the point - cellVtk->GetPointIds()->SetId(indexLocalPointVtk, mapPoints[pPoint]); - // update local index for the points in the cell - ++indexLocalPointVtk; - } - - // insert the new cell - gridVtk->InsertNextCell(cellVtk->GetCellType(), cellVtk->GetPointIds()); - } - - // set the grid points - gridVtk->SetPoints(pointsVtk); - - return gridVtk; - } - - inline auto write(std::string fileName, const vtkSmartPointer & gridVtk) - -> void - { - // write the grid to file - auto writer = vtkSmartPointer::New(); - writer->SetFileName((fileName + ".vtu").c_str()); - -#if VTK_MAJOR_VERSION <= 8 - writer->SetInput(gridVtk); -#else - writer->SetInputData(gridVtk); -#endif - - writer->Write(); - } - - template - requires(cellT::dim == coordT::dim) - auto writer( - std::string fileName, const mito::mesh::mesh_t & mesh, - const geometry::coordinate_system_t & coordinate_system) -> void - { - // create vtk unstructured grid - const auto gridVtk = createVtkUnstructuredGrid(mesh, coordinate_system); - - // write grid to file - write(fileName, gridVtk); - } - - - template - requires(D == coordT::dim) - auto writer( - std::string fileName, const geometry::point_cloud_t & cloud, - const geometry::coordinate_system_t & coordinate_system) -> void - { - // get point cloud compositions - const auto & points = cloud.points(); - - // vtk unstructured grid - auto gridVtk = vtkSmartPointer::New(); - // vtk points and cells - auto pointsVtk = vtkSmartPointer::New(); - - // iterate over the points - for (const auto & point : points) { - insertVtkPoint(coordinate_system.coordinates(point), pointsVtk); - } - - // set the grid points - gridVtk->SetPoints(pointsVtk); - - // write grid to file - write(fileName, gridVtk); - } - - template - requires(cellT::dim == coordT::dim) - auto writer( - std::string fileName, const mito::mesh::mesh_t & mesh, - const geometry::coordinate_system_t & coordinate_system, - const fields::field_c auto & field) -> void - { - // create vtk unstructured grid - const auto gridVtk = createVtkUnstructuredGrid(mesh, coordinate_system); - - auto n_nodes = gridVtk->GetNumberOfPoints(); - - // initialize a vtk array - auto vtkArray = vtkSmartPointer::New(); - vtkArray->SetName("name"); - vtkArray->SetNumberOfComponents(3); - vtkArray->SetNumberOfTuples(n_nodes); - - for (int i = 0; i < gridVtk->GetNumberOfPoints(); ++i) { - auto point = gridVtk->GetPoints()->GetPoint(i); - auto coord = geometry::cartesian_coordinates_t<3>({ point[0], point[1], point[2] }); - - // std::cout << gridVtk->GetNumberOfPoints() << std::endl; - vtkArray->SetTuple(i, field(coord).begin()); - } - - // insert array into output mesh - gridVtk->GetPointData()->AddArray(vtkArray); - - // write grid to file - write(fileName, gridVtk); - } - -} // namespace mito::io::vtk - - -// end of file diff --git a/tests/mito.lib/fem/nodal_field_sphere.cc b/tests/mito.lib/fem/nodal_field_sphere.cc index d9eb8ac9..9dc99c1a 100644 --- a/tests/mito.lib/fem/nodal_field_sphere.cc +++ b/tests/mito.lib/fem/nodal_field_sphere.cc @@ -41,8 +41,12 @@ TEST(Fem, NodalFieldSphere) value = normal_field(coordinates); } - // #ifdef WITH_VTK - // // write mesh to vtk file - // mito::io::vtk::writer("sphere_field", mesh, coord_system, normal_field); - // #endif +#ifdef WITH_VTK + // write mesh to vtk file + auto writer = mito::io::vtk::writer("sphere_field", mesh, coord_system); + // sign {nodal_field} up with the writer + writer.record(nodal_field); + // write output file + writer.write(); +#endif } \ No newline at end of file diff --git a/tests/mito.lib/integration/write_tetra_mesh_to_vtk.cc b/tests/mito.lib/integration/write_tetra_mesh_to_vtk.cc index f42ac9ca..26092b25 100644 --- a/tests/mito.lib/integration/write_tetra_mesh_to_vtk.cc +++ b/tests/mito.lib/integration/write_tetra_mesh_to_vtk.cc @@ -32,5 +32,5 @@ TEST(VtkWriter, WriteTetraMeshToVtk) auto tetra_mesh = mito::mesh::tetra(mesh, coord_system); // write mesh to vtk file - mito::io::vtk::writer("tetra_output", tetra_mesh, coord_system); + mito::io::vtk::writer("tetra_output", tetra_mesh, coord_system).write(); } \ No newline at end of file diff --git a/tests/mito.lib/io/summit_to_summit_mesh_2D.cc b/tests/mito.lib/io/summit_to_summit_mesh_2D.cc index 11b103ed..2cc7c5ae 100644 --- a/tests/mito.lib/io/summit_to_summit_mesh_2D.cc +++ b/tests/mito.lib/io/summit_to_summit_mesh_2D.cc @@ -46,7 +46,7 @@ TEST(SummitToSummit, Mesh2D) #ifdef WITH_VTK // write mesh to vtk file - mito::io::vtk::writer("rectangle_copy", mesh, coord_system); + mito::io::vtk::writer("rectangle_copy", mesh, coord_system).write(); #endif } diff --git a/tests/mito.lib/io/summit_to_vtk_mesh_2D.cc b/tests/mito.lib/io/summit_to_vtk_mesh_2D.cc index 03ec45b5..6a1ddee0 100644 --- a/tests/mito.lib/io/summit_to_vtk_mesh_2D.cc +++ b/tests/mito.lib/io/summit_to_vtk_mesh_2D.cc @@ -22,5 +22,5 @@ TEST(SummitToVTK, Mesh2D) auto mesh = mito::io::summit::reader>(fileStream, coord_system); // write mesh to vtk file - mito::io::vtk::writer("rectangle_output", mesh, coord_system); + mito::io::vtk::writer("rectangle_output", mesh, coord_system).write(); } \ No newline at end of file diff --git a/tests/mito.lib/io/summit_to_vtk_mesh_3D.cc b/tests/mito.lib/io/summit_to_vtk_mesh_3D.cc index 0f3c3d83..551a329a 100644 --- a/tests/mito.lib/io/summit_to_vtk_mesh_3D.cc +++ b/tests/mito.lib/io/summit_to_vtk_mesh_3D.cc @@ -23,5 +23,5 @@ TEST(SummitToVTK, Mesh3D) mito::io::summit::reader>(fileStream, coord_system); // write mesh to vtk file - mito::io::vtk::writer("cube_output", mesh, coord_system); + mito::io::vtk::writer("cube_output", mesh, coord_system).write(); } \ No newline at end of file diff --git a/tests/mito.lib/io/write_point_cloud_vtk.cc b/tests/mito.lib/io/write_point_cloud_vtk.cc index d35b7c06..305cc570 100644 --- a/tests/mito.lib/io/write_point_cloud_vtk.cc +++ b/tests/mito.lib/io/write_point_cloud_vtk.cc @@ -28,6 +28,6 @@ TEST(VtkWriter, WritePointCloudToVtk) const auto point_c = cloud.point(); coord_system.place(point_c, { 0.0, 1.0, 1.0 }); - // print the point cloud - mito::io::vtk::writer("point_cloud_output", cloud, coord_system); + // write point cloud to vtk file + mito::io::vtk::writer("point_cloud_output", cloud, coord_system).write(); } diff --git a/tests/mito.lib/mesh/ball.cc b/tests/mito.lib/mesh/ball.cc index 7bbbd381..aeec1b16 100644 --- a/tests/mito.lib/mesh/ball.cc +++ b/tests/mito.lib/mesh/ball.cc @@ -27,7 +27,7 @@ TEST(Mesh, Ball) #ifdef WITH_VTK // write mesh to vtk file - mito::io::vtk::writer("sphere", boundary_mesh, coord_system); + mito::io::vtk::writer("sphere", boundary_mesh, coord_system).write(); #endif // fetch the boundary of the sphere diff --git a/tests/mito.lib/mesh/filter_ball.cc b/tests/mito.lib/mesh/filter_ball.cc index ea40360e..6a9263e8 100644 --- a/tests/mito.lib/mesh/filter_ball.cc +++ b/tests/mito.lib/mesh/filter_ball.cc @@ -22,7 +22,7 @@ TEST(Mesh, FilterBall) #ifdef WITH_VTK // write mesh to vtk file - mito::io::vtk::writer("ball_internal_faces", mesh_faces, coord_system); + mito::io::vtk::writer("ball_internal_faces", mesh_faces, coord_system).write(); #endif // mesh wireframe @@ -30,7 +30,7 @@ TEST(Mesh, FilterBall) #ifdef WITH_VTK // write mesh to vtk file - mito::io::vtk::writer("ball_wireframe", wireframe, coord_system); + mito::io::vtk::writer("ball_wireframe", wireframe, coord_system).write(); #endif // mesh wireframe @@ -47,6 +47,6 @@ TEST(Mesh, FilterBall) #ifdef WITH_VTK // write mesh to vtk file - mito::io::vtk::writer("sphere_wireframe", boundary_wireframe, coord_system); + mito::io::vtk::writer("sphere_wireframe", boundary_wireframe, coord_system).write(); #endif } diff --git a/tests/mito.lib/mesh/sphere.cc b/tests/mito.lib/mesh/sphere.cc index 8fef3db3..507e3f42 100644 --- a/tests/mito.lib/mesh/sphere.cc +++ b/tests/mito.lib/mesh/sphere.cc @@ -23,7 +23,7 @@ TEST(Mesh, Sphere) #ifdef WITH_VTK // write mesh to vtk file - mito::io::vtk::writer("sphere", mesh, coord_system); + mito::io::vtk::writer("sphere", mesh, coord_system).write(); #endif // fetch the boundary of the sphere From d71cc05ec0e0b3c7dc833b52903f8fd11e51c3cd Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Mon, 27 May 2024 19:08:18 +0200 Subject: [PATCH 24/58] lib/mito/io: move {MeshWriter} and {PointCloudWriter} out of {vtk} directory --- lib/mito/io/{vtk => }/MeshWriter.h | 0 lib/mito/io/{vtk => }/PointCloudWriter.h | 0 lib/mito/io/{vtk => }/api.h | 0 lib/mito/io/{vtk => }/factories.h | 0 lib/mito/io/{vtk => }/forward.h | 0 lib/mito/io/public.h | 13 +++++++++++++ lib/mito/io/vtk/public.h | 14 -------------- 7 files changed, 13 insertions(+), 14 deletions(-) rename lib/mito/io/{vtk => }/MeshWriter.h (100%) rename lib/mito/io/{vtk => }/PointCloudWriter.h (100%) rename lib/mito/io/{vtk => }/api.h (100%) rename lib/mito/io/{vtk => }/factories.h (100%) rename lib/mito/io/{vtk => }/forward.h (100%) diff --git a/lib/mito/io/vtk/MeshWriter.h b/lib/mito/io/MeshWriter.h similarity index 100% rename from lib/mito/io/vtk/MeshWriter.h rename to lib/mito/io/MeshWriter.h diff --git a/lib/mito/io/vtk/PointCloudWriter.h b/lib/mito/io/PointCloudWriter.h similarity index 100% rename from lib/mito/io/vtk/PointCloudWriter.h rename to lib/mito/io/PointCloudWriter.h diff --git a/lib/mito/io/vtk/api.h b/lib/mito/io/api.h similarity index 100% rename from lib/mito/io/vtk/api.h rename to lib/mito/io/api.h diff --git a/lib/mito/io/vtk/factories.h b/lib/mito/io/factories.h similarity index 100% rename from lib/mito/io/vtk/factories.h rename to lib/mito/io/factories.h diff --git a/lib/mito/io/vtk/forward.h b/lib/mito/io/forward.h similarity index 100% rename from lib/mito/io/vtk/forward.h rename to lib/mito/io/forward.h diff --git a/lib/mito/io/public.h b/lib/mito/io/public.h index 4fa55f17..d19a74b1 100644 --- a/lib/mito/io/public.h +++ b/lib/mito/io/public.h @@ -10,11 +10,24 @@ // external packages #include "externals.h" +// get the forward declarations +#include "forward.h" + +// published type factories; this is the file you are looking for... +#include "api.h" + // classes implementation #include "summit/public.h" #ifdef WITH_VTK #include "vtk/public.h" #endif // WITH_VTK +// classes implementation +#include "MeshWriter.h" +#include "PointCloudWriter.h" + +// factories implementation +#include "factories.h" + // end of file diff --git a/lib/mito/io/vtk/public.h b/lib/mito/io/vtk/public.h index c1495e68..7899d66f 100644 --- a/lib/mito/io/vtk/public.h +++ b/lib/mito/io/vtk/public.h @@ -10,22 +10,8 @@ // external packages #include "externals.h" -// get the forward declarations -#include "forward.h" - -// published type factories; this is the file you are looking for... -#include "api.h" - // wrappers #include "vtk_cell.h" #include "vtk_point.h" -// classes implementation -#include "MeshWriter.h" -#include "PointCloudWriter.h" - -// factories implementation -#include "factories.h" - - // end of file From d32905830abbfd65a902e64574df07c91e48cc8b Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Mon, 27 May 2024 19:23:25 +0200 Subject: [PATCH 25/58] lib/mito/io: move external packages needed from {summit} reader/writer to {summit} directory --- lib/mito/io/externals.h | 3 +-- lib/mito/io/summit/externals.h | 15 +++++++++++++++ lib/mito/io/summit/public.h | 3 +++ 3 files changed, 19 insertions(+), 2 deletions(-) create mode 100644 lib/mito/io/summit/externals.h diff --git a/lib/mito/io/externals.h b/lib/mito/io/externals.h index dcf10bf2..e8c3c1ee 100644 --- a/lib/mito/io/externals.h +++ b/lib/mito/io/externals.h @@ -8,8 +8,7 @@ // externals -#include -#include + // support #include "../journal.h" diff --git a/lib/mito/io/summit/externals.h b/lib/mito/io/summit/externals.h new file mode 100644 index 00000000..5d413d92 --- /dev/null +++ b/lib/mito/io/summit/externals.h @@ -0,0 +1,15 @@ +// -*- c++ -*- +// +// Copyright (c) 2020-2024, the MiTo Authors, all rights reserved +// + +// code guard +#pragma once + + +// externals +#include +#include + + +// end of file diff --git a/lib/mito/io/summit/public.h b/lib/mito/io/summit/public.h index 4c9495a6..deeb5525 100644 --- a/lib/mito/io/summit/public.h +++ b/lib/mito/io/summit/public.h @@ -7,6 +7,9 @@ #pragma once +// external packages +#include "externals.h" + // classes implementation #include "reader.h" #include "writer.h" From 4895a80b62919af74067770209ed515a510531cd Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Mon, 27 May 2024 19:24:10 +0200 Subject: [PATCH 26/58] lib/mito/io/summit: add {std} namespace specialization when appropriate --- lib/mito/io/summit/writer.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/mito/io/summit/writer.h b/lib/mito/io/summit/writer.h index 131d30d0..659764e6 100644 --- a/lib/mito/io/summit/writer.h +++ b/lib/mito/io/summit/writer.h @@ -56,7 +56,7 @@ namespace mito::io::summit { for (const auto & cell : mesh.cells()) { outfile << cellT::n_vertices << " "; for (const auto & node : cell.nodes()) { - outfile << distance(points.begin(), points.find(node.point().id())) + 1 << " "; + outfile << std::distance(points.begin(), points.find(node.point())) + 1 << " "; } // TOFIX: material label is always 1 for now outfile << 1 << std::endl; From bf2369a2cba819ea5c93250cc32fd12c65d47ba9 Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Mon, 27 May 2024 19:39:30 +0200 Subject: [PATCH 27/58] lib/mito/io/summit: use a {set} of points for internal book-keeping No need to use a {map}. --- lib/mito/io/summit/externals.h | 2 +- lib/mito/io/summit/writer.h | 12 ++++-------- 2 files changed, 5 insertions(+), 9 deletions(-) diff --git a/lib/mito/io/summit/externals.h b/lib/mito/io/summit/externals.h index 5d413d92..43484377 100644 --- a/lib/mito/io/summit/externals.h +++ b/lib/mito/io/summit/externals.h @@ -9,7 +9,7 @@ // externals #include -#include +#include // end of file diff --git a/lib/mito/io/summit/writer.h b/lib/mito/io/summit/writer.h index 659764e6..e48aa96f 100644 --- a/lib/mito/io/summit/writer.h +++ b/lib/mito/io/summit/writer.h @@ -27,17 +27,13 @@ namespace mito::io::summit { // type of point using point_type = geometry::point_t; - // type of point id - using point_id_type = utilities::index_t; - - // a map between point ids to points - std::unordered_map points; + // a set between of points (to remove duplicates) + std::unordered_set> points; // insert the points corresponding to the mesh nodes for (const auto & cell : mesh.cells()) { for (const auto & node : cell.nodes()) { - const auto & point = node.point(); - points.insert({ point.id(), point }); + points.insert(node.point()); } } @@ -47,7 +43,7 @@ namespace mito::io::summit { outfile << std::size(points) << " " << mesh.nCells() << " " << 1 << std::endl; // write the points to file - for (const auto & [_, point] : points) { + for (const auto & point : points) { const auto & coord = coordinate_system.coordinates(point); outfile << std::setprecision(15) << coord << std::endl; } From 03d4f3e2199355f2fb17db9ff2d1819fddf2dd63 Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Sat, 1 Jun 2024 18:52:09 +0200 Subject: [PATCH 28/58] lib/mito/io: readability improvements --- lib/mito/io/MeshWriter.h | 21 +++++++++++---------- 1 file changed, 11 insertions(+), 10 deletions(-) diff --git a/lib/mito/io/MeshWriter.h b/lib/mito/io/MeshWriter.h index 8603c2ff..826f6d92 100644 --- a/lib/mito/io/MeshWriter.h +++ b/lib/mito/io/MeshWriter.h @@ -22,10 +22,11 @@ namespace mito::io::vtk { static constexpr int D = mesh_type::cell_type::dim; // the type of grid using grid_type = vtkSmartPointer; + using point_type = geometry::point_t; // map mesh points to the index of the vtk points. Points that are shared among // multiple elements have the same index. - using point_to_index_type = std::unordered_map< - geometry::point_t, int, utilities::hash_function>>; + using points_type = + std::unordered_map>; private: auto _create_vtk_grid(const mesh_type & mesh, const coord_system_type & coordinate_system) @@ -53,16 +54,16 @@ namespace mito::io::vtk { // retrieve the corresponding point const auto & point = node.point(); // if the point is not present in the map - if (!_point_to_index.contains(point)) { + if (!_points.contains(point)) { // insert the new vtk point insert_vtk_point(coordinate_system.coordinates(point), pointsVtk); // add the point to the map with its global index - _point_to_index[point] = indexPointVtk; + _points[point] = indexPointVtk; // update global index for the vtk point ++indexPointVtk; } // set the id of the point - cellVtk->GetPointIds()->SetId(indexLocalPointVtk, _point_to_index[point]); + cellVtk->GetPointIds()->SetId(indexLocalPointVtk, _points[point]); // update local index for the points in the cell ++indexLocalPointVtk; } @@ -82,7 +83,7 @@ namespace mito::io::vtk { MeshWriter( std::string filename, const mesh_type & mesh, const coord_system_type & coord_system) : _filename(filename), - _point_to_index(), + _points(), _grid(_create_vtk_grid(mesh, coord_system)) {} @@ -103,8 +104,8 @@ namespace mito::io::vtk { std::cout << n_nodes << "\t" << _grid->GetNumberOfPoints() << std::endl; // assert(n_nodes == _grid->GetNumberOfPoints()); // and equal to the number of indices that we have recorded so far - // assert(n_nodes == int(std::size(_point_to_index))); - std::cout << n_nodes << "\t" << int(std::size(_point_to_index)) << std::endl; + // assert(n_nodes == int(std::size(_points))); + std::cout << n_nodes << "\t" << int(std::size(_points)) << std::endl; // initialize a vtk array auto vtkArray = vtkSmartPointer::New(); @@ -115,7 +116,7 @@ namespace mito::io::vtk { // populate the array with the nodal values for (auto & [node, value] : field) { // get the index corresponding to the current point - auto i = _point_to_index.at(node.point()); + auto i = _points.at(node.point()); vtkArray->SetTuple(i, value.begin()); } @@ -151,7 +152,7 @@ namespace mito::io::vtk { std::string _filename; // a map storing point -> index relation - point_to_index_type _point_to_index; + points_type _points; // the grid grid_type _grid; From 41ecbb8ed7b0f772e3cf83b3012660c5edc15c92 Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Sat, 1 Jun 2024 18:58:27 +0200 Subject: [PATCH 29/58] lib/mito/io/vtk: use STL {set} instead of {map} No need to use a map, as the 'value' (that is the sequential index of the point stored) can be retrieved by the distance of a given point from the beginning of the set memory allocation. --- lib/mito/io/MeshWriter.h | 42 +++++++++++++++++++++++----------------- 1 file changed, 24 insertions(+), 18 deletions(-) diff --git a/lib/mito/io/MeshWriter.h b/lib/mito/io/MeshWriter.h index 826f6d92..7ca377b9 100644 --- a/lib/mito/io/MeshWriter.h +++ b/lib/mito/io/MeshWriter.h @@ -25,20 +25,31 @@ namespace mito::io::vtk { using point_type = geometry::point_t; // map mesh points to the index of the vtk points. Points that are shared among // multiple elements have the same index. - using points_type = - std::unordered_map>; + using points_type = std::unordered_set>; private: auto _create_vtk_grid(const mesh_type & mesh, const coord_system_type & coordinate_system) -> grid_type { - // vtk unstructured grid - auto gridVtk = vtkSmartPointer::New(); + // loop over the cells + for (const auto & cell : mesh.cells()) { + // loop over the nodes of the cell + for (const auto & node : cell.nodes()) { + // add the point to the map with its global index + _points.insert(node.point()); + } + } + // vtk points and cells auto pointsVtk = vtkSmartPointer::New(); - // global index assigned to each vtk point - auto indexPointVtk = 0; + // insert the new vtk point + for (const auto & point : _points) { + insert_vtk_point(coordinate_system.coordinates(point), pointsVtk); + } + + // vtk unstructured grid + auto gridVtk = vtkSmartPointer::New(); // loop over the cells for (const auto & cell : mesh.cells()) { @@ -53,17 +64,12 @@ namespace mito::io::vtk { for (const auto & node : cell.nodes()) { // retrieve the corresponding point const auto & point = node.point(); - // if the point is not present in the map - if (!_points.contains(point)) { - // insert the new vtk point - insert_vtk_point(coordinate_system.coordinates(point), pointsVtk); - // add the point to the map with its global index - _points[point] = indexPointVtk; - // update global index for the vtk point - ++indexPointVtk; - } + // assert that the point is present in the set of points + assert(_points.contains(point)); + // calculate the index of the point + auto index = std::distance(_points.begin(), _points.find(point)); // set the id of the point - cellVtk->GetPointIds()->SetId(indexLocalPointVtk, _points[point]); + cellVtk->GetPointIds()->SetId(indexLocalPointVtk, index); // update local index for the points in the cell ++indexLocalPointVtk; } @@ -116,8 +122,8 @@ namespace mito::io::vtk { // populate the array with the nodal values for (auto & [node, value] : field) { // get the index corresponding to the current point - auto i = _points.at(node.point()); - vtkArray->SetTuple(i, value.begin()); + auto index = std::distance(_points.begin(), _points.find(node.point())); + vtkArray->SetTuple(index, value.begin()); } // insert array into output mesh From 5496474f29788f0684b20d378cbe28990ae80fc3 Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Mon, 3 Jun 2024 12:06:51 +0200 Subject: [PATCH 30/58] lib/mito/io/vtk: move {vtk} writers to {vtk} directories --- lib/mito/io/GridWriter.h | 27 ++++ lib/mito/io/MeshWriter.h | 170 -------------------------- lib/mito/io/PointCloudWriter.h | 145 ---------------------- lib/mito/io/public.h | 14 +-- lib/mito/io/vtk/GridWriterVTK.h | 88 +++++++++++++ lib/mito/io/vtk/MeshWriterVTK.h | 90 ++++++++++++++ lib/mito/io/vtk/PointCloudWriterVTK.h | 59 +++++++++ lib/mito/io/{ => vtk}/api.h | 4 +- lib/mito/io/{ => vtk}/factories.h | 4 +- lib/mito/io/{ => vtk}/forward.h | 4 +- lib/mito/io/vtk/public.h | 15 +++ 11 files changed, 287 insertions(+), 333 deletions(-) create mode 100644 lib/mito/io/GridWriter.h delete mode 100644 lib/mito/io/MeshWriter.h delete mode 100644 lib/mito/io/PointCloudWriter.h create mode 100644 lib/mito/io/vtk/GridWriterVTK.h create mode 100644 lib/mito/io/vtk/MeshWriterVTK.h create mode 100644 lib/mito/io/vtk/PointCloudWriterVTK.h rename lib/mito/io/{ => vtk}/api.h (90%) rename lib/mito/io/{ => vtk}/factories.h (83%) rename lib/mito/io/{ => vtk}/forward.h (90%) diff --git a/lib/mito/io/GridWriter.h b/lib/mito/io/GridWriter.h new file mode 100644 index 00000000..8cb8bd9f --- /dev/null +++ b/lib/mito/io/GridWriter.h @@ -0,0 +1,27 @@ +// -*- c++ -*- +// +// Copyright (c) 2020-2024, the MiTo Authors, all rights reserved +// + +// code guard +#pragma once + + +namespace mito::io { + + class GridWriter { + + public: + GridWriter(std::string filename) : _filename(filename) {} + + virtual auto write() const -> void = 0; + + protected: + // the name of the output file to be written + std::string _filename; + }; + +} // namespace mito::io + + +// end of file diff --git a/lib/mito/io/MeshWriter.h b/lib/mito/io/MeshWriter.h deleted file mode 100644 index 7ca377b9..00000000 --- a/lib/mito/io/MeshWriter.h +++ /dev/null @@ -1,170 +0,0 @@ -// -*- c++ -*- -// -// Copyright (c) 2020-2024, the MiTo Authors, all rights reserved -// - -// code guard -#pragma once - - -namespace mito::io::vtk { - - template - requires(meshT::cell_type::dim == coordSystemT::coordinates_type::dim) - class MeshWriter { - - private: - // the mesh type - using mesh_type = meshT; - // the coordinate system type - using coord_system_type = coordSystemT; - // the dimension of the physical space - static constexpr int D = mesh_type::cell_type::dim; - // the type of grid - using grid_type = vtkSmartPointer; - using point_type = geometry::point_t; - // map mesh points to the index of the vtk points. Points that are shared among - // multiple elements have the same index. - using points_type = std::unordered_set>; - - private: - auto _create_vtk_grid(const mesh_type & mesh, const coord_system_type & coordinate_system) - -> grid_type - { - // loop over the cells - for (const auto & cell : mesh.cells()) { - // loop over the nodes of the cell - for (const auto & node : cell.nodes()) { - // add the point to the map with its global index - _points.insert(node.point()); - } - } - - // vtk points and cells - auto pointsVtk = vtkSmartPointer::New(); - - // insert the new vtk point - for (const auto & point : _points) { - insert_vtk_point(coordinate_system.coordinates(point), pointsVtk); - } - - // vtk unstructured grid - auto gridVtk = vtkSmartPointer::New(); - - // loop over the cells - for (const auto & cell : mesh.cells()) { - - // create vtk cell - auto cellVtk = vtkCellPointer(); - - // local index for the points of the cell - auto indexLocalPointVtk = 0; - - // loop over the nodes of the cell - for (const auto & node : cell.nodes()) { - // retrieve the corresponding point - const auto & point = node.point(); - // assert that the point is present in the set of points - assert(_points.contains(point)); - // calculate the index of the point - auto index = std::distance(_points.begin(), _points.find(point)); - // set the id of the point - cellVtk->GetPointIds()->SetId(indexLocalPointVtk, index); - // update local index for the points in the cell - ++indexLocalPointVtk; - } - - // insert the new cell - gridVtk->InsertNextCell(cellVtk->GetCellType(), cellVtk->GetPointIds()); - } - - // set the grid points - gridVtk->SetPoints(pointsVtk); - - // all done - return gridVtk; - } - - public: - MeshWriter( - std::string filename, const mesh_type & mesh, const coord_system_type & coord_system) : - _filename(filename), - _points(), - _grid(_create_vtk_grid(mesh, coord_system)) - {} - - // TOFIX: use concepts to say that Y is a tensor - template - auto record(const fem::nodal_field_t & field, std::string fieldname = "") -> void - { - // if no name was provided - if (fieldname == "") { - // use the name of the field - fieldname = field.name(); - } - - // get the number of nodes - auto n_nodes = field.n_nodes(); - - // check the number of nodes in the field equals the number of points in the grid - std::cout << n_nodes << "\t" << _grid->GetNumberOfPoints() << std::endl; - // assert(n_nodes == _grid->GetNumberOfPoints()); - // and equal to the number of indices that we have recorded so far - // assert(n_nodes == int(std::size(_points))); - std::cout << n_nodes << "\t" << int(std::size(_points)) << std::endl; - - // initialize a vtk array - auto vtkArray = vtkSmartPointer::New(); - vtkArray->SetName(fieldname.data()); - vtkArray->SetNumberOfComponents(Y::size); - vtkArray->SetNumberOfTuples(field.n_nodes()); - - // populate the array with the nodal values - for (auto & [node, value] : field) { - // get the index corresponding to the current point - auto index = std::distance(_points.begin(), _points.find(node.point())); - vtkArray->SetTuple(index, value.begin()); - } - - // insert array into output mesh - _grid->GetPointData()->AddArray(vtkArray); - - // all done - return; - } - - auto write() const -> void - { - // create a new writer - auto writer = vtkSmartPointer::New(); - // set the name of the output file - writer->SetFileName((_filename + ".vtu").data()); - - // sign the grid up for writing -#if VTK_MAJOR_VERSION <= 8 - writer->SetInput(_grid); -#else - writer->SetInputData(_grid); -#endif - // write the grid to file - writer->Write(); - - // all done - return; - } - - private: - // the name of the output file to be written - std::string _filename; - - // a map storing point -> index relation - points_type _points; - - // the grid - grid_type _grid; - }; - -} // namespace mito::io::vtk - - -// end of file diff --git a/lib/mito/io/PointCloudWriter.h b/lib/mito/io/PointCloudWriter.h deleted file mode 100644 index b654a39a..00000000 --- a/lib/mito/io/PointCloudWriter.h +++ /dev/null @@ -1,145 +0,0 @@ -// -*- c++ -*- -// -// Copyright (c) 2020-2024, the MiTo Authors, all rights reserved -// - -// code guard -#pragma once - - -namespace mito::io::vtk { - - template - requires(cloudT::dim == coordSystemT::coordinates_type::dim) - class PointCloudWriter { - - private: - // the point cloud type - using cloud_type = cloudT; - // the coordinate system type - using coord_system_type = coordSystemT; - // the dimension of the physical space - static constexpr int D = cloudT::dim; - // the type of grid - using grid_type = vtkSmartPointer; - // map mesh points to the index of the vtk points. Points that are shared among - // multiple elements have the same index. - using point_to_index_type = std::unordered_map< - geometry::point_t, int, utilities::hash_function>>; - - private: - auto _create_vtk_grid(const cloud_type & cloud, const coord_system_type & coordinate_system) - -> grid_type - { - // vtk unstructured grid - auto gridVtk = vtkSmartPointer::New(); - // vtk points and cells - auto pointsVtk = vtkSmartPointer::New(); - - // global index assigned to each vtk point - auto indexPointVtk = 0; - - // get point cloud compositions - const auto & points = cloud.points(); - - // iterate over the points - for (const auto & point : points) { - // if the point is not present in the map - if (!_point_to_index.contains(point)) { - insert_vtk_point(coordinate_system.coordinates(point), pointsVtk); - // add the point to the map with its global index - _point_to_index[point] = indexPointVtk; - // update global index for the vtk point - ++indexPointVtk; - } - } - - // set the grid points - gridVtk->SetPoints(pointsVtk); - - // all done - return gridVtk; - } - - public: - PointCloudWriter( - std::string filename, const cloud_type & cloud, - const coord_system_type & coord_system) : - _filename(filename), - _point_to_index(), - _grid(_create_vtk_grid(cloud, coord_system)) - {} - - // TOFIX: use concepts to say that Y is a tensor - template - auto record(const fem::nodal_field_t & field, std::string fieldname = "") -> void - { - // if no name was provided - if (fieldname == "") { - // use the name of the field - fieldname = field.name(); - } - - // get the number of nodes - auto n_nodes = field.n_nodes(); - - // check the number of nodes in the field equals the number of points in the grid - assert(n_nodes == _grid->GetNumberOfPoints()); - // and equal to the number of indices that we have recorded so far - assert(n_nodes == int(std::size(_point_to_index))); - - // initialize a vtk array - auto vtkArray = vtkSmartPointer::New(); - vtkArray->SetName(fieldname.data()); - vtkArray->SetNumberOfComponents(Y::size); - vtkArray->SetNumberOfTuples(field.n_nodes()); - - // populate the array with the nodal values - for (auto & [node, value] : field) { - // get the index corresponding to the current point - auto i = _point_to_index.at(node.point()); - vtkArray->SetTuple(i, value.begin()); - } - - // insert array into output mesh - _grid->GetPointData()->AddArray(vtkArray); - - // all done - return; - } - - auto write() const -> void - { - // create a new writer - auto writer = vtkSmartPointer::New(); - // set the name of the output file - writer->SetFileName((_filename + ".vtu").data()); - - // sign the grid up for writing -#if VTK_MAJOR_VERSION <= 8 - writer->SetInput(_grid); -#else - writer->SetInputData(_grid); -#endif - // write the grid to file - writer->Write(); - - // all done - return; - } - - private: - // the name of the output file to be written - std::string _filename; - - // a map storing point -> index relation - point_to_index_type _point_to_index; - - // the grid - grid_type _grid; - }; - -} // namespace mito::io::vtk - - -// end of file diff --git a/lib/mito/io/public.h b/lib/mito/io/public.h index d19a74b1..09f06291 100644 --- a/lib/mito/io/public.h +++ b/lib/mito/io/public.h @@ -10,11 +10,8 @@ // external packages #include "externals.h" -// get the forward declarations -#include "forward.h" - -// published type factories; this is the file you are looking for... -#include "api.h" +// classes implementation +#include "GridWriter.h" // classes implementation #include "summit/public.h" @@ -22,12 +19,5 @@ #include "vtk/public.h" #endif // WITH_VTK -// classes implementation -#include "MeshWriter.h" -#include "PointCloudWriter.h" - -// factories implementation -#include "factories.h" - // end of file diff --git a/lib/mito/io/vtk/GridWriterVTK.h b/lib/mito/io/vtk/GridWriterVTK.h new file mode 100644 index 00000000..7aef66c8 --- /dev/null +++ b/lib/mito/io/vtk/GridWriterVTK.h @@ -0,0 +1,88 @@ +// -*- c++ -*- +// +// Copyright (c) 2020-2024, the MiTo Authors, all rights reserved +// + +// code guard +#pragma once + + +namespace mito::io::vtk { + + template + class GridWriterVTK : public GridWriter { + + private: + // the type of grid + using grid_type = vtkSmartPointer; + // the type of point + using point_type = geometry::point_t; + // the type of a collection of points + using points_type = std::unordered_set>; + + protected: + GridWriterVTK(std::string filename) : + GridWriter(filename), + _grid(grid_type::New()), + _points() + {} + + public: + auto write() const -> void override + { + // create a new writer + auto writer = vtkSmartPointer::New(); + // set the name of the output file + writer->SetFileName((_filename + ".vtu").data()); + + // sign the grid up for writing +#if VTK_MAJOR_VERSION <= 8 + writer->SetInput(_grid); +#else + writer->SetInputData(_grid); +#endif + // write the grid to file + writer->Write(); + + // all done + return; + } + + template + auto attach_field(const fem::nodal_field_t & field, std::string fieldname) -> void + { + // get the number of nodes + auto n_nodes = field.n_nodes(); + + // check the number of nodes in the field equals the number of points in the grid + assert(n_nodes == _grid->GetNumberOfPoints()); + + // initialize a vtk array + auto vtkArray = vtkSmartPointer::New(); + vtkArray->SetName(fieldname.data()); + vtkArray->SetNumberOfComponents(Y::size); + vtkArray->SetNumberOfTuples(n_nodes); + + // populate the array with the nodal values + for (auto & [node, value] : field) { + // get the index corresponding to the current point + auto index = std::distance(_points.begin(), _points.find(node.point())); + vtkArray->SetTuple(index, value.begin()); + } + + // insert array into output mesh + _grid->GetPointData()->AddArray(vtkArray); + } + + protected: + // the grid + grid_type _grid; + + // a collection of points in the grid + points_type _points; + }; + +} // namespace mito::io::vtk + + +// end of file diff --git a/lib/mito/io/vtk/MeshWriterVTK.h b/lib/mito/io/vtk/MeshWriterVTK.h new file mode 100644 index 00000000..756f2aeb --- /dev/null +++ b/lib/mito/io/vtk/MeshWriterVTK.h @@ -0,0 +1,90 @@ +// -*- c++ -*- +// +// Copyright (c) 2020-2024, the MiTo Authors, all rights reserved +// + +// code guard +#pragma once + + +namespace mito::io::vtk { + + template + requires(meshT::cell_type::dim == coordSystemT::coordinates_type::dim) + class MeshWriterVTK : public GridWriterVTK { + + private: + // the mesh type + using mesh_type = meshT; + // the coordinate system type + using coord_system_type = coordSystemT; + // the dimension of the physical space + static constexpr int D = mesh_type::cell_type::dim; + + private: + auto _create_vtk_grid(const mesh_type & mesh, const coord_system_type & coordinate_system) + { + // loop over the cells + for (const auto & cell : mesh.cells()) { + // loop over the nodes of the cell + for (const auto & node : cell.nodes()) { + // add the point to the collection of points (eliminating duplicates) + _points.insert(node.point()); + } + } + + // vtk points and cells + auto pointsVtk = vtkSmartPointer::New(); + + // insert the new vtk point + for (const auto & point : _points) { + insert_vtk_point(coordinate_system.coordinates(point), pointsVtk); + } + + // loop over the cells + for (const auto & cell : mesh.cells()) { + + // create vtk cell + auto cellVtk = vtkCellPointer(); + + // local index for the points of the cell + auto indexLocalPointVtk = 0; + + // loop over the nodes of the cell + for (const auto & node : cell.nodes()) { + // retrieve the corresponding point + const auto & point = node.point(); + // assert that the point is present in the set of points + assert(_points.contains(point)); + // calculate the index of the point + auto index = std::distance(_points.begin(), _points.find(point)); + // set the id of the point + cellVtk->GetPointIds()->SetId(indexLocalPointVtk, index); + // update local index for the points in the cell + ++indexLocalPointVtk; + } + + // insert the new cell + _grid->InsertNextCell(cellVtk->GetCellType(), cellVtk->GetPointIds()); + } + + // set the grid points + _grid->SetPoints(pointsVtk); + + // all done + return; + } + + public: + MeshWriterVTK( + std::string filename, const mesh_type & mesh, const coord_system_type & coord_system) : + GridWriterVTK(filename) + { + _create_vtk_grid(mesh, coord_system); + } + }; + +} // namespace mito::io::vtk + + +// end of file diff --git a/lib/mito/io/vtk/PointCloudWriterVTK.h b/lib/mito/io/vtk/PointCloudWriterVTK.h new file mode 100644 index 00000000..f81c2f5c --- /dev/null +++ b/lib/mito/io/vtk/PointCloudWriterVTK.h @@ -0,0 +1,59 @@ +// -*- c++ -*- +// +// Copyright (c) 2020-2024, the MiTo Authors, all rights reserved +// + +// code guard +#pragma once + + +namespace mito::io::vtk { + + template + requires(cloudT::dim == coordSystemT::coordinates_type::dim) + class PointCloudWriterVTK : public GridWriterVTK { + + private: + // the point cloud type + using cloud_type = cloudT; + // the coordinate system type + using coord_system_type = coordSystemT; + // the dimension of the physical space + static constexpr int D = cloudT::dim; + + private: + auto _create_vtk_grid(const cloud_type & cloud, const coord_system_type & coordinate_system) + { + // vtk points and cells + auto pointsVtk = vtkSmartPointer::New(); + + // get point cloud compositions + const auto & points = cloud.points(); + + // insert the new vtk point + for (const auto & point : points) { + _points.insert(point); + insert_vtk_point(coordinate_system.coordinates(point), pointsVtk); + } + + // set the grid points + _grid->SetPoints(pointsVtk); + + // all done + return; + } + + public: + PointCloudWriterVTK( + std::string filename, const cloud_type & cloud, + const coord_system_type & coord_system) : + GridWriterVTK(filename) + { + _create_vtk_grid(cloud, coord_system); + } + }; + +} // namespace mito::io::vtk + + +// end of file diff --git a/lib/mito/io/api.h b/lib/mito/io/vtk/api.h similarity index 90% rename from lib/mito/io/api.h rename to lib/mito/io/vtk/api.h index c68cd80d..f90fc160 100644 --- a/lib/mito/io/api.h +++ b/lib/mito/io/vtk/api.h @@ -12,7 +12,7 @@ namespace mito::io::vtk { // mesh writer alias template requires(meshT::cell_type::dim == coordSystemT::coordinates_type::dim) - using mesh_writer_t = MeshWriter; + using mesh_writer_t = MeshWriterVTK; // vtk mesh writer factory template @@ -23,7 +23,7 @@ namespace mito::io::vtk { // point cloud writer alias template requires(cloudT::dim == coordSystemT::coordinates_type::dim) - using cloud_writer_t = PointCloudWriter; + using cloud_writer_t = PointCloudWriterVTK; // point cloud writer factory template diff --git a/lib/mito/io/factories.h b/lib/mito/io/vtk/factories.h similarity index 83% rename from lib/mito/io/factories.h rename to lib/mito/io/vtk/factories.h index 21d48a20..01198882 100644 --- a/lib/mito/io/factories.h +++ b/lib/mito/io/vtk/factories.h @@ -15,7 +15,7 @@ namespace mito::io::vtk { auto writer(std::string filename, const meshT & mesh, const coordSystemT & coord_system) -> mesh_writer_t { - return MeshWriter(filename, mesh, coord_system); + return MeshWriterVTK(filename, mesh, coord_system); } // vtk point cloud writer factory @@ -24,7 +24,7 @@ namespace mito::io::vtk { auto writer(std::string filename, const cloudT & cloud, const coordSystemT & coord_system) -> cloud_writer_t { - return PointCloudWriter(filename, cloud, coord_system); + return PointCloudWriterVTK(filename, cloud, coord_system); } } diff --git a/lib/mito/io/forward.h b/lib/mito/io/vtk/forward.h similarity index 90% rename from lib/mito/io/forward.h rename to lib/mito/io/vtk/forward.h index 8b478f15..58d1152f 100644 --- a/lib/mito/io/forward.h +++ b/lib/mito/io/vtk/forward.h @@ -12,12 +12,12 @@ namespace mito::io::vtk { // class mesh writer template requires(meshT::cell_type::dim == coordSystemT::coordinates_type::dim) - class MeshWriter; + class MeshWriterVTK; // class point cloud writer template requires(cloudT::dim == coordSystemT::coordinates_type::dim) - class PointCloudWriter; + class PointCloudWriterVTK; } diff --git a/lib/mito/io/vtk/public.h b/lib/mito/io/vtk/public.h index 7899d66f..ffdd34a8 100644 --- a/lib/mito/io/vtk/public.h +++ b/lib/mito/io/vtk/public.h @@ -10,8 +10,23 @@ // external packages #include "externals.h" +// get the forward declarations +#include "forward.h" + +// published type factories; this is the file you are looking for... +#include "api.h" + // wrappers #include "vtk_cell.h" #include "vtk_point.h" +// classes +#include "GridWriterVTK.h" +#include "MeshWriterVTK.h" +#include "PointCloudWriterVTK.h" + +// factories implementation +#include "factories.h" + + // end of file From f1f744038feff6869bc245e40f5e96b5cb56f2cc Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Mon, 3 Jun 2024 12:09:58 +0200 Subject: [PATCH 31/58] lib/mito/io/vtk: refer to attributes in parent class with {this->} A bug in gcc does not make the (protected!) attributes of the parent class visible in the derived class scope. --- lib/mito/io/vtk/MeshWriterVTK.h | 12 ++++++------ lib/mito/io/vtk/PointCloudWriterVTK.h | 4 ++-- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/lib/mito/io/vtk/MeshWriterVTK.h b/lib/mito/io/vtk/MeshWriterVTK.h index 756f2aeb..d72459df 100644 --- a/lib/mito/io/vtk/MeshWriterVTK.h +++ b/lib/mito/io/vtk/MeshWriterVTK.h @@ -29,7 +29,7 @@ namespace mito::io::vtk { // loop over the nodes of the cell for (const auto & node : cell.nodes()) { // add the point to the collection of points (eliminating duplicates) - _points.insert(node.point()); + this->_points.insert(node.point()); } } @@ -37,7 +37,7 @@ namespace mito::io::vtk { auto pointsVtk = vtkSmartPointer::New(); // insert the new vtk point - for (const auto & point : _points) { + for (const auto & point : this->_points) { insert_vtk_point(coordinate_system.coordinates(point), pointsVtk); } @@ -55,9 +55,9 @@ namespace mito::io::vtk { // retrieve the corresponding point const auto & point = node.point(); // assert that the point is present in the set of points - assert(_points.contains(point)); + assert(this->_points.contains(point)); // calculate the index of the point - auto index = std::distance(_points.begin(), _points.find(point)); + auto index = std::distance(this->_points.begin(), this->_points.find(point)); // set the id of the point cellVtk->GetPointIds()->SetId(indexLocalPointVtk, index); // update local index for the points in the cell @@ -65,11 +65,11 @@ namespace mito::io::vtk { } // insert the new cell - _grid->InsertNextCell(cellVtk->GetCellType(), cellVtk->GetPointIds()); + this->_grid->InsertNextCell(cellVtk->GetCellType(), cellVtk->GetPointIds()); } // set the grid points - _grid->SetPoints(pointsVtk); + this->_grid->SetPoints(pointsVtk); // all done return; diff --git a/lib/mito/io/vtk/PointCloudWriterVTK.h b/lib/mito/io/vtk/PointCloudWriterVTK.h index f81c2f5c..5f507148 100644 --- a/lib/mito/io/vtk/PointCloudWriterVTK.h +++ b/lib/mito/io/vtk/PointCloudWriterVTK.h @@ -32,12 +32,12 @@ namespace mito::io::vtk { // insert the new vtk point for (const auto & point : points) { - _points.insert(point); + this->_points.insert(point); insert_vtk_point(coordinate_system.coordinates(point), pointsVtk); } // set the grid points - _grid->SetPoints(pointsVtk); + this->_grid->SetPoints(pointsVtk); // all done return; From bec6b7dd12f0dbd2526b9d13e12bf103b6e78665 Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Mon, 3 Jun 2024 12:12:51 +0200 Subject: [PATCH 32/58] lib/mito/io: make {GridWriter} template with respect to the spatial dimension --- lib/mito/io/GridWriter.h | 12 +++++++++++- lib/mito/io/vtk/GridWriterVTK.h | 19 ++++--------------- 2 files changed, 15 insertions(+), 16 deletions(-) diff --git a/lib/mito/io/GridWriter.h b/lib/mito/io/GridWriter.h index 8cb8bd9f..4916c3ce 100644 --- a/lib/mito/io/GridWriter.h +++ b/lib/mito/io/GridWriter.h @@ -9,16 +9,26 @@ namespace mito::io { + template class GridWriter { + private: + // the type of point + using point_type = geometry::point_t; + // the type of a collection of points + using points_type = std::unordered_set>; + public: - GridWriter(std::string filename) : _filename(filename) {} + GridWriter(std::string filename) : _filename(filename), _points() {} virtual auto write() const -> void = 0; protected: // the name of the output file to be written std::string _filename; + + // a collection of points in the grid + points_type _points; }; } // namespace mito::io diff --git a/lib/mito/io/vtk/GridWriterVTK.h b/lib/mito/io/vtk/GridWriterVTK.h index 7aef66c8..9d9b4552 100644 --- a/lib/mito/io/vtk/GridWriterVTK.h +++ b/lib/mito/io/vtk/GridWriterVTK.h @@ -10,22 +10,14 @@ namespace mito::io::vtk { template - class GridWriterVTK : public GridWriter { + class GridWriterVTK : public GridWriter { private: // the type of grid using grid_type = vtkSmartPointer; - // the type of point - using point_type = geometry::point_t; - // the type of a collection of points - using points_type = std::unordered_set>; protected: - GridWriterVTK(std::string filename) : - GridWriter(filename), - _grid(grid_type::New()), - _points() - {} + GridWriterVTK(std::string filename) : GridWriter(filename), _grid(grid_type::New()) {} public: auto write() const -> void override @@ -33,7 +25,7 @@ namespace mito::io::vtk { // create a new writer auto writer = vtkSmartPointer::New(); // set the name of the output file - writer->SetFileName((_filename + ".vtu").data()); + writer->SetFileName((this->_filename + ".vtu").data()); // sign the grid up for writing #if VTK_MAJOR_VERSION <= 8 @@ -66,7 +58,7 @@ namespace mito::io::vtk { // populate the array with the nodal values for (auto & [node, value] : field) { // get the index corresponding to the current point - auto index = std::distance(_points.begin(), _points.find(node.point())); + auto index = std::distance(this->_points.begin(), this->_points.find(node.point())); vtkArray->SetTuple(index, value.begin()); } @@ -77,9 +69,6 @@ namespace mito::io::vtk { protected: // the grid grid_type _grid; - - // a collection of points in the grid - points_type _points; }; } // namespace mito::io::vtk From 62a10034df54867b74f21af4b2c1d24f6b29a77a Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Mon, 3 Jun 2024 14:18:56 +0200 Subject: [PATCH 33/58] lib/mito/io/vtk: add class {FieldWriterVTK} --- lib/mito/io/vtk/FieldWriterVTK.h | 59 ++++++++++++++++++++++++++++++++ lib/mito/io/vtk/api.h | 18 ++++++---- lib/mito/io/vtk/factories.h | 8 ++--- lib/mito/io/vtk/forward.h | 4 +++ lib/mito/io/vtk/public.h | 4 +++ lib/mito/io/vtk/utilities.h | 30 ++++++++++++++++ 6 files changed, 112 insertions(+), 11 deletions(-) create mode 100644 lib/mito/io/vtk/FieldWriterVTK.h create mode 100644 lib/mito/io/vtk/utilities.h diff --git a/lib/mito/io/vtk/FieldWriterVTK.h b/lib/mito/io/vtk/FieldWriterVTK.h new file mode 100644 index 00000000..7660fc6a --- /dev/null +++ b/lib/mito/io/vtk/FieldWriterVTK.h @@ -0,0 +1,59 @@ +// -*- c++ -*- +// +// Copyright (c) 2020-2024, the MiTo Authors, all rights reserved +// + +// code guard +#pragma once + + +namespace mito::io::vtk { + + template + class FieldWriterVTK { + + private: + // the dimension of the physical space + static constexpr int D = coordSystemT::coordinates_type::dim; + // the grid type (e.g. mesh or point cloud) + using grid_type = gridT; + // the coordinate system type + using coord_system_type = coordSystemT; + // the grid writer type + using grid_writer_type = typename grid_writer::type; + + public: + FieldWriterVTK( + std::string filename, const grid_type & grid, const coord_system_type & coord_system) : + _grid_writer(filename, grid, coord_system) + {} + + // TOFIX: use concepts to say that Y is a tensor + template + auto record(const fem::nodal_field_t & field, std::string fieldname = "") -> void + { + // if no name was provided + if (fieldname == "") { + // use the name of the field + fieldname = field.name(); + } + + // delegate to the grid + return _grid_writer.attach_field(field, fieldname); + } + + auto write() const -> void + { + // delegate to the grid + return _grid_writer.write(); + } + + private: + // the grid writer + grid_writer_type _grid_writer; + }; + +} // namespace mito::io::vtk + + +// end of file diff --git a/lib/mito/io/vtk/api.h b/lib/mito/io/vtk/api.h index f90fc160..d9100ae8 100644 --- a/lib/mito/io/vtk/api.h +++ b/lib/mito/io/vtk/api.h @@ -14,22 +14,26 @@ namespace mito::io::vtk { requires(meshT::cell_type::dim == coordSystemT::coordinates_type::dim) using mesh_writer_t = MeshWriterVTK; - // vtk mesh writer factory - template - requires(meshT::cell_type::dim == coordSystemT::coordinates_type::dim) - auto writer(std::string filename, const meshT & mesh, const coordSystemT & coord_system) - -> mesh_writer_t; - // point cloud writer alias template requires(cloudT::dim == coordSystemT::coordinates_type::dim) using cloud_writer_t = PointCloudWriterVTK; + // field writer alias + template + using field_writer_t = FieldWriterVTK; + + // vtk mesh writer factory + template + requires(meshT::cell_type::dim == coordSystemT::coordinates_type::dim) + auto writer(std::string filename, const meshT & mesh, const coordSystemT & coord_system) + -> field_writer_t; + // point cloud writer factory template requires(cloudT::dim == coordSystemT::coordinates_type::dim) auto writer(std::string filename, const cloudT & cloud, const coordSystemT & coord_system) - -> cloud_writer_t; + -> field_writer_t; } diff --git a/lib/mito/io/vtk/factories.h b/lib/mito/io/vtk/factories.h index 01198882..96dde6fe 100644 --- a/lib/mito/io/vtk/factories.h +++ b/lib/mito/io/vtk/factories.h @@ -13,18 +13,18 @@ namespace mito::io::vtk { template requires(meshT::cell_type::dim == coordSystemT::coordinates_type::dim) auto writer(std::string filename, const meshT & mesh, const coordSystemT & coord_system) - -> mesh_writer_t + -> field_writer_t { - return MeshWriterVTK(filename, mesh, coord_system); + return FieldWriterVTK(filename, mesh, coord_system); } // vtk point cloud writer factory template requires(cloudT::dim == coordSystemT::coordinates_type::dim) auto writer(std::string filename, const cloudT & cloud, const coordSystemT & coord_system) - -> cloud_writer_t + -> field_writer_t { - return PointCloudWriterVTK(filename, cloud, coord_system); + return FieldWriterVTK(filename, cloud, coord_system); } } diff --git a/lib/mito/io/vtk/forward.h b/lib/mito/io/vtk/forward.h index 58d1152f..8520009e 100644 --- a/lib/mito/io/vtk/forward.h +++ b/lib/mito/io/vtk/forward.h @@ -18,6 +18,10 @@ namespace mito::io::vtk { template requires(cloudT::dim == coordSystemT::coordinates_type::dim) class PointCloudWriterVTK; + + // class field writer + template + class FieldWriterVTK; } diff --git a/lib/mito/io/vtk/public.h b/lib/mito/io/vtk/public.h index ffdd34a8..c8225d60 100644 --- a/lib/mito/io/vtk/public.h +++ b/lib/mito/io/vtk/public.h @@ -16,6 +16,9 @@ // published type factories; this is the file you are looking for... #include "api.h" +// utilities +#include "utilities.h" + // wrappers #include "vtk_cell.h" #include "vtk_point.h" @@ -24,6 +27,7 @@ #include "GridWriterVTK.h" #include "MeshWriterVTK.h" #include "PointCloudWriterVTK.h" +#include "FieldWriterVTK.h" // factories implementation #include "factories.h" diff --git a/lib/mito/io/vtk/utilities.h b/lib/mito/io/vtk/utilities.h new file mode 100644 index 00000000..10e4389d --- /dev/null +++ b/lib/mito/io/vtk/utilities.h @@ -0,0 +1,30 @@ +// -*- c++ -*- +// +// Copyright (c) 2020-2024, the MiTo Authors, all rights reserved +// + +// code guard +#pragma once + + +namespace mito::io::vtk { + + // the type of the grid writer depending on the type of grid and on the coordinate system + template + struct grid_writer; + + // specialization to {mesh_c} case + template + struct grid_writer { + using type = mesh_writer_t; + }; + + // specialization to {point_cloud_c} case + template + struct grid_writer { + using type = cloud_writer_t; + }; +} + + +// end of file From 3bbfccc15bc98df317b0406387a335408c2f936e Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Mon, 3 Jun 2024 14:38:57 +0200 Subject: [PATCH 34/58] lib/mito/io/summit: reimplement {summit} mesh writer according to new design of class {MeshWriter} --- lib/mito/io/summit/MeshWriterSummit.h | 95 +++++++++++++++++++++++++++ lib/mito/io/summit/api.h | 19 ++++++ lib/mito/io/summit/public.h | 6 ++ lib/mito/io/summit/writer.h | 54 ++------------- 4 files changed, 125 insertions(+), 49 deletions(-) create mode 100644 lib/mito/io/summit/MeshWriterSummit.h create mode 100644 lib/mito/io/summit/api.h diff --git a/lib/mito/io/summit/MeshWriterSummit.h b/lib/mito/io/summit/MeshWriterSummit.h new file mode 100644 index 00000000..31ceb519 --- /dev/null +++ b/lib/mito/io/summit/MeshWriterSummit.h @@ -0,0 +1,95 @@ +// -*- c++ -*- +// +// Copyright (c) 2020-2024, the MiTo Authors, all rights reserved +// + +// code guard +#pragma once + + +namespace mito::io { + + template + requires(meshT::cell_type::dim == coordSystemT::coordinates_type::dim) + class MeshWriterSummit : public GridWriter { + + private: + // the mesh type + using mesh_type = meshT; + // the coordinate system type + using coord_system_type = coordSystemT; + // the dimension of the physical space + static constexpr int D = mesh_type::cell_type::dim; + + public: + // constructor + MeshWriterSummit( + std::string filename, const mesh_type & mesh, const coord_system_type & coord_system) : + GridWriter(filename), + _mesh(mesh), + _coord_system(coord_system) + { + // insert the points corresponding to the mesh nodes + for (const auto & cell : mesh.cells()) { + for (const auto & node : cell.nodes()) { + this->_points.insert(node.point()); + } + } + } + + // destructor + ~MeshWriterSummit() = default; + + public: + // write mesh to file + virtual auto write() const -> void override + { + // assemble the file name + std::string filename = this->_filename + ".summit"; + + // create the output file + std::ofstream outfile(filename); + + // populate the file heading + // TOFIX: number of materials is always 1 for now + outfile << D << std::endl; + outfile << std::size(this->_points) << " " << _mesh.nCells() << " " << 1 << std::endl; + + // write the points to file + for (const auto & point : this->_points) { + const auto & coord = _coord_system.coordinates(point); + outfile << std::setprecision(15) << coord << std::endl; + } + + // write the cells to file + for (const auto & cell : _mesh.cells()) { + outfile << mesh_type::cell_type::n_vertices << " "; + for (const auto & node : cell.nodes()) { + outfile + << std::distance(this->_points.begin(), this->_points.find(node.point())) + + 1 + << " "; + } + // TOFIX: material label is always 1 for now + outfile << 1 << std::endl; + } + + // close the file + outfile.close(); + + // all done + return; + } + + private: + // a const reference to the mesh + const mesh_type & _mesh; + + // a const reference to the coordinate system + const coord_system_type & _coord_system; + }; + +} // namespace mito::io + + +// end of file diff --git a/lib/mito/io/summit/api.h b/lib/mito/io/summit/api.h new file mode 100644 index 00000000..b5176e3a --- /dev/null +++ b/lib/mito/io/summit/api.h @@ -0,0 +1,19 @@ +// -*- c++ -*- +// +// Copyright (c) 2020-2024, the MiTo Authors, all rights reserved +// + +// code guard +#pragma once + + +namespace mito::io::summit { + + // mesh writer alias + template + requires(meshT::cell_type::dim == coordSystemT::coordinates_type::dim) + using mesh_writer_t = MeshWriterSummit; +} + + +// end of file diff --git a/lib/mito/io/summit/public.h b/lib/mito/io/summit/public.h index deeb5525..39949f4b 100644 --- a/lib/mito/io/summit/public.h +++ b/lib/mito/io/summit/public.h @@ -11,6 +11,12 @@ #include "externals.h" // classes implementation +#include "MeshWriterSummit.h" + +// published type factories; this is the file you are looking for... +#include "api.h" + +// read and write functions #include "reader.h" #include "writer.h" diff --git a/lib/mito/io/summit/writer.h b/lib/mito/io/summit/writer.h index e48aa96f..58e07cba 100644 --- a/lib/mito/io/summit/writer.h +++ b/lib/mito/io/summit/writer.h @@ -12,57 +12,13 @@ namespace mito::io::summit { template requires(cellT::dim == coordT::dim) auto writer( - std::string fileName, const mito::mesh::mesh_t & mesh, + std::string filename, const mito::mesh::mesh_t & mesh, const geometry::coordinate_system_t & coordinate_system) -> void { - // the dimension of the physical space - constexpr int D = cellT::dim; - - // append file extension - fileName = fileName + ".summit"; - - // create the output file - std::ofstream outfile(fileName); - - // type of point - using point_type = geometry::point_t; - - // a set between of points (to remove duplicates) - std::unordered_set> points; - - // insert the points corresponding to the mesh nodes - for (const auto & cell : mesh.cells()) { - for (const auto & node : cell.nodes()) { - points.insert(node.point()); - } - } - - // populate the file heading - // TOFIX: number of materials is always 1 for now - outfile << D << std::endl; - outfile << std::size(points) << " " << mesh.nCells() << " " << 1 << std::endl; - - // write the points to file - for (const auto & point : points) { - const auto & coord = coordinate_system.coordinates(point); - outfile << std::setprecision(15) << coord << std::endl; - } - - // write the cells to file - for (const auto & cell : mesh.cells()) { - outfile << cellT::n_vertices << " "; - for (const auto & node : cell.nodes()) { - outfile << std::distance(points.begin(), points.find(node.point())) + 1 << " "; - } - // TOFIX: material label is always 1 for now - outfile << 1 << std::endl; - } - - // close the file - outfile.close(); - - // all done - return; + // create a writer + auto mesh_writer = mesh_writer_t(filename, mesh, coordinate_system); + // write + return mesh_writer.write(); } } // namespace mito::io::summit From b698516ef04736a1c18ea8b25b20d181de2cbe8f Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Mon, 3 Jun 2024 14:46:57 +0200 Subject: [PATCH 35/58] lib/mito/utilities: add forgotten header --- lib/mito/utilities/externals.h | 1 + 1 file changed, 1 insertion(+) diff --git a/lib/mito/utilities/externals.h b/lib/mito/utilities/externals.h index 333d7097..f3d4a1e4 100644 --- a/lib/mito/utilities/externals.h +++ b/lib/mito/utilities/externals.h @@ -14,6 +14,7 @@ #include #include #include +#include // support #include "../journal.h" From efe4b9a41f8c1c959e80178c8d6a63ad46b21d76 Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Mon, 3 Jun 2024 14:47:50 +0200 Subject: [PATCH 36/58] lib/mito/io: move {unordered_set} external dependency to {io/externals.h} header --- lib/mito/io/externals.h | 2 +- lib/mito/io/summit/externals.h | 1 - 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/lib/mito/io/externals.h b/lib/mito/io/externals.h index e8c3c1ee..deb1c886 100644 --- a/lib/mito/io/externals.h +++ b/lib/mito/io/externals.h @@ -8,7 +8,7 @@ // externals - +#include // support #include "../journal.h" diff --git a/lib/mito/io/summit/externals.h b/lib/mito/io/summit/externals.h index 43484377..f73f0c85 100644 --- a/lib/mito/io/summit/externals.h +++ b/lib/mito/io/summit/externals.h @@ -9,7 +9,6 @@ // externals #include -#include // end of file From 7d13dea12235ffff4dcf2f57de8f4a8c8b39f368 Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Mon, 3 Jun 2024 14:49:13 +0200 Subject: [PATCH 37/58] lib/mito/fem: nodal field factory now infers the space dimension from the mesh --- lib/mito/fem/factories.h | 6 +++--- tests/mito.lib/fem/nodal_field_sphere.cc | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/lib/mito/fem/factories.h b/lib/mito/fem/factories.h index c67ee923..d4424d30 100644 --- a/lib/mito/fem/factories.h +++ b/lib/mito/fem/factories.h @@ -17,10 +17,10 @@ namespace mito::fem { } // nodal field factory - template - constexpr auto nodal_field(const mesh::mesh_c auto & mesh, std::string name) + template + constexpr auto nodal_field(const meshT & mesh, std::string name) { - return nodal_field_t(mesh, name); + return nodal_field_t(mesh, name); } } diff --git a/tests/mito.lib/fem/nodal_field_sphere.cc b/tests/mito.lib/fem/nodal_field_sphere.cc index 9dc99c1a..e6e9d8c1 100644 --- a/tests/mito.lib/fem/nodal_field_sphere.cc +++ b/tests/mito.lib/fem/nodal_field_sphere.cc @@ -23,7 +23,7 @@ TEST(Fem, NodalFieldSphere) auto mesh = mito::io::summit::reader>(fileStream, coord_system); // a nodal field on the mesh - auto nodal_field = mito::fem::nodal_field, 3>(mesh, "normal"); + auto nodal_field = mito::fem::nodal_field>(mesh, "normal"); // the normal field to the submanifold constexpr auto normal_field = mito::fields::field([](const coordinates_t & x) -> auto { From c6534c57588f120ad71ccf3dff7e6594cdf48c0f Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Mon, 3 Jun 2024 14:55:05 +0200 Subject: [PATCH 38/58] lib/mito/fem: use {at} instead of {operator[]} for read-only map access --- lib/mito/fem/NodalField.h | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/lib/mito/fem/NodalField.h b/lib/mito/fem/NodalField.h index d81143b5..9fce4f8b 100644 --- a/lib/mito/fem/NodalField.h +++ b/lib/mito/fem/NodalField.h @@ -55,10 +55,13 @@ namespace mito::fem { */ inline auto operator()(const node_type & node) const -> const Y & { - return _map_nodes_to_values[node]; + return _map_nodes_to_values.at(node); } - inline auto operator()(const node_type & node) -> Y & { return _map_nodes_to_values[node]; } + inline auto operator()(const node_type & node) -> Y & + { + return _map_nodes_to_values.at(node); + } /** * accessor for the number of nodes From 4fb03686d1e3ca1a303cfd86916f2c464b237bec Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Wed, 5 Jun 2024 07:52:23 +0200 Subject: [PATCH 39/58] lib/mito/io: readability improvement Directly access {dim} when possible. --- lib/mito/io/summit/MeshWriterSummit.h | 6 +++--- lib/mito/io/summit/api.h | 2 +- lib/mito/io/vtk/FieldWriterVTK.h | 2 +- lib/mito/io/vtk/MeshWriterVTK.h | 6 +++--- lib/mito/io/vtk/PointCloudWriterVTK.h | 2 +- lib/mito/io/vtk/api.h | 8 ++++---- lib/mito/io/vtk/factories.h | 4 ++-- lib/mito/io/vtk/forward.h | 4 ++-- 8 files changed, 17 insertions(+), 17 deletions(-) diff --git a/lib/mito/io/summit/MeshWriterSummit.h b/lib/mito/io/summit/MeshWriterSummit.h index 31ceb519..6313ae1c 100644 --- a/lib/mito/io/summit/MeshWriterSummit.h +++ b/lib/mito/io/summit/MeshWriterSummit.h @@ -10,8 +10,8 @@ namespace mito::io { template - requires(meshT::cell_type::dim == coordSystemT::coordinates_type::dim) - class MeshWriterSummit : public GridWriter { + requires(meshT::dim == coordSystemT::dim) + class MeshWriterSummit : public GridWriter { private: // the mesh type @@ -19,7 +19,7 @@ namespace mito::io { // the coordinate system type using coord_system_type = coordSystemT; // the dimension of the physical space - static constexpr int D = mesh_type::cell_type::dim; + static constexpr int D = mesh_type::dim; public: // constructor diff --git a/lib/mito/io/summit/api.h b/lib/mito/io/summit/api.h index b5176e3a..893ea9bc 100644 --- a/lib/mito/io/summit/api.h +++ b/lib/mito/io/summit/api.h @@ -11,7 +11,7 @@ namespace mito::io::summit { // mesh writer alias template - requires(meshT::cell_type::dim == coordSystemT::coordinates_type::dim) + requires(meshT::dim == coordSystemT::dim) using mesh_writer_t = MeshWriterSummit; } diff --git a/lib/mito/io/vtk/FieldWriterVTK.h b/lib/mito/io/vtk/FieldWriterVTK.h index 7660fc6a..772ae144 100644 --- a/lib/mito/io/vtk/FieldWriterVTK.h +++ b/lib/mito/io/vtk/FieldWriterVTK.h @@ -14,7 +14,7 @@ namespace mito::io::vtk { private: // the dimension of the physical space - static constexpr int D = coordSystemT::coordinates_type::dim; + static constexpr int D = coordSystemT::dim; // the grid type (e.g. mesh or point cloud) using grid_type = gridT; // the coordinate system type diff --git a/lib/mito/io/vtk/MeshWriterVTK.h b/lib/mito/io/vtk/MeshWriterVTK.h index d72459df..93c3e3e0 100644 --- a/lib/mito/io/vtk/MeshWriterVTK.h +++ b/lib/mito/io/vtk/MeshWriterVTK.h @@ -10,8 +10,8 @@ namespace mito::io::vtk { template - requires(meshT::cell_type::dim == coordSystemT::coordinates_type::dim) - class MeshWriterVTK : public GridWriterVTK { + requires(meshT::dim == coordSystemT::dim) + class MeshWriterVTK : public GridWriterVTK { private: // the mesh type @@ -19,7 +19,7 @@ namespace mito::io::vtk { // the coordinate system type using coord_system_type = coordSystemT; // the dimension of the physical space - static constexpr int D = mesh_type::cell_type::dim; + static constexpr int D = mesh_type::dim; private: auto _create_vtk_grid(const mesh_type & mesh, const coord_system_type & coordinate_system) diff --git a/lib/mito/io/vtk/PointCloudWriterVTK.h b/lib/mito/io/vtk/PointCloudWriterVTK.h index 5f507148..062e9025 100644 --- a/lib/mito/io/vtk/PointCloudWriterVTK.h +++ b/lib/mito/io/vtk/PointCloudWriterVTK.h @@ -10,7 +10,7 @@ namespace mito::io::vtk { template - requires(cloudT::dim == coordSystemT::coordinates_type::dim) + requires(cloudT::dim == coordSystemT::dim) class PointCloudWriterVTK : public GridWriterVTK { private: diff --git a/lib/mito/io/vtk/api.h b/lib/mito/io/vtk/api.h index d9100ae8..14d2e782 100644 --- a/lib/mito/io/vtk/api.h +++ b/lib/mito/io/vtk/api.h @@ -11,12 +11,12 @@ namespace mito::io::vtk { // mesh writer alias template - requires(meshT::cell_type::dim == coordSystemT::coordinates_type::dim) + requires(meshT::dim == coordSystemT::dim) using mesh_writer_t = MeshWriterVTK; // point cloud writer alias template - requires(cloudT::dim == coordSystemT::coordinates_type::dim) + requires(cloudT::dim == coordSystemT::dim) using cloud_writer_t = PointCloudWriterVTK; // field writer alias @@ -25,13 +25,13 @@ namespace mito::io::vtk { // vtk mesh writer factory template - requires(meshT::cell_type::dim == coordSystemT::coordinates_type::dim) + requires(meshT::dim == coordSystemT::dim) auto writer(std::string filename, const meshT & mesh, const coordSystemT & coord_system) -> field_writer_t; // point cloud writer factory template - requires(cloudT::dim == coordSystemT::coordinates_type::dim) + requires(cloudT::dim == coordSystemT::dim) auto writer(std::string filename, const cloudT & cloud, const coordSystemT & coord_system) -> field_writer_t; } diff --git a/lib/mito/io/vtk/factories.h b/lib/mito/io/vtk/factories.h index 96dde6fe..f540bbb7 100644 --- a/lib/mito/io/vtk/factories.h +++ b/lib/mito/io/vtk/factories.h @@ -11,7 +11,7 @@ namespace mito::io::vtk { // vtk mesh writer factory template - requires(meshT::cell_type::dim == coordSystemT::coordinates_type::dim) + requires(meshT::dim == coordSystemT::dim) auto writer(std::string filename, const meshT & mesh, const coordSystemT & coord_system) -> field_writer_t { @@ -20,7 +20,7 @@ namespace mito::io::vtk { // vtk point cloud writer factory template - requires(cloudT::dim == coordSystemT::coordinates_type::dim) + requires(cloudT::dim == coordSystemT::dim) auto writer(std::string filename, const cloudT & cloud, const coordSystemT & coord_system) -> field_writer_t { diff --git a/lib/mito/io/vtk/forward.h b/lib/mito/io/vtk/forward.h index 8520009e..adcf0e4a 100644 --- a/lib/mito/io/vtk/forward.h +++ b/lib/mito/io/vtk/forward.h @@ -11,12 +11,12 @@ namespace mito::io::vtk { // class mesh writer template - requires(meshT::cell_type::dim == coordSystemT::coordinates_type::dim) + requires(meshT::dim == coordSystemT::dim) class MeshWriterVTK; // class point cloud writer template - requires(cloudT::dim == coordSystemT::coordinates_type::dim) + requires(cloudT::dim == coordSystemT::dim) class PointCloudWriterVTK; // class field writer From 6f83fef29833a643bea614ed8094a6ff7967be47 Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Wed, 5 Jun 2024 08:13:23 +0200 Subject: [PATCH 40/58] lib/mito: add {concept} of two types having the same dimension --- lib/mito/io/summit/MeshWriterSummit.h | 2 +- lib/mito/io/summit/api.h | 2 +- lib/mito/io/summit/reader.h | 2 +- lib/mito/io/summit/writer.h | 2 +- lib/mito/io/vtk/MeshWriterVTK.h | 2 +- lib/mito/io/vtk/PointCloudWriterVTK.h | 2 +- lib/mito/io/vtk/api.h | 8 ++++---- lib/mito/io/vtk/factories.h | 4 ++-- lib/mito/io/vtk/forward.h | 4 ++-- lib/mito/utilities/api.h | 4 ++++ 10 files changed, 18 insertions(+), 14 deletions(-) diff --git a/lib/mito/io/summit/MeshWriterSummit.h b/lib/mito/io/summit/MeshWriterSummit.h index 6313ae1c..ee2c2cab 100644 --- a/lib/mito/io/summit/MeshWriterSummit.h +++ b/lib/mito/io/summit/MeshWriterSummit.h @@ -10,7 +10,7 @@ namespace mito::io { template - requires(meshT::dim == coordSystemT::dim) + requires(utilities::same_dim_c) class MeshWriterSummit : public GridWriter { private: diff --git a/lib/mito/io/summit/api.h b/lib/mito/io/summit/api.h index 893ea9bc..ccf7e5d2 100644 --- a/lib/mito/io/summit/api.h +++ b/lib/mito/io/summit/api.h @@ -11,7 +11,7 @@ namespace mito::io::summit { // mesh writer alias template - requires(meshT::dim == coordSystemT::dim) + requires(utilities::same_dim_c) using mesh_writer_t = MeshWriterSummit; } diff --git a/lib/mito/io/summit/reader.h b/lib/mito/io/summit/reader.h index 40d45a5c..bfdbea4e 100644 --- a/lib/mito/io/summit/reader.h +++ b/lib/mito/io/summit/reader.h @@ -218,7 +218,7 @@ namespace mito::io::summit { auto reader( std::ifstream & fileStream, geometry::coordinate_system_t & coordinate_system) -> mesh::mesh_t - requires(cellT::dim == coordT::dim) + requires(utilities::same_dim_c) { if (!fileStream.is_open()) { throw std::runtime_error("reader: Mesh file could not be opened"); diff --git a/lib/mito/io/summit/writer.h b/lib/mito/io/summit/writer.h index 58e07cba..ae5d04e3 100644 --- a/lib/mito/io/summit/writer.h +++ b/lib/mito/io/summit/writer.h @@ -10,7 +10,7 @@ namespace mito::io::summit { template - requires(cellT::dim == coordT::dim) + requires(utilities::same_dim_c) auto writer( std::string filename, const mito::mesh::mesh_t & mesh, const geometry::coordinate_system_t & coordinate_system) -> void diff --git a/lib/mito/io/vtk/MeshWriterVTK.h b/lib/mito/io/vtk/MeshWriterVTK.h index 93c3e3e0..9e15811d 100644 --- a/lib/mito/io/vtk/MeshWriterVTK.h +++ b/lib/mito/io/vtk/MeshWriterVTK.h @@ -10,7 +10,7 @@ namespace mito::io::vtk { template - requires(meshT::dim == coordSystemT::dim) + requires(utilities::same_dim_c) class MeshWriterVTK : public GridWriterVTK { private: diff --git a/lib/mito/io/vtk/PointCloudWriterVTK.h b/lib/mito/io/vtk/PointCloudWriterVTK.h index 062e9025..813c5ccc 100644 --- a/lib/mito/io/vtk/PointCloudWriterVTK.h +++ b/lib/mito/io/vtk/PointCloudWriterVTK.h @@ -10,7 +10,7 @@ namespace mito::io::vtk { template - requires(cloudT::dim == coordSystemT::dim) + requires(utilities::same_dim_c) class PointCloudWriterVTK : public GridWriterVTK { private: diff --git a/lib/mito/io/vtk/api.h b/lib/mito/io/vtk/api.h index 14d2e782..ee785e12 100644 --- a/lib/mito/io/vtk/api.h +++ b/lib/mito/io/vtk/api.h @@ -11,12 +11,12 @@ namespace mito::io::vtk { // mesh writer alias template - requires(meshT::dim == coordSystemT::dim) + requires(utilities::same_dim_c) using mesh_writer_t = MeshWriterVTK; // point cloud writer alias template - requires(cloudT::dim == coordSystemT::dim) + requires(utilities::same_dim_c) using cloud_writer_t = PointCloudWriterVTK; // field writer alias @@ -25,13 +25,13 @@ namespace mito::io::vtk { // vtk mesh writer factory template - requires(meshT::dim == coordSystemT::dim) + requires(utilities::same_dim_c) auto writer(std::string filename, const meshT & mesh, const coordSystemT & coord_system) -> field_writer_t; // point cloud writer factory template - requires(cloudT::dim == coordSystemT::dim) + requires(utilities::same_dim_c) auto writer(std::string filename, const cloudT & cloud, const coordSystemT & coord_system) -> field_writer_t; } diff --git a/lib/mito/io/vtk/factories.h b/lib/mito/io/vtk/factories.h index f540bbb7..3eff42e0 100644 --- a/lib/mito/io/vtk/factories.h +++ b/lib/mito/io/vtk/factories.h @@ -11,7 +11,7 @@ namespace mito::io::vtk { // vtk mesh writer factory template - requires(meshT::dim == coordSystemT::dim) + requires(utilities::same_dim_c) auto writer(std::string filename, const meshT & mesh, const coordSystemT & coord_system) -> field_writer_t { @@ -20,7 +20,7 @@ namespace mito::io::vtk { // vtk point cloud writer factory template - requires(cloudT::dim == coordSystemT::dim) + requires(utilities::same_dim_c) auto writer(std::string filename, const cloudT & cloud, const coordSystemT & coord_system) -> field_writer_t { diff --git a/lib/mito/io/vtk/forward.h b/lib/mito/io/vtk/forward.h index adcf0e4a..e66fc5de 100644 --- a/lib/mito/io/vtk/forward.h +++ b/lib/mito/io/vtk/forward.h @@ -11,12 +11,12 @@ namespace mito::io::vtk { // class mesh writer template - requires(meshT::dim == coordSystemT::dim) + requires(utilities::same_dim_c) class MeshWriterVTK; // class point cloud writer template - requires(cloudT::dim == coordSystemT::dim) + requires(utilities::same_dim_c) class PointCloudWriterVTK; // class field writer diff --git a/lib/mito/utilities/api.h b/lib/mito/utilities/api.h index 84855179..1fae0542 100644 --- a/lib/mito/utilities/api.h +++ b/lib/mito/utilities/api.h @@ -28,6 +28,10 @@ namespace mito::utilities { // segmented vector alias template using segmented_vector_t = SegmentedVector; + + // concept of the types having the same dimension + template + concept same_dim_c = requires { F1::dim == F2::dim; }; } From 68df98c97085820e7eeb4e7cd241aaaa5bd1f65d Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Sat, 8 Jun 2024 09:04:09 +0200 Subject: [PATCH 41/58] lib/mito: draft implementation of reference counted nodes This is necessary so as to avoid having the same node being copied into different nodes when inserting a geometric simplex in a mesh. Copying a node now results in copying its shared pointer. The current implementation has some issues, though: 1) While originally the {SharedPointer} class was designed for not owning the resource, now a {SharedPointer} may instantiate a resource. This breaks the design and a mechanism for deletion of the resource should be implemented. 2) Nodes are now 'special' geometric simplices, in that they do reference counting while the other geometric simplices do not. This implies that sometimes they need to be treated as special cases (for example in the extraction of the boundary of a segment). --- lib/mito/fem/NodalField.h | 15 +------------ lib/mito/geometry/GeometricSimplex.h | 10 ++++----- lib/mito/geometry/api.h | 2 +- lib/mito/geometry/utilities.h | 17 +++++++------- lib/mito/io/summit/MeshWriterSummit.h | 4 ++-- lib/mito/io/vtk/GridWriterVTK.h | 3 ++- lib/mito/io/vtk/MeshWriterVTK.h | 4 ++-- lib/mito/manifolds/Manifold.h | 2 +- lib/mito/mesh/metis/Partitioner.icc | 2 +- lib/mito/mesh/tetra.h | 2 +- lib/mito/utilities/SharedPointer.h | 4 ++++ lib/mito/utilities/SharedPointer.icc | 8 +++++++ tests/mito.lib/fem/nodal_field_sphere.cc | 2 +- tests/mito.lib/geometry/cube_volume.cc | 2 +- .../mito.lib/geometry/euclidean_metric_2D.cc | 6 ++--- .../mito.lib/geometry/euclidean_metric_3D.cc | 8 +++---- .../euclidean_submanifold_metric_3D.cc | 6 ++--- .../io/read_mesh_summit_segment_3D.cc | 22 ++++++++++--------- tests/mito.lib/mesh/erase_element.cc | 18 +++++++-------- 19 files changed, 70 insertions(+), 67 deletions(-) diff --git a/lib/mito/fem/NodalField.h b/lib/mito/fem/NodalField.h index 9fce4f8b..e389ce01 100644 --- a/lib/mito/fem/NodalField.h +++ b/lib/mito/fem/NodalField.h @@ -16,21 +16,8 @@ namespace mito::fem { // the node type using node_type = geometry::node_t; - // hash function for nodes - struct hash_function { - auto operator()(const node_type & node) const -> uintptr_t - { - // TOFIX: this hash function is not safe against collisions: either figure out a - // proper hash function or turn nodes into something shareable, so that the hash - // can simply be a check on the address of the underlying resource. We'll figure - // out which one of two routes to pursue when we do discontinuous Galerkin - // build a hash based on the id of the node and that of the point - return node.vertex().id() + node.point().id(); - } - }; - // a map from {key_type} to {Y} values - using map_type = std::unordered_map; + using map_type = std::unordered_map>; public: // constructor diff --git a/lib/mito/geometry/GeometricSimplex.h b/lib/mito/geometry/GeometricSimplex.h index 601c616c..cdc0598f 100644 --- a/lib/mito/geometry/GeometricSimplex.h +++ b/lib/mito/geometry/GeometricSimplex.h @@ -58,7 +58,7 @@ namespace mito::geometry { // of nodes in {_nodes} auto _check_vertices = [this](integer_sequence) -> bool { return ( - math::permutation_sign(_simplex->vertices(), { _nodes[J].vertex()... }) == +1); + math::permutation_sign(_simplex->vertices(), { _nodes[J]->vertex()... }) == +1); }; return _check_vertices(make_integer_sequence{}); @@ -71,7 +71,7 @@ namespace mito::geometry { auto & topology = topology::topology(); // instantiate a simplex with the vertices prescribed by {_nodes} - return topology.simplex({ _nodes[J].vertex()... }); + return topology.simplex({ _nodes[J]->vertex()... }); } public: @@ -131,7 +131,7 @@ namespace mito::geometry { template requires(D > 0) - class GeometricSimplex<0, D> : public utilities::Invalidatable { + class GeometricSimplex<0, D> : public utilities::Invalidatable, public utilities::Shareable { public: // spatial dimension @@ -207,8 +207,8 @@ namespace mito::geometry { constexpr auto operator==( const GeometricSimplex<0, D> & node_a, const GeometricSimplex<0, D> & node_b) -> bool { - // two nodes are the same if their vertex and point are the same - return node_a.vertex() == node_b.vertex() && node_a.point() == node_b.point(); + // two nodes are the same if they have the same id + return node_a.id() == node_b.id(); } } diff --git a/lib/mito/geometry/api.h b/lib/mito/geometry/api.h index 1b06f611..0f427d1a 100644 --- a/lib/mito/geometry/api.h +++ b/lib/mito/geometry/api.h @@ -23,7 +23,7 @@ namespace mito::geometry { // node alias template - using node_t = geometric_simplex_t<0, D>; + using node_t = utilities::shared_ptr>; // segment alias template diff --git a/lib/mito/geometry/utilities.h b/lib/mito/geometry/utilities.h index a0c9d2b6..0886d9fd 100644 --- a/lib/mito/geometry/utilities.h +++ b/lib/mito/geometry/utilities.h @@ -15,14 +15,14 @@ namespace mito::geometry { const cellT & cell, const coordinate_system_t & coordinate_system) -> coordT { // get the coordinates of the first node - auto coord_0 = coordinate_system.coordinates(cell.nodes()[0].point()); + auto coord_0 = coordinate_system.coordinates(cell.nodes()[0]->point()); // the vector going from {coord_0} to the barycenter (initialize with the zero vector) auto result = coord_0 - coord_0; // average the position of each vertex for (const auto & node : cell.nodes()) { - result += coordinate_system.coordinates(node.point()) - coord_0; + result += coordinate_system.coordinates(node->point()) - coord_0; } result /= cellT::n_vertices; @@ -61,11 +61,11 @@ namespace mito::geometry { auto nodes = simplex.nodes(); // get the coordinates of the first node - const auto & p0 = coordinate_system.coordinates(nodes[0].point()); + const auto & p0 = coordinate_system.coordinates(nodes[0]->point()); // compute the director vectors associated with each director edge auto directors = edge_simplex_directors_t{ ( - coordinate_system.coordinates(nodes[J + 1].point()) - p0)... }; + coordinate_system.coordinates(nodes[J + 1]->point()) - p0)... }; // all done return { p0, directors }; @@ -96,14 +96,15 @@ namespace mito::geometry { // given vertex auto has_vertex = [](const vertex_type & vertex) { auto lambda = [&vertex](const node_type & node) { - return node.vertex() == vertex; + return node->vertex() == vertex; }; return lambda; }; // return the geometric simplex return geometric_simplex_type({ node_type( - vertices[K], std::ranges::find_if(nodes, has_vertex(vertices[K]))->point())... }); + vertices[K], + (*std::ranges::find_if(nodes, has_vertex(vertices[K])))->point())... }); }; // build a geometric simplex based on {simplex} with the vertex-point pair as appears in @@ -127,13 +128,13 @@ namespace mito::geometry { // given vertex auto has_vertex = [](const vertex_type & vertex) { auto lambda = [&vertex](const node_type & node) { - return node.vertex() == vertex; + return node->vertex() == vertex; }; return lambda; }; // return the node - return node_type(vertex, std::ranges::find_if(nodes, has_vertex(vertex))->point()); + return node_type(vertex, (*std::ranges::find_if(nodes, has_vertex(vertex)))->point()); } // returns the geometric simplex with opposite orientation to {simplex} diff --git a/lib/mito/io/summit/MeshWriterSummit.h b/lib/mito/io/summit/MeshWriterSummit.h index ee2c2cab..90e5fa8a 100644 --- a/lib/mito/io/summit/MeshWriterSummit.h +++ b/lib/mito/io/summit/MeshWriterSummit.h @@ -32,7 +32,7 @@ namespace mito::io { // insert the points corresponding to the mesh nodes for (const auto & cell : mesh.cells()) { for (const auto & node : cell.nodes()) { - this->_points.insert(node.point()); + this->_points.insert(node->point()); } } } @@ -66,7 +66,7 @@ namespace mito::io { outfile << mesh_type::cell_type::n_vertices << " "; for (const auto & node : cell.nodes()) { outfile - << std::distance(this->_points.begin(), this->_points.find(node.point())) + << std::distance(this->_points.begin(), this->_points.find(node->point())) + 1 << " "; } diff --git a/lib/mito/io/vtk/GridWriterVTK.h b/lib/mito/io/vtk/GridWriterVTK.h index 9d9b4552..5168de7d 100644 --- a/lib/mito/io/vtk/GridWriterVTK.h +++ b/lib/mito/io/vtk/GridWriterVTK.h @@ -58,7 +58,8 @@ namespace mito::io::vtk { // populate the array with the nodal values for (auto & [node, value] : field) { // get the index corresponding to the current point - auto index = std::distance(this->_points.begin(), this->_points.find(node.point())); + auto index = + std::distance(this->_points.begin(), this->_points.find(node->point())); vtkArray->SetTuple(index, value.begin()); } diff --git a/lib/mito/io/vtk/MeshWriterVTK.h b/lib/mito/io/vtk/MeshWriterVTK.h index 9e15811d..80574fe4 100644 --- a/lib/mito/io/vtk/MeshWriterVTK.h +++ b/lib/mito/io/vtk/MeshWriterVTK.h @@ -29,7 +29,7 @@ namespace mito::io::vtk { // loop over the nodes of the cell for (const auto & node : cell.nodes()) { // add the point to the collection of points (eliminating duplicates) - this->_points.insert(node.point()); + this->_points.insert(node->point()); } } @@ -53,7 +53,7 @@ namespace mito::io::vtk { // loop over the nodes of the cell for (const auto & node : cell.nodes()) { // retrieve the corresponding point - const auto & point = node.point(); + const auto & point = node->point(); // assert that the point is present in the set of points assert(this->_points.contains(point)); // calculate the index of the point diff --git a/lib/mito/manifolds/Manifold.h b/lib/mito/manifolds/Manifold.h index d62a5ebd..c4952c73 100644 --- a/lib/mito/manifolds/Manifold.h +++ b/lib/mito/manifolds/Manifold.h @@ -75,7 +75,7 @@ namespace mito::manifolds { constexpr auto coordinates(const node_type & v) const -> const coordinates_type & { // get the coordinates of the point attached to vertex {v} - return _coordinate_system.coordinates(v.point()); + return _coordinate_system.coordinates(v->point()); } constexpr auto parametrization( diff --git a/lib/mito/mesh/metis/Partitioner.icc b/lib/mito/mesh/metis/Partitioner.icc index 01814b90..7bf42cf9 100644 --- a/lib/mito/mesh/metis/Partitioner.icc +++ b/lib/mito/mesh/metis/Partitioner.icc @@ -60,7 +60,7 @@ mito::mesh::metis::Partitioner::_populate_element_connectivity( int v = 0; for (const auto & node : cell.nodes()) { // populate element connectivity with the id of each of its vertices - element_connectivity[n_vertices_per_element * e + v] = vertex_to_id.at(node.vertex()); + element_connectivity[n_vertices_per_element * e + v] = vertex_to_id.at(node->vertex()); // increment ++v; } diff --git a/lib/mito/mesh/tetra.h b/lib/mito/mesh/tetra.h index 16f0f428..f58c7727 100644 --- a/lib/mito/mesh/tetra.h +++ b/lib/mito/mesh/tetra.h @@ -17,7 +17,7 @@ namespace mito::mesh { { // return a new node at the midpoint between {node_a} and {node_b} return mito::geometry::node( - coordinate_system, coordinate_system.midpoint(node_a.point(), node_b.point())); + coordinate_system, coordinate_system.midpoint(node_a->point(), node_b->point())); } template diff --git a/lib/mito/utilities/SharedPointer.h b/lib/mito/utilities/SharedPointer.h index e1e1a9b1..9f963c13 100644 --- a/lib/mito/utilities/SharedPointer.h +++ b/lib/mito/utilities/SharedPointer.h @@ -45,6 +45,10 @@ namespace mito::utilities { // default constructor inline SharedPointer(); + template + requires(std::is_constructible_v) + inline SharedPointer(Args &&... args); + // constructor inline SharedPointer(handle_type); diff --git a/lib/mito/utilities/SharedPointer.icc b/lib/mito/utilities/SharedPointer.icc index b2db4161..e5ed803b 100644 --- a/lib/mito/utilities/SharedPointer.icc +++ b/lib/mito/utilities/SharedPointer.icc @@ -105,6 +105,14 @@ mito::utilities::SharedPointer::SharedPointer(handle_type resource) : _acquire(); } +// TOFIX: delete? Implement a memory-owning version and a non memory owning one +template +template +requires(std::is_constructible_v) +mito::utilities::SharedPointer::SharedPointer(Args &&... args) : + _handle(new resourceT(args...)) +{} + // the copy constructor template mito::utilities::SharedPointer::SharedPointer( diff --git a/tests/mito.lib/fem/nodal_field_sphere.cc b/tests/mito.lib/fem/nodal_field_sphere.cc index e6e9d8c1..42a9f175 100644 --- a/tests/mito.lib/fem/nodal_field_sphere.cc +++ b/tests/mito.lib/fem/nodal_field_sphere.cc @@ -36,7 +36,7 @@ TEST(Fem, NodalFieldSphere) // fill information in nodal field for (auto & [node, value] : nodal_field) { // get the coordinates of the node - auto & coordinates = coord_system.coordinates(node.point()); + auto & coordinates = coord_system.coordinates(node->point()); // compute the value of the normal field at those coordinates value = normal_field(coordinates); } diff --git a/tests/mito.lib/geometry/cube_volume.cc b/tests/mito.lib/geometry/cube_volume.cc index e13cce4e..0993d452 100644 --- a/tests/mito.lib/geometry/cube_volume.cc +++ b/tests/mito.lib/geometry/cube_volume.cc @@ -58,7 +58,7 @@ volume_determinant( for (const auto & node : element_nodes) { // fill up pointsTensor container for (int d = 0; d < D; ++d) { - pointsTensor[v * V + d] = coord_system.coordinates(node.point())[d]; + pointsTensor[v * V + d] = coord_system.coordinates(node->point())[d]; } pointsTensor[v * V + D] = 1.0; // update element vertices counter diff --git a/tests/mito.lib/geometry/euclidean_metric_2D.cc b/tests/mito.lib/geometry/euclidean_metric_2D.cc index d48ba02b..d82a4e5f 100644 --- a/tests/mito.lib/geometry/euclidean_metric_2D.cc +++ b/tests/mito.lib/geometry/euclidean_metric_2D.cc @@ -22,9 +22,9 @@ area( const mito::geometry::node_t<2> & v2) -> mito::scalar_t { // get vertex coordinates - auto x0 = coordinate_system.coordinates(v0.point()); - auto x1 = coordinate_system.coordinates(v1.point()); - auto x2 = coordinate_system.coordinates(v2.point()); + auto x0 = coordinate_system.coordinates(v0->point()); + auto x1 = coordinate_system.coordinates(v1->point()); + auto x2 = coordinate_system.coordinates(v2->point()); // build director vectors auto director0 = x1 - x0; diff --git a/tests/mito.lib/geometry/euclidean_metric_3D.cc b/tests/mito.lib/geometry/euclidean_metric_3D.cc index 3e4fa8b9..a253b7d6 100644 --- a/tests/mito.lib/geometry/euclidean_metric_3D.cc +++ b/tests/mito.lib/geometry/euclidean_metric_3D.cc @@ -23,10 +23,10 @@ volume( const mito::geometry::node_t<3> & v2, const mito::geometry::node_t<3> & v3) -> mito::scalar_t { // build director segments - auto x0 = coordinate_system.coordinates(v0.point()); - auto x1 = coordinate_system.coordinates(v1.point()); - auto x2 = coordinate_system.coordinates(v2.point()); - auto x3 = coordinate_system.coordinates(v3.point()); + auto x0 = coordinate_system.coordinates(v0->point()); + auto x1 = coordinate_system.coordinates(v1->point()); + auto x2 = coordinate_system.coordinates(v2->point()); + auto x3 = coordinate_system.coordinates(v3->point()); // build director vectors auto director0 = x1 - x0; diff --git a/tests/mito.lib/geometry/euclidean_submanifold_metric_3D.cc b/tests/mito.lib/geometry/euclidean_submanifold_metric_3D.cc index 3921203e..b8e2fd69 100644 --- a/tests/mito.lib/geometry/euclidean_submanifold_metric_3D.cc +++ b/tests/mito.lib/geometry/euclidean_submanifold_metric_3D.cc @@ -27,9 +27,9 @@ area( const mito::geometry::node_t<3> & v2) -> mito::scalar_t { // get vertex coordinates - auto x0 = coordinate_system.coordinates(v0.point()); - auto x1 = coordinate_system.coordinates(v1.point()); - auto x2 = coordinate_system.coordinates(v2.point()); + auto x0 = coordinate_system.coordinates(v0->point()); + auto x1 = coordinate_system.coordinates(v1->point()); + auto x2 = coordinate_system.coordinates(v2->point()); // build director vectors auto director0 = x1 - x0; diff --git a/tests/mito.lib/io/read_mesh_summit_segment_3D.cc b/tests/mito.lib/io/read_mesh_summit_segment_3D.cc index 5507863a..6ae9136a 100644 --- a/tests/mito.lib/io/read_mesh_summit_segment_3D.cc +++ b/tests/mito.lib/io/read_mesh_summit_segment_3D.cc @@ -23,11 +23,12 @@ TEST(SummitReader, LoadSummitSegmentsMesh3D) // assert you read 10 cells EXPECT_EQ(mesh.nCells(), 10); - // assert you found 2 nodes on the boundary - { - auto boundary_mesh = mito::mesh::boundary(mesh); - EXPECT_EQ(boundary_mesh.nCells(), 2); - } + // TOFIX + // // assert you found 2 nodes on the boundary + // { + // auto boundary_mesh = mito::mesh::boundary(mesh); + // EXPECT_EQ(boundary_mesh.nCells(), 2); + // } // get a cell auto & cell = mesh.cells()[5]; @@ -38,9 +39,10 @@ TEST(SummitReader, LoadSummitSegmentsMesh3D) // assert you read 9 cells EXPECT_EQ(mesh.nCells(), 9); - // assert you found 4 nodes on the boundary - { - auto boundary_mesh = mito::mesh::boundary(mesh); - EXPECT_EQ(boundary_mesh.nCells(), 4); - } + // TOFIX + // // assert you found 4 nodes on the boundary + // { + // auto boundary_mesh = mito::mesh::boundary(mesh); + // EXPECT_EQ(boundary_mesh.nCells(), 4); + // } } \ No newline at end of file diff --git a/tests/mito.lib/mesh/erase_element.cc b/tests/mito.lib/mesh/erase_element.cc index 9fdeeec0..111a5888 100644 --- a/tests/mito.lib/mesh/erase_element.cc +++ b/tests/mito.lib/mesh/erase_element.cc @@ -40,11 +40,11 @@ TEST(Mesh, EraseElement) EXPECT_EQ(mito::mesh::boundary_size(mesh), 4); // assert that there exists a segment connecting vertex 0 and 1 - EXPECT_TRUE(topology.exists({ node_0.vertex(), node_1.vertex() })); + EXPECT_TRUE(topology.exists({ node_0->vertex(), node_1->vertex() })); // assert that there exists a segment connecting vertex 1 and 3 - EXPECT_TRUE(topology.exists({ node_1.vertex(), node_3.vertex() })); + EXPECT_TRUE(topology.exists({ node_1->vertex(), node_3->vertex() })); // assert that there exists a segment connecting vertex 3 and 0 - EXPECT_TRUE(topology.exists({ node_3.vertex(), node_0.vertex() })); + EXPECT_TRUE(topology.exists({ node_3->vertex(), node_0->vertex() })); // erase a simplex mesh.erase(geom_cell0); @@ -57,11 +57,11 @@ TEST(Mesh, EraseElement) EXPECT_EQ(mito::mesh::boundary_size(mesh), 5); // assert that a segment connecting vertex 0 and 1 no longer exists in the topology - EXPECT_FALSE(topology.exists({ node_0.vertex(), node_1.vertex() })); + EXPECT_FALSE(topology.exists({ node_0->vertex(), node_1->vertex() })); // assert that a segment connecting vertex 1 and 3 no longer exists in the topology - EXPECT_FALSE(topology.exists({ node_1.vertex(), node_3.vertex() })); + EXPECT_FALSE(topology.exists({ node_1->vertex(), node_3->vertex() })); // assert that a segment connecting vertex 3 and 0 no longer exists in the topology - EXPECT_FALSE(topology.exists({ node_3.vertex(), node_0.vertex() })); + EXPECT_FALSE(topology.exists({ node_3->vertex(), node_0->vertex() })); // check that erasing a cell twice does not result in an error mesh.erase(geom_cell0); @@ -77,9 +77,9 @@ TEST(Mesh, EraseElement) EXPECT_EQ(mito::mesh::boundary_size(mesh), 4); // assert that a segment connecting vertex 1 and 2 no longer exists in the topology - EXPECT_FALSE(topology.exists({ node_1.vertex(), node_2.vertex() })); + EXPECT_FALSE(topology.exists({ node_1->vertex(), node_2->vertex() })); // assert that a segment connecting vertex 2 and 3 no longer exists in the topology - EXPECT_FALSE(topology.exists({ node_2.vertex(), node_3.vertex() })); + EXPECT_FALSE(topology.exists({ node_2->vertex(), node_3->vertex() })); // assert that a segment connecting vertex 3 and 1 no longer exists in the topology - EXPECT_FALSE(topology.exists({ node_3.vertex(), node_1.vertex() })); + EXPECT_FALSE(topology.exists({ node_3->vertex(), node_1->vertex() })); } From a31e1aaa183a10bced9a7a250277f9d69e554fcd Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Sat, 8 Jun 2024 09:46:26 +0200 Subject: [PATCH 42/58] lib/mito/utilities: add wrapper class for {std::shared_ptr} This class extends the {std::shared_ptr} in order to provide a constructor with a suitable argument list, as well as an id, which enables hashing and <==> operators. --- lib/mito/utilities/SharedPointer.h | 4 -- lib/mito/utilities/SharedPointer.icc | 8 --- lib/mito/utilities/StdSharedPointer.h | 72 +++++++++++++++++++++++++ lib/mito/utilities/StdSharedPointer.icc | 27 ++++++++++ lib/mito/utilities/api.h | 4 ++ lib/mito/utilities/forward.h | 4 ++ lib/mito/utilities/public.h | 1 + 7 files changed, 108 insertions(+), 12 deletions(-) create mode 100644 lib/mito/utilities/StdSharedPointer.h create mode 100644 lib/mito/utilities/StdSharedPointer.icc diff --git a/lib/mito/utilities/SharedPointer.h b/lib/mito/utilities/SharedPointer.h index 9f963c13..e1e1a9b1 100644 --- a/lib/mito/utilities/SharedPointer.h +++ b/lib/mito/utilities/SharedPointer.h @@ -45,10 +45,6 @@ namespace mito::utilities { // default constructor inline SharedPointer(); - template - requires(std::is_constructible_v) - inline SharedPointer(Args &&... args); - // constructor inline SharedPointer(handle_type); diff --git a/lib/mito/utilities/SharedPointer.icc b/lib/mito/utilities/SharedPointer.icc index e5ed803b..b2db4161 100644 --- a/lib/mito/utilities/SharedPointer.icc +++ b/lib/mito/utilities/SharedPointer.icc @@ -105,14 +105,6 @@ mito::utilities::SharedPointer::SharedPointer(handle_type resource) : _acquire(); } -// TOFIX: delete? Implement a memory-owning version and a non memory owning one -template -template -requires(std::is_constructible_v) -mito::utilities::SharedPointer::SharedPointer(Args &&... args) : - _handle(new resourceT(args...)) -{} - // the copy constructor template mito::utilities::SharedPointer::SharedPointer( diff --git a/lib/mito/utilities/StdSharedPointer.h b/lib/mito/utilities/StdSharedPointer.h new file mode 100644 index 00000000..d36e8650 --- /dev/null +++ b/lib/mito/utilities/StdSharedPointer.h @@ -0,0 +1,72 @@ +// -*- c++ -*- +// +// Copyright (c) 2020-2024, the MiTo Authors, all rights reserved +// + + +// code guard +#pragma once + + +namespace mito::utilities { + + // declaration + template + class StdSharedPointer : public std::shared_ptr { + // types + public: + using shared_ptr_type = StdSharedPointer; + using resource_type = resourceT; + + // interface + public: + // returns the id of this (oriented) simplex + inline auto id() const -> index_t; + + // meta methods + public: + // constructor + template + requires(std::is_constructible_v) + inline StdSharedPointer(Args &&... args); + + // destructor + inline ~StdSharedPointer() = default; + + // copy constructor + inline StdSharedPointer(const shared_ptr_type &) = default; + + // move constructor + inline StdSharedPointer(shared_ptr_type &&) noexcept = default; + + // assignment operator + inline StdSharedPointer & operator=(const shared_ptr_type &) = default; + + // move assignment operator + inline StdSharedPointer & operator=(shared_ptr_type &&) noexcept = default; + }; + + template + inline bool operator==( + const StdSharedPointer & lhs, const StdSharedPointer & rhs) + { + return lhs.id() == rhs.id(); + } + + template + inline auto operator<=>( + const StdSharedPointer & lhs, const StdSharedPointer & rhs) + { + // delegate to ids + return lhs.id() <=> rhs.id(); + } +} + + +// get the inline definitions +#define mito_utilities_StdSharedPointer_icc +#include "StdSharedPointer.icc" +#undef mito_utilities_StdSharedPointer_icc + + +// end of file diff --git a/lib/mito/utilities/StdSharedPointer.icc b/lib/mito/utilities/StdSharedPointer.icc new file mode 100644 index 00000000..58037996 --- /dev/null +++ b/lib/mito/utilities/StdSharedPointer.icc @@ -0,0 +1,27 @@ +// -*- c++ -*- +// +// Copyright (c) 2020-2024, the MiTo Authors, all rights reserved +// + + +#if !defined(mito_utilities_StdSharedPointer_icc) +#error This header file contains implementation details of class mito::utilities::StdSharedPointer +#else + +template +auto +mito::utilities::StdSharedPointer::id() const -> index_t +{ + // the id is the (immutable) address of this object + return reinterpret_cast>(this->get()); +} + +template +template +requires(std::is_constructible_v) +mito::utilities::StdSharedPointer::StdSharedPointer(Args &&... args) : + std::shared_ptr(std::make_shared(args...)) +{} + +#endif +// end of file diff --git a/lib/mito/utilities/api.h b/lib/mito/utilities/api.h index 1fae0542..a4304a99 100644 --- a/lib/mito/utilities/api.h +++ b/lib/mito/utilities/api.h @@ -13,6 +13,10 @@ namespace mito::utilities { template using shared_ptr = SharedPointer; + // wrapper of std shared pointer alias + template + using std_shared_ptr = StdSharedPointer; + // index type alias template using index_t = std::uintptr_t; diff --git a/lib/mito/utilities/forward.h b/lib/mito/utilities/forward.h index adf1622d..28a199be 100644 --- a/lib/mito/utilities/forward.h +++ b/lib/mito/utilities/forward.h @@ -28,6 +28,10 @@ namespace mito::utilities { // requires reference_countable_c class SharedPointer; + // wrapper class for std shared pointer + template + class StdSharedPointer; + // class segmented allocator template class SegmentedAllocator; diff --git a/lib/mito/utilities/public.h b/lib/mito/utilities/public.h index 3cb21a03..af95c6dc 100644 --- a/lib/mito/utilities/public.h +++ b/lib/mito/utilities/public.h @@ -21,6 +21,7 @@ #include "Shareable.h" #include "SharedPointer.h" #include "Invalidatable.h" +#include "StdSharedPointer.h" #include "SegmentedAllocator.h" #include "SegmentedAllocatorIterator.h" #include "SegmentedVector.h" From 655108da1cc7b001ce39d88ee00e857a07be5e82 Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Sat, 8 Jun 2024 09:47:25 +0200 Subject: [PATCH 43/58] lib/mito/geometry: nodes are now wrapped into std shared pointers --- lib/mito/geometry/GeometricSimplex.h | 2 +- lib/mito/geometry/api.h | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/lib/mito/geometry/GeometricSimplex.h b/lib/mito/geometry/GeometricSimplex.h index cdc0598f..ebf48053 100644 --- a/lib/mito/geometry/GeometricSimplex.h +++ b/lib/mito/geometry/GeometricSimplex.h @@ -131,7 +131,7 @@ namespace mito::geometry { template requires(D > 0) - class GeometricSimplex<0, D> : public utilities::Invalidatable, public utilities::Shareable { + class GeometricSimplex<0, D> : public utilities::Invalidatable { public: // spatial dimension diff --git a/lib/mito/geometry/api.h b/lib/mito/geometry/api.h index 0f427d1a..69be0bdd 100644 --- a/lib/mito/geometry/api.h +++ b/lib/mito/geometry/api.h @@ -23,7 +23,7 @@ namespace mito::geometry { // node alias template - using node_t = utilities::shared_ptr>; + using node_t = utilities::std_shared_ptr>; // segment alias template From 1aabd15b39424da1d28a8568daa64e61b47613f8 Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Sat, 8 Jun 2024 12:28:31 +0200 Subject: [PATCH 44/58] lib/mito/mesh: remove unused {typedef} --- lib/mito/mesh/Boundary.h | 2 -- 1 file changed, 2 deletions(-) diff --git a/lib/mito/mesh/Boundary.h b/lib/mito/mesh/Boundary.h index b20bf957..b8b55b7d 100644 --- a/lib/mito/mesh/Boundary.h +++ b/lib/mito/mesh/Boundary.h @@ -16,8 +16,6 @@ namespace mito::mesh { static constexpr int D = cell_type::dim; // get the order of the cell static constexpr int N = cell_type::order; - // a collection of nodes - using nodes_type = cell_type::nodes_type; // get the family this cell type belongs to (e.g. geometric simplices) template using cell_family_type = typename cell_type::cell_family_type; From 7f7f4423f0cbea14648cca765ea6d8b3e596f595 Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Sat, 8 Jun 2024 15:43:09 +0200 Subject: [PATCH 45/58] lib/mito: the boundary of a mesh of segments returns now a collection of nodes --- lib/mito/geometry/utilities.h | 2 +- lib/mito/mesh/Boundary.h | 10 ++++++- lib/mito/mesh/Boundary.icc | 28 ++++++++++++++++++- lib/mito/mesh/api.h | 2 +- .../io/read_mesh_summit_segment_3D.cc | 22 +++++++-------- 5 files changed, 48 insertions(+), 16 deletions(-) diff --git a/lib/mito/geometry/utilities.h b/lib/mito/geometry/utilities.h index 0886d9fd..734e4f36 100644 --- a/lib/mito/geometry/utilities.h +++ b/lib/mito/geometry/utilities.h @@ -115,7 +115,7 @@ namespace mito::geometry { // builds the node based on vertex {vertex} with the vertex-point pairing as appears in {nodes} template auto geometric_simplex(const topology::simplex_t & simplex, const nodesT & nodes) - -> geometric_simplex_t + -> node_t requires(N == 0) { using vertex_type = topology::vertex_t; diff --git a/lib/mito/mesh/Boundary.h b/lib/mito/mesh/Boundary.h index b8b55b7d..f6cfda99 100644 --- a/lib/mito/mesh/Boundary.h +++ b/lib/mito/mesh/Boundary.h @@ -21,11 +21,19 @@ namespace mito::mesh { using cell_family_type = typename cell_type::cell_family_type; // the boundary mesh type using mesh_boundary_type = mesh_t>; + // the node type + using node_type = cell_type::node_type; + // a cloud of nodes (the boundary of a 1D mesh) + using node_cloud_type = std::unordered_set>; public: // returns the boundary mesh of {mesh} static inline auto boundary(const mesh_type & mesh) -> mesh_boundary_type - requires(N > 0); + requires(N > 1); + + // returns the boundary mesh of {mesh} + static inline auto boundary(const mesh_type & mesh) -> node_cloud_type + requires(N == 1); // returns the size of the boundary of {mesh} static inline auto boundary_size(const mesh_type & mesh) -> int diff --git a/lib/mito/mesh/Boundary.icc b/lib/mito/mesh/Boundary.icc index fde3fa15..56dc3939 100644 --- a/lib/mito/mesh/Boundary.icc +++ b/lib/mito/mesh/Boundary.icc @@ -34,7 +34,7 @@ requires(N > 0) template auto mito::mesh::Boundary::boundary(const mesh_type & mesh) -> mesh_boundary_type -requires(N > 0) +requires(N > 1) { // instantiate a new mesh for the boundary elements mesh_boundary_type boundary_mesh; @@ -57,6 +57,32 @@ requires(N > 0) return boundary_mesh; } +template +auto +mito::mesh::Boundary::boundary(const mesh_type & mesh) -> node_cloud_type +requires(N == 1) +{ + // instantiate a new node cloud for the boundary nodes + node_cloud_type node_cloud; + + // loop on the mesh cells + for (const auto & cell : mesh.cells()) { + // create view of {subcell}s on the boundary of the mesh + auto subcellsOnBoundary = + cell.simplex()->composition() | std::views::filter([&mesh](const auto & subcell) { + return mesh.isOnBoundary(subcell); + }); + // loop on the {subcell} on the boundary + for (const auto & subcell : subcellsOnBoundary) { + // add {cell} to the boundary mesh + node_cloud.insert(mito::geometry::geometric_simplex(subcell, cell.nodes())); + } + } + + // return the cloud of boundary nodes + return node_cloud; +} + #endif // end of file diff --git a/lib/mito/mesh/api.h b/lib/mito/mesh/api.h index 269ec5ad..fce65f65 100644 --- a/lib/mito/mesh/api.h +++ b/lib/mito/mesh/api.h @@ -15,7 +15,7 @@ namespace mito::mesh { // assemble boundary mesh of {mesh} template class cellT> - auto boundary(const mesh_t> & mesh) -> mesh_t> + auto boundary(const mesh_t> & mesh) { return Boundary>>::boundary(mesh); } diff --git a/tests/mito.lib/io/read_mesh_summit_segment_3D.cc b/tests/mito.lib/io/read_mesh_summit_segment_3D.cc index 6ae9136a..f1121c0d 100644 --- a/tests/mito.lib/io/read_mesh_summit_segment_3D.cc +++ b/tests/mito.lib/io/read_mesh_summit_segment_3D.cc @@ -23,12 +23,11 @@ TEST(SummitReader, LoadSummitSegmentsMesh3D) // assert you read 10 cells EXPECT_EQ(mesh.nCells(), 10); - // TOFIX - // // assert you found 2 nodes on the boundary - // { - // auto boundary_mesh = mito::mesh::boundary(mesh); - // EXPECT_EQ(boundary_mesh.nCells(), 2); - // } + // assert you found 2 nodes on the boundary + { + auto boundary_nodes = mito::mesh::boundary(mesh); + EXPECT_EQ(std::size(boundary_nodes), 2); + } // get a cell auto & cell = mesh.cells()[5]; @@ -39,10 +38,9 @@ TEST(SummitReader, LoadSummitSegmentsMesh3D) // assert you read 9 cells EXPECT_EQ(mesh.nCells(), 9); - // TOFIX - // // assert you found 4 nodes on the boundary - // { - // auto boundary_mesh = mito::mesh::boundary(mesh); - // EXPECT_EQ(boundary_mesh.nCells(), 4); - // } + // assert you found 4 nodes on the boundary + { + auto boundary_nodes = mito::mesh::boundary(mesh); + EXPECT_EQ(std::size(boundary_nodes), 4); + } } \ No newline at end of file From cdf291ffb1452f38931f84e19d1ebf8ac78b340f Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Sat, 8 Jun 2024 17:04:24 +0200 Subject: [PATCH 46/58] lib/mito/geometry: remove unused type traits for 0-order {GeometricSimplex} --- lib/mito/geometry/GeometricSimplex.h | 16 ---------------- 1 file changed, 16 deletions(-) diff --git a/lib/mito/geometry/GeometricSimplex.h b/lib/mito/geometry/GeometricSimplex.h index ebf48053..7bb57d33 100644 --- a/lib/mito/geometry/GeometricSimplex.h +++ b/lib/mito/geometry/GeometricSimplex.h @@ -143,25 +143,9 @@ namespace mito::geometry { // the type of the vertices using vertex_type = topology::vertex_t; - // the type of the topological simplex - using simplex_type = topology::oriented_simplex_t<0>; - - // number of vertices of simplex - static constexpr int n_vertices = 0; - - // the node type - using node_type = node_t; - - // a collection of nodes - using nodes_type = std::array; - // the point type using point_type = point_t; - // typedef for the topological family type (simplicial) - template - using cell_topological_family_type = typename topology::cell_family; - public: // get the coordinates of the point GeometricSimplex(const vertex_type & vertex, const point_type & point) : From 091e28a94359ff7cc8c1ebb8cf18e375fd319760 Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Sat, 8 Jun 2024 17:05:21 +0200 Subject: [PATCH 47/58] lib/mito/geometry: add concept of a geometric simplex and of a geometric segment --- lib/mito/geometry/forward.h | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/lib/mito/geometry/forward.h b/lib/mito/geometry/forward.h index aed99cc6..51ae759d 100644 --- a/lib/mito/geometry/forward.h +++ b/lib/mito/geometry/forward.h @@ -34,6 +34,18 @@ namespace mito::geometry { requires((N >= 0) && (N <= D) && (D > 0)) class GeometricSimplex; + // concept of a geometric simplex + template + concept geometric_simplex_c = requires(F c) { + // require that F only binds to {GeometricSimplex} specializations + [](const GeometricSimplex &) { + }(c); + }; + + // concept of a geometric segment (geometric simplex with order 1) + template + concept geometric_segment_c = geometric_simplex_c && F::order == 1; + // class point cloud template class PointCloud; From 3297c01c9ed180e401cc54d0d55c84f7d82cb860 Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Sat, 8 Jun 2024 17:07:26 +0200 Subject: [PATCH 48/58] lib/mito/mesh: readability improvement in class {Boundary} Instead of implementing two almost identical methods for fetching the boundary of a mesh of segments and of any other mesh, devise the suitable type of the boundary (either a cloud of nodes or a proper boundary mesh) at compile time. --- lib/mito/mesh/Boundary.h | 54 ++++++++++++++++++++++++-------------- lib/mito/mesh/Boundary.icc | 35 +++--------------------- 2 files changed, 39 insertions(+), 50 deletions(-) diff --git a/lib/mito/mesh/Boundary.h b/lib/mito/mesh/Boundary.h index f6cfda99..c613ea5a 100644 --- a/lib/mito/mesh/Boundary.h +++ b/lib/mito/mesh/Boundary.h @@ -4,14 +4,15 @@ namespace mito::mesh { - template - class Boundary { + // type for a boundary mesh of a mesh of cells {cellT} + template + struct boundary_mesh; - private: + // type for a boundary mesh of a mesh of cells {cellT} + template + struct boundary_mesh { // typedef my template parameter - using mesh_type = meshT; - // typedef cell type - using cell_type = typename mesh_type::cell_type; + using cell_type = cellT; // the dimension of physical space static constexpr int D = cell_type::dim; // get the order of the cell @@ -19,25 +20,40 @@ namespace mito::mesh { // get the family this cell type belongs to (e.g. geometric simplices) template using cell_family_type = typename cell_type::cell_family_type; - // the boundary mesh type - using mesh_boundary_type = mesh_t>; + // the boundary type: a boundary mesh + using type = mesh_t>; + }; + + // specialization for geometric segments + template + struct boundary_mesh { + // typedef my template parameter + using cell_type = cellT; // the node type using node_type = cell_type::node_type; - // a cloud of nodes (the boundary of a 1D mesh) - using node_cloud_type = std::unordered_set>; + // the boundary type: a cloud of nodes + using type = std::unordered_set>; + }; - public: - // returns the boundary mesh of {mesh} - static inline auto boundary(const mesh_type & mesh) -> mesh_boundary_type - requires(N > 1); + template + class Boundary { + + private: + // typedef my template parameter + using mesh_type = meshT; + // the cell type + using cell_type = typename mesh_type::cell_type; + // the dimension of physical space + static constexpr int D = cell_type::dim; + // boundary type: either a mesh or a collection of nodes (boundary of a mesh of segments) + using boundary_mesh_type = typename boundary_mesh::type; - // returns the boundary mesh of {mesh} - static inline auto boundary(const mesh_type & mesh) -> node_cloud_type - requires(N == 1); + public: + // returns the boundary of {mesh} + static inline auto boundary(const mesh_type & mesh) -> boundary_mesh_type; // returns the size of the boundary of {mesh} - static inline auto boundary_size(const mesh_type & mesh) -> int - requires(N > 0); + static inline auto boundary_size(const mesh_type & mesh) -> int; }; } diff --git a/lib/mito/mesh/Boundary.icc b/lib/mito/mesh/Boundary.icc index 56dc3939..2ac1495d 100644 --- a/lib/mito/mesh/Boundary.icc +++ b/lib/mito/mesh/Boundary.icc @@ -12,7 +12,6 @@ template auto mito::mesh::Boundary::boundary_size(const mesh_type & mesh) -> int -requires(N > 0) { // number of boundary cells int count = 0; @@ -33,11 +32,10 @@ requires(N > 0) template auto -mito::mesh::Boundary::boundary(const mesh_type & mesh) -> mesh_boundary_type -requires(N > 1) +mito::mesh::Boundary::boundary(const mesh_type & mesh) -> boundary_mesh_type { // instantiate a new mesh for the boundary elements - mesh_boundary_type boundary_mesh; + boundary_mesh_type boundary; // loop on the mesh cells for (const auto & cell : mesh.cells()) { @@ -49,39 +47,14 @@ requires(N > 1) // loop on the {subcell} on the boundary for (const auto & subcell : subcellsOnBoundary) { // add {cell} to the boundary mesh - boundary_mesh.insert(mito::geometry::geometric_simplex(subcell, cell.nodes())); + boundary.insert(mito::geometry::geometric_simplex(subcell, cell.nodes())); } } // return the boundary mesh - return boundary_mesh; + return boundary; } -template -auto -mito::mesh::Boundary::boundary(const mesh_type & mesh) -> node_cloud_type -requires(N == 1) -{ - // instantiate a new node cloud for the boundary nodes - node_cloud_type node_cloud; - - // loop on the mesh cells - for (const auto & cell : mesh.cells()) { - // create view of {subcell}s on the boundary of the mesh - auto subcellsOnBoundary = - cell.simplex()->composition() | std::views::filter([&mesh](const auto & subcell) { - return mesh.isOnBoundary(subcell); - }); - // loop on the {subcell} on the boundary - for (const auto & subcell : subcellsOnBoundary) { - // add {cell} to the boundary mesh - node_cloud.insert(mito::geometry::geometric_simplex(subcell, cell.nodes())); - } - } - - // return the cloud of boundary nodes - return node_cloud; -} #endif From 089a4da24cb7960ee7060dc183d01ae4adf46026 Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Sat, 8 Jun 2024 17:11:35 +0200 Subject: [PATCH 49/58] lib/mito: readability improvement Remove {requires} from {concept} definition when possible. --- lib/mito/fields/forward.h | 3 +-- lib/mito/tensor/api.h | 7 +++---- lib/mito/utilities/api.h | 2 +- 3 files changed, 5 insertions(+), 7 deletions(-) diff --git a/lib/mito/fields/forward.h b/lib/mito/fields/forward.h index 2e097e54..f8107fe6 100644 --- a/lib/mito/fields/forward.h +++ b/lib/mito/fields/forward.h @@ -62,9 +62,8 @@ namespace mito::fields { // concept of two fields being compatible with each other (i.e. defined on the same coordinates) template // a scalar field is a field returning a scalar - concept compatible_fields_c = requires { + concept compatible_fields_c = std::is_same_v; - }; } diff --git a/lib/mito/tensor/api.h b/lib/mito/tensor/api.h index 28a709dd..cbb64f61 100644 --- a/lib/mito/tensor/api.h +++ b/lib/mito/tensor/api.h @@ -57,16 +57,15 @@ namespace mito { // concept of a matrix template - concept matrix_c = tensor_c and requires { F::rank == 2; }; + concept matrix_c = tensor_c and F::rank == 2; // concept of a diagonal matrix template - concept diagonal_matrix_c = matrix_c and requires { F::diagonal; }; + concept diagonal_matrix_c = matrix_c and F::diagonal; // concept of a symmetric matrix template - concept symmetric_matrix_c = - (matrix_c and requires { F::symmetric; }) or diagonal_matrix_c; + concept symmetric_matrix_c = (matrix_c and F::symmetric) or diagonal_matrix_c; // I-th basis vector in dimension D template diff --git a/lib/mito/utilities/api.h b/lib/mito/utilities/api.h index a4304a99..4002a598 100644 --- a/lib/mito/utilities/api.h +++ b/lib/mito/utilities/api.h @@ -35,7 +35,7 @@ namespace mito::utilities { // concept of the types having the same dimension template - concept same_dim_c = requires { F1::dim == F2::dim; }; + concept same_dim_c = F1::dim == F2::dim; } From c0bb3f86f18298a0c6fce0654553a8f066eda91a Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Sun, 9 Jun 2024 06:40:39 +0200 Subject: [PATCH 50/58] lib/mito/io/summit: in mesh reader use templates to avoid repeating the same code for segments, triangles and tets --- lib/mito/io/summit/reader.h | 189 ++++++++---------------------------- 1 file changed, 39 insertions(+), 150 deletions(-) diff --git a/lib/mito/io/summit/reader.h b/lib/mito/io/summit/reader.h index bfdbea4e..52f425c3 100644 --- a/lib/mito/io/summit/reader.h +++ b/lib/mito/io/summit/reader.h @@ -37,94 +37,37 @@ namespace mito::io::summit { return; } - template - auto readSegment( - std::ifstream & fileStream, mesh::mesh_t> & mesh, - const std::vector> & nodes) -> void + template + auto readElement( + std::ifstream & fileStream, meshT & mesh, + const std::vector> & nodes) -> void { - int index0 = 0; - fileStream >> index0; - --index0; - - int index1 = 0; - fileStream >> index1; - --index1; - - const auto & node_0 = nodes[index0]; - const auto & node_1 = nodes[index1]; - - mesh.insert({ node_0, node_1 }); - - // QUESTION: Can the label be more than one? - // read label for cell - // TOFIX: Ignored for now - std::string cell_label; - fileStream >> cell_label; - - // all done - return; - } - - template - auto readTriangle( - std::ifstream & fileStream, mesh::mesh_t> & mesh, - const std::vector> & nodes) -> void - { - int index0 = 0; - fileStream >> index0; - --index0; - - int index1 = 0; - fileStream >> index1; - --index1; - - int index2 = 0; - fileStream >> index2; - --index2; - - const auto & node_0 = nodes[index0]; - const auto & node_1 = nodes[index1]; - const auto & node_2 = nodes[index2]; - - mesh.insert({ node_0, node_1, node_2 }); - - // QUESTION: Can the label be more than one? - // read label for cell - // TOFIX: Ignored for now - std::string cell_label; - fileStream >> cell_label; - - // all done - return; - } - - template - auto readTetrahedron( - std::ifstream & fileStream, mesh::mesh_t> & mesh, - const std::vector> & nodes) -> void - { - int index0 = 0; - fileStream >> index0; - --index0; - - int index1 = 0; - fileStream >> index1; - --index1; - - int index2 = 0; - fileStream >> index2; - --index2; - - int index3 = 0; - fileStream >> index3; - --index3; + // get the number of vertices + constexpr int N = meshT::cell_type::n_vertices; + + // the element connectivity + std::array index; + // for each node + for (int i = 0; i < N; ++i) { + // read from file the id of the node + int id = 0; + fileStream >> id; + // start from zero (the summit format starts counting from one) + --id; + // take a note of it + index[i] = id; + } - const auto & node_0 = nodes[index0]; - const auto & node_1 = nodes[index1]; - const auto & node_2 = nodes[index2]; - const auto & node_3 = nodes[index3]; + // helper function to call the mesh insert method with {N} nodes + constexpr auto _insert = []( + meshT & mesh, const std::array index, + const std::vector> & nodes, + std::index_sequence) { + mesh.insert({ nodes[index[I]]... }); + }; - mesh.insert({ node_0, node_1, node_2, node_3 }); + // insert the nodes in the mesh + _insert(mesh, index, nodes, std::make_index_sequence{}); // QUESTION: Can the label be more than one? // read label for cell @@ -136,53 +79,23 @@ namespace mito::io::summit { return; } - template - auto readSegments( - std::ifstream & fileStream, mesh::mesh_t> & mesh, int N_cells, - const std::vector> & nodes) -> void - { - for (int i = 0; i < N_cells; ++i) { - int cell_type = 0; - fileStream >> cell_type; - - if (cell_type == 2) { - readSegment(fileStream, mesh, nodes); - } - } - - // all done - return; - } - - template - auto readTriangles( - std::ifstream & fileStream, mesh::mesh_t> & mesh, int N_cells, - const std::vector> & nodes) -> void + template + auto readElements( + std::ifstream & fileStream, mesh::mesh_t & mesh, int N_cells, + const std::vector> & nodes) -> void { - for (int i = 0; i < N_cells; ++i) { - int cell_type = 0; - fileStream >> cell_type; + // get the number of vertices + constexpr int N = cellT::n_vertices; - if (cell_type == 3) { - readTriangle(fileStream, mesh, nodes); - } - } - - // all done - return; - } - - template - auto readTetrahedra( - std::ifstream & fileStream, mesh::mesh_t> & mesh, int N_cells, - const std::vector> & nodes) -> void - { + // for each element for (int i = 0; i < N_cells; ++i) { + // read the cell type from file int cell_type = 0; fileStream >> cell_type; - if (cell_type == 4) { - readTetrahedron(fileStream, mesh, nodes); + // TOFIX: should be cell_type == summit::cell::type + if (cell_type == N) { + readElement(fileStream, mesh, nodes); } } @@ -190,30 +103,6 @@ namespace mito::io::summit { return; } - template - auto readElements( - std::ifstream & fileStream, mesh::mesh_t> & mesh, int N_cells, - const std::vector> & nodes) -> void - { - return readSegments(fileStream, mesh, N_cells, nodes); - } - - template - auto readElements( - std::ifstream & fileStream, mesh::mesh_t> & mesh, int N_cells, - const std::vector> & nodes) -> void - { - return readTetrahedra(fileStream, mesh, N_cells, nodes); - } - - template - auto readElements( - std::ifstream & fileStream, mesh::mesh_t> & mesh, int N_cells, - const std::vector> & nodes) -> void - { - return readTriangles(fileStream, mesh, N_cells, nodes); - } - template auto reader( std::ifstream & fileStream, geometry::coordinate_system_t & coordinate_system) From 8fbcd524dfe360aa3c25634db07f74eabba91d49 Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Sun, 9 Jun 2024 06:45:57 +0200 Subject: [PATCH 51/58] lib/mito/geometry: add concept for all geometric simplices --- lib/mito/geometry/forward.h | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/lib/mito/geometry/forward.h b/lib/mito/geometry/forward.h index 51ae759d..be6ef700 100644 --- a/lib/mito/geometry/forward.h +++ b/lib/mito/geometry/forward.h @@ -46,6 +46,14 @@ namespace mito::geometry { template concept geometric_segment_c = geometric_simplex_c && F::order == 1; + // concept of a geometric triangle (geometric simplex with order 2) + template + concept geometric_triangle_c = geometric_simplex_c && F::order == 2; + + // concept of a geometric tetrahedron (geometric simplex with order 3) + template + concept geometric_tetrahedron_c = geometric_simplex_c && F::order == 3; + // class point cloud template class PointCloud; From 5c771a9ec4b2a0903ec3c567aad5279d31874c5b Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Sun, 9 Jun 2024 06:53:28 +0200 Subject: [PATCH 52/58] lib/mito/io/summit: use type traits on the mesh cells This will make it easier to generalize the reader / writer to the case of nonsimplicial elements. --- lib/mito/io/summit/MeshWriterSummit.h | 4 ++- lib/mito/io/summit/public.h | 1 + lib/mito/io/summit/reader.h | 20 +++++++-------- lib/mito/io/summit/summit_cell.h | 37 +++++++++++++++++++++++++++ 4 files changed, 50 insertions(+), 12 deletions(-) create mode 100644 lib/mito/io/summit/summit_cell.h diff --git a/lib/mito/io/summit/MeshWriterSummit.h b/lib/mito/io/summit/MeshWriterSummit.h index 90e5fa8a..fe16d69e 100644 --- a/lib/mito/io/summit/MeshWriterSummit.h +++ b/lib/mito/io/summit/MeshWriterSummit.h @@ -16,6 +16,8 @@ namespace mito::io { private: // the mesh type using mesh_type = meshT; + // the cell type + using cell_type = typename mesh_type::cell_type; // the coordinate system type using coord_system_type = coordSystemT; // the dimension of the physical space @@ -63,7 +65,7 @@ namespace mito::io { // write the cells to file for (const auto & cell : _mesh.cells()) { - outfile << mesh_type::cell_type::n_vertices << " "; + outfile << summit::cell::type << " "; for (const auto & node : cell.nodes()) { outfile << std::distance(this->_points.begin(), this->_points.find(node->point())) diff --git a/lib/mito/io/summit/public.h b/lib/mito/io/summit/public.h index 39949f4b..93297af2 100644 --- a/lib/mito/io/summit/public.h +++ b/lib/mito/io/summit/public.h @@ -11,6 +11,7 @@ #include "externals.h" // classes implementation +#include "summit_cell.h" #include "MeshWriterSummit.h" // published type factories; this is the file you are looking for... diff --git a/lib/mito/io/summit/reader.h b/lib/mito/io/summit/reader.h index 52f425c3..14bb400f 100644 --- a/lib/mito/io/summit/reader.h +++ b/lib/mito/io/summit/reader.h @@ -37,13 +37,13 @@ namespace mito::io::summit { return; } - template + template auto readElement( - std::ifstream & fileStream, meshT & mesh, - const std::vector> & nodes) -> void + std::ifstream & fileStream, mesh::mesh_t & mesh, + const std::vector> & nodes) -> void { // get the number of vertices - constexpr int N = meshT::cell_type::n_vertices; + constexpr int N = cellT::n_vertices; // the element connectivity std::array index; @@ -60,8 +60,8 @@ namespace mito::io::summit { // helper function to call the mesh insert method with {N} nodes constexpr auto _insert = []( - meshT & mesh, const std::array index, - const std::vector> & nodes, + mesh::mesh_t & mesh, const std::array index, + const std::vector> & nodes, std::index_sequence) { mesh.insert({ nodes[index[I]]... }); }; @@ -84,17 +84,15 @@ namespace mito::io::summit { std::ifstream & fileStream, mesh::mesh_t & mesh, int N_cells, const std::vector> & nodes) -> void { - // get the number of vertices - constexpr int N = cellT::n_vertices; - // for each element for (int i = 0; i < N_cells; ++i) { // read the cell type from file int cell_type = 0; fileStream >> cell_type; - // TOFIX: should be cell_type == summit::cell::type - if (cell_type == N) { + // if the cell type read from file matches with the cell of the mesh to be populated + if (cell_type == summit::cell::type) { + // read the element and insert it in the mesh readElement(fileStream, mesh, nodes); } } diff --git a/lib/mito/io/summit/summit_cell.h b/lib/mito/io/summit/summit_cell.h new file mode 100644 index 00000000..68c189d8 --- /dev/null +++ b/lib/mito/io/summit/summit_cell.h @@ -0,0 +1,37 @@ +// -*- c++ -*- +// +// Copyright (c) 2020-2024, the MiTo Authors, all rights reserved +// + +// code guard +#pragma once + + +namespace mito::io::summit { + + // the summit type of a geometric simplex + template + struct cell; + + // specialization for a geometric segment + template + struct cell { + static constexpr int type = 2; + }; + + // specialization for a geometric triangle + template + struct cell { + static constexpr int type = 3; + }; + + // specialization for a geometric tetrahedron + template + struct cell { + static constexpr int type = 4; + }; + +} // namespace io::summit + + +// end of file From fb0b55c6df2e8d4dd191b1390c9ca753e4e2517e Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Sun, 9 Jun 2024 17:06:17 +0200 Subject: [PATCH 53/58] lib/mito/geometry: publish {point_type} trait in class {CoordinateSystem} --- lib/mito/geometry/CoordinateSystem.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/lib/mito/geometry/CoordinateSystem.h b/lib/mito/geometry/CoordinateSystem.h index feb0a661..61e1b694 100644 --- a/lib/mito/geometry/CoordinateSystem.h +++ b/lib/mito/geometry/CoordinateSystem.h @@ -16,12 +16,12 @@ namespace mito::geometry { using coordinates_type = coordT; // the dimension of the physical space static constexpr int dim = coordinates_type::dim; + // a point + using point_type = point_t; private: // alias {dim} for internal purposes static constexpr int D = dim; - // a point - using point_type = point_t; // id type of point using point_id_type = utilities::index_t; // a map between points and their coordinates From 7366d74adf09b00f1385d045e36fabb1e006e611 Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Sun, 9 Jun 2024 17:45:48 +0200 Subject: [PATCH 54/58] lib/mito/io: redesign of summit and vtk writers The new design is informed of the difference between points and nodes. This prepares the ground for supporting both CG and DG meshes. --- lib/mito/io/GridWriter.h | 37 ------------------- lib/mito/io/Writer.h | 30 ++++++++++++++++ lib/mito/io/public.h | 2 +- lib/mito/io/summit/MeshWriterSummit.h | 29 +++++++++------ lib/mito/io/vtk/GridWriterVTK.h | 33 +++-------------- lib/mito/io/vtk/MeshWriterVTK.h | 52 +++++++++++++++++++++------ lib/mito/io/vtk/PointCloudWriterVTK.h | 10 +++++- 7 files changed, 104 insertions(+), 89 deletions(-) delete mode 100644 lib/mito/io/GridWriter.h create mode 100644 lib/mito/io/Writer.h diff --git a/lib/mito/io/GridWriter.h b/lib/mito/io/GridWriter.h deleted file mode 100644 index 4916c3ce..00000000 --- a/lib/mito/io/GridWriter.h +++ /dev/null @@ -1,37 +0,0 @@ -// -*- c++ -*- -// -// Copyright (c) 2020-2024, the MiTo Authors, all rights reserved -// - -// code guard -#pragma once - - -namespace mito::io { - - template - class GridWriter { - - private: - // the type of point - using point_type = geometry::point_t; - // the type of a collection of points - using points_type = std::unordered_set>; - - public: - GridWriter(std::string filename) : _filename(filename), _points() {} - - virtual auto write() const -> void = 0; - - protected: - // the name of the output file to be written - std::string _filename; - - // a collection of points in the grid - points_type _points; - }; - -} // namespace mito::io - - -// end of file diff --git a/lib/mito/io/Writer.h b/lib/mito/io/Writer.h new file mode 100644 index 00000000..632ed5ba --- /dev/null +++ b/lib/mito/io/Writer.h @@ -0,0 +1,30 @@ +// -*- c++ -*- +// +// Copyright (c) 2020-2024, the MiTo Authors, all rights reserved +// + +// code guard +#pragma once + + +namespace mito::io { + + // a class that writes to file + class Writer { + + public: + // constructor + Writer(std::string filename) : _filename(filename) {} + + // the {write} method + virtual auto write() const -> void = 0; + + protected: + // the name of the output file to be written + std::string _filename; + }; + +} // namespace mito::io + + +// end of file diff --git a/lib/mito/io/public.h b/lib/mito/io/public.h index 09f06291..63fd4234 100644 --- a/lib/mito/io/public.h +++ b/lib/mito/io/public.h @@ -11,7 +11,7 @@ #include "externals.h" // classes implementation -#include "GridWriter.h" +#include "Writer.h" // classes implementation #include "summit/public.h" diff --git a/lib/mito/io/summit/MeshWriterSummit.h b/lib/mito/io/summit/MeshWriterSummit.h index fe16d69e..78f6cbb4 100644 --- a/lib/mito/io/summit/MeshWriterSummit.h +++ b/lib/mito/io/summit/MeshWriterSummit.h @@ -9,9 +9,10 @@ namespace mito::io { + // class that writes a mesh to file in {summit} format template requires(utilities::same_dim_c) - class MeshWriterSummit : public GridWriter { + class MeshWriterSummit : public Writer { private: // the mesh type @@ -22,19 +23,24 @@ namespace mito::io { using coord_system_type = coordSystemT; // the dimension of the physical space static constexpr int D = mesh_type::dim; + // the type of point + using point_type = typename coord_system_type::point_type; + // the type of a collection of points + using points_type = std::unordered_set>; public: // constructor MeshWriterSummit( std::string filename, const mesh_type & mesh, const coord_system_type & coord_system) : - GridWriter(filename), + Writer(filename), _mesh(mesh), - _coord_system(coord_system) + _coord_system(coord_system), + _points() { - // insert the points corresponding to the mesh nodes + // insert the points corresponding to the mesh nodes (eliminating duplicates) for (const auto & cell : mesh.cells()) { for (const auto & node : cell.nodes()) { - this->_points.insert(node->point()); + _points.insert(node->point()); } } } @@ -55,10 +61,10 @@ namespace mito::io { // populate the file heading // TOFIX: number of materials is always 1 for now outfile << D << std::endl; - outfile << std::size(this->_points) << " " << _mesh.nCells() << " " << 1 << std::endl; + outfile << std::size(_points) << " " << _mesh.nCells() << " " << 1 << std::endl; // write the points to file - for (const auto & point : this->_points) { + for (const auto & point : _points) { const auto & coord = _coord_system.coordinates(point); outfile << std::setprecision(15) << coord << std::endl; } @@ -67,10 +73,8 @@ namespace mito::io { for (const auto & cell : _mesh.cells()) { outfile << summit::cell::type << " "; for (const auto & node : cell.nodes()) { - outfile - << std::distance(this->_points.begin(), this->_points.find(node->point())) - + 1 - << " "; + outfile << std::distance(_points.begin(), _points.find(node->point())) + 1 + << " "; } // TOFIX: material label is always 1 for now outfile << 1 << std::endl; @@ -89,6 +93,9 @@ namespace mito::io { // a const reference to the coordinate system const coord_system_type & _coord_system; + + // a collection of points in the mesh + points_type _points; }; } // namespace mito::io diff --git a/lib/mito/io/vtk/GridWriterVTK.h b/lib/mito/io/vtk/GridWriterVTK.h index 5168de7d..3d493b84 100644 --- a/lib/mito/io/vtk/GridWriterVTK.h +++ b/lib/mito/io/vtk/GridWriterVTK.h @@ -10,14 +10,16 @@ namespace mito::io::vtk { template - class GridWriterVTK : public GridWriter { + class GridWriterVTK : public Writer { private: // the type of grid using grid_type = vtkSmartPointer; protected: - GridWriterVTK(std::string filename) : GridWriter(filename), _grid(grid_type::New()) {} + // constructor + // (protected so this class cannot be instantiated unless by the derived classes) + GridWriterVTK(std::string filename) : Writer(filename), _grid(grid_type::New()) {} public: auto write() const -> void override @@ -40,33 +42,6 @@ namespace mito::io::vtk { return; } - template - auto attach_field(const fem::nodal_field_t & field, std::string fieldname) -> void - { - // get the number of nodes - auto n_nodes = field.n_nodes(); - - // check the number of nodes in the field equals the number of points in the grid - assert(n_nodes == _grid->GetNumberOfPoints()); - - // initialize a vtk array - auto vtkArray = vtkSmartPointer::New(); - vtkArray->SetName(fieldname.data()); - vtkArray->SetNumberOfComponents(Y::size); - vtkArray->SetNumberOfTuples(n_nodes); - - // populate the array with the nodal values - for (auto & [node, value] : field) { - // get the index corresponding to the current point - auto index = - std::distance(this->_points.begin(), this->_points.find(node->point())); - vtkArray->SetTuple(index, value.begin()); - } - - // insert array into output mesh - _grid->GetPointData()->AddArray(vtkArray); - } - protected: // the grid grid_type _grid; diff --git a/lib/mito/io/vtk/MeshWriterVTK.h b/lib/mito/io/vtk/MeshWriterVTK.h index 80574fe4..371dc901 100644 --- a/lib/mito/io/vtk/MeshWriterVTK.h +++ b/lib/mito/io/vtk/MeshWriterVTK.h @@ -20,6 +20,10 @@ namespace mito::io::vtk { using coord_system_type = coordSystemT; // the dimension of the physical space static constexpr int D = mesh_type::dim; + // the type of node + using node_type = typename mesh_type::cell_type::node_type; + // the type of a collection of nodes + using nodes_type = std::unordered_set>; private: auto _create_vtk_grid(const mesh_type & mesh, const coord_system_type & coordinate_system) @@ -28,8 +32,8 @@ namespace mito::io::vtk { for (const auto & cell : mesh.cells()) { // loop over the nodes of the cell for (const auto & node : cell.nodes()) { - // add the point to the collection of points (eliminating duplicates) - this->_points.insert(node->point()); + // add the node to the collection of nodes + _nodes.insert(node); } } @@ -37,8 +41,8 @@ namespace mito::io::vtk { auto pointsVtk = vtkSmartPointer::New(); // insert the new vtk point - for (const auto & point : this->_points) { - insert_vtk_point(coordinate_system.coordinates(point), pointsVtk); + for (const auto & node : _nodes) { + insert_vtk_point(coordinate_system.coordinates(node->point()), pointsVtk); } // loop over the cells @@ -52,12 +56,10 @@ namespace mito::io::vtk { // loop over the nodes of the cell for (const auto & node : cell.nodes()) { - // retrieve the corresponding point - const auto & point = node->point(); - // assert that the point is present in the set of points - assert(this->_points.contains(point)); - // calculate the index of the point - auto index = std::distance(this->_points.begin(), this->_points.find(point)); + // assert that the node is present in the collection of nodes + assert(_nodes.contains(node)); + // calculate the index of the node + auto index = std::distance(_nodes.begin(), _nodes.find(node)); // set the id of the point cellVtk->GetPointIds()->SetId(indexLocalPointVtk, index); // update local index for the points in the cell @@ -82,6 +84,36 @@ namespace mito::io::vtk { { _create_vtk_grid(mesh, coord_system); } + + template + auto attach_field(const fem::nodal_field_t & field, std::string fieldname) -> void + { + // get the number of nodes + auto n_nodes = field.n_nodes(); + + // check the number of nodes in the field equals the number of points in the grid + assert(n_nodes == this->_grid->GetNumberOfPoints()); + + // initialize a vtk array + auto vtkArray = vtkSmartPointer::New(); + vtkArray->SetName(fieldname.data()); + vtkArray->SetNumberOfComponents(Y::size); + vtkArray->SetNumberOfTuples(n_nodes); + + // populate the array with the nodal values + for (auto & [node, value] : field) { + // get the index corresponding to the current node + auto index = std::distance(_nodes.begin(), _nodes.find(node)); + vtkArray->SetTuple(index, value.begin()); + } + + // insert array into output mesh + this->_grid->GetPointData()->AddArray(vtkArray); + } + + private: + // a collection of nodes in the grid + nodes_type _nodes; }; } // namespace mito::io::vtk diff --git a/lib/mito/io/vtk/PointCloudWriterVTK.h b/lib/mito/io/vtk/PointCloudWriterVTK.h index 813c5ccc..91737e3a 100644 --- a/lib/mito/io/vtk/PointCloudWriterVTK.h +++ b/lib/mito/io/vtk/PointCloudWriterVTK.h @@ -20,6 +20,10 @@ namespace mito::io::vtk { using coord_system_type = coordSystemT; // the dimension of the physical space static constexpr int D = cloudT::dim; + // the type of point + using point_type = typename coord_system_type::point_type; + // the type of a collection of points + using points_type = std::unordered_set>; private: auto _create_vtk_grid(const cloud_type & cloud, const coord_system_type & coordinate_system) @@ -32,7 +36,7 @@ namespace mito::io::vtk { // insert the new vtk point for (const auto & point : points) { - this->_points.insert(point); + _points.insert(point); insert_vtk_point(coordinate_system.coordinates(point), pointsVtk); } @@ -51,6 +55,10 @@ namespace mito::io::vtk { { _create_vtk_grid(cloud, coord_system); } + + private: + // a collection of points in the mesh + points_type _points; }; } // namespace mito::io::vtk From 76ae72e34ca1c3ec601b068ec9febaf4633c678b Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Sun, 9 Jun 2024 17:48:22 +0200 Subject: [PATCH 55/58] tests/mito.lib/io: add check that the number of nodes is conserved in the original and reread mesh --- tests/mito.lib/io/summit_to_summit_mesh_2D.cc | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/tests/mito.lib/io/summit_to_summit_mesh_2D.cc b/tests/mito.lib/io/summit_to_summit_mesh_2D.cc index 2cc7c5ae..b21c50f1 100644 --- a/tests/mito.lib/io/summit_to_summit_mesh_2D.cc +++ b/tests/mito.lib/io/summit_to_summit_mesh_2D.cc @@ -16,6 +16,8 @@ TEST(SummitToSummit, Mesh2D) { int original_mesh_cells = 0; int reread_mesh_cells = 0; + int original_mesh_nodes = 0; + int reread_mesh_nodes = 0; { // the coordinate system @@ -29,6 +31,9 @@ TEST(SummitToSummit, Mesh2D) // get the original number of mesh cells original_mesh_cells = mesh.nCells(); + // get the original number of mesh nodes by counting the nodes of a nodal field built on it + original_mesh_nodes = mito::fem::nodal_field(mesh, "field").n_nodes(); + // write summit mesh mito::io::summit::writer("rectangle_copy", mesh, coord_system); } @@ -44,6 +49,9 @@ TEST(SummitToSummit, Mesh2D) // get the reread number of mesh cells reread_mesh_cells = mesh.nCells(); + // get the reread number of mesh nodes by counting the nodes of a nodal field built on it + reread_mesh_nodes = mito::fem::nodal_field(mesh, "field").n_nodes(); + #ifdef WITH_VTK // write mesh to vtk file mito::io::vtk::writer("rectangle_copy", mesh, coord_system).write(); @@ -52,4 +60,7 @@ TEST(SummitToSummit, Mesh2D) // check that the number of mesh cells is the same EXPECT_TRUE(original_mesh_cells == reread_mesh_cells); + + // check that the number of mesh nodes is the same + EXPECT_TRUE(original_mesh_nodes == reread_mesh_nodes); } \ No newline at end of file From 3589a5dc824b05fc5c4d3b3a069740c2bc920a0f Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Sun, 9 Jun 2024 18:12:25 +0200 Subject: [PATCH 56/58] lib/mito/io: add support for reading a mesh in {summit} format and returning a discontinuous Galerkin mesh --- lib/mito/io/summit/api.h | 3 ++ lib/mito/io/summit/reader.h | 53 ++++++++++++++----- tests/mito.lib/io/summit_to_summit_mesh_2D.cc | 22 ++++++-- 3 files changed, 59 insertions(+), 19 deletions(-) diff --git a/lib/mito/io/summit/api.h b/lib/mito/io/summit/api.h index ccf7e5d2..1963ee3d 100644 --- a/lib/mito/io/summit/api.h +++ b/lib/mito/io/summit/api.h @@ -13,6 +13,9 @@ namespace mito::io::summit { template requires(utilities::same_dim_c) using mesh_writer_t = MeshWriterSummit; + + // the type of mesh (either CG or DG) + enum GalerkinMeshType {CG, DG}; } diff --git a/lib/mito/io/summit/reader.h b/lib/mito/io/summit/reader.h index 14bb400f..5f373f19 100644 --- a/lib/mito/io/summit/reader.h +++ b/lib/mito/io/summit/reader.h @@ -37,7 +37,7 @@ namespace mito::io::summit { return; } - template + template auto readElement( std::ifstream & fileStream, mesh::mesh_t & mesh, const std::vector> & nodes) -> void @@ -58,16 +58,41 @@ namespace mito::io::summit { index[i] = id; } - // helper function to call the mesh insert method with {N} nodes - constexpr auto _insert = []( - mesh::mesh_t & mesh, const std::array index, - const std::vector> & nodes, - std::index_sequence) { - mesh.insert({ nodes[index[I]]... }); - }; + // if it is a continuous Galerkin mesh + if constexpr (galerkinT == CG) { - // insert the nodes in the mesh - _insert(mesh, index, nodes, std::make_index_sequence{}); + // helper function to call the mesh insert method with {N} nodes + constexpr auto _insert = []( + mesh::mesh_t & mesh, const std::array index, + const std::vector> & nodes, + std::index_sequence) { + // insert in the mesh a geometric simplex with these nodes + mesh.insert({ nodes[index[I]]... }); + }; + + // insert the nodes in the mesh + _insert(mesh, index, nodes, std::make_index_sequence{}); + + } + // otherwise + else { + // assert that it is then a discontinuous Galerkin mesh + static_assert(galerkinT == DG); + + // helper function to call the mesh insert method with {N} nodes + constexpr auto _insert = []( + mesh::mesh_t & mesh, const std::array index, + const std::vector> & nodes, + std::index_sequence) { + // insert in the mesh a geometric simplex with a new instance of the nodes riding on + // same vertex and same point + mesh.insert({ geometry::node_t( + nodes[index[I]]->vertex(), nodes[index[I]]->point())... }); + }; + + // insert the nodes in the mesh + _insert(mesh, index, nodes, std::make_index_sequence{}); + } // QUESTION: Can the label be more than one? // read label for cell @@ -79,7 +104,7 @@ namespace mito::io::summit { return; } - template + template auto readElements( std::ifstream & fileStream, mesh::mesh_t & mesh, int N_cells, const std::vector> & nodes) -> void @@ -93,7 +118,7 @@ namespace mito::io::summit { // if the cell type read from file matches with the cell of the mesh to be populated if (cell_type == summit::cell::type) { // read the element and insert it in the mesh - readElement(fileStream, mesh, nodes); + readElement(fileStream, mesh, nodes); } } @@ -101,7 +126,7 @@ namespace mito::io::summit { return; } - template + template auto reader( std::ifstream & fileStream, geometry::coordinate_system_t & coordinate_system) -> mesh::mesh_t @@ -152,7 +177,7 @@ namespace mito::io::summit { readVertices(fileStream, coordinate_system, N_vertices, nodes); // read the cells - readElements(fileStream, mesh, N_cells, nodes); + readElements(fileStream, mesh, N_cells, nodes); // sanity check: the number of nodes in the map is N_vertices nodes.shrink_to_fit(); diff --git a/tests/mito.lib/io/summit_to_summit_mesh_2D.cc b/tests/mito.lib/io/summit_to_summit_mesh_2D.cc index b21c50f1..869a367c 100644 --- a/tests/mito.lib/io/summit_to_summit_mesh_2D.cc +++ b/tests/mito.lib/io/summit_to_summit_mesh_2D.cc @@ -12,7 +12,9 @@ using coordinates_t = mito::geometry::coordinates_t<2, mito::geometry::CARTESIAN>; -TEST(SummitToSummit, Mesh2D) +template +void +test() { int original_mesh_cells = 0; int reread_mesh_cells = 0; @@ -25,8 +27,8 @@ TEST(SummitToSummit, Mesh2D) // read summit mesh std::ifstream fileStream("rectangle.summit"); - auto mesh = - mito::io::summit::reader>(fileStream, coord_system); + auto mesh = mito::io::summit::reader, galerkinT>( + fileStream, coord_system); // get the original number of mesh cells original_mesh_cells = mesh.nCells(); @@ -43,8 +45,8 @@ TEST(SummitToSummit, Mesh2D) // read the written summit mesh std::ifstream fileStream("rectangle_copy.summit"); - auto mesh = - mito::io::summit::reader>(fileStream, coord_system); + auto mesh = mito::io::summit::reader, galerkinT>( + fileStream, coord_system); // get the reread number of mesh cells reread_mesh_cells = mesh.nCells(); @@ -63,4 +65,14 @@ TEST(SummitToSummit, Mesh2D) // check that the number of mesh nodes is the same EXPECT_TRUE(original_mesh_nodes == reread_mesh_nodes); +} + +TEST(SummitToSummit, Mesh2D_CG) +{ + test(); +} + +TEST(SummitToSummit, Mesh2D_DG) +{ + test(); } \ No newline at end of file From 2a308905ea1d8147b9846af3681b82f591c9f32b Mon Sep 17 00:00:00 2001 From: Anwar Koshakji Date: Fri, 5 Jul 2024 13:12:46 +0200 Subject: [PATCH 57/58] .github/workflows: temporarily restricted {numpy} under {2.0} until its relocated headers can be handled gracefully --- .github/workflows/cmake-macos.yaml | 2 +- .github/workflows/cmake-ubuntu.yaml | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/cmake-macos.yaml b/.github/workflows/cmake-macos.yaml index eec851cb..a756e981 100644 --- a/.github/workflows/cmake-macos.yaml +++ b/.github/workflows/cmake-macos.yaml @@ -49,7 +49,7 @@ jobs: - name: install dependencies run: | python3 -m pip install --upgrade pip - pip3 install distro numpy pybind11 pytest vtk + pip3 install distro 'numpy<2.0' pybind11 pytest vtk - name: checkout gtest uses: actions/checkout@v3 diff --git a/.github/workflows/cmake-ubuntu.yaml b/.github/workflows/cmake-ubuntu.yaml index c85536e0..3055477b 100644 --- a/.github/workflows/cmake-ubuntu.yaml +++ b/.github/workflows/cmake-ubuntu.yaml @@ -53,7 +53,7 @@ jobs: - name: install dependencies run: | python3 -m pip install --upgrade pip - pip3 install distro numpy pybind11 pytest vtk + pip3 install distro 'numpy<2.0' pybind11 pytest vtk - name: checkout pyre uses: actions/checkout@v3 From ed60c8a89ad6f15efd5cda346d3fe00fc16fe80b Mon Sep 17 00:00:00 2001 From: Bianca Giovanardi Date: Mon, 19 Aug 2024 13:07:39 +0200 Subject: [PATCH 58/58] lib/mito/utilities: minor redesign of class {StdSharedPointer} Deriving {StdSharedPointer} from {std::shared_ptr} leads to memory leaks (in this case, memory not being freed). --- lib/mito/utilities/StdSharedPointer.h | 10 +++++++++- lib/mito/utilities/StdSharedPointer.icc | 14 ++++++++++++-- 2 files changed, 21 insertions(+), 3 deletions(-) diff --git a/lib/mito/utilities/StdSharedPointer.h b/lib/mito/utilities/StdSharedPointer.h index d36e8650..97b3666a 100644 --- a/lib/mito/utilities/StdSharedPointer.h +++ b/lib/mito/utilities/StdSharedPointer.h @@ -12,17 +12,21 @@ namespace mito::utilities { // declaration template - class StdSharedPointer : public std::shared_ptr { + class StdSharedPointer { // types public: using shared_ptr_type = StdSharedPointer; using resource_type = resourceT; + using handle_type = resourceT *; // interface public: // returns the id of this (oriented) simplex inline auto id() const -> index_t; + // operator-> + auto operator->() const noexcept -> handle_type; + // meta methods public: // constructor @@ -44,6 +48,10 @@ namespace mito::utilities { // move assignment operator inline StdSharedPointer & operator=(shared_ptr_type &&) noexcept = default; + + // data members + private: + std::shared_ptr _ptr; }; template diff --git a/lib/mito/utilities/StdSharedPointer.icc b/lib/mito/utilities/StdSharedPointer.icc index 58037996..ddc764e9 100644 --- a/lib/mito/utilities/StdSharedPointer.icc +++ b/lib/mito/utilities/StdSharedPointer.icc @@ -13,15 +13,25 @@ auto mito::utilities::StdSharedPointer::id() const -> index_t { // the id is the (immutable) address of this object - return reinterpret_cast>(this->get()); + return reinterpret_cast>(_ptr.get()); } template template requires(std::is_constructible_v) mito::utilities::StdSharedPointer::StdSharedPointer(Args &&... args) : - std::shared_ptr(std::make_shared(args...)) + _ptr(std::make_shared(args...)) {} +// operator-> +template +auto +mito::utilities::StdSharedPointer::operator->() const noexcept + -> mito::utilities::StdSharedPointer::handle_type +{ + // return the handle + return _ptr.operator->(); +} + #endif // end of file