Nodal fields redesign

.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_sphere.cc)
 
 # io
 mito_test_driver(tests/mito.lib/io/read_mesh_summit_2D.cc)

.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

.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

lib/mito/fem/NodalField.h

@@ -9,103 +9,81 @@
 
 namespace mito::fem {
 
-    // TOFIX: order of template parameters
-    template <typename T, int D>
+    template <int D, class Y>
     class NodalField {
+
+      private:
+        // the node type
+        using node_type = geometry::node_t<D>;
+
+        // a map from {key_type} to {Y} values
+        using map_type = std::unordered_map<node_type, Y, utilities::hash_function<node_type>>;
+
       public:
-        NodalField(int nodes, std::string name = "") :
-            _nodes(nodes),
-            _name(name),
-            _nodalField(_nodes * D, 0.0)
-        {}
+        // constructor
+        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();
+                }
+            }
+        }
 
-        ~NodalField() {}
+        // destructor
+        ~NodalField() = default;
 
+      public:
         /**
          * Operator()
          */
-        const T & operator()(const int a, const int i) const
+        inline auto operator()(const node_type & node) const -> const Y &
         {
-            assert(i < D);
-            return _nodalField[D * a + i];
+            return _map_nodes_to_values.at(node);
         }
 
-        T & operator()(const int a, const int i)
+        inline auto operator()(const node_type & node) -> Y &
         {
-            assert(i < D);
-            return _nodalField[D * a + i];
+            return _map_nodes_to_values.at(node);
         }
 
         /**
-         * Iterators
+         * accessor for the number of nodes
          */
-        typename std::vector<T>::iterator begin() { return std::begin(_nodalField); }
-
-        typename std::vector<T>::iterator end() { return std::end(_nodalField); }
-
-        typename std::vector<T>::const_iterator begin() const { return std::cbegin(_nodalField); }
-
-        typename std::vector<T>::const_iterator end() const { return std::cend(_nodalField); }
+        inline auto n_nodes() const -> int { return _map_nodes_to_values.size(); }
 
         /**
-         * 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(_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:
-        /**
-         * number of nodes
-         */
-        int _nodes;
+        // the underlying mapping of nodes to nodal values
+        map_type _map_nodes_to_values;
 
-        /**
-         * name of the nodal field
-         */
+        // the name of the field
         std::string _name;
-
-        /**
-         * nodal field
-         */
-        std::vector<T> _nodalField;
     };
 
-    template <typename T, int D>
-    std::ostream & operator<<(std::ostream & os, const NodalField<T, D> & nodalField)
+    template <int D, class Y>
+    std::ostream & operator<<(std::ostream & os, const NodalField<D, Y> & 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);
+        for (auto i = 0; i < nodalField.n_nodes(); ++i) {
+            os << "Node " << i << ": " << nodalField(i) << std::endl;
         }
-        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 << ")";
-        }
-
-        os << "]";
 
         return os;
     }

lib/mito/fem/QuadratureField.h

@@ -14,32 +14,32 @@ 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<Y>;
+        using storage_type = pyre::memory::heap_t<Y>;
         // putting it all together
-        using grid_t = pyre::grid::grid_t<pack_t, storage_t>;
+        using grid_type = pyre::grid::grid_t<pack_type, storage_type>;
         // index
-        using index_t = pack_t::index_type;
+        using index_type = pack_type::index_type;
 
       public:
         /**
          * 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_type{ { nElements, Q } }, name)
         {}
 
       private:
-        inline QuadratureField(int /*nElements*/, const pack_t & packing, std::string name = "") :
+        inline QuadratureField(const pack_type && packing, std::string name) :
             _grid{ packing, packing.cells() },
             _name(name)
         {}
 
       public:
         // destructor
-        ~QuadratureField() {}
+        ~QuadratureField() = default;
 
       public:
         /**
@@ -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_t & index)
+        inline auto operator[](const index_type & index) -> Y &
         {
             // all done
             return _grid[index];
         }
 
-        inline const Y & operator[](const index_t & index) const
+        inline auto operator[](const index_type & index) const -> const Y &
         {
             // all done
             return _grid[index];
@@ -82,27 +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; }
+        constexpr auto n_quad_points() const noexcept -> int { return Q; }
 
         /**
-         * const accessor for name
+         * accessor for name
          */
-        inline std::string name() const noexcept { return _name; }
-
-        /**
-         * setter method for name
-         */
-        inline void name(std::string name) noexcept
-        {
-            _name = name;
-            return;
-        }
+        inline auto name() const noexcept -> std::string { return _name; }
 
         // support for ranged for loops (wrapping grid)
         inline auto begin() const { return std::cbegin(_grid); }
@@ -111,10 +102,10 @@ namespace mito::fem {
         inline auto end() { return std::end(_grid); }
 
       private:
-        // instantiate the grid
-        grid_t _grid;
+        // the underlying grid
+        grid_type _grid;
 
-        // name of the field
+        // the name of the field
         std::string _name;
     };
 

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 <typename T, int D>
-    constexpr auto nodal_field(int nodes, std::string name = "");
-
+    template <class Y, int D>
+    constexpr auto nodal_field(const mesh::mesh_c auto & mesh, std::string name);
 }
 
 

lib/mito/fem/externals.h

@@ -8,13 +8,15 @@
 
 
 // externals
+#define HAVE_TENSOR
+#define HAVE_COMPACT_PACKINGS
 #include <pyre/grid.h>
 #include <string>
-#include <cassert>
 
 // support
 #include "../journal.h"
 #include "../tensor.h"
+#include "../mesh.h"
 
 
 // end of file

lib/mito/fem/factories.h

@@ -17,10 +17,10 @@ namespace mito::fem {
     }
 
     // nodal field factory
-    template <typename T, int D>
-    constexpr auto nodal_field(int nodes, std::string name)
+    template <class Y, mesh::mesh_c meshT>
+    constexpr auto nodal_field(const meshT & mesh, std::string name)
     {
-        return nodal_field_t<T, D>(nodes, name);
+        return nodal_field_t<meshT::dim, Y>(mesh, name);
     }
 }
 

lib/mito/fem/forward.h

@@ -10,12 +10,12 @@
 namespace mito::fem {
 
     // class nodal field
-    template <typename T, int D>
+    template <int D, class Y>
     class NodalField;
 
     // nodal field alias
-    template <typename T, int D>
-    using nodal_field_t = NodalField<T, D>;
+    template <int D, class Y>
+    using nodal_field_t = NodalField<D, Y>;
 
     // class quadrature field
     template <int Q, class Y>

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 <class FIELD1, class FIELD2>
     // a scalar field  is a field returning a scalar
-    concept compatible_fields_c = requires {
+    concept compatible_fields_c =
         std::is_same_v<typename FIELD1::coordinates_type, typename FIELD2::coordinates_type>;
-    };
 }
 
 

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<dim>;
 
       private:
         // alias {dim} for internal purposes
         static constexpr int D = dim;
-        // a point
-        using point_type = point_t<D>;
         // id type of point
         using point_id_type = utilities::index_t<point_type>;
         // a map between points and their coordinates