mitadel · Pawel024 · Jan 15, 2026 · Jan 17, 2026 · Jan 20, 2026 · Jan 20, 2026
diff --git a/.cmake/mito_tests_mito_lib.cmake b/.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,15 +57,21 @@ 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/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/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)

diff --git a/.github/workflows/cmake-macos.yaml b/.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

diff --git a/.github/workflows/cmake-ubuntu.yaml b/.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

diff --git a/lib/mito/fem/blocks/GradGradBlock.h b/lib/mito/fem/blocks/GradGradBlock.h
@@ -50,7 +50,7 @@ namespace mito::fem::blocks {
                 constexpr auto w = measure * quadrature_rule.weight(q);
 
                 // precompute the common factor
-                auto factor = w * tensor::determinant(element.jacobian()(xi));
+                auto factor = w * element.volume_element()(xi);
 
                 // loop on the nodes of the element
                 tensor::constexpr_for_1<n_nodes>([&]<int a>() {

diff --git a/lib/mito/fem/blocks/L2NormBlock.h b/lib/mito/fem/blocks/L2NormBlock.h
@@ -57,7 +57,7 @@ namespace mito::fem::blocks {
                 constexpr auto w = measure * quadrature_rule.weight(q);
 
                 // precompute the common factor
-                auto factor = w * tensor::determinant(element.jacobian()(xi));
+                auto factor = w * element.volume_element()(xi);
 
                 // populate the elementary contribution to the matrix
                 norm += factor * _function(xi) * _function(xi);

diff --git a/lib/mito/fem/blocks/MassBlock.h b/lib/mito/fem/blocks/MassBlock.h
@@ -50,7 +50,7 @@ namespace mito::fem::blocks {
                 constexpr auto w = measure * quadrature_rule.weight(q);
 
                 // precompute the common factor
-                auto factor = w * tensor::determinant(element.jacobian()(xi));
+                auto factor = w * element.volume_element()(xi);
 
                 // loop on the nodes of the element
                 tensor::constexpr_for_1<n_nodes>([&]<int a>() {

diff --git a/lib/mito/fem/blocks/SourceTermBlock.h b/lib/mito/fem/blocks/SourceTermBlock.h
@@ -62,7 +62,7 @@ namespace mito::fem::blocks {
                 constexpr auto w = measure * quadrature_rule.weight(q);
 
                 // precompute the common factor
-                auto factor = w * tensor::determinant(element.jacobian()(xi));
+                auto factor = w * element.volume_element()(xi);
 
                 // loop on the nodes of the element
                 tensor::constexpr_for_1<n_nodes>([&]<int a>() {

diff --git a/lib/mito/fem/elements/forward.h b/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;
 
     // concept of {parametrizedElementT} being compatible with finite elements of type

diff --git a/lib/mito/fem/elements/seg1/DiscretizerCG.h b/lib/mito/fem/elements/seg1/DiscretizerCG.h
@@ -9,9 +9,9 @@
 
 namespace mito::fem {
 
-    // discretizer specialization for a segment in 1D with order 1 continuous Galerkin
-    template <>
-    struct Discretizer<finite_element_family<geometry::segment_t<1>, 1>, discretization_t::CG> {
+    // discretizer specialization for a segment in D dimensions with order 1 continuous Galerkin
+    template <int D>
+    struct Discretizer<finite_element_family<geometry::segment_t<D>, 1>, discretization_t::CG> {
         template <
             typename manifoldT, typename constraintsT, typename connectivity_table_type,
             typename map_type, typename constrained_nodes_type>
@@ -21,7 +21,7 @@ namespace mito::fem {
             constrained_nodes_type & constrained_nodes)
         {
             // the finite element type
-            using finite_element_type = finite_element_family<geometry::segment_t<1>, 1>;
+            using finite_element_type = finite_element_family<geometry::segment_t<D>, 1>;
 
             // the discretization node type
             using discretization_node_type = typename finite_element_type::discretization_node_type;

diff --git a/lib/mito/fem/elements/seg1/IsoparametricSegmentP1.h b/lib/mito/fem/elements/seg1/IsoparametricSegmentP1.h
@@ -8,8 +8,8 @@
 
 
 // DESIGN NOTES
-// Class {IsoparametricSegmentP1} represents a first order simplex (segment) living in 1D cartesian
-// space, equipped with linear shape functions defined in the parametric space.
+// Class {IsoparametricSegmentP1} represents a 1-simplex (segment) embedded in a D-dimensional
+// physical space, equipped with linear shape functions defined in the parametric space.
 
 
 namespace mito::fem {
@@ -22,6 +22,14 @@ namespace mito::fem {
         using parametrized_element_type = parametrizedElementT;
         // the underlying mesh cell type
         using mesh_cell_type = typename parametrized_element_type::cell_type;
+        // the metric volume form type
+        using metric_volume_form_type = typename parametrized_element_type::metric_volume_form_type;
+        // the coordinates type of the physical space
+        using coordinates_type = typename metric_volume_form_type::input_type;
+        // the dimension of the physical space
+        static constexpr int dim = coordinates_type::dim;
+        // the (euclidean) metric of the physical space in {coordinates_type} coordinates
+        using ambient_metric_type = geometry::euclidean_metric<coordinates_type>;
 
         // the degree of the finite element
         static constexpr int degree = 1;
@@ -99,43 +107,53 @@ namespace mito::fem {
         {
             // assemble the jacobian as a function of parametric coordinates
             auto jacobian_function = functions::function(
-                [&](const parametric_coordinates_type & xi) -> tensor::matrix_t<1> {
+                [&](const parametric_coordinates_type & xi) -> tensor::matrix_t<dim, 1> {
                     // get the shape functions derivatives
                     constexpr auto dphi_0 = shape_functions.dshape<0>();
                     constexpr auto dphi_1 = shape_functions.dshape<1>();
 
-                    // store the coordinates of the vertices of the triangle in physical space
+                    // store the coordinates of the vertices of the segment in physical space
                     auto x0 = _element.parametrization()({ 0.0 });
                     auto x1 = _element.parametrization()({ 1.0 });
 
                     // compute the jacobian of the isoparametric mapping: dx/dxi
-                    auto dx_dxi = x0 * dphi_0(xi) + x1 * dphi_1(xi);
-                    // wrap the result in a 1x1 matrix
-                    return tensor::matrix_t<1>{ dx_dxi };
+                    return tensor::dyadic(x0, dphi_0(xi)) + tensor::dyadic(x1, dphi_1(xi));
                 });
 
             // and return it
             return jacobian_function;
         }
 
+        // get the volume element of the isoparametric mapping, i.e. the metric volume form
+        // contracted with the tangent vector of the element
+        constexpr auto volume_element() const
+        {
+            // assemble the volume element as a function of parametric coordinates
+            auto volume_element_function = functions::function(
+                [&](const parametric_coordinates_type & xi) -> tensor::scalar_t {
+                    // the metric volume form at the physical location of {xi}
+                    auto w = _element.metric_volume_form()(_element.parametrization())(xi);
+                    // the tangent vector to the element at {xi}
+                    auto tangent = tensor::column<0>(jacobian()(xi));
+                    // contract the metric volume form with the tangent vector
+                    return w(tangent);
+                });
+            // and return it
+            return volume_element_function;
+        }
+
         // get the gradient of the a-th shape function as a function of parametric coordinates
         template <int a>
         requires(a >= 0 && a < n_nodes)
         constexpr auto gradient() const
         {
-            // assemble the gradient as a function of parametric coordinates
-            auto gradient_function = functions::function(
-                [&](const parametric_coordinates_type & xi) -> tensor::vector_t<1> {
-                    // the jacobian of the mapping from the reference element to the physical
-                    // element evaluated at {xi}
-                    auto J = jacobian()(xi);
-                    // the derivative of the coordinates with respect to the parametric coordinates
-                    auto J_inv = tensor::inverse(J);
-                    // return the spatial gradients of the shape functions evaluated at {xi}
-                    return shape_functions.dshape<a>()(xi) * J_inv;
-                });
-            // and return it
-            return gradient_function;
+            // the pullback of the ambient metric along the element parametrization
+            auto g =
+                geometry::pullback_metric<ambient_metric_type>::field(_element.parametrization());
+            // the intrinsic gradient of the a-th shape function on the parametric space
+            auto dphi = operators::gradient(shape_functions.shape<a>(), g);
+            // push the gradient forward to the physical space with the jacobian
+            return jacobian() * dphi;
         }
 
       private:

diff --git a/lib/mito/fem/elements/seg1/api.h b/lib/mito/fem/elements/seg1/api.h
@@ -10,11 +10,11 @@
 namespace mito::fem {
 
     // specialization of {finite_element_family} for first order segments
-    template <>
-    struct finite_element_family<geometry::segment_t<1>, 1> {
+    template <int D>
+    struct finite_element_family<geometry::segment_t<D>, 1> {
 
         // the dimension of the physical space
-        static constexpr int dim = 1;
+        static constexpr int dim = D;
         // the underlying mesh cell type
         using mesh_cell_type = geometry::segment_t<dim>;
         // the degree of the finite element

diff --git a/lib/mito/fem/elements/tri1/DiscretizerCG.h b/lib/mito/fem/elements/tri1/DiscretizerCG.h
@@ -10,8 +10,8 @@
 namespace mito::fem {
 
     // discretizer specialization for {IsoparametricTriangleP1} with continuous Galerkin
-    template <>
-    struct Discretizer<finite_element_family<geometry::triangle_t<2>, 1>, discretization_t::CG> {
+    template <int D>
+    struct Discretizer<finite_element_family<geometry::triangle_t<D>, 1>, discretization_t::CG> {
         template <
             typename manifoldT, typename constraintsT, typename connectivity_table_type,
             typename map_type, typename constrained_nodes_type>
@@ -21,7 +21,7 @@ namespace mito::fem {
             constrained_nodes_type & constrained_nodes)
         {
             // the finite element type
-            using finite_element_type = finite_element_family<geometry::triangle_t<2>, 1>;
+            using finite_element_type = finite_element_family<geometry::triangle_t<D>, 1>;
 
             // the discretization node type
             using discretization_node_type = typename finite_element_type::discretization_node_type;

diff --git a/lib/mito/fem/elements/tri1/IsoparametricTriangleP1.h b/lib/mito/fem/elements/tri1/IsoparametricTriangleP1.h
@@ -8,8 +8,8 @@
 
 
 // DESIGN NOTES
-// Class {IsoparametricTriangleP1} represents a second order simplex (triangle) living in 2D
-// cartesian space, equipped with linear shape functions defined in the parametric space.
+// Class {IsoparametricTriangleP1} represents a 2-simplex (triangle) embedded in a D-dimensional
+// physical space, equipped with linear shape functions defined in the parametric space.
 
 
 namespace mito::fem {
@@ -22,6 +22,14 @@ namespace mito::fem {
         using parametrized_element_type = parametrizedElementT;
         // the underlying mesh cell type
         using mesh_cell_type = typename parametrized_element_type::cell_type;
+        // the metric volume form type
+        using metric_volume_form_type = typename parametrized_element_type::metric_volume_form_type;
+        // the coordinates type of the physical space
+        using coordinates_type = typename metric_volume_form_type::input_type;
+        // the dimension of the physical space
+        static constexpr int dim = coordinates_type::dim;
+        // the (euclidean) metric of the physical space in {coordinates_type} coordinates
+        using ambient_metric_type = geometry::euclidean_metric<coordinates_type>;
 
         // the degree of the finite element
         static constexpr int degree = 1;
@@ -99,7 +107,7 @@ namespace mito::fem {
         {
             // assemble the jacobian as a function of parametric coordinates
             auto jacobian_function = functions::function(
-                [&](const parametric_coordinates_type & xi) -> tensor::matrix_t<2> {
+                [&](const parametric_coordinates_type & xi) -> tensor::matrix_t<dim, 2> {
                     // get the shape functions derivatives
                     constexpr auto dphi_0 = shape_functions.dshape<0>();
                     constexpr auto dphi_1 = shape_functions.dshape<1>();
@@ -120,25 +128,36 @@ namespace mito::fem {
             return jacobian_function;
         }
 
+        // get the volume element of the isoparametric mapping, i.e. the metric volume form
+        // contracted with the tangent vectors of the element
+        constexpr auto volume_element() const
+        {
+            // assemble the volume element as a function of parametric coordinates
+            auto volume_element_function = functions::function(
+                [&](const parametric_coordinates_type & xi) -> tensor::scalar_t {
+                    // the metric volume form at the physical location of {xi}
+                    auto w = _element.metric_volume_form()(_element.parametrization())(xi);
+                    // the tangent vectors to the element at {xi}
+                    auto J = jacobian()(xi);
+                    // contract the metric volume form with the tangent vectors
+                    return w(tensor::column<0>(J), tensor::column<1>(J));
+                });
+            // and return it
+            return volume_element_function;
+        }
+
         // get the gradient of the a-th shape function as a function of parametric coordinates
         template <int a>
         requires(a >= 0 && a < n_nodes)
         constexpr auto gradient() const
         {
-            // assemble the gradient as a function of parametric coordinates
-            auto gradient_function = functions::function(
-                [&](const parametric_coordinates_type & xi) -> tensor::vector_t<2> {
-                    // the jacobian of the mapping from the reference element to the physical
-                    // element evaluated at {xi}
-                    auto J = jacobian()(xi);
-                    // the derivative of the coordinates with respect to the parametric coordinates
-                    auto J_inv = tensor::inverse(J);
-                    // return the spatial gradients of the shape functions evaluated at {xi}
-                    return shape_functions.dshape<a>()(xi) * J_inv;
-                });
-
-            // and return it
-            return gradient_function;
+            // the pullback of the ambient metric along the element parametrization
+            auto g =
+                geometry::pullback_metric<ambient_metric_type>::field(_element.parametrization());
+            // the intrinsic gradient of the a-th shape function on the parametric space
+            auto dphi = operators::gradient(shape_functions.shape<a>(), g);
+            // push the gradient forward to the physical space with the jacobian
+            return jacobian() * dphi;
         }
 
       private:

diff --git a/lib/mito/fem/elements/tri1/api.h b/lib/mito/fem/elements/tri1/api.h
@@ -9,11 +9,12 @@
 
 namespace mito::fem {
 
-    template <>
-    struct finite_element_family<geometry::triangle_t<2>, 1> {
+    // specialization of {finite_element_family} for first order triangles
+    template <int D>
+    struct finite_element_family<geometry::triangle_t<D>, 1> {
 
         // the dimension of the physical space
-        static constexpr int dim = 2;
+        static constexpr int dim = D;
         // the underlying mesh cell type
         using mesh_cell_type = geometry::triangle_t<dim>;
         // the degree of the finite element

diff --git a/lib/mito/fem/elements/tri2/DiscretizerCG.h b/lib/mito/fem/elements/tri2/DiscretizerCG.h
@@ -10,8 +10,8 @@
 namespace mito::fem {
 
     // discretizer specialization for {IsoparametricTriangleP2} with continuous Galerkin
-    template <>
-    struct Discretizer<finite_element_family<geometry::triangle_t<2>, 2>, discretization_t::CG> {
+    template <int D>
+    struct Discretizer<finite_element_family<geometry::triangle_t<D>, 2>, discretization_t::CG> {
         template <
             typename manifoldT, typename constraintsT, typename connectivity_table_type,
             typename map_type, typename constrained_nodes_type>
@@ -21,7 +21,7 @@ namespace mito::fem {
             constrained_nodes_type & constrained_nodes)
         {
             // the finite element type
-            using finite_element_type = finite_element_family<geometry::triangle_t<2>, 2>;
+            using finite_element_type = finite_element_family<geometry::triangle_t<D>, 2>;
 
             // assemble the mesh node type
             using mesh_node_type = typename manifoldT::mesh_type::node_type;