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

@@ -30,6 +30,7 @@ mito_test_driver(tests/mito.lib/geometry/cell_directors.cc)
 mito_test_driver(tests/mito.lib/geometry/point.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/cube_volume.cc)
 mito_test_driver(tests/mito.lib/geometry/metric.cc)
@@ -47,8 +48,10 @@ 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_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/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)
@@ -57,6 +60,7 @@ mito_test_driver(tests/mito.lib/fem/localize_field.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)
 
 # io
 mito_test_driver(tests/mito.lib/io/summit_mesh_reader_2D.cc)

benchmarks/mito.lib/pdes/poisson.cc

@@ -13,8 +13,6 @@ using coordinates_t = mito::geometry::coordinates_t<2, mito::geometry::CARTESIAN
 using cell_t = mito::geometry::triangle_t<2>;
 // second degree finite elements
 constexpr int degree = 2;
-// assemble the finite element type
-using finite_element_t = mito::fem::isoparametric_simplex_t<degree, cell_t>;
 
 // the reference simplex
 using reference_simplex_t = mito::geometry::reference_triangle_t;
@@ -59,6 +57,9 @@ main()
     // create the body manifold
     auto manifold = mito::manifolds::manifold(mesh, coord_system);
 
+    // assemble the finite element type
+    using finite_element_t = mito::fem::isoparametric_simplex_t<degree, decltype(manifold)>;
+
     // get the boundary mesh
     auto boundary_mesh = mito::mesh::boundary(mesh);
 

lib/mito/fem/blocks/GradGradBlock.h

@@ -44,8 +44,8 @@ namespace mito::fem::blocks {
                 constexpr auto w =
                     element_type::canonical_element_type::area * quadrature_rule.weight(q);
 
-                // precompute the common factor
-                auto factor = w * tensor::determinant(element.jacobian()(xi));
+                // precompute the common factor using the element's volume element
+                auto factor = w * element.volume_element()(xi);
 
                 // loop on the nodes of the element
                 tensor::constexpr_for_1<n_nodes>([&]<int a>() {

lib/mito/fem/blocks/L2NormBlock.h

@@ -51,8 +51,8 @@ namespace mito::fem::blocks {
                 constexpr auto w =
                     element_type::canonical_element_type::area * quadrature_rule.weight(q);
 
-                // precompute the common factor
-                auto factor = w * tensor::determinant(element.jacobian()(xi));
+                // precompute the common factor using the element's volume element
+                auto factor = w * element.volume_element()(xi);
 
                 // populate the elementary contribution to the matrix
                 elementary_contribution += factor * _function(xi) * _function(xi);

lib/mito/fem/blocks/MassBlock.h

@@ -44,8 +44,8 @@ namespace mito::fem::blocks {
                 constexpr auto w =
                     element_type::canonical_element_type::area * quadrature_rule.weight(q);
 
-                // precompute the common factor
-                auto factor = w * tensor::determinant(element.jacobian()(xi));
+                // precompute the common factor using the element's volume element
+                auto factor = w * element.volume_element()(xi);
 
                 // loop on the nodes of the element
                 tensor::constexpr_for_1<n_nodes>([&]<int a>() {

lib/mito/fem/blocks/SourceTermBlock.h

@@ -56,8 +56,8 @@ namespace mito::fem::blocks {
                 constexpr auto w =
                     element_type::canonical_element_type::area * quadrature_rule.weight(q);
 
-                // precompute the common factor
-                auto factor = w * tensor::determinant(element.jacobian()(xi));
+                // precompute the common factor using the element's volume element
+                auto factor = w * element.volume_element()(xi);
 
                 // loop on the nodes of the element
                 tensor::constexpr_for_1<n_nodes>([&]<int a>() {

lib/mito/fem/elements/IsoparametricSegment.h

@@ -8,34 +8,52 @@
 
 
 // DESIGN NOTES
-// Class {IsoparametricSegment} represents a first order simplex (segment) equipped with parametric
-// coordinates.
+// Class {IsoparametricSegment} represents a 1-simplex equipped with parametric
+// coordinates. Template parameters:
+//   - coordsT: The coordinate type (determines ambient dimension D and coordinate system)
+//   - VolumeFormT: The type of the volume form (N-form for integration)
+// This design supports segments in any coordinate system and any embedding (N=D
+// for non-embedded, N<D for embedded manifolds).
 
 
 namespace mito::fem {
 
+    template <geometry::coordinates_c coordsT, class VolumeFormT>
     class IsoparametricSegment : public utilities::Invalidatable {
       public:
-        // the dimension of the physical space
-        static constexpr int dim = 1;
+        // dimension of the parametric space
+        static constexpr int N = 1;
+        // ambient dimension from coordinate type
+        static constexpr int D = coordsT::dim;
+
+        // the Jacobian is a D×N matrix (for a segment: D×1, a column vector)
+        using jacobian_type = tensor::matrix_t<D, N>;
+
+        // the metric space provides g, g_inv for this coordinate system
+        using metric_space_type = geometry::metric_space<coordsT>;
+
         // the discretization node type
         using discretization_node_type = discrete::discretization_node_t;
         // the underlying cell type
-        using cell_type = geometry::segment_t<dim>;
+        using cell_type = geometry::segment_t<D>;
 
       protected:
-        // cartesian coordinates in 1D
-        using coordinates_type = geometry::coordinates_t<1, geometry::CARTESIAN>;
-        // the coordinate system type
+        // coordinate types
+        using coordinates_type = coordsT;
         using coordinate_system_type = geometry::coordinate_system_t<coordinates_type>;
         // the vector type
-        using vector_type = tensor::vector_t<1>;
+        using vector_type = tensor::vector_t<D>;
+        // volume form type
+        using volume_form_type = VolumeFormT;
 
       public:
         // the default constructor
         constexpr IsoparametricSegment(
-            const cell_type & cell, const coordinate_system_type & coord_system) :
+            const cell_type & cell, const coordinate_system_type & coord_system,
+            const volume_form_type & volume_form) :
             _cell(cell),
+            _coord_system(coord_system),
+            _volume_form(volume_form),
             _x0{ coord_system.coordinates(cell.nodes()[0]->point()) - coordinates_type{} },
             _x1{ coord_system.coordinates(cell.nodes()[1]->point()) - coordinates_type{} }
         {}
@@ -59,11 +77,26 @@ namespace mito::fem {
         // get the geometric simplex
         constexpr auto cell() const noexcept -> const cell_type & { return _cell; }
 
+      protected:
+        // access to metric tensor at a point (for gradient computation)
+        static constexpr auto g(const coordinates_type & x) { return metric_space_type::g(x); }
+
+        static constexpr auto g_inv(const coordinates_type & x)
+        {
+            return metric_space_type::g_inv(x);
+        }
+
       protected:
         // QUESTION: do we need to maintain a reference to the geometric simplex?
         // a const reference to the geometric simplex
         const cell_type & _cell;
 
+        // a const reference to the coordinate system
+        const coordinate_system_type & _coord_system;
+
+        // the volume form (received from manifold/discretizer)
+        const volume_form_type & _volume_form;
+
         // the coordinates of the discretization nodes of the segment
         const vector_type _x0;
         const vector_type _x1;

lib/mito/fem/elements/IsoparametricTriangle.h

@@ -8,34 +8,52 @@
 
 
 // DESIGN NOTES
-// Class {IsoparametricTriangle} represents a second order simplex equipped barycentric coordinates.
+// Class {IsoparametricTriangle} represents a 2-simplex equipped with barycentric
+// coordinates. Template parameters:
+//   - coordsT: The coordinate type (determines ambient dimension D and coordinate system)
+//   - VolumeFormT: The type of the volume form (N-form for integration)
+// This design supports triangles in any coordinate system and any
+// embedding (N=D for non-embedded, N<D for embedded manifolds).
 
 
 namespace mito::fem {
 
+    template <geometry::coordinates_c coordsT, class VolumeFormT>
     class IsoparametricTriangle : public utilities::Invalidatable {
       public:
-        // the dimension of the physical space
-        static constexpr int dim = 2;
+        // dimension of the parametric space
+        static constexpr int N = 2;
+        // ambient dimension from coordinate type
+        static constexpr int D = coordsT::dim;
+
+        // the Jacobian is a D×N matrix (for a triangle: D×2)
+        using jacobian_type = tensor::matrix_t<D, N>;
+
+        // the metric space provides g, g_inv for this coordinate system
+        using metric_space_type = geometry::metric_space<coordsT>;
+
         // the discretization node type
         using discretization_node_type = discrete::discretization_node_t;
         // the underlying cell type
-        using cell_type = geometry::triangle_t<dim>;
+        using cell_type = geometry::triangle_t<D>;
 
       protected:
-        // cartesian coordinates in 2D
-        using coordinates_type = geometry::coordinates_t<2, geometry::CARTESIAN>;
-        // the coordinate system type
+        // coordinate types
+        using coordinates_type = coordsT;
         using coordinate_system_type = geometry::coordinate_system_t<coordinates_type>;
         // the vector type
-        using vector_type = tensor::vector_t<2>;
+        using vector_type = tensor::vector_t<D>;
+        // volume form type
+        using volume_form_type = VolumeFormT;
 
       public:
         // the default constructor
         constexpr IsoparametricTriangle(
-            const cell_type & cell, const coordinate_system_type & coord_system) :
+            const cell_type & cell, const coordinate_system_type & coord_system,
+            const volume_form_type & volume_form) :
             _cell(cell),
             _coord_system(coord_system),
+            _volume_form(volume_form),
             _x0{ coord_system.coordinates(cell.nodes()[0]->point()) - coordinates_type{} },
             _x1{ coord_system.coordinates(cell.nodes()[1]->point()) - coordinates_type{} },
             _x2{ coord_system.coordinates(cell.nodes()[2]->point()) - coordinates_type{} }
@@ -63,13 +81,25 @@ namespace mito::fem {
         // get the mapping from parametric coordinates to physical coordinates
         constexpr auto parametrization() const { return _cell.parametrization(_coord_system); }
 
+      protected:
+        // access to metric tensor at a point (for gradient computation)
+        static constexpr auto g(const coordinates_type & x) { return metric_space_type::g(x); }
+
+        static constexpr auto g_inv(const coordinates_type & x)
+        {
+            return metric_space_type::g_inv(x);
+        }
+
       protected:
         // a const reference to the geometric simplex
         const cell_type & _cell;
 
         // a const reference to the coordinate system
         const coordinate_system_type & _coord_system;
 
+        // the volume form (received from manifold/discretizer)
+        const volume_form_type & _volume_form;
+
         // the coordinates of the discretization nodes of the triangle
         const vector_type _x0;
         const vector_type _x1;

lib/mito/fem/elements/api.h

@@ -8,11 +8,7 @@
 
 
 namespace mito::fem {
-
-    // type alias for convenient access to the isoparametric simplex type
-    template <int degree, geometry::geometric_simplex_c geometricSimplexT>
-    using isoparametric_simplex_t = typename isoparametric_simplex<degree, geometricSimplexT>::type;
-
+    // no specializations needed here - see elements_library.h for the unified API
 }
 
 

lib/mito/fem/elements/elements_library.h

@@ -7,9 +7,64 @@
 #pragma once
 
 
+// DESIGN NOTES:
+// This file provides the main user-facing API for selecting finite elements.
+// Users specify only the polynomial degree and manifold type:
+//
+//   using element_t = mito::fem::isoparametric_simplex_t<degree, ManifoldT>;
+//
+// The type system automatically extracts:
+//   - Coordinate type (coordsT) from ManifoldT::coordinates_type
+//   - Volume form type (VolumeFormT) from ManifoldT::volume_form_type
+//   - Cell dimensionality to select segment vs triangle
+
+
 #include "seg1/public.h"
 #include "tri1/public.h"
 #include "tri2/public.h"
 
 
+namespace mito::fem {
+
+    template <int degree, class ManifoldT>
+    struct isoparametric_simplex;
+
+    // specialization for degree 1, 1D manifolds (segments)
+    template <class ManifoldT>
+    requires(ManifoldT::cell_type::order == 1)    // 1D cell (segment)
+    struct isoparametric_simplex<1, ManifoldT> {
+        using coordinates_type = typename ManifoldT::coordinates_type;
+        using volume_form_type = typename ManifoldT::volume_form_type;
+
+        using type = IsoparametricSegmentP1<coordinates_type, volume_form_type>;
+    };
+
+    // specialization for degree 1, 2D manifolds (triangles)
+    template <class ManifoldT>
+    requires(ManifoldT::cell_type::order == 2)    // 2D cell (triangle)
+    struct isoparametric_simplex<1, ManifoldT> {
+        using coordinates_type = typename ManifoldT::coordinates_type;
+        using volume_form_type = typename ManifoldT::volume_form_type;
+
+        using type = IsoparametricTriangleP1<coordinates_type, volume_form_type>;
+    };
+
+    // specialization for degree 2, 2D manifolds (triangles)
+    template <class ManifoldT>
+    requires(ManifoldT::cell_type::order == 2)    // 2D cell (triangle)
+    struct isoparametric_simplex<2, ManifoldT> {
+        using coordinates_type = typename ManifoldT::coordinates_type;
+        using volume_form_type = typename ManifoldT::volume_form_type;
+
+        using type = IsoparametricTriangleP2<coordinates_type, volume_form_type>;
+    };
+
+    // convenience alias - primary way users specify element types
+    // usage: isoparametric_simplex_t<degree, decltype(manifold)>
+    template <int degree, class ManifoldT>
+    using isoparametric_simplex_t = typename isoparametric_simplex<degree, ManifoldT>::type;
+
+}    // namespace mito::fem
+
+
 // end of file

lib/mito/fem/elements/forward.h

@@ -13,9 +13,10 @@ namespace mito::fem {
     template <class elementT, discretization_t discretizationT>
     struct Discretizer;
 
-    // struct storing the type of an isoparametric simplex of polynomial degree {degree} on a
-    // geometric simplex of type {geometricSimplexT}
-    template <int degree, geometry::geometric_simplex_c geometricSimplexT>
+    // struct storing the type of an isoparametric simplex of polynomial degree {degree}
+    // on a manifold of type {ManifoldT}
+    // the manifold type encapsulates the coordinate system, volume form, and cell type
+    template <int degree, class ManifoldT>
     struct isoparametric_simplex;
 
 }