diff --git a/.cmake/mito_tests_mito_lib.cmake b/.cmake/mito_tests_mito_lib.cmake
index f872aade..3f1c30d2 100644
--- 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,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/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)
@@ -65,6 +70,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)
diff --git a/.github/workflows/cmake-macos.yaml b/.github/workflows/cmake-macos.yaml
index 5ea6c67b..6d7e92f2 100644
--- 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
index 448a6e81..38edd138 100644
--- 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
index 7bcbcd1e..9123a0cf 100644
--- 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
index 4e0be6dc..7ad2109d 100644
--- 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
index b57199b0..91f2546f 100644
--- 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
index 9c9c565e..a7362f74 100644
--- 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
index 9ecbc31a..4fe226be 100644
--- 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
index c77b391f..de49226e 100644
--- 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
index 09760dcf..fe3ec073 100644
--- 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
index b61c5c27..4894e473 100644
--- 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
index 92782591..d33a9530 100644
--- 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
index c6c56e35..7e1d1999 100644
--- 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
index d2efc5f9..688b638b 100644
--- 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
index b6322be3..d14e1bf8 100644
--- 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;
diff --git a/lib/mito/fem/elements/tri2/IsoparametricTriangleP2.h b/lib/mito/fem/elements/tri2/IsoparametricTriangleP2.h
index 5f6e41aa..e8e2b9cf 100644
--- a/lib/mito/fem/elements/tri2/IsoparametricTriangleP2.h
+++ b/lib/mito/fem/elements/tri2/IsoparametricTriangleP2.h
@@ -8,8 +8,8 @@
 
 
 // DESIGN NOTES
-// Class {IsoparametricTriangleP2} represents a second order simplex (triangle) living in 2D
-// cartesian space, equipped with quadratic shape functions defined in the parametric space.
+// Class {IsoparametricTriangleP2} represents a 2-simplex (triangle) embedded in a D-dimensional
+// physical space, equipped with quadratic 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 = 2;
@@ -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> {
                     // store the coordinates of the vertices of the triangle in physical space
                     auto x0 = _element.parametrization()({ 0.0, 0.0 });
                     auto x1 = _element.parametrization()({ 1.0, 0.0 });
@@ -127,24 +135,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/tri2/api.h b/lib/mito/fem/elements/tri2/api.h
index 396b5a08..02a8e13e 100644
--- a/lib/mito/fem/elements/tri2/api.h
+++ b/lib/mito/fem/elements/tri2/api.h
@@ -9,11 +9,12 @@
 
 namespace mito::fem {
 
-    template <>
-    struct finite_element_family<geometry::triangle_t<2>, 2> {
+    // specialization of {finite_element_family} for second order triangles
+    template <int D>
+    struct finite_element_family<geometry::triangle_t<D>, 2> {
 
         // 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/operators/differential.h b/lib/mito/operators/differential.h
index 063957bf..d89645bb 100644
--- a/lib/mito/operators/differential.h
+++ b/lib/mito/operators/differential.h
@@ -64,6 +64,23 @@ namespace mito::operators {
         return _grad(field, std::make_index_sequence<D>{});
     }
 
+    // function to compute the contravariant gradient of a scalar field with respect to a metric
+    template <fields::scalar_field_c F, class G>
+    requires(
+        (fields::tensor_field_c<G> or fields::scalar_field_c<G>)
+        and fields::compatible_fields_c<F, G>)
+    constexpr auto gradient(const F & field, const G & metric_field)
+    {
+        // if the metric collapses to a scalar
+        if constexpr (fields::scalar_field_c<G>) {
+            // divide by the metric directly
+            return gradient(field) / metric_field;
+        } else {
+            // multiply with the inverse of the metric
+            return functions::inverse(metric_field) * gradient(field);
+        }
+    }
+
     // function to compute the divergence of a vector field
     template <fields::vector_field_c F>
     constexpr auto divergence(const F & field)
diff --git a/tests/mito.lib/fem/block_grad_grad_embedded_segment.cc b/tests/mito.lib/fem/block_grad_grad_embedded_segment.cc
new file mode 100644
index 00000000..aed050bc
--- /dev/null
+++ b/tests/mito.lib/fem/block_grad_grad_embedded_segment.cc
@@ -0,0 +1,97 @@
+// -*- c++ -*-
+//
+// Copyright (c) 2020-2026, the MiTo Authors, all rights reserved
+//
+
+#include <gtest/gtest.h>
+#include <mito.h>
+
+
+// the type of coordinates (2D physical space)
+using coordinates_t = mito::geometry::coordinates_t<2, mito::geometry::CARTESIAN>;
+// the metric space type
+using metric_space_t = mito::geometry::euclidean_metric_space<coordinates_t>;
+// the type of discretization node
+using discretization_node_t = mito::discrete::discretization_node_t;
+// the type of cell (segment embedded in 2D)
+using cell_t = mito::geometry::segment_t<2>;
+// the reference simplex
+using reference_simplex_t = cell_t::reference_simplex_type;
+// Gauss quadrature on segments with degree of exactness 2
+using quadrature_rule_t =
+    mito::quadrature::quadrature_rule_t<mito::quadrature::GAUSS, reference_simplex_t, 2>;
+
+// instantiate the quadrature rule
+constexpr auto quadrature_rule = quadrature_rule_t();
+
+
+TEST(Fem, BlockGradGradEmbeddedSegment)
+{
+    // the coordinate system
+    auto coord_system = mito::geometry::coordinate_system_t<coordinates_t>();
+
+    // an atlas under the coordinate system
+    auto atlas = mito::manifolds::atlas<cell_t>(coord_system);
+
+    // build nodes (unit-length diagonal segment embedded in 2D)
+    constexpr auto inv_sqrt2 = 1.0 / std::sqrt(2.0);
+    auto node_0 = mito::geometry::node(coord_system, { 0.0, 0.0 });
+    auto node_1 = mito::geometry::node(coord_system, { inv_sqrt2, inv_sqrt2 });
+
+    // make a geometric simplex
+    auto segment = mito::geometry::segment(node_0, node_1);
+
+    // the normal field to the segment (perpendicular to the segment direction)
+    // diagonal direction is (inv_sqrt2, inv_sqrt2), so normal is (inv_sqrt2, -inv_sqrt2)
+    // (rotated 90° clockwise to get positive orientation with w(normal, tangent) > 0)
+    auto normal_field = mito::functions::constant<coordinates_t>(
+        mito::tensor::vector_t<2>{ inv_sqrt2, -inv_sqrt2 });
+
+    // strip the namespace from the placeholder for forms contractions
+    using mito::tensor::_;
+
+    // the ambient metric volume form
+    constexpr auto w = metric_space_t::w;
+
+    // the restriction of the metric volume form to the segment
+    auto wS = mito::functions::function(
+        [w, normal_field](const coordinates_t & x) { return w(x)(normal_field(x), _); });
+
+    // make a manifold element from the segment
+    auto element =
+        mito::manifolds::parametrized_element(segment, atlas.parametrization(segment), wS);
+
+    {
+        // build the discretization nodes
+        auto discretization_node_0 = discretization_node_t();
+        auto discretization_node_1 = discretization_node_t();
+
+        // the degree of the finite element
+        constexpr int degree = 1;
+        // assemble the finite element type
+        using finite_element_t = mito::fem::finite_element_family<cell_t, degree>;
+
+        // a finite element on the embedded segment
+        auto element_p1 = mito::fem::finite_element<finite_element_t>(
+            element, { discretization_node_0, discretization_node_1 });
+
+        // a grad-grad matrix block
+        auto grad_grad_block =
+            mito::fem::blocks::grad_grad_block<finite_element_t, quadrature_rule_t>();
+
+        // the analytical elementary stiffness matrix (same as 1D for unit-length segment)
+        auto analytical_block = mito::tensor::matrix_t<2>{ 1.0, -1.0, -1.0, 1.0 };
+
+        // compute the elementary contribution of the block
+        auto computed_block = grad_grad_block.compute(element_p1);
+
+        // compute the error
+        auto error = mito::tensor::norm(computed_block - analytical_block);
+
+        // check the error is zero to machine precision
+        EXPECT_DOUBLE_EQ(0.0, error);
+    }
+
+    // all done
+    return;
+}
diff --git a/tests/mito.lib/fem/block_grad_grad_embedded_triangle.cc b/tests/mito.lib/fem/block_grad_grad_embedded_triangle.cc
new file mode 100644
index 00000000..a96298f3
--- /dev/null
+++ b/tests/mito.lib/fem/block_grad_grad_embedded_triangle.cc
@@ -0,0 +1,100 @@
+// -*- c++ -*-
+//
+// Copyright (c) 2020-2026, the MiTo Authors, all rights reserved
+//
+
+#include <gtest/gtest.h>
+#include <mito.h>
+
+
+// the type of coordinates (3D physical space)
+using coordinates_t = mito::geometry::coordinates_t<3, mito::geometry::CARTESIAN>;
+// the metric space type
+using metric_space_t = mito::geometry::euclidean_metric_space<coordinates_t>;
+// the type of discretization node
+using discretization_node_t = mito::discrete::discretization_node_t;
+// the type of cell (triangle embedded in 3D)
+using cell_t = mito::geometry::triangle_t<3>;
+// the reference simplex
+using reference_simplex_t = cell_t::reference_simplex_type;
+// Gauss quadrature on triangles with degree of exactness 2
+using quadrature_rule_t =
+    mito::quadrature::quadrature_rule_t<mito::quadrature::GAUSS, reference_simplex_t, 2>;
+
+// instantiate the quadrature rule
+constexpr auto quadrature_rule = quadrature_rule_t();
+
+
+TEST(Fem, BlockGradGradEmbeddedTriangle)
+{
+    // the coordinate system
+    auto coord_system = mito::geometry::coordinate_system_t<coordinates_t>();
+
+    // an atlas under the coordinate system
+    auto atlas = mito::manifolds::atlas<cell_t>(coord_system);
+
+    // build nodes for a unit right triangle rigidly rotated out of the xy-plane
+    // (the legs {t_1} = (0.6, 0.0, 0.8) and {t_2} = (0.0, 1.0, 0.0) are unit and orthogonal)
+    auto node_0 = mito::geometry::node(coord_system, { 0.0, 0.0, 0.0 });
+    auto node_1 = mito::geometry::node(coord_system, { 0.6, 0.0, 0.8 });
+    auto node_2 = mito::geometry::node(coord_system, { 0.0, 1.0, 0.0 });
+
+    // make a geometric simplex
+    auto triangle = mito::geometry::triangle(node_0, node_1, node_2);
+
+    // the normal field to the triangle: {n} = {t_1} x {t_2} normalized,
+    // oriented such that w(n, t_1, t_2) > 0
+    auto normal_field =
+        mito::functions::constant<coordinates_t>(mito::tensor::vector_t<3>{ -0.8, 0.0, 0.6 });
+
+    // strip the namespace from the placeholder for forms contractions
+    using mito::tensor::_;
+
+    // the ambient metric volume form
+    constexpr auto w = metric_space_t::w;
+
+    // the restriction of the metric volume form to the triangle
+    auto wS = mito::functions::function(
+        [w, normal_field](const coordinates_t & x) { return w(x)(normal_field(x), _, _); });
+
+    // make a manifold element from the triangle
+    auto element =
+        mito::manifolds::parametrized_element(triangle, atlas.parametrization(triangle), wS);
+
+    {
+        // build the discretization nodes
+        auto discretization_node_0 = discretization_node_t();
+        auto discretization_node_1 = discretization_node_t();
+        auto discretization_node_2 = discretization_node_t();
+
+        // the degree of the finite element
+        constexpr int degree = 1;
+        // assemble the finite element type
+        using finite_element_t = mito::fem::finite_element_family<cell_t, degree>;
+
+        // a finite element on the embedded triangle
+        auto element_p1 = mito::fem::finite_element<finite_element_t>(
+            element, { discretization_node_0, discretization_node_1, discretization_node_2 });
+
+        // a grad-grad matrix block
+        auto grad_grad_block =
+            mito::fem::blocks::grad_grad_block<finite_element_t, quadrature_rule_t>();
+
+        // the analytical elementary stiffness matrix (the stiffness matrix is invariant under
+        // rigid rotations, so it is the same as for the unit right triangle in 2D)
+        auto analytical_block = 1.0 / 2.0 * mito::tensor::matrix_t<3>{ 2.0, -1.0, -1.0, -1.0, 1.0,
+                                                                       0.0, -1.0, 0.0,  1.0 };
+
+        // compute the elementary contribution of the block
+        auto computed_block = grad_grad_block.compute(element_p1);
+
+        // compute the error
+        auto error = mito::tensor::norm(computed_block - analytical_block);
+
+        // check the error is zero to machine precision
+        EXPECT_NEAR(0.0, error, 1.0e-15);
+    }
+
+    // all done
+    return;
+}
diff --git a/tests/mito.lib/fem/block_mass_embedded_segment.cc b/tests/mito.lib/fem/block_mass_embedded_segment.cc
new file mode 100644
index 00000000..4c846b40
--- /dev/null
+++ b/tests/mito.lib/fem/block_mass_embedded_segment.cc
@@ -0,0 +1,96 @@
+// -*- c++ -*-
+//
+// Copyright (c) 2020-2026, the MiTo Authors, all rights reserved
+//
+
+#include <gtest/gtest.h>
+#include <mito.h>
+
+
+// the type of coordinates (2D physical space)
+using coordinates_t = mito::geometry::coordinates_t<2, mito::geometry::CARTESIAN>;
+// the metric space type
+using metric_space_t = mito::geometry::euclidean_metric_space<coordinates_t>;
+// the type of discretization node
+using discretization_node_t = mito::discrete::discretization_node_t;
+// the type of cell (segment embedded in 2D)
+using cell_t = mito::geometry::segment_t<2>;
+// the reference simplex
+using reference_simplex_t = cell_t::reference_simplex_type;
+// Gauss quadrature on segments with degree of exactness 2
+using quadrature_rule_t =
+    mito::quadrature::quadrature_rule_t<mito::quadrature::GAUSS, reference_simplex_t, 2>;
+
+// instantiate the quadrature rule
+constexpr auto quadrature_rule = quadrature_rule_t();
+
+
+TEST(Fem, BlockMassEmbeddedSegment)
+{
+    // the coordinate system
+    auto coord_system = mito::geometry::coordinate_system_t<coordinates_t>();
+
+    // an atlas under the coordinate system
+    auto atlas = mito::manifolds::atlas<cell_t>(coord_system);
+
+    // build nodes (unit-length diagonal segment embedded in 2D)
+    constexpr auto inv_sqrt2 = 1.0 / std::sqrt(2.0);
+    auto node_0 = mito::geometry::node(coord_system, { 0.0, 0.0 });
+    auto node_1 = mito::geometry::node(coord_system, { inv_sqrt2, inv_sqrt2 });
+
+    // make a geometric simplex
+    auto segment = mito::geometry::segment(node_0, node_1);
+
+    // the normal field to the segment (perpendicular to the segment direction)
+    // diagonal direction is (inv_sqrt2, inv_sqrt2), so normal is (inv_sqrt2, -inv_sqrt2)
+    // (rotated 90° clockwise to get positive orientation with w(normal, tangent) > 0)
+    auto normal_field = mito::functions::constant<coordinates_t>(
+        mito::tensor::vector_t<2>{ inv_sqrt2, -inv_sqrt2 });
+
+    // strip the namespace from the placeholder for forms contractions
+    using mito::tensor::_;
+
+    // the ambient metric volume form
+    constexpr auto w = metric_space_t::w;
+
+    // the restriction of the metric volume form to the segment
+    auto wS = mito::functions::function(
+        [w, normal_field](const coordinates_t & x) { return w(x)(normal_field(x), _); });
+
+    // make a manifold element from the segment
+    auto element =
+        mito::manifolds::parametrized_element(segment, atlas.parametrization(segment), wS);
+
+    {
+        // build the discretization nodes
+        auto discretization_node_0 = discretization_node_t();
+        auto discretization_node_1 = discretization_node_t();
+
+        // the degree of the finite element
+        constexpr int degree = 1;
+        // assemble the finite element type
+        using finite_element_t = mito::fem::finite_element_family<cell_t, degree>;
+
+        // a finite element on the embedded segment
+        auto element_p1 = mito::fem::finite_element<finite_element_t>(
+            element, { discretization_node_0, discretization_node_1 });
+
+        // a mass matrix block
+        auto mass_block = mito::fem::blocks::mass_block<finite_element_t, quadrature_rule_t>();
+
+        // the analytical elementary mass matrix (same as 1D for unit-length segment)
+        auto analytical_block = 1.0 / 6.0 * mito::tensor::matrix_t<2>{ 2.0, 1.0, 1.0, 2.0 };
+
+        // compute the elementary contribution of the block
+        auto computed_block = mass_block.compute(element_p1);
+
+        // compute the error
+        auto error = mito::tensor::norm(computed_block - analytical_block);
+
+        // check the error is reasonable
+        EXPECT_NEAR(0.0, error, 1.5e-16);
+    }
+
+    // all done
+    return;
+}
diff --git a/tests/mito.lib/fem/block_mass_embedded_triangle.cc b/tests/mito.lib/fem/block_mass_embedded_triangle.cc
new file mode 100644
index 00000000..cce42fce
--- /dev/null
+++ b/tests/mito.lib/fem/block_mass_embedded_triangle.cc
@@ -0,0 +1,99 @@
+// -*- c++ -*-
+//
+// Copyright (c) 2020-2026, the MiTo Authors, all rights reserved
+//
+
+#include <gtest/gtest.h>
+#include <mito.h>
+
+
+// the type of coordinates (3D physical space)
+using coordinates_t = mito::geometry::coordinates_t<3, mito::geometry::CARTESIAN>;
+// the metric space type
+using metric_space_t = mito::geometry::euclidean_metric_space<coordinates_t>;
+// the type of discretization node
+using discretization_node_t = mito::discrete::discretization_node_t;
+// the type of cell (triangle embedded in 3D)
+using cell_t = mito::geometry::triangle_t<3>;
+// the reference simplex
+using reference_simplex_t = cell_t::reference_simplex_type;
+// Gauss quadrature on triangles with degree of exactness 2
+using quadrature_rule_t =
+    mito::quadrature::quadrature_rule_t<mito::quadrature::GAUSS, reference_simplex_t, 2>;
+
+// instantiate the quadrature rule
+constexpr auto quadrature_rule = quadrature_rule_t();
+
+
+TEST(Fem, BlockMassEmbeddedTriangle)
+{
+    // the coordinate system
+    auto coord_system = mito::geometry::coordinate_system_t<coordinates_t>();
+
+    // an atlas under the coordinate system
+    auto atlas = mito::manifolds::atlas<cell_t>(coord_system);
+
+    // build nodes for a unit right triangle rigidly rotated out of the xy-plane
+    // (the legs {t_1} = (0.6, 0.0, 0.8) and {t_2} = (0.0, 1.0, 0.0) are unit and orthogonal)
+    auto node_0 = mito::geometry::node(coord_system, { 0.0, 0.0, 0.0 });
+    auto node_1 = mito::geometry::node(coord_system, { 0.6, 0.0, 0.8 });
+    auto node_2 = mito::geometry::node(coord_system, { 0.0, 1.0, 0.0 });
+
+    // make a geometric simplex
+    auto triangle = mito::geometry::triangle(node_0, node_1, node_2);
+
+    // the normal field to the triangle: {n} = {t_1} x {t_2} normalized,
+    // oriented such that w(n, t_1, t_2) > 0
+    auto normal_field =
+        mito::functions::constant<coordinates_t>(mito::tensor::vector_t<3>{ -0.8, 0.0, 0.6 });
+
+    // strip the namespace from the placeholder for forms contractions
+    using mito::tensor::_;
+
+    // the ambient metric volume form
+    constexpr auto w = metric_space_t::w;
+
+    // the restriction of the metric volume form to the triangle
+    auto wS = mito::functions::function(
+        [w, normal_field](const coordinates_t & x) { return w(x)(normal_field(x), _, _); });
+
+    // make a manifold element from the triangle
+    auto element =
+        mito::manifolds::parametrized_element(triangle, atlas.parametrization(triangle), wS);
+
+    {
+        // build the discretization nodes
+        auto discretization_node_0 = discretization_node_t();
+        auto discretization_node_1 = discretization_node_t();
+        auto discretization_node_2 = discretization_node_t();
+
+        // the degree of the finite element
+        constexpr int degree = 1;
+        // assemble the finite element type
+        using finite_element_t = mito::fem::finite_element_family<cell_t, degree>;
+
+        // a finite element on the embedded triangle
+        auto element_p1 = mito::fem::finite_element<finite_element_t>(
+            element, { discretization_node_0, discretization_node_1, discretization_node_2 });
+
+        // a mass matrix block
+        auto mass_block = mito::fem::blocks::mass_block<finite_element_t, quadrature_rule_t>();
+
+        // the analytical elementary mass matrix (the mass matrix is invariant under rigid
+        // rotations, so it is the same as for the unit right triangle in 2D)
+        auto analytical_block =
+            1.0 / 24.0 * mito::tensor::matrix_t<3>{ 2.0, 1.0, 1.0, 1.0, 2.0, 1.0, 1.0, 1.0, 2.0 };
+
+        // compute the elementary contribution of the block
+        auto computed_block = mass_block.compute(element_p1);
+
+        // compute the error
+        auto error = mito::tensor::norm(computed_block - analytical_block);
+
+        // check the error is zero to machine precision
+        EXPECT_NEAR(0.0, error, 1.0e-15);
+    }
+
+    // all done
+    return;
+}
diff --git a/tests/mito.lib/fem/isoparametric_embedded_segment.cc b/tests/mito.lib/fem/isoparametric_embedded_segment.cc
new file mode 100644
index 00000000..38e7872a
--- /dev/null
+++ b/tests/mito.lib/fem/isoparametric_embedded_segment.cc
@@ -0,0 +1,191 @@
+// -*- c++ -*-
+//
+// Copyright (c) 2020-2026, the MiTo Authors, all rights reserved
+//
+
+#include <gtest/gtest.h>
+#include <mito.h>
+
+
+// the type of coordinates (2D physical space)
+using coordinates_t = mito::geometry::coordinates_t<2, mito::geometry::CARTESIAN>;
+// the metric space type
+using metric_space_t = mito::geometry::euclidean_metric_space<coordinates_t>;
+// the type of discretization node
+using discretization_node_t = mito::discrete::discretization_node_t;
+// the type of cell (segment embedded in 2D)
+using cell_t = mito::geometry::segment_t<2>;
+// the reference simplex
+using reference_simplex_t = cell_t::reference_simplex_type;
+// Gauss quadrature on segments with degree of exactness 2
+using quadrature_rule_t =
+    mito::quadrature::quadrature_rule_t<mito::quadrature::GAUSS, reference_simplex_t, 2>;
+
+
+// instantiate the quadrature rule
+constexpr auto quadrature_rule = quadrature_rule_t();
+
+
+// test that all shape functions sum to 1.0 at any quadrature point
+auto
+test_partition_of_unity(const auto & element)
+{
+    // the number of quadrature points per element
+    constexpr int n_quads = quadrature_rule_t::npoints;
+
+    // the number of nodes per element
+    constexpr int n_nodes = mito::utilities::base_type<decltype(element)>::n_nodes;
+
+    // loop on the quadrature points
+    mito::tensor::constexpr_for_1<n_quads>([&]<int q>() {
+        // the parametric coordinates of the quadrature point
+        constexpr auto xi = quadrature_rule.point(q);
+
+        // compute the sum of the shape functions at {xi} for all nodes
+        constexpr auto sum =
+            ([]<int... a>(
+                 const auto & element, const auto & xi, mito::tensor::integer_sequence<a...>) {
+                return ((element.template shape<a>()(xi)) + ...);
+            })(element, xi, mito::tensor::make_integer_sequence<n_nodes>{});
+
+        // check the sum of the shape functions
+        static_assert(1.0 == sum);
+    });
+
+    // all done
+    return;
+}
+
+// test that the gradients of all shape functions sum to 0.0 at any quadrature point
+auto
+test_gradient_consistency(const auto & element)
+{
+    // the number of quadrature points per element
+    constexpr int n_quads = quadrature_rule_t::npoints;
+
+    // the number of nodes per element
+    constexpr int n_nodes = mito::utilities::base_type<decltype(element)>::n_nodes;
+
+    // loop on the quadrature points
+    mito::tensor::constexpr_for_1<n_quads>([&]<int q>() {
+        // the parametric coordinates of the quadrature point
+        constexpr auto xi = quadrature_rule.point(q);
+
+        // compute the sum of the shape functions gradients at {xi} for all nodes
+        auto sum =
+            ([]<int... a>(
+                 const auto & element, const auto & xi, mito::tensor::integer_sequence<a...>) {
+                return ((element.template gradient<a>()(xi)) + ...);
+            })(element, xi, mito::tensor::make_integer_sequence<n_nodes>{});
+
+        // check that the sum of the shape functions gradients is the zero vector
+        EXPECT_NEAR(0.0, mito::tensor::norm(sum), 3.0e-16);
+    });
+
+    // all done
+    return;
+}
+
+// test that the arc length computed via the volume element matches the expected value
+auto
+test_arc_length(const auto & element, double expected_length)
+{
+    // the element type
+    using element_t = mito::utilities::base_type<decltype(element)>;
+
+    // the number of quadrature points per element
+    constexpr int n_quads = quadrature_rule_t::npoints;
+
+    // compute the arc length by integrating 1 over the segment using the volume element
+    double arc_length = 0.0;
+    mito::tensor::constexpr_for_1<n_quads>([&]<int q>() {
+        constexpr auto xi = quadrature_rule.point(q);
+        constexpr auto w = element_t::canonical_element_type::measure * quadrature_rule.weight(q);
+        arc_length += w * element.volume_element()(xi);
+    });
+
+    // check the arc length
+    EXPECT_NEAR(expected_length, arc_length, 1.0e-15);
+
+    // all done
+    return;
+}
+
+TEST(Fem, IsoparametricEmbeddedSegment)
+{
+    // the coordinate system
+    auto coord_system = mito::geometry::coordinate_system_t<coordinates_t>();
+
+    // an atlas under the coordinate system
+    auto atlas = mito::manifolds::atlas<cell_t>(coord_system);
+
+    // build nodes for a diagonal segment from (0,0) to (3,4) - length 5
+    auto node_0 = mito::geometry::node(coord_system, { 0.0, 0.0 });
+    auto node_1 = mito::geometry::node(coord_system, { 3.0, 4.0 });
+
+    // make a geometric simplex
+    auto segment = mito::geometry::segment(node_0, node_1);
+
+    // the normal field to the segment (perpendicular to the segment direction)
+    // segment direction is (3,4)/5 = (0.6, 0.8), so normal is (0.8, -0.6)
+    // (rotated 90° clockwise to get positive orientation with w(normal, tangent) > 0)
+    auto normal_field =
+        mito::functions::constant<coordinates_t>(mito::tensor::vector_t<2>{ 0.8, -0.6 });
+
+    // strip the namespace from the placeholder for forms contractions
+    using mito::tensor::_;
+
+    // the ambient metric volume form
+    constexpr auto w = metric_space_t::w;
+
+    // the restriction of the metric volume form to the segment
+    auto wS = mito::functions::function(
+        [w, normal_field](const coordinates_t & x) { return w(x)(normal_field(x), _); });
+
+    // make a manifold element from the segment
+    auto element =
+        mito::manifolds::parametrized_element(segment, atlas.parametrization(segment), wS);
+
+    {
+        // build the discretization nodes
+        auto discretization_node_0 = discretization_node_t();
+        auto discretization_node_1 = discretization_node_t();
+
+        // the degree of the finite element
+        constexpr int degree = 1;
+        // assemble the finite element type
+        using finite_element_t = mito::fem::finite_element_family<cell_t, degree>;
+
+        // first order isoparametric finite element on the embedded segment
+        auto element_p1 = mito::fem::finite_element<finite_element_t>(
+            element, { discretization_node_0, discretization_node_1 });
+
+        // check that first order shape functions are a partition of unity
+        test_partition_of_unity(element_p1);
+
+        // check that the gradients of first order shape functions sum to 0.0
+        test_gradient_consistency(element_p1);
+
+        // check that the arc length is sqrt(3^2 + 4^2) = 5
+        test_arc_length(element_p1, 5.0);
+
+        // check the gradient values at the midpoint for the (0,0) to (3,4) segment
+        {
+            auto xi = reference_simplex_t::parametric_coordinates_type{ 0.5 };
+            auto grad_0 = element_p1.gradient<0>()(xi);
+            auto grad_1 = element_p1.gradient<1>()(xi);
+
+            // for linear shape functions, the gradients should be constant and opposite
+            // the gradient should be in the direction of the segment: (3/5, 4/5) / 5 = (3/25, 4/25)
+            // phi_0 decreases from 1 to 0, so grad_0 = (-3/25, -4/25)
+            // phi_1 increases from 0 to 1, so grad_1 = (3/25, 4/25)
+            EXPECT_NEAR(-3.0 / 25.0, grad_0[0], 1.0e-15);
+            EXPECT_NEAR(-4.0 / 25.0, grad_0[1], 1.0e-15);
+            EXPECT_NEAR(3.0 / 25.0, grad_1[0], 1.0e-15);
+            EXPECT_NEAR(4.0 / 25.0, grad_1[1], 1.0e-15);
+        }
+    }
+
+    // all done
+    return;
+}
diff --git a/tests/mito.lib/fem/isoparametric_embedded_triangle.cc b/tests/mito.lib/fem/isoparametric_embedded_triangle.cc
new file mode 100644
index 00000000..8ecafcb0
--- /dev/null
+++ b/tests/mito.lib/fem/isoparametric_embedded_triangle.cc
@@ -0,0 +1,203 @@
+// -*- c++ -*-
+//
+// Copyright (c) 2020-2026, the MiTo Authors, all rights reserved
+//
+
+#include <gtest/gtest.h>
+#include <mito.h>
+
+
+// the type of coordinates (3D physical space)
+using coordinates_t = mito::geometry::coordinates_t<3, mito::geometry::CARTESIAN>;
+// the metric space type
+using metric_space_t = mito::geometry::euclidean_metric_space<coordinates_t>;
+// the type of discretization node
+using discretization_node_t = mito::discrete::discretization_node_t;
+// the type of cell (triangle embedded in 3D)
+using cell_t = mito::geometry::triangle_t<3>;
+// the reference simplex
+using reference_simplex_t = cell_t::reference_simplex_type;
+// Gauss quadrature on triangles with degree of exactness 2
+using quadrature_rule_t =
+    mito::quadrature::quadrature_rule_t<mito::quadrature::GAUSS, reference_simplex_t, 2>;
+
+
+// instantiate the quadrature rule
+constexpr auto quadrature_rule = quadrature_rule_t();
+
+
+// test that all shape functions sum to 1.0 at any quadrature point
+auto
+test_partition_of_unity(const auto & element)
+{
+    // the number of quadrature points per element
+    constexpr int n_quads = quadrature_rule_t::npoints;
+
+    // the number of nodes per element
+    constexpr int n_nodes = mito::utilities::base_type<decltype(element)>::n_nodes;
+
+    // loop on the quadrature points
+    mito::tensor::constexpr_for_1<n_quads>([&]<int q>() {
+        // the parametric coordinates of the quadrature point
+        constexpr auto xi = quadrature_rule.point(q);
+
+        // compute the sum of the shape functions at {xi} for all nodes
+        constexpr auto sum =
+            ([]<int... a>(
+                 const auto & element, const auto & xi, mito::tensor::integer_sequence<a...>) {
+                return ((element.template shape<a>()(xi)) + ...);
+            })(element, xi, mito::tensor::make_integer_sequence<n_nodes>{});
+
+        // check the sum of the shape functions
+        static_assert(1.0 == sum);
+    });
+
+    // all done
+    return;
+}
+
+// test that the gradients of all shape functions sum to 0.0 at any quadrature point
+auto
+test_gradient_consistency(const auto & element)
+{
+    // the number of quadrature points per element
+    constexpr int n_quads = quadrature_rule_t::npoints;
+
+    // the number of nodes per element
+    constexpr int n_nodes = mito::utilities::base_type<decltype(element)>::n_nodes;
+
+    // loop on the quadrature points
+    mito::tensor::constexpr_for_1<n_quads>([&]<int q>() {
+        // the parametric coordinates of the quadrature point
+        constexpr auto xi = quadrature_rule.point(q);
+
+        // compute the sum of the shape functions gradients at {xi} for all nodes
+        auto sum =
+            ([]<int... a>(
+                 const auto & element, const auto & xi, mito::tensor::integer_sequence<a...>) {
+                return ((element.template gradient<a>()(xi)) + ...);
+            })(element, xi, mito::tensor::make_integer_sequence<n_nodes>{});
+
+        // check that the sum of the shape functions gradients is the zero vector
+        EXPECT_NEAR(0.0, mito::tensor::norm(sum), 3.0e-16);
+    });
+
+    // all done
+    return;
+}
+
+// test that the area computed via the volume element matches the expected value
+auto
+test_area(const auto & element, double expected_area)
+{
+    // the element type
+    using element_t = mito::utilities::base_type<decltype(element)>;
+
+    // the number of quadrature points per element
+    constexpr int n_quads = quadrature_rule_t::npoints;
+
+    // compute the area by integrating 1 over the triangle using the volume element
+    double area = 0.0;
+    mito::tensor::constexpr_for_1<n_quads>([&]<int q>() {
+        constexpr auto xi = quadrature_rule.point(q);
+        constexpr auto w = element_t::canonical_element_type::measure * quadrature_rule.weight(q);
+        area += w * element.volume_element()(xi);
+    });
+
+    // check the area
+    EXPECT_NEAR(expected_area, area, 1.0e-15);
+
+    // all done
+    return;
+}
+
+TEST(Fem, IsoparametricEmbeddedTriangle)
+{
+    // the coordinate system
+    auto coord_system = mito::geometry::coordinate_system_t<coordinates_t>();
+
+    // an atlas under the coordinate system
+    auto atlas = mito::manifolds::atlas<cell_t>(coord_system);
+
+    // build nodes for a unit right triangle rigidly rotated out of the xy-plane
+    // (the legs {t_1} = (0.6, 0.0, 0.8) and {t_2} = (0.0, 1.0, 0.0) are unit and orthogonal)
+    auto node_0 = mito::geometry::node(coord_system, { 0.0, 0.0, 0.0 });
+    auto node_1 = mito::geometry::node(coord_system, { 0.6, 0.0, 0.8 });
+    auto node_2 = mito::geometry::node(coord_system, { 0.0, 1.0, 0.0 });
+
+    // make a geometric simplex
+    auto triangle = mito::geometry::triangle(node_0, node_1, node_2);
+
+    // the normal field to the triangle: {n} = {t_1} x {t_2} normalized,
+    // oriented such that w(n, t_1, t_2) > 0
+    auto normal_field =
+        mito::functions::constant<coordinates_t>(mito::tensor::vector_t<3>{ -0.8, 0.0, 0.6 });
+
+    // strip the namespace from the placeholder for forms contractions
+    using mito::tensor::_;
+
+    // the ambient metric volume form
+    constexpr auto w = metric_space_t::w;
+
+    // the restriction of the metric volume form to the triangle
+    auto wS = mito::functions::function(
+        [w, normal_field](const coordinates_t & x) { return w(x)(normal_field(x), _, _); });
+
+    // make a manifold element from the triangle
+    auto element =
+        mito::manifolds::parametrized_element(triangle, atlas.parametrization(triangle), wS);
+
+    {
+        // build the discretization nodes
+        auto discretization_node_0 = discretization_node_t();
+        auto discretization_node_1 = discretization_node_t();
+        auto discretization_node_2 = discretization_node_t();
+
+        // the degree of the finite element
+        constexpr int degree = 1;
+        // assemble the finite element type
+        using finite_element_t = mito::fem::finite_element_family<cell_t, degree>;
+
+        // first order isoparametric finite element on the embedded triangle
+        auto element_p1 = mito::fem::finite_element<finite_element_t>(
+            element, { discretization_node_0, discretization_node_1, discretization_node_2 });
+
+        // check that first order shape functions are a partition of unity
+        test_partition_of_unity(element_p1);
+
+        // check that the gradients of first order shape functions sum to 0.0
+        test_gradient_consistency(element_p1);
+
+        // check that the area is that of the unit right triangle
+        test_area(element_p1, 0.5);
+
+        // check the gradient values at the barycenter
+        {
+            auto xi = reference_simplex_t::parametric_coordinates_type{ 1.0 / 3.0, 1.0 / 3.0 };
+            auto grad_0 = element_p1.gradient<0>()(xi);
+            auto grad_1 = element_p1.gradient<1>()(xi);
+            auto grad_2 = element_p1.gradient<2>()(xi);
+
+            // the tangent vectors (the triangle legs) and the normal vector
+            constexpr auto t_1 = mito::tensor::vector_t<3>{ 0.6, 0.0, 0.8 };
+            constexpr auto t_2 = mito::tensor::vector_t<3>{ 0.0, 1.0, 0.0 };
+            constexpr auto n = mito::tensor::vector_t<3>{ -0.8, 0.0, 0.6 };
+
+            // the gradients lie in the plane of the triangle
+            EXPECT_NEAR(0.0, grad_0 * n, 1.0e-15);
+            EXPECT_NEAR(0.0, grad_1 * n, 1.0e-15);
+            EXPECT_NEAR(0.0, grad_2 * n, 1.0e-15);
+
+            // the directional derivatives along the legs recover the parametric derivatives
+            EXPECT_NEAR(-1.0, grad_0 * t_1, 1.0e-15);
+            EXPECT_NEAR(-1.0, grad_0 * t_2, 1.0e-15);
+            EXPECT_NEAR(1.0, grad_1 * t_1, 1.0e-15);
+            EXPECT_NEAR(0.0, grad_1 * t_2, 1.0e-15);
+            EXPECT_NEAR(0.0, grad_2 * t_1, 1.0e-15);
+            EXPECT_NEAR(1.0, grad_2 * t_2, 1.0e-15);
+        }
+    }
+
+    // all done
+    return;
+}
diff --git a/tests/mito.lib/geometry/euclidean_submanifold_metric_2D.cc b/tests/mito.lib/geometry/euclidean_submanifold_metric_2D.cc
new file mode 100644
index 00000000..52d8257a
--- /dev/null
+++ b/tests/mito.lib/geometry/euclidean_submanifold_metric_2D.cc
@@ -0,0 +1,180 @@
+// -*- c++ -*-
+//
+// Copyright (c) 2020-2026, the MiTo Authors, all rights reserved
+//
+
+#include <gtest/gtest.h>
+#include <mito/geometry.h>
+#include <cmath>
+#include <vector>
+
+
+// cartesian coordinates in 2D
+using coordinates_t = mito::geometry::coordinates_t<2, mito::geometry::CARTESIAN>;
+
+// the basis for vectors
+static constexpr auto e_x = mito::tensor::e_0<2>;
+static constexpr auto e_y = mito::tensor::e_1<2>;
+
+
+// the placeholder for empty slots in contractions
+using mito::tensor::_;
+
+
+auto
+length(
+    const auto & w, const mito::geometry::coordinate_system_t<coordinates_t> & coordinate_system,
+    const mito::geometry::node_t<2> & v0, const mito::geometry::node_t<2> & v1)
+    -> mito::tensor::scalar_t
+{
+    // get vertex coordinates
+    auto x0 = coordinate_system.coordinates(v0->point());
+    auto x1 = coordinate_system.coordinates(v1->point());
+
+    // build director vector
+    auto director0 = x1 - x0;
+
+    // compute signed length of segment
+    auto length = w(director0);
+
+    // all done
+    return length;
+}
+
+
+TEST(Tensor, EuclideanSubmanifoldMetric2D)
+{
+    // the basis one-forms
+    constexpr auto dx = mito::tensor::one_form(e_x);
+    constexpr auto dy = mito::tensor::one_form(e_y);
+
+    // the coordinate system
+    auto coord_system = mito::geometry::coordinate_system<coordinates_t>();
+
+    // pick two sets of coordinates
+    constexpr auto x_0 = mito::geometry::cartesian::coordinates({ 0.0, 0.0 });
+    constexpr auto x_1 = mito::geometry::cartesian::coordinates({ 3.0, 4.0 });
+
+    // assemble the normal vector to the submanifold from tangent
+    constexpr auto tangent_vector = x_1 - x_0;
+    constexpr auto normal_vector =
+        mito::tensor::vector_t<2>{ tangent_vector[1], -tangent_vector[0] }
+        / mito::tensor::norm(tangent_vector);
+
+    // the 2D metric volume element
+    constexpr auto w = mito::tensor::wedge(dx, dy);
+
+    // the 1D restriction of the 2D metric volume element
+    constexpr auto wS = w(normal_vector, _);
+
+    // build nodes of a line segment
+    auto node_0 = mito::geometry::node(coord_system, x_0);
+    auto node_1 = mito::geometry::node(coord_system, x_1);
+
+    // check that the length of the line segment is correct
+    EXPECT_DOUBLE_EQ(5.0, length(wS, coord_system, node_0, node_1));
+
+    // check that the odd permutation of the vertices gives a negative length
+    EXPECT_DOUBLE_EQ(-5.0, length(wS, coord_system, node_1, node_0));
+
+    // the normal form to the submanifold
+    constexpr auto normal_form = mito::tensor::one_form(normal_vector);
+
+    // rebuild the volume form as wV = wedge(normal_form, wS)
+    constexpr auto wV = mito::tensor::wedge(normal_form, wS);
+
+    // check that wV coincides with w
+    EXPECT_DOUBLE_EQ(w(e_x, e_y), wV(e_x, e_y));
+}
+
+
+// helper to test a single segment: computes the restricted volume form and verifies length
+auto
+test_segment(
+    mito::geometry::coordinate_system_t<coordinates_t> & coord_system, const auto & w,
+    const coordinates_t & x0, const coordinates_t & x1, double expected_length) -> void
+{
+    // compute tangent and normal vectors
+    auto tangent_vector = x1 - x0;
+    auto tangent_norm = mito::tensor::norm(tangent_vector);
+
+    // skip degenerate segments
+    if (tangent_norm < 1.0e-14) {
+        return;
+    }
+
+    // normal is tangent rotated 90° clockwise: (tx, ty) -> (ty, -tx)
+    // this gives positive length for w(normal, tangent)
+    auto normal_vector =
+        mito::tensor::vector_t<2>{ tangent_vector[1], -tangent_vector[0] } / tangent_norm;
+
+    // the 1D restriction of the 2D metric volume element
+    auto wS = w(normal_vector, _);
+
+    // build nodes
+    auto node_0 = mito::geometry::node(coord_system, x0);
+    auto node_1 = mito::geometry::node(coord_system, x1);
+
+    // check length (positive orientation)
+    EXPECT_NEAR(expected_length, length(wS, coord_system, node_0, node_1), 1.0e-14);
+
+    // check length (negative orientation)
+    EXPECT_NEAR(-expected_length, length(wS, coord_system, node_1, node_0), 1.0e-14);
+}
+
+
+TEST(Tensor, EuclideanSubmanifoldMetric2DEdgeCases)
+{
+    // the basis one-forms
+    constexpr auto dx = mito::tensor::one_form(e_x);
+    constexpr auto dy = mito::tensor::one_form(e_y);
+
+    // the 2D metric volume element
+    constexpr auto w = mito::tensor::wedge(dx, dy);
+
+    // the coordinate system
+    auto coord_system = mito::geometry::coordinate_system<coordinates_t>();
+
+    // define test cases: (x0, y0), (x1, y1), expected_length
+    // covering various orientations and edge cases
+    struct SegmentTestCase {
+        double x0, y0, x1, y1;
+        double expected_length;
+        const char * description;
+    };
+
+    std::vector<SegmentTestCase> test_cases = {
+        // axis-aligned segments
+        { 0.0, 0.0, 1.0, 0.0, 1.0, "horizontal unit segment" },
+        { 0.0, 0.0, 0.0, 1.0, 1.0, "vertical unit segment" },
+
+        // 45-degree diagonal
+        { 0.0, 0.0, 1.0, 1.0, std::sqrt(2.0), "diagonal" },
+
+        // Pythagorean triple segments
+        { 0.0, 0.0, 3.0, 4.0, 5.0, "3-4-5 segment" },
+        { 0.0, 0.0, 4.0, 3.0, 5.0, "4-3-5 segment (swapped)" },
+
+        // steep and shallow angles
+        { 0.0, 0.0, 0.1, 10.0, std::sqrt(0.01 + 100.0), "nearly vertical" },
+        { 0.0, 0.0, 10.0, 0.1, std::sqrt(100.0 + 0.01), "nearly horizontal" },
+
+        // non-origin segment
+        { 1.0, 2.0, 4.0, 6.0, 5.0, "3-4-5 offset from origin" },
+
+        // numerical stability at different scales
+        { 0.0, 0.0, 1.0e-6, 1.0e-6, std::sqrt(2.0) * 1.0e-6, "tiny segment" },
+        { 0.0, 0.0, 300.0, 400.0, 500.0, "large segment" },
+    };
+
+    // run all test cases
+    for (const auto & tc : test_cases) {
+        SCOPED_TRACE(tc.description);
+        auto x0 = mito::geometry::cartesian::coordinates({ tc.x0, tc.y0 });
+        auto x1 = mito::geometry::cartesian::coordinates({ tc.x1, tc.y1 });
+        test_segment(coord_system, w, x0, x1, tc.expected_length);
+    }
+}
+
+
+// end of file