Add support for FEM isoparametric elements (segments and triangles) generalized to an arbitrary dimension D of the physical space

.cmake/mito_tests_mito_lib.cmake

@@ -35,6 +35,7 @@ mito_test_driver(tests/mito.lib/geometry/point.cc)
 mito_test_driver(tests/mito.lib/geometry/euclidean_metric_1D.cc)
 mito_test_driver(tests/mito.lib/geometry/euclidean_metric_2D.cc)
 mito_test_driver(tests/mito.lib/geometry/euclidean_metric_3D.cc)
+mito_test_driver(tests/mito.lib/geometry/euclidean_submanifold_metric_2D.cc)
 mito_test_driver(tests/mito.lib/geometry/euclidean_submanifold_metric_3D.cc)
 mito_test_driver(tests/mito.lib/geometry/induced_metric_segment.cc)
 mito_test_driver(tests/mito.lib/geometry/induced_metric_triangle.cc)
@@ -56,8 +57,12 @@ mito_test_driver(tests/mito.lib/discrete/mesh_field.cc)
 # fem
 mito_test_driver(tests/mito.lib/fem/block_grad_grad.cc)
 mito_test_driver(tests/mito.lib/fem/block_grad_grad_segment.cc)
+mito_test_driver(tests/mito.lib/fem/block_grad_grad_embedded_segment.cc)
+mito_test_driver(tests/mito.lib/fem/block_grad_grad_embedded_triangle.cc)
 mito_test_driver(tests/mito.lib/fem/block_mass.cc)
 mito_test_driver(tests/mito.lib/fem/block_mass_segment.cc)
+mito_test_driver(tests/mito.lib/fem/block_mass_embedded_segment.cc)
+mito_test_driver(tests/mito.lib/fem/block_mass_embedded_triangle.cc)
 mito_test_driver(tests/mito.lib/fem/block_algebra.cc)
 mito_test_driver(tests/mito.lib/fem/shape_functions_triangle_construction.cc)
 mito_test_driver(tests/mito.lib/fem/shape_functions_triangle_p1.cc)
@@ -66,6 +71,8 @@ mito_test_driver(tests/mito.lib/fem/isoparametric_triangle.cc)
 mito_test_driver(tests/mito.lib/fem/fem_field.cc)
 mito_test_driver(tests/mito.lib/fem/shape_functions_segment_p1.cc)
 mito_test_driver(tests/mito.lib/fem/isoparametric_segment.cc)
+mito_test_driver(tests/mito.lib/fem/isoparametric_embedded_segment.cc)
+mito_test_driver(tests/mito.lib/fem/isoparametric_embedded_triangle.cc)
 
 # io
 mito_test_driver(tests/mito.lib/io/summit_mesh_reader_2D.cc)
@@ -189,8 +196,10 @@ mito_test_driver(tests/mito.lib/quadrature/quadrature_parametric_segment.cc)
 mito_test_driver(tests/mito.lib/quadrature/quadrature_parametric_triangle.cc)
 mito_test_driver(tests/mito.lib/quadrature/quadrature_parametric_tetrahedron.cc)
 mito_test_driver(tests/mito.lib/quadrature/quadrature_segment_1D.cc)
+mito_test_driver(tests/mito.lib/quadrature/quadrature_segment_2D.cc)
 mito_test_driver(tests/mito.lib/quadrature/quadrature_segment_3D.cc)
 mito_test_driver(tests/mito.lib/quadrature/quadrature_triangle_2D.cc)
+mito_test_driver(tests/mito.lib/quadrature/quadrature_triangle_3D.cc)
 mito_test_driver(tests/mito.lib/quadrature/quadrature_square_3D.cc)
 
 # integration

.github/workflows/cmake-macos.yaml

@@ -82,7 +82,7 @@ jobs:
           repository: pyre/pyre
           fetch-depth: 0
           path: ${{github.workspace}}/pyre
-          ref: tensor
+          ref: main
 
       - name: retrieve pyre commit SHA
         id: pyre-sha

.github/workflows/cmake-ubuntu.yaml

@@ -57,7 +57,7 @@ jobs:
           repository: pyre/pyre
           fetch-depth: 0
           path: ${{github.workspace}}/pyre
-          ref: tensor
+          ref: main
 
       - name: retrieve pyre commit SHA
         id: pyre-sha

lib/mito/fem/blocks/GradientGradientBlock.h

@@ -23,10 +23,6 @@ namespace mito::fem::blocks {
         // the type of the coefficient field
         using coefficient_field_type = coefficientFieldT;
 
-      public:
-        // instantiate the quadrature rule
-        static constexpr auto quadrature_rule = quadrature_rule_type();
-
       public:
         // constructor
         GradientGradientBlock(const coefficient_field_type & coefficient) :
@@ -39,53 +35,40 @@ namespace mito::fem::blocks {
         requires element_of_type_c<elementT, element_type>
         auto compute(const elementT & element) const -> elementary_shape
         {
-            // the number of nodes per element
-            constexpr int n_nodes = element_type::n_nodes;
-
-            // the number of quadrature points per element
-            constexpr int n_quads = quadrature_rule_type::npoints;
+            // the parametric coordinates type
+            using parametric_coordinates_type = typename elementT::parametric_coordinates_type;
 
             // the elementary matrix
-            elementary_shape elementary_matrix{};
-
-            // loop on the quadrature points
-            tensor::constexpr_for_1<n_quads>([&]<int q>() {
-                // the parametric coordinates of the quadrature point
-                constexpr auto xi = quadrature_rule.point(q);
-
-                // the coordinates of the quadrature point
-                auto x = element.parametrization()(xi);
-
-                // evaluate the coefficient field at the quadrature point
-                auto coefficient = _coefficient(x);
+            return manifolds::cell_integrator<quadrature_rule_type>(element.element())
+                .integrate(mito::functions::function([&](const parametric_coordinates_type & xi) {
+                    // the elementary contribution at quadrature point {xi}
+                    elementary_shape elementary_matrix{};
 
-                // the measure of the canonical simplex
-                constexpr auto measure =
-                    element_type::mesh_cell_type::reference_simplex_type::measure;
+                    // the number of nodes per element
+                    constexpr int n_nodes = element_type::n_nodes;
 
-                // the quadrature weight at this point scaled with the area of the canonical simplex
-                constexpr auto w = measure * quadrature_rule.weight(q);
+                    // the coordinates of the quadrature point
+                    auto x = element.parametrization()(xi);
 
-                // precompute the common factor
-                auto factor = w * tensor::determinant(element.jacobian()(xi));
+                    // evaluate the coefficient at the quadrature point
+                    auto coefficient = _coefficient(x);
 
-                // loop on the nodes of the element
-                tensor::constexpr_for_1<n_nodes>([&]<int a>() {
-                    // evaluate the spatial gradient of the element's a-th shape function at {xi}
-                    auto dphi_a = element.template gradient<a>()(xi);
                     // loop on the nodes of the element
-                    tensor::constexpr_for_1<n_nodes>([&]<int b>() {
-                        // evaluate the spatial gradient of the element's b-th shape function at
-                        // {xi}
-                        auto dphi_b = element.template gradient<b>()(xi);
-                        // populate the elementary contribution to the matrix
-                        elementary_matrix[{ a, b }] += factor * dphi_a * (coefficient * dphi_b);
+                    tensor::constexpr_for_1<n_nodes>([&]<int a>() {
+                        // evaluate the spatial gradient of the a-th shape function at {xi}
+                        auto dphi_a = element.template gradient<a>()(xi);
+                        // loop on the nodes of the element
+                        tensor::constexpr_for_1<n_nodes>([&]<int b>() {
+                            // evaluate the spatial gradient of the b-th shape function at {xi}
+                            auto dphi_b = element.template gradient<b>()(xi);
+                            // populate the elementary contribution to the matrix
+                            elementary_matrix[{ a, b }] = dphi_a * (coefficient * dphi_b);
+                        });
                     });
-                });
-            });
 
-            // all done
-            return elementary_matrix;
+                    // all done
+                    return elementary_matrix;
+                }));
         }
 
       private:

lib/mito/fem/blocks/L2NormBlock.h

@@ -16,18 +16,16 @@ namespace mito::fem::blocks {
     class L2NormBlock {
 
       public:
-        // my template parameters
+        // my finite element type
         using element_type = finiteElementT;
-        using elementary_shape = tensor::scalar_t;
+        // my quadrature rule
         using quadrature_rule_type = quadratureRuleT;
+        // my elementary shape
+        using elementary_shape = tensor::scalar_t;
 
         // the type of the function to compute the L2 norm of
         using function_type = functionT;
 
-      public:
-        // instantiate the quadrature rule
-        static constexpr auto quadrature_rule = quadrature_rule_type();
-
       public:
         // constructor
         L2NormBlock(const function_type & function) : _function(function) {}
@@ -38,33 +36,24 @@ namespace mito::fem::blocks {
         requires element_of_type_c<elementT, element_type>
         auto compute(const elementT & element) const -> elementary_shape
         {
-            // the number of quadrature points per element
-            constexpr int n_quads = quadrature_rule_type::npoints;
-
-            // the elementary contribution to the L2 norm
-            elementary_shape norm{};
-
-            // loop on the quadrature points
-            tensor::constexpr_for_1<n_quads>([&]<int q>() {
-                // the parametric coordinates of the quadrature point
-                constexpr auto xi = quadrature_rule.point(q);
-
-                // the measure of the canonical simplex
-                constexpr auto measure =
-                    element_type::mesh_cell_type::reference_simplex_type::measure;
+            // the parametric coordinates type
+            using parametric_coordinates_type = typename elementT::parametric_coordinates_type;
 
-                // the quadrature weight at this point scaled with the area of the canonical simplex
-                constexpr auto w = measure * quadrature_rule.weight(q);
+            // the elementary matrix
+            return manifolds::cell_integrator<quadrature_rule_type>(element.element())
+                .integrate(mito::functions::function([&](const parametric_coordinates_type & xi) {
+                    // the elementary contribution at quadrature point {xi}
+                    elementary_shape norm{};
 
-                // precompute the common factor
-                auto factor = w * tensor::determinant(element.jacobian()(xi));
+                    // evaluate the function at the quadrature point
+                    auto fx = _function(xi);
 
-                // populate the elementary contribution to the matrix
-                norm += factor * _function(xi) * _function(xi);
-            });
+                    // assemble the elementary contribution
+                    norm = fx * fx;
 
-            // all done
-            return norm;
+                    // all done
+                    return norm;
+                }));
         }
 
       private:

lib/mito/fem/blocks/ValueBlock.h

@@ -25,10 +25,6 @@ namespace mito::fem::blocks {
         // the type of the coefficient field
         using coefficient_field_type = coefficientFieldT;
 
-      public:
-        // instantiate the quadrature rule
-        static constexpr auto quadrature_rule = quadrature_rule_type();
-
       public:
         // constructor
         ValueBlock(const coefficient_field_type & coefficient) : _coefficient(coefficient) {}
@@ -39,44 +35,35 @@ namespace mito::fem::blocks {
         requires(element_of_type_c<elementT, element_type>)
         auto compute(const elementT & element) const -> elementary_shape
         {
-            // the number of nodes per element
-            constexpr int n_nodes = element_type::n_nodes;
-
-            // the number of quadrature points per element
-            constexpr int n_quads = quadrature_rule_type::npoints;
+            // the parametric coordinates type
+            using parametric_coordinates_type = typename elementT::parametric_coordinates_type;
 
             // the elementary vector
-            elementary_shape elementary_vector{};
-
-            // loop on the quadrature points
-            tensor::constexpr_for_1<n_quads>([&]<int q>() {
-                // the parametric coordinates of the quadrature point
-                constexpr auto xi = quadrature_rule.point(q);
-
-                // the coordinates of the quadrature point
-                auto coord = element.parametrization()(xi);
-
-                // the measure of the canonical simplex
-                constexpr auto measure =
-                    element_type::mesh_cell_type::reference_simplex_type::measure;
-
-                // the quadrature weight at this point scaled with the area of the canonical simplex
-                constexpr auto w = measure * quadrature_rule.weight(q);
-
-                // precompute the common factor
-                auto factor = w * _coefficient(coord) * tensor::determinant(element.jacobian()(xi));
-
-                // loop on the nodes of the element
-                tensor::constexpr_for_1<n_nodes>([&]<int a>() {
-                    // evaluate the a-th shape function at {xi}
-                    auto phi_a = element.template shape<a>()(xi);
-                    // populate the elementary contribution to the vector
-                    elementary_vector[{ a }] += factor * phi_a;
-                });
-            });
-
-            // all done
-            return elementary_vector;
+            return manifolds::cell_integrator<quadrature_rule_type>(element.element())
+                .integrate(mito::functions::function([&](const parametric_coordinates_type & xi) {
+                    // the elementary contribution at quadrature point {xi}
+                    elementary_shape elementary_vector{};
+
+                    // the number of nodes per element
+                    constexpr int n_nodes = element_type::n_nodes;
+
+                    // the coordinates of the quadrature point
+                    auto x = element.parametrization()(xi);
+
+                    // evaluate the coefficient at the quadrature point
+                    auto coefficient = _coefficient(x);
+
+                    // loop on the nodes of the element
+                    tensor::constexpr_for_1<n_nodes>([&]<int a>() {
+                        // evaluate the element's a-th shape function at {xi}
+                        const auto phi_a = element.template shape<a>()(xi);
+                        // populate the elementary contribution to the vector
+                        elementary_vector[{ a }] = coefficient * phi_a;
+                    });
+
+                    // all done
+                    return elementary_vector;
+                }));
         }
 
       private:

lib/mito/fem/blocks/ValueGradientBlock.h

@@ -23,10 +23,6 @@ namespace mito::fem::blocks {
         // the type of the coefficient field
         using coefficient_field_type = coefficientFieldT;
 
-      public:
-        // instantiate the quadrature rule
-        static constexpr auto quadrature_rule = quadrature_rule_type();
-
       public:
         // constructor
         ValueGradientBlock(const coefficient_field_type & coefficient) : _coefficient(coefficient)
@@ -38,53 +34,40 @@ namespace mito::fem::blocks {
         requires element_of_type_c<elementT, element_type>
         auto compute(const elementT & element) const -> elementary_shape
         {
-            // the number of nodes per element
-            constexpr int n_nodes = element_type::n_nodes;
-
-            // the number of quadrature points per element
-            constexpr int n_quads = quadrature_rule_type::npoints;
+            // the parametric coordinates type
+            using parametric_coordinates_type = typename elementT::parametric_coordinates_type;
 
             // the elementary matrix
-            elementary_shape elementary_matrix{};
-
-            // loop on the quadrature points
-            tensor::constexpr_for_1<n_quads>([&]<int q>() {
-                // the parametric coordinates of the quadrature point
-                constexpr auto xi = quadrature_rule.point(q);
-
-                // the coordinates of the quadrature point
-                auto x = element.parametrization()(xi);
-
-                // evaluate the coefficient at the quadrature point
-                auto coefficient = _coefficient(x);
+            return manifolds::cell_integrator<quadrature_rule_type>(element.element())
+                .integrate(mito::functions::function([&](const parametric_coordinates_type & xi) {
+                    // the elementary contribution at quadrature point {xi}
+                    elementary_shape elementary_matrix{};
 
-                // the measure of the canonical simplex
-                constexpr auto measure =
-                    element_type::mesh_cell_type::reference_simplex_type::measure;
+                    // the number of nodes per element
+                    constexpr int n_nodes = element_type::n_nodes;
 
-                // the quadrature weight at this point scaled with the area of the canonical simplex
-                constexpr auto w = measure * quadrature_rule.weight(q);
+                    // the coordinates of the quadrature point
+                    auto x = element.parametrization()(xi);
 
-                // precompute the common factor
-                auto factor = w * tensor::determinant(element.jacobian()(xi));
+                    // evaluate the coefficient at the quadrature point
+                    auto coefficient = _coefficient(x);
 
-                // loop on the nodes of the element
-                tensor::constexpr_for_1<n_nodes>([&]<int a>() {
-                    // evaluate the element's a-th shape function at {xi}
-                    auto phi_a = element.template shape<a>()(xi);
                     // loop on the nodes of the element
-                    tensor::constexpr_for_1<n_nodes>([&]<int b>() {
-                        // evaluate the spatial gradient of the element's b-th shape function at
-                        // {xi}
-                        auto dphi_b = element.template gradient<b>()(xi);
-                        // populate the elementary contribution to the matrix
-                        elementary_matrix[{ a, b }] += factor * phi_a * coefficient * dphi_b;
+                    tensor::constexpr_for_1<n_nodes>([&]<int a>() {
+                        // evaluate the element's a-th shape function at {xi}
+                        const auto phi_a = element.template shape<a>()(xi);
+                        // loop on the nodes of the element
+                        tensor::constexpr_for_1<n_nodes>([&]<int b>() {
+                            // evaluate the spatial gradient of the b-th shape function at {xi}
+                            auto dphi_b = element.template gradient<b>()(xi);
+                            // populate the elementary contribution to the matrix
+                            elementary_matrix[{ a, b }] = phi_a * coefficient * dphi_b;
+                        });
                     });
-                });
-            });
 
-            // all done
-            return elementary_matrix;
+                    // all done
+                    return elementary_matrix;
+                }));
         }
 
       private: