diff --git a/.cmake/mito_tests_mito_lib.cmake b/.cmake/mito_tests_mito_lib.cmake
index c4d6f232..32ba8770 100644
--- a/.cmake/mito_tests_mito_lib.cmake
+++ b/.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)
diff --git a/benchmarks/mito.lib/pdes/poisson.cc b/benchmarks/mito.lib/pdes/poisson.cc
index 69845671..f4e35f8b 100644
--- a/benchmarks/mito.lib/pdes/poisson.cc
+++ b/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);
 
diff --git a/lib/mito/fem/blocks/GradGradBlock.h b/lib/mito/fem/blocks/GradGradBlock.h
index 65d8eb2d..1acb3a11 100644
--- a/lib/mito/fem/blocks/GradGradBlock.h
+++ b/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>() {
diff --git a/lib/mito/fem/blocks/L2NormBlock.h b/lib/mito/fem/blocks/L2NormBlock.h
index 0d0a00f4..90ac3be2 100644
--- a/lib/mito/fem/blocks/L2NormBlock.h
+++ b/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);
diff --git a/lib/mito/fem/blocks/MassBlock.h b/lib/mito/fem/blocks/MassBlock.h
index d8300760..98db27c6 100644
--- a/lib/mito/fem/blocks/MassBlock.h
+++ b/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>() {
diff --git a/lib/mito/fem/blocks/SourceTermBlock.h b/lib/mito/fem/blocks/SourceTermBlock.h
index ed9ecb8e..2c451302 100644
--- a/lib/mito/fem/blocks/SourceTermBlock.h
+++ b/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>() {
diff --git a/lib/mito/fem/elements/IsoparametricSegment.h b/lib/mito/fem/elements/IsoparametricSegment.h
index cc3b384b..84e5eff4 100644
--- a/lib/mito/fem/elements/IsoparametricSegment.h
+++ b/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;
diff --git a/lib/mito/fem/elements/IsoparametricTriangle.h b/lib/mito/fem/elements/IsoparametricTriangle.h
index 5426f2ad..76cf08b0 100644
--- a/lib/mito/fem/elements/IsoparametricTriangle.h
+++ b/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,6 +81,15 @@ 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;
@@ -70,6 +97,9 @@ namespace mito::fem {
         // 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;
diff --git a/lib/mito/fem/elements/api.h b/lib/mito/fem/elements/api.h
index 89cc64b5..37ec1a84 100644
--- a/lib/mito/fem/elements/api.h
+++ b/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
 }
 
 
diff --git a/lib/mito/fem/elements/elements_library.h b/lib/mito/fem/elements/elements_library.h
index 60e59338..3355bc54 100644
--- a/lib/mito/fem/elements/elements_library.h
+++ b/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
diff --git a/lib/mito/fem/elements/forward.h b/lib/mito/fem/elements/forward.h
index 36e247d2..a0a68fec 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;
 
 }
diff --git a/lib/mito/fem/elements/seg1/DiscretizerCG.h b/lib/mito/fem/elements/seg1/DiscretizerCG.h
index 0043129d..5a2a01df 100644
--- a/lib/mito/fem/elements/seg1/DiscretizerCG.h
+++ b/lib/mito/fem/elements/seg1/DiscretizerCG.h
@@ -10,8 +10,9 @@
 namespace mito::fem {
 
     // discretizer specialization for {IsoparametricSegmentP1} with continuous Galerkin
-    template <>
-    struct Discretizer<IsoparametricSegmentP1, discretization_t::CG> {
+    // agnostic of the coordinate type and volume form
+    template <geometry::coordinates_c coordsT, class VolumeFormT>
+    struct Discretizer<IsoparametricSegmentP1<coordsT, VolumeFormT>, discretization_t::CG> {
         template <
             typename manifoldT, typename constraintsT, typename elements_type, typename map_type,
             typename constrained_nodes_type>
@@ -19,17 +20,21 @@ namespace mito::fem {
             const manifoldT & manifold, const constraintsT & constraints, elements_type & elements,
             map_type & node_map, constrained_nodes_type & constrained_nodes)
         {
+            // the element type
+            using element_type = IsoparametricSegmentP1<coordsT, VolumeFormT>;
 
             // the discretization node type
-            using discretization_node_type =
-                typename IsoparametricSegmentP1::discretization_node_type;
+            using discretization_node_type = typename element_type::discretization_node_type;
 
             // the connectivity type
-            using connectivity_type = typename IsoparametricSegmentP1::connectivity_type;
+            using connectivity_type = typename element_type::connectivity_type;
 
             // get the coordinate system of the manifold
             const auto & coord_system = manifold.coordinate_system();
 
+            // get the volume form from the manifold
+            const auto & volume_form = manifold.volume_form();
+
             // loop on the cells of the mesh
             for (const auto & cell : manifold.elements()) {
 
@@ -43,8 +48,10 @@ namespace mito::fem {
                 auto node_1 =
                     node_map.insert({ nodes[1], discretization_node_type() }).first->second;
 
-                // create a finite element for each cell and add it to the pile
-                elements.emplace(cell, coord_system, connectivity_type{ node_0, node_1 });
+                // create a finite element for each cell and add it to the pile, passing the volume
+                // form to the element
+                elements.emplace(
+                    cell, coord_system, connectivity_type{ node_0, node_1 }, volume_form);
             }
 
             // populate the constrained nodes
diff --git a/lib/mito/fem/elements/seg1/IsoparametricSegmentP1.h b/lib/mito/fem/elements/seg1/IsoparametricSegmentP1.h
index eee35058..ed24a93a 100644
--- a/lib/mito/fem/elements/seg1/IsoparametricSegmentP1.h
+++ b/lib/mito/fem/elements/seg1/IsoparametricSegmentP1.h
@@ -8,15 +8,27 @@
 
 
 // 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 equipped with linear
+// shape functions defined in the parametric space. This unified implementation works for:
+//   - Segments in any coordinate system (Cartesian, polar, spherical, etc.)
+//   - Segments in any embedding (1D in 1D, 1D in 2D, 1D in 3D, etc.)
 
 
 namespace mito::fem {
 
-    class IsoparametricSegmentP1 : public IsoparametricSegment {
+    template <geometry::coordinates_c coordsT, class VolumeFormT>
+    class IsoparametricSegmentP1 : public IsoparametricSegment<coordsT, VolumeFormT> {
+
+      private:
+        using base_type = IsoparametricSegment<coordsT, VolumeFormT>;
+        using metric_space_type = typename base_type::metric_space_type;
 
       public:
+        static constexpr int D = base_type::D;    // ambient dimension
+        static constexpr int N =
+            base_type::N;                // dimension of the parametric space (=1 for segment)
+        static constexpr int dim = D;    // expected by FunctionSpace
+
         // the degree of the finite element
         static constexpr int degree = 1;
         // the type of shape functions
@@ -31,14 +43,21 @@ namespace mito::fem {
         // the number of discretization nodes
         static constexpr int n_nodes = shape_functions_type::N;
         // a collection of discretization nodes
-        using connectivity_type = std::array<discretization_node_type, n_nodes>;
+        using connectivity_type = std::array<typename base_type::discretization_node_type, n_nodes>;
+
+        // import types from base
+        using typename base_type::cell_type;
+        using typename base_type::coordinate_system_type;
+        using typename base_type::coordinates_type;
+        using typename base_type::vector_type;
+        using typename base_type::volume_form_type;
 
       public:
         // the default constructor
         inline IsoparametricSegmentP1(
             const cell_type & geometric_simplex, const coordinate_system_type & coord_system,
-            const connectivity_type & connectivity) :
-            IsoparametricSegment(geometric_simplex, coord_system),
+            const connectivity_type & connectivity, const volume_form_type & volume_form) :
+            base_type(geometric_simplex, coord_system, volume_form),
             _connectivity(connectivity)
         {}
 
@@ -72,7 +91,8 @@ namespace mito::fem {
             constexpr auto phi_1 = shape_functions.shape<1>();
 
             // return the isoparametric mapping from parametric to physical coordinates
-            return mito::functions::linear_combination(std::array{ _x0, _x1 }, phi_0, phi_1);
+            return mito::functions::linear_combination(
+                std::array{ this->_x0, this->_x1 }, phi_0, phi_1);
         }
 
         // get the shape function associated with local node {a}
@@ -85,40 +105,71 @@ namespace mito::fem {
         }
 
         // get the jacobian of the isoparametric mapping from parametric to actual coordinates
+        // returns a D×1 matrix (the tangent/director vector as a column)
         constexpr auto jacobian() const
         {
             // 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<D, 1> {
                     // get the shape functions derivatives
                     constexpr auto dphi_0 = shape_functions.dshape<0>();
                     constexpr auto dphi_1 = shape_functions.dshape<1>();
 
-                    // 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 };
+                    // compute the tangent vector: dx/dxi
+                    auto dx_dxi = this->_x0 * dphi_0(xi) + this->_x1 * dphi_1(xi);
+
+                    // wrap the tangent vector as a D×1 matrix
+                    return tensor::as_column_matrix<D>(dx_dxi);
                 });
 
             // and return it
             return jacobian_function;
         }
 
+        // volume element: contract the volume form with the tangent vector
+        // this follows the same pattern as Manifold::_volume, but for one tangent vector
+        constexpr auto volume_element() const
+        {
+            return functions::function([&](const parametric_coordinates_type & xi) {
+                // physical point
+                auto x = parametrization()(xi);
+
+                // get the tangent vector from the Jacobian
+                auto J = jacobian()(xi);
+                // extract the tangent vector from the D×1 matrix
+                auto tangent = pyre::tensor::col<0>(J);
+
+                // contract the volume form with the tangent vector
+                // for N=1, no factorial needed (1/1! = 1)
+                return this->_volume_form(x)(tangent);
+            });
+        }
+
         // get the gradient of the a-th shape function as a function of parametric coordinates
+        // returns a contravariant vector (raised index)
         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 gradient_function =
+                functions::function([&](const parametric_coordinates_type & xi) -> vector_type {
+                    // physical point
+                    auto x = parametrization()(xi);
+
+                    // jacobian matrix (D×1 for curves)
                     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;
+
+                    // metric tensors at the physical point
+                    auto g = base_type::g(x);
+                    auto g_inv = base_type::g_inv(x);
+
+                    // parametric derivative (scalar for N=1)
+                    constexpr auto dphi_dxi = shape_functions.dshape<a>();
+                    auto dphi_dxi_val = dphi_dxi(xi);
+
+                    // compute gradient (handles both square and embedded cases)
+                    return fem::compute_gradient<D, N>(J, g, g_inv, dphi_dxi_val);
                 });
             // and return it
             return gradient_function;
diff --git a/lib/mito/fem/elements/seg1/api.h b/lib/mito/fem/elements/seg1/api.h
index cbdd9a0e..37ec1a84 100644
--- a/lib/mito/fem/elements/seg1/api.h
+++ b/lib/mito/fem/elements/seg1/api.h
@@ -8,13 +8,7 @@
 
 
 namespace mito::fem {
-
-    // specialization for linear shape functions on segments in 1D
-    template <>
-    struct isoparametric_simplex<1, geometry::segment_t<1>> {
-        using type = IsoparametricSegmentP1;
-    };
-
+    // no specializations needed here - see elements_library.h for the unified API
 }
 
 
diff --git a/lib/mito/fem/elements/tri1/DiscretizerCG.h b/lib/mito/fem/elements/tri1/DiscretizerCG.h
index 2c4f2bc0..90bd58de 100644
--- a/lib/mito/fem/elements/tri1/DiscretizerCG.h
+++ b/lib/mito/fem/elements/tri1/DiscretizerCG.h
@@ -10,8 +10,9 @@
 namespace mito::fem {
 
     // discretizer specialization for {IsoparametricTriangleP1} with continuous Galerkin
-    template <>
-    struct Discretizer<IsoparametricTriangleP1, discretization_t::CG> {
+    // agnostic of the coordinate type and volume form
+    template <geometry::coordinates_c coordsT, class VolumeFormT>
+    struct Discretizer<IsoparametricTriangleP1<coordsT, VolumeFormT>, discretization_t::CG> {
         template <
             typename manifoldT, typename constraintsT, typename elements_type, typename map_type,
             typename constrained_nodes_type>
@@ -19,17 +20,21 @@ namespace mito::fem {
             const manifoldT & manifold, const constraintsT & constraints, elements_type & elements,
             map_type & node_map, constrained_nodes_type & constrained_nodes)
         {
+            // the element type
+            using element_type = IsoparametricTriangleP1<coordsT, VolumeFormT>;
 
             // the discretization node type
-            using discretization_node_type =
-                typename IsoparametricTriangleP1::discretization_node_type;
+            using discretization_node_type = typename element_type::discretization_node_type;
 
             // the connectivity type
-            using connectivity_type = typename IsoparametricTriangleP1::connectivity_type;
+            using connectivity_type = typename element_type::connectivity_type;
 
             // get the coordinate system of the manifold
             const auto & coord_system = manifold.coordinate_system();
 
+            // get the volume form from the manifold
+            const auto & volume_form = manifold.volume_form();
+
             // loop on the cells of the mesh
             for (const auto & cell : manifold.elements()) {
 
@@ -45,8 +50,10 @@ namespace mito::fem {
                 auto node_2 =
                     node_map.insert({ nodes[2], discretization_node_type() }).first->second;
 
-                // create a finite element for each cell and add it to the pile
-                elements.emplace(cell, coord_system, connectivity_type{ node_0, node_1, node_2 });
+                // create a finite element for each cell and add it to the pile, passing the volume
+                // form to the element
+                elements.emplace(
+                    cell, coord_system, connectivity_type{ node_0, node_1, node_2 }, volume_form);
             }
 
             // populate the constrained nodes
diff --git a/lib/mito/fem/elements/tri1/IsoparametricTriangleP1.h b/lib/mito/fem/elements/tri1/IsoparametricTriangleP1.h
index 6364a68e..552df9ec 100644
--- a/lib/mito/fem/elements/tri1/IsoparametricTriangleP1.h
+++ b/lib/mito/fem/elements/tri1/IsoparametricTriangleP1.h
@@ -8,15 +8,27 @@
 
 
 // DESIGN NOTES
-// Class {IsoparametricTriangleP1} represents a second order simplex living in 2D cartesian space,
-// equipped with linear shape functions defined in the parametric space.
+// Class {IsoparametricTriangleP1} represents a 2-simplex equipped with linear shape
+// functions defined in the parametric space. This unified implementation works for:
+//   - Triangles in any coordinate system (Cartesian, polar, spherical, etc.)
+//   - Triangles in any embedding (2D in 2D, 2D in 3D, etc.)
 
 
 namespace mito::fem {
 
-    class IsoparametricTriangleP1 : public IsoparametricTriangle {
+    template <geometry::coordinates_c coordsT, class VolumeFormT>
+    class IsoparametricTriangleP1 : public IsoparametricTriangle<coordsT, VolumeFormT> {
+
+      private:
+        using base_type = IsoparametricTriangle<coordsT, VolumeFormT>;
+        using metric_space_type = typename base_type::metric_space_type;
 
       public:
+        static constexpr int D = base_type::D;    // ambient dimension
+        static constexpr int N =
+            base_type::N;                // dimension of the parametric space (=2 for triangle)
+        static constexpr int dim = D;    // expected by FunctionSpace
+
         // the degree of the finite element
         static constexpr int degree = 1;
         // the type of shape functions
@@ -31,14 +43,21 @@ namespace mito::fem {
         // the number of discretization discretization nodes
         static constexpr int n_nodes = shape_functions_type::N;
         // a collection of discretization discretization nodes
-        using connectivity_type = std::array<discretization_node_type, n_nodes>;
+        using connectivity_type = std::array<typename base_type::discretization_node_type, n_nodes>;
+
+        // import types from base
+        using typename base_type::cell_type;
+        using typename base_type::coordinate_system_type;
+        using typename base_type::coordinates_type;
+        using typename base_type::vector_type;
+        using typename base_type::volume_form_type;
 
       public:
         // the default constructor
         inline IsoparametricTriangleP1(
             const cell_type & geometric_simplex, const coordinate_system_type & coord_system,
-            const connectivity_type & connectivity) :
-            IsoparametricTriangle(geometric_simplex, coord_system),
+            const connectivity_type & connectivity, const volume_form_type & volume_form) :
+            base_type(geometric_simplex, coord_system, volume_form),
             _connectivity(connectivity)
         {}
 
@@ -74,41 +93,76 @@ namespace mito::fem {
         }
 
         // get the jacobian of the isoparametric mapping from barycentric to actual coordinates
+        // returns a D×2 matrix (columns are the two tangent/director vectors)
         constexpr auto jacobian() const
         {
             // assemble the jacobian as a function of barycentric coordinates
             auto jacobian_function = functions::function(
-                [&](const parametric_coordinates_type & xi) -> tensor::matrix_t<2> {
-                    // get the shape functions derivatives
+                [&](const parametric_coordinates_type & xi) -> tensor::matrix_t<D, 2> {
+                    // get the shape functions derivatives (these are 2D vectors in parametric
+                    // space)
                     constexpr auto dphi_0 = shape_functions.dshape<0>();
                     constexpr auto dphi_1 = shape_functions.dshape<1>();
                     constexpr auto dphi_2 = shape_functions.dshape<2>();
 
                     // compute the gradient of the isoparametric mapping
+                    // this is the outer product of position vectors with gradient vectors
                     return (
-                        tensor::dyadic(_x0, dphi_0(xi)) + tensor::dyadic(_x1, dphi_1(xi))
-                        + tensor::dyadic(_x2, dphi_2(xi)));
+                        tensor::dyadic(this->_x0, dphi_0(xi))
+                        + tensor::dyadic(this->_x1, dphi_1(xi))
+                        + tensor::dyadic(this->_x2, dphi_2(xi)));
                 });
 
             // and return it
             return jacobian_function;
         }
 
+        // volume element: contract the volume form with the two tangent vectors
+        // we don't include the 1/N! factorial here because it's already included
+        // in the canonical_element_type::area used by the assembly blocks
+        constexpr auto volume_element() const
+        {
+            return functions::function([&](const parametric_coordinates_type & xi) {
+                // physical point
+                auto x = this->parametrization()(xi);
+
+                // get the two tangent vectors from the Jacobian
+                auto J = jacobian()(xi);
+
+                // extract tangent vectors from the D×2 matrix
+                auto tangent_0 = pyre::tensor::col<0>(J);
+                auto tangent_1 = pyre::tensor::col<1>(J);
+
+                // contract the volume form with both tangent vectors
+                return this->_volume_form(x)(tangent_0, tangent_1);
+            });
+        }
+
         // get the gradient of the a-th shape function as a function of barycentric coordinates
+        // returns a contravariant vector (raised index)
         template <int a>
         requires(a >= 0 && a < n_nodes)
         constexpr auto gradient() const
         {
             // assemble the gradient as a function of barycentric 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 gradient_function =
+                functions::function([&](const parametric_coordinates_type & xi) -> vector_type {
+                    // physical point
+                    auto x = this->parametrization()(xi);
+
+                    // jacobian matrix (D×2 for surfaces)
                     auto J = jacobian()(xi);
-                    // the derivative of the coordinates with respect to the barycentric 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;
+
+                    // metric tensors at the physical point
+                    auto g = base_type::g(x);
+                    auto g_inv = base_type::g_inv(x);
+
+                    // parametric derivative (2D vector for N=2)
+                    constexpr auto dphi_dxi = shape_functions.dshape<a>();
+                    auto dphi_dxi_val = dphi_dxi(xi);
+
+                    // compute gradient (handles both square and embedded cases)
+                    return fem::compute_gradient<D, N>(J, g, g_inv, dphi_dxi_val);
                 });
             // and return it
             return gradient_function;
diff --git a/lib/mito/fem/elements/tri1/api.h b/lib/mito/fem/elements/tri1/api.h
index 388c39f6..58ed3454 100644
--- a/lib/mito/fem/elements/tri1/api.h
+++ b/lib/mito/fem/elements/tri1/api.h
@@ -8,13 +8,7 @@
 
 
 namespace mito::fem {
-
-    // specialization for linear shape functions on triangles in 2D
-    template <>
-    struct isoparametric_simplex<1, geometry::triangle_t<2>> {
-        using type = IsoparametricTriangleP1;
-    };
-
+    // no element specializations needed here - see elements_library.h for the unified API
 }
 
 
diff --git a/lib/mito/fem/elements/tri1/public.h b/lib/mito/fem/elements/tri1/public.h
index d5dc5c0d..8e6ca05c 100644
--- a/lib/mito/fem/elements/tri1/public.h
+++ b/lib/mito/fem/elements/tri1/public.h
@@ -12,8 +12,5 @@
 #include "IsoparametricTriangleP1.h"
 #include "DiscretizerCG.h"
 
-// published types and factories
-#include "api.h"
-
 
 // end of file
diff --git a/lib/mito/fem/elements/tri2/DiscretizerCG.h b/lib/mito/fem/elements/tri2/DiscretizerCG.h
index 6ab25111..2f091675 100644
--- a/lib/mito/fem/elements/tri2/DiscretizerCG.h
+++ b/lib/mito/fem/elements/tri2/DiscretizerCG.h
@@ -10,8 +10,9 @@
 namespace mito::fem {
 
     // discretizer specialization for {IsoparametricTriangleP2} with continuous Galerkin
-    template <>
-    struct Discretizer<IsoparametricTriangleP2, discretization_t::CG> {
+    // agnostic of the coordinate type and volume form
+    template <geometry::coordinates_c coordsT, class VolumeFormT>
+    struct Discretizer<IsoparametricTriangleP2<coordsT, VolumeFormT>, discretization_t::CG> {
         template <
             typename manifoldT, typename constraintsT, typename elements_type, typename map_type,
             typename constrained_nodes_type>
@@ -19,18 +20,20 @@ namespace mito::fem {
             const manifoldT & manifold, const constraintsT & constraints, elements_type & elements,
             map_type & node_map, constrained_nodes_type & constrained_nodes)
         {
+            // the element type
+            using element_type = IsoparametricTriangleP2<coordsT, VolumeFormT>;
+
             // the dimension of the physical space
-            constexpr int dim = IsoparametricTriangleP2::dim;
+            constexpr int D = element_type::D;
 
             // assemble the mesh node type
-            using mesh_node_type = geometry::node_t<dim>;
+            using mesh_node_type = geometry::node_t<D>;
 
             // the discretization node type
-            using discretization_node_type =
-                typename IsoparametricTriangleP2::discretization_node_type;
+            using discretization_node_type = typename element_type::discretization_node_type;
 
             // the connectivity type
-            using connectivity_type = typename IsoparametricTriangleP2::connectivity_type;
+            using connectivity_type = typename element_type::connectivity_type;
 
             // id type of mesh nodes
             using mesh_node_id_t = utilities::index_t<mesh_node_type>;
@@ -45,6 +48,9 @@ namespace mito::fem {
             // get the coordinate system of the manifold
             const auto & coord_system = manifold.coordinate_system();
 
+            // get the volume form from the manifold
+            const auto & volume_form = manifold.volume_form();
+
             // loop on the cells of the mesh
             for (const auto & cell : manifold.elements()) {
 
@@ -76,10 +82,12 @@ namespace mito::fem {
                 auto node_5 = mid_nodes_map.insert({ ordered_nodes_5, discretization_node_type() })
                                   .first->second;
 
-                // create a finite element for each cell and add it to the pile
+                // create a finite element for each cell and add it to the pile, passing the volume
+                // form to the element
                 elements.emplace(
                     cell, coord_system,
-                    connectivity_type{ node_0, node_1, node_2, node_3, node_4, node_5 });
+                    connectivity_type{ node_0, node_1, node_2, node_3, node_4, node_5 },
+                    volume_form);
             }
 
             // populate the constrained nodes
diff --git a/lib/mito/fem/elements/tri2/IsoparametricTriangleP2.h b/lib/mito/fem/elements/tri2/IsoparametricTriangleP2.h
index 6a13b1fc..8a2e72f6 100644
--- a/lib/mito/fem/elements/tri2/IsoparametricTriangleP2.h
+++ b/lib/mito/fem/elements/tri2/IsoparametricTriangleP2.h
@@ -8,15 +8,27 @@
 
 
 // DESIGN NOTES
-// Class {IsoparametricTriangleP2} represents a second order simplex living in 2D cartesian space
-// equipped with quadratic shape functions defined in the parametric space.
+// Class {IsoparametricTriangleP2} represents a 2-simplex equipped with quadratic shape
+// functions defined in the parametric space. This unified implementation works for:
+//   - Triangles in any coordinate system (Cartesian, polar, spherical, etc.)
+//   - Triangles in any embedding (2D in 2D, 2D in 3D, etc.)
 
 
 namespace mito::fem {
 
-    class IsoparametricTriangleP2 : public IsoparametricTriangle {
+    template <geometry::coordinates_c coordsT, class VolumeFormT>
+    class IsoparametricTriangleP2 : public IsoparametricTriangle<coordsT, VolumeFormT> {
+
+      private:
+        using base_type = IsoparametricTriangle<coordsT, VolumeFormT>;
+        using metric_space_type = typename base_type::metric_space_type;
 
       public:
+        static constexpr int D = base_type::D;    // ambient dimension
+        static constexpr int N =
+            base_type::N;                // dimension of the parametric space (=2 for triangle)
+        static constexpr int dim = D;    // expected by FunctionSpace
+
         // the degree of the finite element
         static constexpr int degree = 2;
         // the type of shape functions
@@ -26,19 +38,26 @@ namespace mito::fem {
         // type of a point in parametric coordinates
         using parametric_coordinates_type =
             typename canonical_element_type::parametric_coordinates_type;
-        // the linear shape functions
+        // the quadratic shape functions
         static constexpr auto shape_functions = shape_functions_type();
         // the number of discretization nodes
         static constexpr int n_nodes = shape_functions_type::N;
         // a collection of discretization nodes
-        using connectivity_type = std::array<discretization_node_type, n_nodes>;
+        using connectivity_type = std::array<typename base_type::discretization_node_type, n_nodes>;
+
+        // import types from base
+        using typename base_type::cell_type;
+        using typename base_type::coordinate_system_type;
+        using typename base_type::coordinates_type;
+        using typename base_type::vector_type;
+        using typename base_type::volume_form_type;
 
       public:
         // the default constructor
         inline IsoparametricTriangleP2(
             const cell_type & geometric_simplex, const coordinate_system_type & coord_system,
-            const connectivity_type & connectivity) :
-            IsoparametricTriangle(geometric_simplex, coord_system),
+            const connectivity_type & connectivity, const volume_form_type & volume_form) :
+            base_type(geometric_simplex, coord_system, volume_form),
             _connectivity(connectivity)
         {}
 
@@ -74,14 +93,15 @@ namespace mito::fem {
         }
 
         // get the jacobian of the isoparametric mapping from barycentric to actual coordinates
+        // returns a D×2 matrix (columns are the two tangent/director vectors)
         constexpr auto jacobian() const
         {
             // assemble the jacobian as a function of barycentric coordinates
             auto jacobian_function = functions::function(
-                [&](const parametric_coordinates_type & xi) -> tensor::matrix_t<2> {
-                    auto x3 = 0.5 * (_x0 + _x1);
-                    auto x4 = 0.5 * (_x1 + _x2);
-                    auto x5 = 0.5 * (_x2 + _x0);
+                [&](const parametric_coordinates_type & xi) -> tensor::matrix_t<D, 2> {
+                    auto x3 = 0.5 * (this->_x0 + this->_x1);
+                    auto x4 = 0.5 * (this->_x1 + this->_x2);
+                    auto x5 = 0.5 * (this->_x2 + this->_x0);
 
                     // get the shape functions derivatives
                     constexpr auto dphi_0 = shape_functions.dshape<0>();
@@ -93,8 +113,9 @@ namespace mito::fem {
 
                     // compute the gradient of the isoparametric mapping
                     return (
-                        tensor::dyadic(_x0, dphi_0(xi)) + tensor::dyadic(_x1, dphi_1(xi))
-                        + tensor::dyadic(_x2, dphi_2(xi)) + tensor::dyadic(x3, dphi_3(xi))
+                        tensor::dyadic(this->_x0, dphi_0(xi))
+                        + tensor::dyadic(this->_x1, dphi_1(xi))
+                        + tensor::dyadic(this->_x2, dphi_2(xi)) + tensor::dyadic(x3, dphi_3(xi))
                         + tensor::dyadic(x4, dphi_4(xi)) + tensor::dyadic(x5, dphi_5(xi)));
                 });
 
@@ -102,21 +123,52 @@ namespace mito::fem {
             return jacobian_function;
         }
 
+        // volume element: contract the volume form with the two tangent vectors
+        // we don't include the 1/N! factorial here because it's already included
+        // in the canonical_element_type::area used by the assembly blocks
+        constexpr auto volume_element() const
+        {
+            return functions::function([&](const parametric_coordinates_type & xi) {
+                // physical point
+                auto x = this->parametrization()(xi);
+
+                // get the two tangent vectors from the Jacobian
+                auto J = jacobian()(xi);
+
+                // extract tangent vectors from the D×2 matrix
+                auto tangent_0 = pyre::tensor::col<0>(J);
+                auto tangent_1 = pyre::tensor::col<1>(J);
+
+                // contract the volume form with both tangent vectors
+                return this->_volume_form(x)(tangent_0, tangent_1);
+            });
+        }
+
         // get the gradient of the a-th shape function as a function of barycentric coordinates
+        // returns a contravariant vector (raised index)
         template <int a>
         requires(a >= 0 && a < n_nodes)
         constexpr auto gradient() const
         {
             // assemble the gradient as a function of barycentric 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 gradient_function =
+                functions::function([&](const parametric_coordinates_type & xi) -> vector_type {
+                    // physical point
+                    auto x = this->parametrization()(xi);
+
+                    // jacobian matrix (D×2 for surfaces)
                     auto J = jacobian()(xi);
-                    // the derivative of the coordinates with respect to the barycentric 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;
+
+                    // metric tensors at the physical point
+                    auto g = base_type::g(x);
+                    auto g_inv = base_type::g_inv(x);
+
+                    // parametric derivative (2D vector for N=2)
+                    constexpr auto dphi_dxi = shape_functions.dshape<a>();
+                    auto dphi_dxi_val = dphi_dxi(xi);
+
+                    // compute gradient (handles both square and embedded cases)
+                    return fem::compute_gradient<D, N>(J, g, g_inv, dphi_dxi_val);
                 });
             // and return it
             return gradient_function;
diff --git a/lib/mito/fem/elements/tri2/api.h b/lib/mito/fem/elements/tri2/api.h
index b7045d0b..58ed3454 100644
--- a/lib/mito/fem/elements/tri2/api.h
+++ b/lib/mito/fem/elements/tri2/api.h
@@ -8,13 +8,7 @@
 
 
 namespace mito::fem {
-
-    // specialization for quadratic shape functions on triangles in 2D
-    template <>
-    struct isoparametric_simplex<2, geometry::triangle_t<2>> {
-        using type = IsoparametricTriangleP2;
-    };
-
+    // no element specializations needed here - see elements_library.h for the unified API
 }
 
 
diff --git a/lib/mito/fem/gradient.h b/lib/mito/fem/gradient.h
new file mode 100644
index 00000000..6223003b
--- /dev/null
+++ b/lib/mito/fem/gradient.h
@@ -0,0 +1,48 @@
+// -*- c++ -*-
+//
+// Copyright (c) 2020-2026, the MiTo Authors, all rights reserved
+//
+
+// code guard
+#pragma once
+
+// DESIGN NOTES
+// This function computes the contravariant gradient of a shape function using the Riemannian
+// gradient formula:
+//   - N<D (embedded): grad(phi) = J * (J^T g J)^{-1} * dphi/dxi
+//   - N=D (square): the above simplifies to grad(phi) = g^{-1} * J^{-T} * dphi/dxi
+
+namespace mito::fem {
+
+    template <int D, int N, class JacobianT, class MetricT, class InvMetricT, class DerivativeT>
+    requires(D >= N && N > 0)
+    constexpr auto compute_gradient(
+        const JacobianT & J, const MetricT & g, const InvMetricT & g_inv,
+        const DerivativeT & dphi_dxi)
+    {
+        if constexpr (N == D) {
+            // square case: no embedding (triangle in 2D, segment in 1D, etc.)
+            auto J_inv = tensor::inverse(J);
+            auto J_invT = tensor::transpose(J_inv);
+            return g_inv * (J_invT * dphi_dxi);
+        } else {
+            // embedded case: N < D (curves in 2D/3D, surfaces in 3D, etc.)
+
+            auto g_induced = tensor::transpose(J) * g * J;    // induced metric
+            auto g_induced_inv = tensor::inverse(g_induced);
+
+            // for N=1, convert scalar to 1-element vector for type compatibility
+            auto dphi = [&]() {
+                if constexpr (N == 1)
+                    return tensor::vector_t<1>{ dphi_dxi };
+                else
+                    return dphi_dxi;
+            }();
+            return J * (g_induced_inv * dphi);
+        }
+    }
+
+}    // namespace mito::fem
+
+
+// end of file
diff --git a/lib/mito/fem/public.h b/lib/mito/fem/public.h
index 8e1a6216..4dbaf1ea 100644
--- a/lib/mito/fem/public.h
+++ b/lib/mito/fem/public.h
@@ -19,6 +19,9 @@
 // blocks implementation
 #include "blocks.h"
 
+// gradient computation utilities
+#include "gradient.h"
+
 // utilities implementation
 #include "utilities.h"
 
diff --git a/lib/mito/manifolds/Manifold.h b/lib/mito/manifolds/Manifold.h
index 55fb6bee..cd1631e7 100644
--- a/lib/mito/manifolds/Manifold.h
+++ b/lib/mito/manifolds/Manifold.h
@@ -16,8 +16,6 @@ namespace mito::manifolds {
       private:
         // typedef for node
         using node_type = cellT::node_type;
-        // the volume form type
-        using volume_form_type = volumeFormT;
         // the physical dimension of the manifold (that is that of the cell)
         static constexpr int D = cellT::dim;
         // the dimension of the manifold (that is that of the cell)
@@ -34,6 +32,8 @@ namespace mito::manifolds {
         using coordinates_type = coordsT;
         // typedef for a coordinates system
         using coordinate_system_type = geometry::coordinate_system_t<coordinates_type>;
+        // the volume form type (public for element type extraction)
+        using volume_form_type = volumeFormT;
 
       public:
         constexpr Manifold(
@@ -73,6 +73,12 @@ namespace mito::manifolds {
             return _coordinate_system;
         }
 
+        // accessor for the volume form
+        constexpr auto volume_form() const noexcept -> const volume_form_type &
+        {
+            return _volume_form;
+        }
+
         constexpr auto elements() const noexcept -> const cells_type & { return _mesh.cells(); }
 
         constexpr auto nElements() const noexcept -> int { return std::size(_mesh.cells()); }
diff --git a/lib/mito/tensor/factories.h b/lib/mito/tensor/factories.h
index 008b86ea..2eded014 100644
--- a/lib/mito/tensor/factories.h
+++ b/lib/mito/tensor/factories.h
@@ -47,6 +47,16 @@ namespace mito::tensor {
             return matrix * vector * v;
         });
     }
+
+    // wrap a vector as a D×1 column matrix
+    // useful e.g. for representing tangent vectors as Jacobian matrices for 1D elements
+    template <int D>
+    constexpr auto as_column_matrix(const vector_t<D> & vec) -> matrix_t<D, 1>
+    {
+        return []<std::size_t... Is>(const auto & v, std::index_sequence<Is...>) {
+            return matrix_t<D, 1>{ v[Is]... };
+        }(vec, std::make_index_sequence<D>{});
+    }
 }
 
 
diff --git a/tests/mito.lib/fem/block_grad_grad.cc b/tests/mito.lib/fem/block_grad_grad.cc
index 1abeab8f..1aff81c4 100644
--- a/tests/mito.lib/fem/block_grad_grad.cc
+++ b/tests/mito.lib/fem/block_grad_grad.cc
@@ -38,9 +38,16 @@ TEST(Fem, IsoparametricTriangle)
     // make a geometric simplex
     auto geometric_simplex = mito::geometry::triangle<2>({ node_0, node_1, node_2 });
 
+    // create a mesh with a single triangle
+    auto mesh = mito::mesh::mesh<cell_t>();
+    mesh.insert({ node_0, node_1, node_2 });
+
+    // create the manifold
+    auto manifold = mito::manifolds::manifold(mesh, coord_system);
+
     {
         // first order isoparametric triangle
-        using element_p1_t = mito::fem::isoparametric_simplex_t<1, cell_t>;
+        using element_p1_t = mito::fem::isoparametric_simplex_t<1, decltype(manifold)>;
 
         // build the discretization nodes
         auto discretization_node_0 = discretization_node_t();
@@ -50,7 +57,8 @@ TEST(Fem, IsoparametricTriangle)
         // a finite element
         auto element_p1 = element_p1_t(
             geometric_simplex, coord_system,
-            { discretization_node_0, discretization_node_1, discretization_node_2 });
+            { discretization_node_0, discretization_node_1, discretization_node_2 },
+            manifold.volume_form());
 
         // a grad-grad matrix block
         auto grad_grad_block =
@@ -72,7 +80,7 @@ TEST(Fem, IsoparametricTriangle)
 
     {
         // second order isoparametric triangle
-        using element_p2_t = mito::fem::isoparametric_simplex_t<2, cell_t>;
+        using element_p2_t = mito::fem::isoparametric_simplex_t<2, decltype(manifold)>;
 
         // build the discretization nodes
         auto discretization_node_0 = discretization_node_t();
@@ -86,7 +94,8 @@ TEST(Fem, IsoparametricTriangle)
         auto element_p2 = element_p2_t(
             geometric_simplex, coord_system,
             { discretization_node_0, discretization_node_1, discretization_node_2,
-              discretization_node_3, discretization_node_4, discretization_node_5 });
+              discretization_node_3, discretization_node_4, discretization_node_5 },
+            manifold.volume_form());
 
         // a grad-grad matrix block
         auto grad_grad_block =
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..d45c8d90
--- /dev/null
+++ b/tests/mito.lib/fem/block_grad_grad_embedded_segment.cc
@@ -0,0 +1,86 @@
+// -*- c++ -*-
+//
+// Copyright (c) 2020-2024, the MiTo Authors, all rights reserved
+//
+
+#include <gtest/gtest.h>
+#include <mito.h>
+
+
+// the type of coordinates (2D physical space for embedded segments)
+using coordinates_t = mito::geometry::coordinates_t<2, mito::geometry::CARTESIAN>;
+// the type of coordinate system
+using coord_system_t = mito::geometry::coordinate_system_t<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 = mito::geometry::reference_segment_t;
+// 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 = coord_system_t();
+
+    // 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 geometric_simplex = mito::geometry::segment<2>({ node_0, node_1 });
+
+    // create a mesh with a single segment
+    auto mesh = mito::mesh::mesh<cell_t>();
+    mesh.insert({ node_0, node_1 });
+
+    // create normal field for the submanifold (perpendicular to diagonal)
+    // 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 });
+
+    // create the submanifold
+    auto manifold = mito::manifolds::submanifold(mesh, coord_system, normal_field);
+
+    {
+        // first order isoparametric embedded segment
+        using element_p1_t = mito::fem::isoparametric_simplex_t<1, decltype(manifold)>;
+
+        // build the discretization nodes
+        auto discretization_node_0 = discretization_node_t();
+        auto discretization_node_1 = discretization_node_t();
+
+        // a finite element
+        auto element_p1 = element_p1_t(
+            geometric_simplex, coord_system, { discretization_node_0, discretization_node_1 },
+            manifold.volume_form());
+
+        // a grad-grad matrix block
+        auto grad_grad_block =
+            mito::fem::blocks::grad_grad_block<element_p1_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;
+}
\ No newline at end of file
diff --git a/tests/mito.lib/fem/block_grad_grad_segment.cc b/tests/mito.lib/fem/block_grad_grad_segment.cc
index 1d6a9553..edf67e04 100644
--- a/tests/mito.lib/fem/block_grad_grad_segment.cc
+++ b/tests/mito.lib/fem/block_grad_grad_segment.cc
@@ -37,9 +37,16 @@ TEST(Fem, BlockGradGradSegment)
     // make a geometric simplex
     auto geometric_simplex = mito::geometry::segment<1>({ node_0, node_1 });
 
+    // create a mesh with a single segment
+    auto mesh = mito::mesh::mesh<cell_t>();
+    mesh.insert({ node_0, node_1 });
+
+    // create the manifold
+    auto manifold = mito::manifolds::manifold(mesh, coord_system);
+
     {
         // first order isoparametric segment
-        using element_p1_t = mito::fem::isoparametric_simplex_t<1, cell_t>;
+        using element_p1_t = mito::fem::isoparametric_simplex_t<1, decltype(manifold)>;
 
         // build the discretization nodes
         auto discretization_node_0 = discretization_node_t();
@@ -47,7 +54,8 @@ TEST(Fem, BlockGradGradSegment)
 
         // a finite element
         auto element_p1 = element_p1_t(
-            geometric_simplex, coord_system, { discretization_node_0, discretization_node_1 });
+            geometric_simplex, coord_system, { discretization_node_0, discretization_node_1 },
+            manifold.volume_form());
 
         // a grad-grad matrix block
         auto grad_grad_block =
diff --git a/tests/mito.lib/fem/block_mass.cc b/tests/mito.lib/fem/block_mass.cc
index 06d19987..e4543c25 100644
--- a/tests/mito.lib/fem/block_mass.cc
+++ b/tests/mito.lib/fem/block_mass.cc
@@ -38,9 +38,16 @@ TEST(Fem, IsoparametricTriangle)
     // make a geometric simplex
     auto geometric_simplex = mito::geometry::triangle<2>({ node_0, node_1, node_2 });
 
+    // create a mesh with a single triangle
+    auto mesh = mito::mesh::mesh<cell_t>();
+    mesh.insert({ node_0, node_1, node_2 });
+
+    // create the manifold
+    auto manifold = mito::manifolds::manifold(mesh, coord_system);
+
     {
         // first order isoparametric triangle
-        using element_p1_t = mito::fem::isoparametric_simplex_t<1, cell_t>;
+        using element_p1_t = mito::fem::isoparametric_simplex_t<1, decltype(manifold)>;
 
         // build the discretization nodes
         auto discretization_node_0 = discretization_node_t();
@@ -50,7 +57,8 @@ TEST(Fem, IsoparametricTriangle)
         // a finite element
         auto element_p1 = element_p1_t(
             geometric_simplex, coord_system,
-            { discretization_node_0, discretization_node_1, discretization_node_2 });
+            { discretization_node_0, discretization_node_1, discretization_node_2 },
+            manifold.volume_form());
 
         // a mass matrix block
         auto mass_block = mito::fem::blocks::mass_block<element_p1_t, quadrature_rule_t>();
@@ -71,7 +79,7 @@ TEST(Fem, IsoparametricTriangle)
 
     {
         // second order isoparametric triangle
-        using element_p2_t = mito::fem::isoparametric_simplex_t<2, cell_t>;
+        using element_p2_t = mito::fem::isoparametric_simplex_t<2, decltype(manifold)>;
 
         // build the discretization nodes
         auto discretization_node_0 = discretization_node_t();
@@ -85,7 +93,8 @@ TEST(Fem, IsoparametricTriangle)
         auto element_p2 = element_p2_t(
             geometric_simplex, coord_system,
             { discretization_node_0, discretization_node_1, discretization_node_2,
-              discretization_node_3, discretization_node_4, discretization_node_5 });
+              discretization_node_3, discretization_node_4, discretization_node_5 },
+            manifold.volume_form());
 
         // a mass matrix block
         auto mass_block = mito::fem::blocks::mass_block<element_p2_t, quadrature_rule_t>();
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..c5f4082f
--- /dev/null
+++ b/tests/mito.lib/fem/block_mass_embedded_segment.cc
@@ -0,0 +1,83 @@
+// -*- c++ -*-
+//
+// Copyright (c) 2020-2024, the MiTo Authors, all rights reserved
+//
+
+#include <gtest/gtest.h>
+#include <mito.h>
+
+
+// the type of coordinates (2D physical space for embedded segments)
+using coordinates_t = mito::geometry::coordinates_t<2, mito::geometry::CARTESIAN>;
+// the type of coordinate system
+using coord_system_t = mito::geometry::coordinate_system_t<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 = mito::geometry::reference_segment_t;
+// 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 = coord_system_t();
+
+    // 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 geometric_simplex = mito::geometry::segment<2>({ node_0, node_1 });
+
+    // create a mesh with the single segment
+    auto mesh = mito::mesh::mesh<cell_t>();
+    mesh.insert({ node_0, node_1 });
+
+    // create a submanifold (1D embedded in 2D requires normal field)
+    // 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 });
+    auto manifold = mito::manifolds::submanifold(mesh, coord_system, normal_field);
+
+    {
+        // first order isoparametric embedded segment
+        using element_p1_t = mito::fem::isoparametric_simplex_t<1, decltype(manifold)>;
+
+        // build the discretization nodes
+        auto discretization_node_0 = discretization_node_t();
+        auto discretization_node_1 = discretization_node_t();
+
+        // a finite element
+        auto element_p1 = element_p1_t(
+            geometric_simplex, coord_system, { discretization_node_0, discretization_node_1 },
+            manifold.volume_form());
+
+        // a mass matrix block
+        auto mass_block = mito::fem::blocks::mass_block<element_p1_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_segment.cc b/tests/mito.lib/fem/block_mass_segment.cc
index 64602992..88e4dff9 100644
--- a/tests/mito.lib/fem/block_mass_segment.cc
+++ b/tests/mito.lib/fem/block_mass_segment.cc
@@ -37,9 +37,16 @@ TEST(Fem, BlockMassSegment)
     // make a geometric simplex
     auto geometric_simplex = mito::geometry::segment<1>({ node_0, node_1 });
 
+    // create a mesh with the single segment
+    auto mesh = mito::mesh::mesh<cell_t>();
+    mesh.insert({ node_0, node_1 });
+
+    // create a manifold
+    auto manifold = mito::manifolds::manifold(mesh, coord_system);
+
     {
         // first order isoparametric segment
-        using element_p1_t = mito::fem::isoparametric_simplex_t<1, cell_t>;
+        using element_p1_t = mito::fem::isoparametric_simplex_t<1, decltype(manifold)>;
 
         // build the discretization nodes
         auto discretization_node_0 = discretization_node_t();
@@ -47,7 +54,8 @@ TEST(Fem, BlockMassSegment)
 
         // a finite element
         auto element_p1 = element_p1_t(
-            geometric_simplex, coord_system, { discretization_node_0, discretization_node_1 });
+            geometric_simplex, coord_system, { discretization_node_0, discretization_node_1 },
+            manifold.volume_form());
 
         // a mass matrix block
         auto mass_block = mito::fem::blocks::mass_block<element_p1_t, quadrature_rule_t>();
diff --git a/tests/mito.lib/fem/fem_field.cc b/tests/mito.lib/fem/fem_field.cc
index 99121425..8703e219 100644
--- a/tests/mito.lib/fem/fem_field.cc
+++ b/tests/mito.lib/fem/fem_field.cc
@@ -12,10 +12,6 @@ using coordinates_t = mito::geometry::coordinates_t<2, mito::geometry::CARTESIAN
 
 // simplicial cells in 2D
 using cell_t = mito::geometry::triangle_t<2>;
-// first degree finite elements
-constexpr int degree = 1;
-// assemble the finite element type
-using finite_element_t = mito::fem::isoparametric_simplex_t<degree, cell_t>;
 // the x scalar field in 2D
 constexpr auto x = mito::functions::component<coordinates_t, 0>;
 // the y scalar field in 2D
@@ -34,6 +30,11 @@ TEST(Fem, FemField)
     // create the body manifold
     auto manifold = mito::manifolds::manifold(mesh, coord_system);
 
+    // first degree finite elements
+    constexpr int degree = 1;
+    // assemble the finite element type
+    using finite_element_t = mito::fem::isoparametric_simplex_t<degree, decltype(manifold)>;
+
     // TOFIX: it should not be mandatory to set constraints to create a function space, let's remove
     // this bit once we implement constraints properly
     // get the boundary mesh
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..ae3a22d1
--- /dev/null
+++ b/tests/mito.lib/fem/isoparametric_embedded_segment.cc
@@ -0,0 +1,179 @@
+// -*- c++ -*-
+//
+// Copyright (c) 2020-2024, the MiTo Authors, all rights reserved
+//
+
+#include <gtest/gtest.h>
+#include <mito.h>
+
+
+// the type of coordinates (2D physical space for embedded segments)
+using coordinates_t = mito::geometry::coordinates_t<2, mito::geometry::CARTESIAN>;
+// the type of coordinate system
+using coord_system_t = mito::geometry::coordinate_system_t<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 = mito::geometry::reference_segment_t;
+// 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 the sum of the shape functions gradients
+        EXPECT_NEAR(0.0, sum[0], 3.0e-16);
+        EXPECT_NEAR(0.0, sum[1], 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::area * 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 = coord_system_t();
+
+    // 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 geometric_simplex = mito::geometry::segment<2>({ node_0, node_1 });
+
+    // create a mesh with a single segment
+    auto mesh = mito::mesh::mesh<cell_t>();
+    mesh.insert({ node_0, node_1 });
+
+    // create normal field for the submanifold (perpendicular to 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 });
+
+    // create the submanifold
+    auto manifold = mito::manifolds::submanifold(mesh, coord_system, normal_field);
+
+    // first order isoparametric embedded segment
+    using element_p1_t = mito::fem::isoparametric_simplex_t<1, decltype(manifold)>;
+
+    // build the discretization nodes
+    auto discretization_node_0 = discretization_node_t();
+    auto discretization_node_1 = discretization_node_t();
+
+    // a finite element
+    auto element_p1 = element_p1_t(
+        geometric_simplex, coord_system, { discretization_node_0, discretization_node_1 },
+        manifold.volume_form());
+
+    // 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 = mito::geometry::reference_segment_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_segment.cc b/tests/mito.lib/fem/isoparametric_segment.cc
index 257f3456..a0080a83 100644
--- a/tests/mito.lib/fem/isoparametric_segment.cc
+++ b/tests/mito.lib/fem/isoparametric_segment.cc
@@ -99,9 +99,16 @@ TEST(Fem, IsoparametricSegment)
     // make a geometric simplex
     auto geometric_simplex = mito::geometry::segment<1>({ node_0, node_1 });
 
+    // create a mesh with a single segment
+    auto mesh = mito::mesh::mesh<cell_t>();
+    mesh.insert({ node_0, node_1 });
+
+    // create the manifold
+    auto manifold = mito::manifolds::manifold(mesh, coord_system);
+
     {
         // first order isoparametric segment
-        using element_p1_t = mito::fem::isoparametric_simplex_t<1, cell_t>;
+        using element_p1_t = mito::fem::isoparametric_simplex_t<1, decltype(manifold)>;
 
         // build the discretization nodes
         auto discretization_node_0 = discretization_node_t();
@@ -109,7 +116,8 @@ TEST(Fem, IsoparametricSegment)
 
         // a finite element
         auto element_p1 = element_p1_t(
-            geometric_simplex, coord_system, { discretization_node_0, discretization_node_1 });
+            geometric_simplex, coord_system, { discretization_node_0, discretization_node_1 },
+            manifold.volume_form());
 
         // check that first order shape functions are a partition of unity
         test_partition_of_unity(element_p1);
diff --git a/tests/mito.lib/fem/isoparametric_triangle.cc b/tests/mito.lib/fem/isoparametric_triangle.cc
index ce16d413..818fd426 100644
--- a/tests/mito.lib/fem/isoparametric_triangle.cc
+++ b/tests/mito.lib/fem/isoparametric_triangle.cc
@@ -101,9 +101,16 @@ TEST(Fem, IsoparametricTriangle)
     // make a geometric simplex
     auto geometric_simplex = mito::geometry::triangle<2>({ node_0, node_1, node_2 });
 
+    // create a mesh with a single triangle
+    auto mesh = mito::mesh::mesh<cell_t>();
+    mesh.insert({ node_0, node_1, node_2 });
+
+    // create the manifold
+    auto manifold = mito::manifolds::manifold(mesh, coord_system);
+
     {
         // first order isoparametric triangle
-        using element_p1_t = mito::fem::isoparametric_simplex_t<1, cell_t>;
+        using element_p1_t = mito::fem::isoparametric_simplex_t<1, decltype(manifold)>;
 
         // build the discretization nodes
         auto discretization_node_0 = discretization_node_t();
@@ -113,7 +120,8 @@ TEST(Fem, IsoparametricTriangle)
         // a finite element
         auto element_p1 = element_p1_t(
             geometric_simplex, coord_system,
-            { discretization_node_0, discretization_node_1, discretization_node_2 });
+            { discretization_node_0, discretization_node_1, discretization_node_2 },
+            manifold.volume_form());
 
         // check that first order shape functions are a partition of unity
         test_partition_of_unity(element_p1);
@@ -124,7 +132,7 @@ TEST(Fem, IsoparametricTriangle)
 
     {
         // second order isoparametric triangle
-        using element_p2_t = mito::fem::isoparametric_simplex_t<2, cell_t>;
+        using element_p2_t = mito::fem::isoparametric_simplex_t<2, decltype(manifold)>;
 
         // build the discretization nodes
         auto discretization_node_0 = discretization_node_t();
@@ -138,7 +146,8 @@ TEST(Fem, IsoparametricTriangle)
         auto element_p2 = element_p2_t(
             geometric_simplex, coord_system,
             { discretization_node_0, discretization_node_1, discretization_node_2,
-              discretization_node_3, discretization_node_4, discretization_node_5 });
+              discretization_node_3, discretization_node_4, discretization_node_5 },
+            manifold.volume_form());
 
         // check that second order shape functions are a partition of unity
         test_partition_of_unity(element_p2);
diff --git a/tests/mito.lib/fem/localize_field.cc b/tests/mito.lib/fem/localize_field.cc
index e1321821..44037af1 100644
--- a/tests/mito.lib/fem/localize_field.cc
+++ b/tests/mito.lib/fem/localize_field.cc
@@ -11,6 +11,8 @@
 using coordinates_t = mito::geometry::coordinates_t<2, mito::geometry::CARTESIAN>;
 // the type of coordinate system
 using coord_system_t = mito::geometry::coordinate_system_t<coordinates_t>;
+// the type of cell
+using cell_t = mito::geometry::triangle_t<2>;
 // the x scalar field in 2D
 constexpr auto x = mito::functions::component<coordinates_t, 0>;
 // the y scalar field in 2D
@@ -36,8 +38,16 @@ TEST(Fem, LocalizeField)
     // create a geometric simplex
     auto geometric_simplex = mito::geometry::triangle<2>({ node_0, node_1, node_2 });
 
+    // create a mesh with a single triangle
+    auto mesh = mito::mesh::mesh<cell_t>();
+    mesh.insert({ node_0, node_1, node_2 });
+
+    // create the manifold
+    auto manifold = mito::manifolds::manifold(mesh, coord_system);
+
     // an isoparametric triangle
-    auto element = mito::fem::IsoparametricTriangle(geometric_simplex, coord_system);
+    auto element =
+        mito::fem::IsoparametricTriangle(geometric_simplex, coord_system, manifold.volume_form());
 
     // localize the field on the simplex
     auto localized_field = mito::fem::localize(field, element);
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..36c8b979
--- /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 vectors
+    auto director0 = x1 - x0;
+
+    // compute volume of triangle
+    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 });
+
+    // the normal vector to the submanifold
+    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