Redesign of continuous and discrete fields

.cmake/mito_tests_mito_lib.cmake

@@ -41,7 +41,7 @@ mito_test_driver(tests/mito.lib/geometry/spherical_metric_space.cc)
 mito_test_driver(tests/mito.lib/constraints/dirichlet.cc)
 
 # discrete
-mito_test_driver(tests/mito.lib/discrete/quadrature_field.cc)
+mito_test_driver(tests/mito.lib/discrete/point_field.cc)
 mito_test_driver(tests/mito.lib/discrete/mesh_field.cc)
 
 # fem
@@ -51,6 +51,8 @@ mito_test_driver(tests/mito.lib/fem/shape_functions_triangle_construction.cc)
 mito_test_driver(tests/mito.lib/fem/shape_functions_triangle_p1.cc)
 mito_test_driver(tests/mito.lib/fem/shape_functions_triangle_p2.cc)
 mito_test_driver(tests/mito.lib/fem/isoparametric_triangle.cc)
+mito_test_driver(tests/mito.lib/fem/domain_field.cc)
+mito_test_driver(tests/mito.lib/fem/fem_field.cc)
 
 # io
 mito_test_driver(tests/mito.lib/io/summit_mesh_reader_2D.cc)

benchmarks/mito.lib/fields/laplacian.cc

@@ -47,7 +47,7 @@ laplacian_mito(const coordinates_t & x)
     constexpr auto x1 = mito::functions::component<coordinates_t, 1>;
 
     // a scalar field
-    constexpr auto f = mito::fields::field(mito::functions::pow<4>(x0 * x1));
+    constexpr auto f = mito::functions::pow<4>(x0 * x1);
 
     // the gradient of {f}
     constexpr auto gradient = mito::fields::gradient(f);

benchmarks/mito.lib/integration/integration.cc

@@ -70,7 +70,7 @@ main()
     auto manifold = mito::manifolds::manifold(tetra_mesh, coord_system);
 
     // instantiate a scalar field
-    auto f = mito::fields::field(mito::functions::cos(x_0 * x_1));
+    auto f = mito::functions::cos(x_0 * x_1);
 
     // instantiate a GAUSS integrator with degree of exactness equal to 2
     auto integrator = mito::quadrature::integrator<mito::quadrature::GAUSS, 2>(manifold);

benchmarks/mito.lib/pdes/poisson.cc

@@ -63,7 +63,7 @@ main()
     auto boundary_mesh = mito::mesh::boundary(mesh);
 
     // the zero field
-    auto zero = mito::fields::field(mito::functions::zero<coordinates_t>);
+    auto zero = mito::functions::zero<coordinates_t>;
 
     // set homogeneous Dirichlet boundary condition
     auto constraints = mito::constraints::dirichlet_bc(boundary_mesh, zero);
@@ -75,9 +75,8 @@ main()
     auto fem_lhs_block = mito::fem::blocks::grad_grad_block<finite_element_t, quadrature_rule_t>();
 
     // the right hand side
-    auto f = mito::fields::field(
-        2.0 * std::numbers::pi * std::numbers::pi * mito::functions::sin(std::numbers::pi * x)
-        * mito::functions::sin(std::numbers::pi * y));
+    auto f = 2.0 * std::numbers::pi * std::numbers::pi * mito::functions::sin(std::numbers::pi * x)
+           * mito::functions::sin(std::numbers::pi * y);
     // channel << "Right hand side: " << f(coordinates_t{ 0.5, 0.5 }) << journal::endl;
 
     // a source term block
@@ -105,17 +104,22 @@ main()
     // free the solver
     solver.destroy();
 
+    // get the solution field
+    const auto & solution = discrete_system.solution();
+
     // the exact solution field
-    auto u_ex = mito::fields::field(
-        mito::functions::sin(std::numbers::pi * x) * mito::functions::sin(std::numbers::pi * y));
+    auto u_ex =
+        mito::functions::sin(std::numbers::pi * x) * mito::functions::sin(std::numbers::pi * y);
 
     // compute the L2 error
-    auto error_L2 = discrete_system.compute_l2_error<quadrature_rule_t>(u_ex);
+    auto error_L2 = mito::fem::compute_l2_norm<quadrature_rule_t>(
+        function_space, solution, mito::fem::domain_field(u_ex));
     // report
     channel << "L2 error: " << error_L2 << journal::endl;
 
     // compute the H1 error
-    auto error_H1 = discrete_system.compute_h1_error<quadrature_rule_t>(u_ex);
+    auto error_H1 = mito::fem::compute_h1_norm<quadrature_rule_t>(
+        function_space, solution, mito::fem::domain_field(u_ex));
     // report
     channel << "H1 error: " << error_H1 << journal::endl;
 
@@ -133,8 +137,6 @@ main()
     // write output file
     writer.write();
 
-    // get the solution field
-    const auto & solution = discrete_system.solution();
     // write mesh to vtk file
     auto writer_solution =
         mito::io::vtk::field_writer("poisson_square_solution", function_space, coord_system);

extensions/mito/mito.cc

@@ -73,9 +73,7 @@ PYBIND11_MODULE(mito, m)
 
     // the mito scalar field 2D
     using scalar_function_2D_t = std::function<mito::tensor::scalar_t(const coordinates_2D_t &)>;
-    using scalar_function_functor_2D_t = mito::functions::FunctionFromFunctor<scalar_function_2D_t>;
-    using scalar_field_2D_t =
-        decltype(mito::fields::field(std::declval<scalar_function_functor_2D_t>()));
+    using scalar_field_2D_t = mito::functions::FunctionFromFunctor<scalar_function_2D_t>;
     mito::py::class_<scalar_field_2D_t>(m, "ScalarField2D")
         // the constructor
         .def(mito::py::init<scalar_function_2D_t>())
@@ -90,9 +88,7 @@ PYBIND11_MODULE(mito, m)
 
     // the mito scalar field 3D
     using scalar_function_3D_t = std::function<mito::tensor::scalar_t(const coordinates_3D_t &)>;
-    using scalar_function_functor_3D_t = mito::functions::FunctionFromFunctor<scalar_function_3D_t>;
-    using scalar_field_3D_t =
-        decltype(mito::fields::field(std::declval<scalar_function_functor_3D_t>()));
+    using scalar_field_3D_t = mito::functions::FunctionFromFunctor<scalar_function_3D_t>;
     mito::py::class_<scalar_field_3D_t>(m, "ScalarField3D")
         // the constructor
         .def(mito::py::init<scalar_function_3D_t>())

lib/mito/discrete/DiscreteField.h

@@ -12,15 +12,16 @@ namespace mito::discrete {
     template <class keyT, class valueT>
     class DiscreteField {
 
-      private:
-        // the key type
-        using key_type = keyT;
-        // the value type
+      public:
+        // the input type of the field
+        using input_type = keyT;
+        // the value type returned by the field
         using value_type = valueT;
 
-        // a map from {key_type} to {value_type} values
+      private:
+        // a map from {input_type} to {value_type} values
         using map_type =
-            std::unordered_map<key_type, value_type, utilities::hash_function<key_type>>;
+            std::unordered_map<input_type, value_type, utilities::hash_function<input_type>>;
 
       public:
         // constructor
@@ -55,26 +56,30 @@ namespace mito::discrete {
         /**
          * accessor for the value of a given entry
          */
-        inline auto operator()(const key_type & key) const -> const value_type &
+        inline auto operator()(const input_type & key) const -> const value_type &
         {
             return _map_entry_to_values.at(key);
         }
 
         /**
          * mutator for the value of a given entry
          */
-        inline auto operator()(const key_type & key) -> value_type &
+        inline auto operator()(const input_type & key) -> value_type &
         {
             return _map_entry_to_values.at(key);
         }
 
         /**
          * insert a new entry to the field
          */
-        inline auto insert(const key_type & key, const value_type & entry = value_type()) -> void
+        inline auto insert(const input_type & key, const value_type & entry = value_type()) -> void
         {
+            // check that the entry is not already present
+            // (emplace does not overwrite existing entries)
+            assert(_map_entry_to_values.find(key) == std::cend(_map_entry_to_values));
+
             // add a new entry to the field
-            _map_entry_to_values[key] = entry;
+            _map_entry_to_values.emplace(key, entry);
 
             // all done
             return;
@@ -99,7 +104,7 @@ namespace mito::discrete {
         inline auto cend() const { return std::cend(_map_entry_to_values); }
 
       private:
-        // the underlying mapping of entries to nodal values
+        // the underlying mapping of entries to values
         map_type _map_entry_to_values;
 
         // the name of the field

lib/mito/discrete/externals.h

@@ -9,6 +9,7 @@
 
 // externals
 #include <string>
+#include <cassert>
 
 // support
 #include "../mesh.h"

lib/mito/fem/DiscreteSystem.h

@@ -36,8 +36,8 @@ namespace mito::fem {
         // solutions
         // the solution field type
         using solution_field_type = tensor::scalar_t;
-        // the nodal field type
-        using nodal_field_type = discrete::nodal_field_t<solution_field_type>;
+        // the fem field type
+        using fem_field_type = fem_field_t<solution_field_type>;
         // the number of nodes per element
         static constexpr int n_element_nodes = element_type::n_nodes;
 
@@ -49,7 +49,8 @@ namespace mito::fem {
             _function_space(function_space),
             _weakform(weakform),
             _equation_map(),
-            _solution_field(nodal_field<solution_field_type>(function_space, label + ".solution")),
+            _solution_field(
+                function_space.template fem_field<solution_field_type>(label + ".solution")),
             _linear_system(label)
         {
             // make a channel
@@ -198,126 +199,27 @@ namespace mito::fem {
 
             // get the node map from the function space
             auto node_map = _function_space.node_map();
-            // fill information in nodal field
-            for (auto & [node, value] : _solution_field) {
-                // get the equation number of {node}
-                int eq = _equation_map.at(node);
+            // fill information in finite element field
+            for (auto & [node, eq] : _equation_map) {
                 if (eq != -1) {
-                    // read the solution at {eq}
-                    value = u[eq];
+                    // note the solution on the solution field
+                    _solution_field(node) = u[eq];
                 }
             }
 
             // all done
             return;
         }
 
-        // accessor to the solution nodal field
-        constexpr auto solution() const noexcept -> const nodal_field_type &
+        // accessor to the solution finite element field
+        constexpr auto solution() const noexcept -> const fem_field_type &
         {
             return _solution_field;
         }
 
         // accessor to the number of equations
         constexpr auto n_equations() const noexcept -> int { return _n_equations; }
 
-        // compute the L2 norm of the solution
-        template <class quadratureRuleT>
-        constexpr auto compute_l2_error(const fields::field_c auto & u_exact) const
-            -> tensor::scalar_t
-        {
-            // initialize the norm
-            auto norm = tensor::scalar_t{ 0.0 };
-
-            // helper lambda to assemble the solution on a given element
-            constexpr auto _assemble_element_solution = []<int... a>(
-                                                            const auto & element,
-                                                            const auto & _solution_field,
-                                                            tensor::integer_sequence<a...>) {
-                // assemble the solution field from the shape functions
-                return (
-                    (element.template shape<a>() * _solution_field(element.connectivity()[a]))
-                    + ...);
-            };
-
-            // loop on all the cells of the mesh
-            for (const auto & element : _function_space.elements()) {
-                // assemble the solution field on this element
-                auto u_numerical = _assemble_element_solution(
-                    element, _solution_field, tensor::make_integer_sequence<n_element_nodes>());
-                // assemble the elementary error field as a function of the barycentric coordinates
-                auto u_error =
-                    u_numerical - u_exact.function()(element.parametrization().function());
-
-                // compute the elementary contribution to the L2 norm
-                norm += blocks::l2_norm_block<element_type, quadratureRuleT>(u_error.function())
-                            .compute(element);
-            }
-
-            return std::sqrt(norm);
-        }
-
-        // compute the H1 norm of the solution
-        template <class quadratureRuleT>
-        constexpr auto compute_h1_error(const fields::field_c auto & u_exact) const
-            -> tensor::scalar_t
-        {
-            // initialize the norm
-            auto norm = tensor::scalar_t{ 0.0 };
-
-            // QUESTION: is there a proper place to put this method? Perhaps this belongs to the
-            // isoparametric element class?
-            // helper lambda to assemble the solution on a given element
-            constexpr auto _assemble_element_solution = []<int... a>(
-                                                            const auto & element,
-                                                            const auto & _solution_field,
-                                                            tensor::integer_sequence<a...>) {
-                // assemble the solution field from the shape functions
-                return (
-                    (element.template shape<a>() * _solution_field(element.connectivity()[a]))
-                    + ...);
-            };
-
-            // TOFIX: we should deduce this from {_assemble_element_solution}
-            // helper lambda to assemble the solution on a given element
-            constexpr auto _assemble_element_solution_gradient =
-                []<int... a>(
-                    const auto & element, const auto & _solution_field,
-                    tensor::integer_sequence<a...>) {
-                    // assemble the solution field from the shape functions
-                    return (
-                        (element.template gradient<a>()
-                         * _solution_field(element.connectivity()[a]))
-                        + ...);
-                };
-
-            // loop on all the cells of the mesh
-            for (const auto & element : _function_space.elements()) {
-                // assemble the solution field on this element
-                auto u_numerical = _assemble_element_solution(
-                    element, _solution_field, tensor::make_integer_sequence<n_element_nodes>());
-                // assemble the elementary error field as a function of the barycentric coordinates
-                auto u_error =
-                    u_numerical - u_exact.function()(element.parametrization().function());
-                //
-                auto u_numerical_gradient = _assemble_element_solution_gradient(
-                    element, _solution_field, tensor::make_integer_sequence<n_element_nodes>());
-                // assemble the elementary error gradient as a function of the barycentric
-                // coordinates
-                auto u_error_gradient =
-                    u_numerical_gradient
-                    - fields::gradient(u_exact).function()(element.parametrization().function());
-                // compute the elementary contributions to the H1 norm
-                norm += blocks::l2_norm_block<element_type, quadratureRuleT>(u_error.function())
-                            .compute(element)
-                      + blocks::l2_norm_block<element_type, quadratureRuleT>(
-                            u_error_gradient.function())
-                            .compute(element);
-            }
-
-            return std::sqrt(norm);
-        }
-
       private:
         // a const reference to the function space
         const function_space_type & _function_space;
@@ -328,8 +230,8 @@ namespace mito::fem {
         // the equation map
         equation_map_type _equation_map;
 
-        // the solution nodal field
-        nodal_field_type _solution_field;
+        // the solution finite element field
+        fem_field_type _solution_field;
 
         // the linear system of equations
         linear_system_type _linear_system;