From fd8cf0481873b2c970bda353c4a14df458229aef Mon Sep 17 00:00:00 2001 From: Pawel024 Date: Thu, 15 Jan 2026 18:10:54 +0100 Subject: [PATCH 01/24] fem: generalized jacobian measure using the area formula --- lib/mito/fem/blocks/GradGradBlock.h | 2 +- lib/mito/fem/blocks/L2NormBlock.h | 2 +- lib/mito/fem/blocks/MassBlock.h | 2 +- lib/mito/fem/blocks/SourceTermBlock.h | 2 +- lib/mito/fem/blocks/public.h | 3 +++ lib/mito/fem/blocks/utilities.h | 36 +++++++++++++++++++++++++++ 6 files changed, 43 insertions(+), 4 deletions(-) create mode 100644 lib/mito/fem/blocks/utilities.h diff --git a/lib/mito/fem/blocks/GradGradBlock.h b/lib/mito/fem/blocks/GradGradBlock.h index 65d8eb2d..277e2610 100644 --- a/lib/mito/fem/blocks/GradGradBlock.h +++ b/lib/mito/fem/blocks/GradGradBlock.h @@ -45,7 +45,7 @@ namespace mito::fem::blocks { element_type::canonical_element_type::area * quadrature_rule.weight(q); // precompute the common factor - auto factor = w * tensor::determinant(element.jacobian()(xi)); + auto factor = w * jacobian_measure(element.jacobian()(xi)); // loop on the nodes of the element tensor::constexpr_for_1([&]() { diff --git a/lib/mito/fem/blocks/L2NormBlock.h b/lib/mito/fem/blocks/L2NormBlock.h index 0d0a00f4..b5f86732 100644 --- a/lib/mito/fem/blocks/L2NormBlock.h +++ b/lib/mito/fem/blocks/L2NormBlock.h @@ -52,7 +52,7 @@ namespace mito::fem::blocks { element_type::canonical_element_type::area * quadrature_rule.weight(q); // precompute the common factor - auto factor = w * tensor::determinant(element.jacobian()(xi)); + auto factor = w * jacobian_measure(element.jacobian()(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..8643efef 100644 --- a/lib/mito/fem/blocks/MassBlock.h +++ b/lib/mito/fem/blocks/MassBlock.h @@ -45,7 +45,7 @@ namespace mito::fem::blocks { element_type::canonical_element_type::area * quadrature_rule.weight(q); // precompute the common factor - auto factor = w * tensor::determinant(element.jacobian()(xi)); + auto factor = w * jacobian_measure(element.jacobian()(xi)); // loop on the nodes of the element tensor::constexpr_for_1([&]() { diff --git a/lib/mito/fem/blocks/SourceTermBlock.h b/lib/mito/fem/blocks/SourceTermBlock.h index ed9ecb8e..13914a96 100644 --- a/lib/mito/fem/blocks/SourceTermBlock.h +++ b/lib/mito/fem/blocks/SourceTermBlock.h @@ -57,7 +57,7 @@ namespace mito::fem::blocks { element_type::canonical_element_type::area * quadrature_rule.weight(q); // precompute the common factor - auto factor = w * tensor::determinant(element.jacobian()(xi)); + auto factor = w * jacobian_measure(element.jacobian()(xi)); // loop on the nodes of the element tensor::constexpr_for_1([&]() { diff --git a/lib/mito/fem/blocks/public.h b/lib/mito/fem/blocks/public.h index 20b5b6e9..d202940c 100644 --- a/lib/mito/fem/blocks/public.h +++ b/lib/mito/fem/blocks/public.h @@ -16,6 +16,9 @@ // published types factories; this is the file you are looking for... #include "api.h" +// utilities (must be included before blocks that use jacobian_measure) +#include "utilities.h" + // classes implementation #include "AssemblyBlock.h" #include "GradGradBlock.h" diff --git a/lib/mito/fem/blocks/utilities.h b/lib/mito/fem/blocks/utilities.h new file mode 100644 index 00000000..c46ede26 --- /dev/null +++ b/lib/mito/fem/blocks/utilities.h @@ -0,0 +1,36 @@ +// -*- c++ -*- +// +// Copyright (c) 2020-2024, the MiTo Authors, all rights reserved +// + +// code guard +#pragma once + + +namespace mito::fem::blocks { + + // Helper function to compute the Jacobian measure (integration scaling factor) + // for different Jacobian matrix shapes. + + // For square Jacobians (NxN): return |det(J)| + template + requires(tensor::square_matrix_c) + constexpr auto jacobian_measure(const matrixT & J) + { + return std::abs(tensor::determinant(J)); + } + + // For 2x1 Jacobians (embedded 1D in 2D): return sqrt(det(J^T * J)) + // This is the square root of the Gram determinant as used in the area formula + template + requires(tensor::matrix_c && matrixT::dims[0] == 2 && matrixT::dims[1] == 1) + constexpr auto jacobian_measure(const matrixT & J) + { + auto JtJ = tensor::transpose(J) * J; // 1x1 matrix + return std::sqrt(JtJ[{ 0, 0 }]); + } + +} // namespace mito + + +// end of file From 3995dab7b814276bf1b525464f3f98a5d82f71d3 Mon Sep 17 00:00:00 2001 From: Pawel024 Date: Sat, 17 Jan 2026 22:42:50 +0100 Subject: [PATCH 02/24] fem: completed base implementation of segments embedded in 2d --- .../elements/IsoparametricEmbeddedSegment.h | 75 ++++++++++ lib/mito/fem/elements/elements_library.h | 1 + lib/mito/fem/elements/public.h | 1 + .../elements/seg1_embedded/DiscretizerCG.h | 65 ++++++++ .../IsoparametricEmbeddedSegmentP1.h | 140 ++++++++++++++++++ lib/mito/fem/elements/seg1_embedded/api.h | 21 +++ lib/mito/fem/elements/seg1_embedded/public.h | 18 +++ 7 files changed, 321 insertions(+) create mode 100644 lib/mito/fem/elements/IsoparametricEmbeddedSegment.h create mode 100644 lib/mito/fem/elements/seg1_embedded/DiscretizerCG.h create mode 100644 lib/mito/fem/elements/seg1_embedded/IsoparametricEmbeddedSegmentP1.h create mode 100644 lib/mito/fem/elements/seg1_embedded/api.h create mode 100644 lib/mito/fem/elements/seg1_embedded/public.h diff --git a/lib/mito/fem/elements/IsoparametricEmbeddedSegment.h b/lib/mito/fem/elements/IsoparametricEmbeddedSegment.h new file mode 100644 index 00000000..0c7bbce7 --- /dev/null +++ b/lib/mito/fem/elements/IsoparametricEmbeddedSegment.h @@ -0,0 +1,75 @@ +// -*- c++ -*- +// +// Copyright (c) 2020-2024, the MiTo Authors, all rights reserved +// + +// code guard +#pragma once + + +// DESIGN NOTES +// Class {IsoparametricEmbeddedSegment} represents a first order simplex (segment) embedded in 2D +// space, equipped with parametric coordinates. + + +namespace mito::fem { + + class IsoparametricEmbeddedSegment : public utilities::Invalidatable { + public: + // the dimension of the physical space + static constexpr int dim = 2; + // the discretization node type + using discretization_node_type = discrete::discretization_node_t; + // the underlying cell type + using cell_type = geometry::segment_t; + + protected: + // cartesian coordinates in 2D + using coordinates_type = geometry::coordinates_t<2, geometry::CARTESIAN>; + // the coordinate system type + using coordinate_system_type = geometry::coordinate_system_t; + // the vector type + using vector_type = tensor::vector_t<2>; + + public: + // the default constructor + constexpr IsoparametricEmbeddedSegment( + const cell_type & cell, const coordinate_system_type & coord_system) : + _cell(cell), + _x0{ coord_system.coordinates(cell.nodes()[0]->point()) - coordinates_type{} }, + _x1{ coord_system.coordinates(cell.nodes()[1]->point()) - coordinates_type{} } + {} + + // destructor + constexpr ~IsoparametricEmbeddedSegment() = default; + + // delete move constructor + constexpr IsoparametricEmbeddedSegment(IsoparametricEmbeddedSegment &&) noexcept = delete; + + // delete copy constructor + constexpr IsoparametricEmbeddedSegment(const IsoparametricEmbeddedSegment &) = delete; + + // delete assignment operator + constexpr IsoparametricEmbeddedSegment & operator=(const IsoparametricEmbeddedSegment &) = delete; + + // delete move assignment operator + constexpr IsoparametricEmbeddedSegment & operator=(IsoparametricEmbeddedSegment &&) noexcept = delete; + + public: + // get the geometric simplex + constexpr auto cell() const noexcept -> const cell_type & { return _cell; } + + protected: + // QUESTION: do we need to maintain a reference to the geometric simplex? + // a const reference to the geometric simplex + const cell_type & _cell; + + // the coordinates of the discretization nodes of the segment + const vector_type _x0; + const vector_type _x1; + }; + +} // namespace mito + + +// end of file diff --git a/lib/mito/fem/elements/elements_library.h b/lib/mito/fem/elements/elements_library.h index 60e59338..077f1a98 100644 --- a/lib/mito/fem/elements/elements_library.h +++ b/lib/mito/fem/elements/elements_library.h @@ -8,6 +8,7 @@ #include "seg1/public.h" +#include "seg1_embedded/public.h" #include "tri1/public.h" #include "tri2/public.h" diff --git a/lib/mito/fem/elements/public.h b/lib/mito/fem/elements/public.h index 7ed2002a..633f46f0 100644 --- a/lib/mito/fem/elements/public.h +++ b/lib/mito/fem/elements/public.h @@ -19,6 +19,7 @@ // classes implementation #include "IsoparametricSegment.h" #include "IsoparametricTriangle.h" +#include "IsoparametricEmbeddedSegment.h" #include "Discretizer.h" // library of finite elements diff --git a/lib/mito/fem/elements/seg1_embedded/DiscretizerCG.h b/lib/mito/fem/elements/seg1_embedded/DiscretizerCG.h new file mode 100644 index 00000000..fdceca3c --- /dev/null +++ b/lib/mito/fem/elements/seg1_embedded/DiscretizerCG.h @@ -0,0 +1,65 @@ +// -*- c++ -*- +// +// Copyright (c) 2020-2024, the MiTo Authors, all rights reserved +// + +// code guard +#pragma once + + +namespace mito::fem { + + // discretizer specialization for {IsoparametricEmbeddedSegmentP1} with continuous Galerkin + template <> + struct Discretizer { + template < + typename manifoldT, typename constraintsT, typename elements_type, typename map_type, + typename constrained_nodes_type> + static void apply( + const manifoldT & manifold, const constraintsT & constraints, elements_type & elements, + map_type & node_map, constrained_nodes_type & constrained_nodes) + { + + // the discretization node type + using discretization_node_type = + typename IsoparametricEmbeddedSegmentP1::discretization_node_type; + + // the connectivity type + using connectivity_type = typename IsoparametricEmbeddedSegmentP1::connectivity_type; + + // get the coordinate system of the manifold + const auto & coord_system = manifold.coordinate_system(); + + // loop on the cells of the mesh + for (const auto & cell : manifold.elements()) { + + // get the nodes of the cell + const auto & nodes = cell.nodes(); + + // add the nodes to the map (if the mesh node is already present in the map, + // then the present discretization node is used) + auto node_0 = + node_map.insert({ nodes[0], discretization_node_type() }).first->second; + 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 }); + } + + // populate the constrained nodes + // In 1D, constraints.domain() is a set of nodes, not a mesh with cells, so we loop on the nodes directly + for (const auto & node : constraints.domain()) { + // get the discretization node associated with the mesh node from the map + auto it = node_map.find(node); + // add the node to the constrained nodes + constrained_nodes.insert(it->second); + } + + // all done + return; + } + }; +} + +// end of file diff --git a/lib/mito/fem/elements/seg1_embedded/IsoparametricEmbeddedSegmentP1.h b/lib/mito/fem/elements/seg1_embedded/IsoparametricEmbeddedSegmentP1.h new file mode 100644 index 00000000..71529598 --- /dev/null +++ b/lib/mito/fem/elements/seg1_embedded/IsoparametricEmbeddedSegmentP1.h @@ -0,0 +1,140 @@ +// -*- c++ -*- +// +// Copyright (c) 2020-2024, the MiTo Authors, all rights reserved +// + +// code guard +#pragma once + + +// DESIGN NOTES +// Class {IsoparametricEmbeddedSegmentP1} represents a first order simplex (segment) living in 2D +// cartesian space, equipped with linear shape functions defined in the 1D parametric space. +// The key difference from IsoparametricSegmentP1 is the non-square Jacobian (2x1 matrix) and +// the gradient computation that accounts for the embedded geometry. + + +namespace mito::fem { + + class IsoparametricEmbeddedSegmentP1 : public IsoparametricEmbeddedSegment { + + public: + // the degree of the finite element + static constexpr int degree = 1; + // the type of shape functions (reuse from seg1, same 1D parametric shape functions) + using shape_functions_type = ShapeSegmentP1; + // the canonical element type + using canonical_element_type = typename shape_functions_type::reference_element_type; + // the parametric coordinates type + using parametric_coordinates_type = + typename canonical_element_type::parametric_coordinates_type; + // the linear 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; + + public: + // the default constructor + inline IsoparametricEmbeddedSegmentP1( + const cell_type & geometric_simplex, const coordinate_system_type & coord_system, + const connectivity_type & connectivity) : + IsoparametricEmbeddedSegment(geometric_simplex, coord_system), + _connectivity(connectivity) + {} + + // destructor + inline ~IsoparametricEmbeddedSegmentP1() = default; + + // delete move constructor + constexpr IsoparametricEmbeddedSegmentP1(IsoparametricEmbeddedSegmentP1 &&) noexcept = delete; + + // delete copy constructor + constexpr IsoparametricEmbeddedSegmentP1(const IsoparametricEmbeddedSegmentP1 &) = delete; + + // delete assignment operator + constexpr IsoparametricEmbeddedSegmentP1 & operator=(const IsoparametricEmbeddedSegmentP1 &) = delete; + + // delete move assignment operator + constexpr IsoparametricEmbeddedSegmentP1 & operator=(IsoparametricEmbeddedSegmentP1 &&) noexcept = delete; + + public: + // get the discretization nodes + constexpr auto connectivity() const noexcept -> const connectivity_type & + { + return _connectivity; + } + + // get the isoparametric mapping from parametric coordinates to physical coordinates + constexpr auto parametrization() const + { + // get the shape functions + constexpr auto phi_0 = shape_functions.shape<0>(); + constexpr auto phi_1 = shape_functions.shape<1>(); + + // return the isoparametric mapping from parametric to physical coordinates + return mito::fields::linear_combination( + std::array{ _x0, _x1 }, phi_0, phi_1); + } + + // get the shape function associated with local node {a} + template + requires(a >= 0 && a < n_nodes) + constexpr auto shape() const + { + // return the shape functions + return shape_functions.shape(); + } + + // get the jacobian of the isoparametric mapping from parametric to actual coordinates + 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<2, 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 2x1 matrix + return tensor::matrix_t<2, 1>{ dx_dxi[0], dx_dxi[1] }; + }); + + // and return it + return jacobian_function; + } + + // get the gradient of the a-th shape function as a function of parametric coordinates + template + requires(a >= 0 && a < n_nodes) + constexpr auto gradient() const + { + // assemble the gradient as a function of parametric coordinates + auto gradient_function = + fields::field([&](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); + // for the embedded case, the gradient is: grad(phi) = dphi/dxi * J / |J|^2 + // this follows from: dphi/dx = grad(phi)^T * J, with grad(phi) parallel to J + auto J_norm_squared = J[{0, 0}] * J[{0, 0}] + J[{1, 0}] * J[{1, 0}]; + auto J_scaled = tensor::vector_t<2>{ J[{0, 0}] / J_norm_squared, J[{1, 0}] / J_norm_squared }; + // return the spatial gradients of the shape functions evaluated at {xi} + return shape_functions.dshape()(xi) * J_scaled; + }); + // and return it + return gradient_function; + } + + private: + // the discretization nodes of the simplex + const connectivity_type _connectivity; + }; + +} // namespace mito + + +// end of file diff --git a/lib/mito/fem/elements/seg1_embedded/api.h b/lib/mito/fem/elements/seg1_embedded/api.h new file mode 100644 index 00000000..30edb575 --- /dev/null +++ b/lib/mito/fem/elements/seg1_embedded/api.h @@ -0,0 +1,21 @@ +// -*- c++ -*- +// +// Copyright (c) 2020-2024, the MiTo Authors, all rights reserved +// + +// code guard +#pragma once + + +namespace mito::fem { + + // specialization for linear shape functions on segments embedded in 2D + template <> + struct isoparametric_simplex<1, geometry::segment_t<2>> { + using type = IsoparametricEmbeddedSegmentP1; + }; + +} + + +// end of file diff --git a/lib/mito/fem/elements/seg1_embedded/public.h b/lib/mito/fem/elements/seg1_embedded/public.h new file mode 100644 index 00000000..1a64d4cd --- /dev/null +++ b/lib/mito/fem/elements/seg1_embedded/public.h @@ -0,0 +1,18 @@ +// -*- c++ -*- +// +// Copyright (c) 2020-2024, the MiTo Authors, all rights reserved +// + +// code guard +#pragma once + + +// classes implementation (reuses ShapeSegmentP1 from seg1) +#include "IsoparametricEmbeddedSegmentP1.h" +#include "DiscretizerCG.h" + +// published types and factories +#include "api.h" + + +// end of file From 5ae76d1e32d71e6f39e70321d7feede5771cbb86 Mon Sep 17 00:00:00 2001 From: Pawel024 Date: Tue, 20 Jan 2026 13:47:34 +0100 Subject: [PATCH 03/24] fem: copy block_mass and block_grad_grad tests for embedded segments --- .cmake/mito_tests_mito_lib.cmake | 2 + .../fem/block_grad_grad_embedded_segment.cc | 72 +++++++++++++++++++ .../fem/block_mass_embedded_segment.cc | 72 +++++++++++++++++++ 3 files changed, 146 insertions(+) create mode 100644 tests/mito.lib/fem/block_grad_grad_embedded_segment.cc create mode 100644 tests/mito.lib/fem/block_mass_embedded_segment.cc diff --git a/.cmake/mito_tests_mito_lib.cmake b/.cmake/mito_tests_mito_lib.cmake index 70da8119..dd40ae94 100644 --- a/.cmake/mito_tests_mito_lib.cmake +++ b/.cmake/mito_tests_mito_lib.cmake @@ -47,8 +47,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) 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..e811827e --- /dev/null +++ b/tests/mito.lib/fem/block_grad_grad_embedded_segment.cc @@ -0,0 +1,72 @@ +// -*- c++ -*- +// +// Copyright (c) 2020-2024, the MiTo Authors, all rights reserved +// + +#include +#include + + +// 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; +// 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; + +// 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 horizontal segment embedded in 2D) + auto node_0 = mito::geometry::node(coord_system, { 0.0, 0.0 }); + auto node_1 = mito::geometry::node(coord_system, { 1.0, 0.0 }); + + // make a geometric simplex + auto geometric_simplex = mito::geometry::segment<2>({ node_0, node_1 }); + + { + // first order isoparametric embedded segment + using element_p1_t = mito::fem::isoparametric_simplex_t<1, cell_t>; + + // 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 }); + + // a grad-grad matrix block + auto grad_grad_block = + mito::fem::blocks::grad_grad_block(); + + // 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_mass_embedded_segment.cc b/tests/mito.lib/fem/block_mass_embedded_segment.cc new file mode 100644 index 00000000..0b7da468 --- /dev/null +++ b/tests/mito.lib/fem/block_mass_embedded_segment.cc @@ -0,0 +1,72 @@ +// -*- c++ -*- +// +// Copyright (c) 2020-2024, the MiTo Authors, all rights reserved +// + +#include +#include + + +// 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; +// 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; + +// 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 horizontal segment embedded in 2D) + auto node_0 = mito::geometry::node(coord_system, { 0.0, 0.0 }); + auto node_1 = mito::geometry::node(coord_system, { 1.0, 0.0 }); + + // make a geometric simplex + auto geometric_simplex = mito::geometry::segment<2>({ node_0, node_1 }); + + { + // first order isoparametric embedded segment + using element_p1_t = mito::fem::isoparametric_simplex_t<1, cell_t>; + + // 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 }); + + // a mass matrix block + auto mass_block = mito::fem::blocks::mass_block(); + + // 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; +} From e30e749a9088f355004148c31e843ec4bf9548b5 Mon Sep 17 00:00:00 2001 From: Pawel024 Date: Tue, 20 Jan 2026 13:53:35 +0100 Subject: [PATCH 04/24] fem: copy tests from isoparametric_segment.cc for embedded segments --- .cmake/mito_tests_mito_lib.cmake | 1 + .../fem/isoparametric_embedded_segment.cc | 125 ++++++++++++++++++ 2 files changed, 126 insertions(+) create mode 100644 tests/mito.lib/fem/isoparametric_embedded_segment.cc diff --git a/.cmake/mito_tests_mito_lib.cmake b/.cmake/mito_tests_mito_lib.cmake index dd40ae94..3e019667 100644 --- a/.cmake/mito_tests_mito_lib.cmake +++ b/.cmake/mito_tests_mito_lib.cmake @@ -59,6 +59,7 @@ mito_test_driver(tests/mito.lib/fem/domain_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/tests/mito.lib/fem/isoparametric_embedded_segment.cc b/tests/mito.lib/fem/isoparametric_embedded_segment.cc new file mode 100644 index 00000000..2d849726 --- /dev/null +++ b/tests/mito.lib/fem/isoparametric_embedded_segment.cc @@ -0,0 +1,125 @@ +// -*- c++ -*- +// +// Copyright (c) 2020-2024, the MiTo Authors, all rights reserved +// + +#include +#include + + +// the type of coordinates +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; +// the type of discretization node +using discretization_node_t = mito::discrete::discretization_node_t; +// the type of cell +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; + + +// 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::n_nodes; + + // loop on the quadrature points + mito::tensor::constexpr_for_1([&]() { + // 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 = + ([]( + const auto & element, const auto & xi, mito::tensor::integer_sequence) { + return ((element.template shape()(xi)) + ...); + })(element, xi, mito::tensor::make_integer_sequence{}); + + // 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::n_nodes; + + // loop on the quadrature points + mito::tensor::constexpr_for_1([&]() { + // 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 = + ([]( + const auto & element, const auto & xi, mito::tensor::integer_sequence) { + return ((element.template gradient()(xi)) + ...); + })(element, xi, mito::tensor::make_integer_sequence{}); + + // 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(Fem, IsoparametricEmbeddedSegment) +{ + // the coordinate system + auto coord_system = coord_system_t(); + + // build nodes + auto node_0 = mito::geometry::node(coord_system, { 0.0, 0.0 }); + auto node_1 = mito::geometry::node(coord_system, { 1.0, 0.0 }); + + // make a geometric simplex + auto geometric_simplex = mito::geometry::segment<2>({ node_0, node_1 }); + + { + // first order isoparametric segment + using element_p1_t = mito::fem::isoparametric_simplex_t<1, cell_t>; + + // 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 }); + + // 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); + } + + // all done + return; +} From b9433657793d3cdce0b1aa1df09dedf4d9baa9f3 Mon Sep 17 00:00:00 2001 From: Pawel024 Date: Tue, 20 Jan 2026 13:56:50 +0100 Subject: [PATCH 05/24] fem: add tests for gradient values and arc length to isoparametric_embedded_segment.cc --- .../fem/isoparametric_embedded_segment.cc | 76 ++++++++++++++----- 1 file changed, 58 insertions(+), 18 deletions(-) diff --git a/tests/mito.lib/fem/isoparametric_embedded_segment.cc b/tests/mito.lib/fem/isoparametric_embedded_segment.cc index 2d849726..d808e7c6 100644 --- a/tests/mito.lib/fem/isoparametric_embedded_segment.cc +++ b/tests/mito.lib/fem/isoparametric_embedded_segment.cc @@ -7,13 +7,13 @@ #include -// the type of coordinates +// 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; // the type of discretization node using discretization_node_t = mito::discrete::discretization_node_t; -// the type of cell +// 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; @@ -87,37 +87,77 @@ test_gradient_consistency(const auto & element) return; } +// test that the arc length computed via the Jacobian measure matches the expected value +auto +test_arc_length(const auto & element, double expected_length) +{ + // the element type + using element_t = mito::utilities::base_type; + + // 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 Jacobian measure + double arc_length = 0.0; + mito::tensor::constexpr_for_1([&]() { + constexpr auto xi = quadrature_rule.point(q); + constexpr auto w = element_t::canonical_element_type::area * quadrature_rule.weight(q); + arc_length += w * mito::fem::blocks::jacobian_measure(element.jacobian()(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 + // 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, { 1.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 }); - { - // first order isoparametric segment - using element_p1_t = mito::fem::isoparametric_simplex_t<1, cell_t>; + // first order isoparametric embedded segment + using element_p1_t = mito::fem::isoparametric_simplex_t<1, cell_t>; + + // build the discretization nodes + auto discretization_node_0 = discretization_node_t(); + auto discretization_node_1 = discretization_node_t(); - // 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 }); - // a finite element - auto element_p1 = element_p1_t( - geometric_simplex, coord_system, - { 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 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 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 From d408726a4dd802ee6e39d4368703384b912618a8 Mon Sep 17 00:00:00 2001 From: Pawel024 Date: Tue, 20 Jan 2026 14:06:57 +0100 Subject: [PATCH 06/24] fem: update block tests for embedded segments to use a diagonal segment --- tests/mito.lib/fem/block_grad_grad_embedded_segment.cc | 5 +++-- tests/mito.lib/fem/block_mass_embedded_segment.cc | 5 +++-- 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/tests/mito.lib/fem/block_grad_grad_embedded_segment.cc b/tests/mito.lib/fem/block_grad_grad_embedded_segment.cc index e811827e..9f8aa387 100644 --- a/tests/mito.lib/fem/block_grad_grad_embedded_segment.cc +++ b/tests/mito.lib/fem/block_grad_grad_embedded_segment.cc @@ -30,9 +30,10 @@ TEST(Fem, BlockGradGradEmbeddedSegment) // the coordinate system auto coord_system = coord_system_t(); - // build nodes (unit-length horizontal segment embedded in 2D) + // 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, { 1.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 }); diff --git a/tests/mito.lib/fem/block_mass_embedded_segment.cc b/tests/mito.lib/fem/block_mass_embedded_segment.cc index 0b7da468..094d3c1c 100644 --- a/tests/mito.lib/fem/block_mass_embedded_segment.cc +++ b/tests/mito.lib/fem/block_mass_embedded_segment.cc @@ -30,9 +30,10 @@ TEST(Fem, BlockMassEmbeddedSegment) // the coordinate system auto coord_system = coord_system_t(); - // build nodes (unit-length horizontal segment embedded in 2D) + // 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, { 1.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 }); From 01741b6d67ce7e273e7285d482dd2b9d973dbe44 Mon Sep 17 00:00:00 2001 From: Pawel024 Date: Fri, 23 Jan 2026 12:07:35 +0100 Subject: [PATCH 07/24] fem: update paths to match fields redesign --- .../elements/seg1_embedded/IsoparametricEmbeddedSegmentP1.h | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/lib/mito/fem/elements/seg1_embedded/IsoparametricEmbeddedSegmentP1.h b/lib/mito/fem/elements/seg1_embedded/IsoparametricEmbeddedSegmentP1.h index 71529598..397563a2 100644 --- a/lib/mito/fem/elements/seg1_embedded/IsoparametricEmbeddedSegmentP1.h +++ b/lib/mito/fem/elements/seg1_embedded/IsoparametricEmbeddedSegmentP1.h @@ -74,8 +74,7 @@ namespace mito::fem { constexpr auto phi_1 = shape_functions.shape<1>(); // return the isoparametric mapping from parametric to physical coordinates - return mito::fields::linear_combination( - std::array{ _x0, _x1 }, phi_0, phi_1); + return mito::functions::linear_combination(std::array{ _x0, _x1 }, phi_0, phi_1); } // get the shape function associated with local node {a} @@ -114,7 +113,7 @@ namespace mito::fem { { // assemble the gradient as a function of parametric coordinates auto gradient_function = - fields::field([&](const parametric_coordinates_type & xi) -> tensor::vector_t<2> { + 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); From 4fe483216b622356bbb45e9946bf4edcb4250d72 Mon Sep 17 00:00:00 2001 From: Pawel024 Date: Fri, 23 Jan 2026 12:10:19 +0100 Subject: [PATCH 08/24] fem: formatting edits --- .../elements/IsoparametricEmbeddedSegment.h | 6 ++++-- .../fem/elements/seg1_embedded/DiscretizerCG.h | 3 ++- .../IsoparametricEmbeddedSegmentP1.h | 18 +++++++++++------- 3 files changed, 17 insertions(+), 10 deletions(-) diff --git a/lib/mito/fem/elements/IsoparametricEmbeddedSegment.h b/lib/mito/fem/elements/IsoparametricEmbeddedSegment.h index 0c7bbce7..a9a67b85 100644 --- a/lib/mito/fem/elements/IsoparametricEmbeddedSegment.h +++ b/lib/mito/fem/elements/IsoparametricEmbeddedSegment.h @@ -50,10 +50,12 @@ namespace mito::fem { constexpr IsoparametricEmbeddedSegment(const IsoparametricEmbeddedSegment &) = delete; // delete assignment operator - constexpr IsoparametricEmbeddedSegment & operator=(const IsoparametricEmbeddedSegment &) = delete; + constexpr IsoparametricEmbeddedSegment & operator=(const IsoparametricEmbeddedSegment &) = + delete; // delete move assignment operator - constexpr IsoparametricEmbeddedSegment & operator=(IsoparametricEmbeddedSegment &&) noexcept = delete; + constexpr IsoparametricEmbeddedSegment & operator=( + IsoparametricEmbeddedSegment &&) noexcept = delete; public: // get the geometric simplex diff --git a/lib/mito/fem/elements/seg1_embedded/DiscretizerCG.h b/lib/mito/fem/elements/seg1_embedded/DiscretizerCG.h index fdceca3c..9292b69b 100644 --- a/lib/mito/fem/elements/seg1_embedded/DiscretizerCG.h +++ b/lib/mito/fem/elements/seg1_embedded/DiscretizerCG.h @@ -48,7 +48,8 @@ namespace mito::fem { } // populate the constrained nodes - // In 1D, constraints.domain() is a set of nodes, not a mesh with cells, so we loop on the nodes directly + // In 1D, constraints.domain() is a set of nodes, not a mesh with cells, so we loop on + // the nodes directly for (const auto & node : constraints.domain()) { // get the discretization node associated with the mesh node from the map auto it = node_map.find(node); diff --git a/lib/mito/fem/elements/seg1_embedded/IsoparametricEmbeddedSegmentP1.h b/lib/mito/fem/elements/seg1_embedded/IsoparametricEmbeddedSegmentP1.h index 397563a2..d147f23c 100644 --- a/lib/mito/fem/elements/seg1_embedded/IsoparametricEmbeddedSegmentP1.h +++ b/lib/mito/fem/elements/seg1_embedded/IsoparametricEmbeddedSegmentP1.h @@ -48,16 +48,19 @@ namespace mito::fem { inline ~IsoparametricEmbeddedSegmentP1() = default; // delete move constructor - constexpr IsoparametricEmbeddedSegmentP1(IsoparametricEmbeddedSegmentP1 &&) noexcept = delete; + constexpr IsoparametricEmbeddedSegmentP1(IsoparametricEmbeddedSegmentP1 &&) noexcept = + delete; // delete copy constructor constexpr IsoparametricEmbeddedSegmentP1(const IsoparametricEmbeddedSegmentP1 &) = delete; // delete assignment operator - constexpr IsoparametricEmbeddedSegmentP1 & operator=(const IsoparametricEmbeddedSegmentP1 &) = delete; + constexpr IsoparametricEmbeddedSegmentP1 & operator=( + const IsoparametricEmbeddedSegmentP1 &) = delete; // delete move assignment operator - constexpr IsoparametricEmbeddedSegmentP1 & operator=(IsoparametricEmbeddedSegmentP1 &&) noexcept = delete; + constexpr IsoparametricEmbeddedSegmentP1 & operator=( + IsoparametricEmbeddedSegmentP1 &&) noexcept = delete; public: // get the discretization nodes @@ -112,15 +115,16 @@ namespace mito::fem { 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> { + 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); // for the embedded case, the gradient is: grad(phi) = dphi/dxi * J / |J|^2 // this follows from: dphi/dx = grad(phi)^T * J, with grad(phi) parallel to J - auto J_norm_squared = J[{0, 0}] * J[{0, 0}] + J[{1, 0}] * J[{1, 0}]; - auto J_scaled = tensor::vector_t<2>{ J[{0, 0}] / J_norm_squared, J[{1, 0}] / J_norm_squared }; + auto J_norm_squared = J[{ 0, 0 }] * J[{ 0, 0 }] + J[{ 1, 0 }] * J[{ 1, 0 }]; + auto J_scaled = tensor::vector_t<2>{ J[{ 0, 0 }] / J_norm_squared, + J[{ 1, 0 }] / J_norm_squared }; // return the spatial gradients of the shape functions evaluated at {xi} return shape_functions.dshape()(xi) * J_scaled; }); From c1e17e275692a2e65cb2975c399b91fde327858e Mon Sep 17 00:00:00 2001 From: Pawel024 Date: Fri, 23 Jan 2026 12:12:59 +0100 Subject: [PATCH 09/24] tests/fem: formatting edits --- tests/mito.lib/fem/block_grad_grad_embedded_segment.cc | 3 +-- tests/mito.lib/fem/block_mass_embedded_segment.cc | 6 ++---- 2 files changed, 3 insertions(+), 6 deletions(-) diff --git a/tests/mito.lib/fem/block_grad_grad_embedded_segment.cc b/tests/mito.lib/fem/block_grad_grad_embedded_segment.cc index 9f8aa387..e7f9ec99 100644 --- a/tests/mito.lib/fem/block_grad_grad_embedded_segment.cc +++ b/tests/mito.lib/fem/block_grad_grad_embedded_segment.cc @@ -48,8 +48,7 @@ TEST(Fem, BlockGradGradEmbeddedSegment) // 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 }); // a grad-grad matrix block auto grad_grad_block = diff --git a/tests/mito.lib/fem/block_mass_embedded_segment.cc b/tests/mito.lib/fem/block_mass_embedded_segment.cc index 094d3c1c..03320623 100644 --- a/tests/mito.lib/fem/block_mass_embedded_segment.cc +++ b/tests/mito.lib/fem/block_mass_embedded_segment.cc @@ -48,15 +48,13 @@ TEST(Fem, BlockMassEmbeddedSegment) // 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 }); // a mass matrix block auto mass_block = mito::fem::blocks::mass_block(); // 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 }; + 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); From b313d85b2b56d337c8361ca1daf726a9fa87614b Mon Sep 17 00:00:00 2001 From: Pawel024 Date: Fri, 23 Jan 2026 12:16:10 +0100 Subject: [PATCH 10/24] fem: replace hardcoded dimension values with the variable dim --- lib/mito/fem/elements/IsoparametricEmbeddedSegment.h | 4 ++-- lib/mito/fem/elements/IsoparametricSegment.h | 4 ++-- lib/mito/fem/elements/IsoparametricTriangle.h | 4 ++-- 3 files changed, 6 insertions(+), 6 deletions(-) diff --git a/lib/mito/fem/elements/IsoparametricEmbeddedSegment.h b/lib/mito/fem/elements/IsoparametricEmbeddedSegment.h index a9a67b85..5d9fa921 100644 --- a/lib/mito/fem/elements/IsoparametricEmbeddedSegment.h +++ b/lib/mito/fem/elements/IsoparametricEmbeddedSegment.h @@ -25,11 +25,11 @@ namespace mito::fem { protected: // cartesian coordinates in 2D - using coordinates_type = geometry::coordinates_t<2, geometry::CARTESIAN>; + using coordinates_type = geometry::coordinates_t; // the coordinate system type using coordinate_system_type = geometry::coordinate_system_t; // the vector type - using vector_type = tensor::vector_t<2>; + using vector_type = tensor::vector_t; public: // the default constructor diff --git a/lib/mito/fem/elements/IsoparametricSegment.h b/lib/mito/fem/elements/IsoparametricSegment.h index cc3b384b..c7fb5dcf 100644 --- a/lib/mito/fem/elements/IsoparametricSegment.h +++ b/lib/mito/fem/elements/IsoparametricSegment.h @@ -25,11 +25,11 @@ namespace mito::fem { protected: // cartesian coordinates in 1D - using coordinates_type = geometry::coordinates_t<1, geometry::CARTESIAN>; + using coordinates_type = geometry::coordinates_t; // the coordinate system type using coordinate_system_type = geometry::coordinate_system_t; // the vector type - using vector_type = tensor::vector_t<1>; + using vector_type = tensor::vector_t; public: // the default constructor diff --git a/lib/mito/fem/elements/IsoparametricTriangle.h b/lib/mito/fem/elements/IsoparametricTriangle.h index 5426f2ad..a216e57b 100644 --- a/lib/mito/fem/elements/IsoparametricTriangle.h +++ b/lib/mito/fem/elements/IsoparametricTriangle.h @@ -24,11 +24,11 @@ namespace mito::fem { protected: // cartesian coordinates in 2D - using coordinates_type = geometry::coordinates_t<2, geometry::CARTESIAN>; + using coordinates_type = geometry::coordinates_t; // the coordinate system type using coordinate_system_type = geometry::coordinate_system_t; // the vector type - using vector_type = tensor::vector_t<2>; + using vector_type = tensor::vector_t; public: // the default constructor From 2b9e098023bbce3a7c8a85835ad4bc108c29198b Mon Sep 17 00:00:00 2001 From: Pawel024 Date: Fri, 23 Jan 2026 12:24:55 +0100 Subject: [PATCH 11/24] fem: make base isoparametric segment class a template of the dimension --- .../elements/IsoparametricEmbeddedSegment.h | 77 ------------------- lib/mito/fem/elements/IsoparametricSegment.h | 8 +- lib/mito/fem/elements/public.h | 1 - .../elements/seg1/IsoparametricSegmentP1.h | 4 +- .../IsoparametricEmbeddedSegmentP1.h | 6 +- 5 files changed, 9 insertions(+), 87 deletions(-) delete mode 100644 lib/mito/fem/elements/IsoparametricEmbeddedSegment.h diff --git a/lib/mito/fem/elements/IsoparametricEmbeddedSegment.h b/lib/mito/fem/elements/IsoparametricEmbeddedSegment.h deleted file mode 100644 index 5d9fa921..00000000 --- a/lib/mito/fem/elements/IsoparametricEmbeddedSegment.h +++ /dev/null @@ -1,77 +0,0 @@ -// -*- c++ -*- -// -// Copyright (c) 2020-2024, the MiTo Authors, all rights reserved -// - -// code guard -#pragma once - - -// DESIGN NOTES -// Class {IsoparametricEmbeddedSegment} represents a first order simplex (segment) embedded in 2D -// space, equipped with parametric coordinates. - - -namespace mito::fem { - - class IsoparametricEmbeddedSegment : public utilities::Invalidatable { - public: - // the dimension of the physical space - static constexpr int dim = 2; - // the discretization node type - using discretization_node_type = discrete::discretization_node_t; - // the underlying cell type - using cell_type = geometry::segment_t; - - protected: - // cartesian coordinates in 2D - using coordinates_type = geometry::coordinates_t; - // the coordinate system type - using coordinate_system_type = geometry::coordinate_system_t; - // the vector type - using vector_type = tensor::vector_t; - - public: - // the default constructor - constexpr IsoparametricEmbeddedSegment( - const cell_type & cell, const coordinate_system_type & coord_system) : - _cell(cell), - _x0{ coord_system.coordinates(cell.nodes()[0]->point()) - coordinates_type{} }, - _x1{ coord_system.coordinates(cell.nodes()[1]->point()) - coordinates_type{} } - {} - - // destructor - constexpr ~IsoparametricEmbeddedSegment() = default; - - // delete move constructor - constexpr IsoparametricEmbeddedSegment(IsoparametricEmbeddedSegment &&) noexcept = delete; - - // delete copy constructor - constexpr IsoparametricEmbeddedSegment(const IsoparametricEmbeddedSegment &) = delete; - - // delete assignment operator - constexpr IsoparametricEmbeddedSegment & operator=(const IsoparametricEmbeddedSegment &) = - delete; - - // delete move assignment operator - constexpr IsoparametricEmbeddedSegment & operator=( - IsoparametricEmbeddedSegment &&) noexcept = delete; - - public: - // get the geometric simplex - constexpr auto cell() const noexcept -> const cell_type & { return _cell; } - - protected: - // QUESTION: do we need to maintain a reference to the geometric simplex? - // a const reference to the geometric simplex - const cell_type & _cell; - - // the coordinates of the discretization nodes of the segment - const vector_type _x0; - const vector_type _x1; - }; - -} // namespace mito - - -// end of file diff --git a/lib/mito/fem/elements/IsoparametricSegment.h b/lib/mito/fem/elements/IsoparametricSegment.h index c7fb5dcf..d0f8be28 100644 --- a/lib/mito/fem/elements/IsoparametricSegment.h +++ b/lib/mito/fem/elements/IsoparametricSegment.h @@ -9,22 +9,22 @@ // DESIGN NOTES // Class {IsoparametricSegment} represents a first order simplex (segment) equipped with parametric -// coordinates. +// coordinates. The template parameter D specifies the dimension of the physical space in which +// the segment lives (D=1 for a segment in 1D, D=2 for a segment embedded in 2D, etc.). namespace mito::fem { + template class IsoparametricSegment : public utilities::Invalidatable { public: - // the dimension of the physical space - static constexpr int dim = 1; // the discretization node type using discretization_node_type = discrete::discretization_node_t; // the underlying cell type using cell_type = geometry::segment_t; protected: - // cartesian coordinates in 1D + // cartesian coordinates in dim dimensions using coordinates_type = geometry::coordinates_t; // the coordinate system type using coordinate_system_type = geometry::coordinate_system_t; diff --git a/lib/mito/fem/elements/public.h b/lib/mito/fem/elements/public.h index 633f46f0..7ed2002a 100644 --- a/lib/mito/fem/elements/public.h +++ b/lib/mito/fem/elements/public.h @@ -19,7 +19,6 @@ // classes implementation #include "IsoparametricSegment.h" #include "IsoparametricTriangle.h" -#include "IsoparametricEmbeddedSegment.h" #include "Discretizer.h" // library of finite elements diff --git a/lib/mito/fem/elements/seg1/IsoparametricSegmentP1.h b/lib/mito/fem/elements/seg1/IsoparametricSegmentP1.h index eee35058..876d8cec 100644 --- a/lib/mito/fem/elements/seg1/IsoparametricSegmentP1.h +++ b/lib/mito/fem/elements/seg1/IsoparametricSegmentP1.h @@ -14,7 +14,7 @@ namespace mito::fem { - class IsoparametricSegmentP1 : public IsoparametricSegment { + class IsoparametricSegmentP1 : public IsoparametricSegment<1> { public: // the degree of the finite element @@ -38,7 +38,7 @@ namespace mito::fem { inline IsoparametricSegmentP1( const cell_type & geometric_simplex, const coordinate_system_type & coord_system, const connectivity_type & connectivity) : - IsoparametricSegment(geometric_simplex, coord_system), + IsoparametricSegment<1>(geometric_simplex, coord_system), _connectivity(connectivity) {} diff --git a/lib/mito/fem/elements/seg1_embedded/IsoparametricEmbeddedSegmentP1.h b/lib/mito/fem/elements/seg1_embedded/IsoparametricEmbeddedSegmentP1.h index d147f23c..c9a6e0cb 100644 --- a/lib/mito/fem/elements/seg1_embedded/IsoparametricEmbeddedSegmentP1.h +++ b/lib/mito/fem/elements/seg1_embedded/IsoparametricEmbeddedSegmentP1.h @@ -8,7 +8,7 @@ // DESIGN NOTES -// Class {IsoparametricEmbeddedSegmentP1} represents a first order simplex (segment) living in 2D +// Class {IsoparametricEmbeddedSegmentP1} represents a first order simplex (segment) embedded in 2D // cartesian space, equipped with linear shape functions defined in the 1D parametric space. // The key difference from IsoparametricSegmentP1 is the non-square Jacobian (2x1 matrix) and // the gradient computation that accounts for the embedded geometry. @@ -16,7 +16,7 @@ namespace mito::fem { - class IsoparametricEmbeddedSegmentP1 : public IsoparametricEmbeddedSegment { + class IsoparametricEmbeddedSegmentP1 : public IsoparametricSegment<2> { public: // the degree of the finite element @@ -40,7 +40,7 @@ namespace mito::fem { inline IsoparametricEmbeddedSegmentP1( const cell_type & geometric_simplex, const coordinate_system_type & coord_system, const connectivity_type & connectivity) : - IsoparametricEmbeddedSegment(geometric_simplex, coord_system), + IsoparametricSegment<2>(geometric_simplex, coord_system), _connectivity(connectivity) {} From f3f345125c9fdc8775f268e7212e5e133a07f7e1 Mon Sep 17 00:00:00 2001 From: Pawel024 Date: Tue, 27 Jan 2026 13:50:25 +0100 Subject: [PATCH 12/24] fem: update elements api to select the right simplex based on degree and manifold type. Delete the seg1_embedded folder, we'll unify it with seg1 --- lib/mito/fem/elements/api.h | 6 +- lib/mito/fem/elements/elements_library.h | 56 ++++++- lib/mito/fem/elements/forward.h | 7 +- lib/mito/fem/elements/seg1/api.h | 8 +- .../elements/seg1_embedded/DiscretizerCG.h | 66 -------- .../IsoparametricEmbeddedSegmentP1.h | 143 ------------------ lib/mito/fem/elements/seg1_embedded/api.h | 21 --- lib/mito/fem/elements/seg1_embedded/public.h | 18 --- lib/mito/fem/elements/tri1/api.h | 8 +- lib/mito/fem/elements/tri1/public.h | 3 - lib/mito/fem/elements/tri2/api.h | 8 +- 11 files changed, 63 insertions(+), 281 deletions(-) delete mode 100644 lib/mito/fem/elements/seg1_embedded/DiscretizerCG.h delete mode 100644 lib/mito/fem/elements/seg1_embedded/IsoparametricEmbeddedSegmentP1.h delete mode 100644 lib/mito/fem/elements/seg1_embedded/api.h delete mode 100644 lib/mito/fem/elements/seg1_embedded/public.h 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 - using isoparametric_simplex_t = typename isoparametric_simplex::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 077f1a98..3355bc54 100644 --- a/lib/mito/fem/elements/elements_library.h +++ b/lib/mito/fem/elements/elements_library.h @@ -7,10 +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; +// +// 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 "seg1_embedded/public.h" #include "tri1/public.h" #include "tri2/public.h" +namespace mito::fem { + + template + struct isoparametric_simplex; + + // specialization for degree 1, 1D manifolds (segments) + template + 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; + }; + + // specialization for degree 1, 2D manifolds (triangles) + template + 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; + }; + + // specialization for degree 2, 2D manifolds (triangles) + template + 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; + }; + + // convenience alias - primary way users specify element types + // usage: isoparametric_simplex_t + template + using isoparametric_simplex_t = typename isoparametric_simplex::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 struct Discretizer; - // struct storing the type of an isoparametric simplex of polynomial degree {degree} on a - // geometric simplex of type {geometricSimplexT} - template + // 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 struct isoparametric_simplex; } 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/seg1_embedded/DiscretizerCG.h b/lib/mito/fem/elements/seg1_embedded/DiscretizerCG.h deleted file mode 100644 index 9292b69b..00000000 --- a/lib/mito/fem/elements/seg1_embedded/DiscretizerCG.h +++ /dev/null @@ -1,66 +0,0 @@ -// -*- c++ -*- -// -// Copyright (c) 2020-2024, the MiTo Authors, all rights reserved -// - -// code guard -#pragma once - - -namespace mito::fem { - - // discretizer specialization for {IsoparametricEmbeddedSegmentP1} with continuous Galerkin - template <> - struct Discretizer { - template < - typename manifoldT, typename constraintsT, typename elements_type, typename map_type, - typename constrained_nodes_type> - static void apply( - const manifoldT & manifold, const constraintsT & constraints, elements_type & elements, - map_type & node_map, constrained_nodes_type & constrained_nodes) - { - - // the discretization node type - using discretization_node_type = - typename IsoparametricEmbeddedSegmentP1::discretization_node_type; - - // the connectivity type - using connectivity_type = typename IsoparametricEmbeddedSegmentP1::connectivity_type; - - // get the coordinate system of the manifold - const auto & coord_system = manifold.coordinate_system(); - - // loop on the cells of the mesh - for (const auto & cell : manifold.elements()) { - - // get the nodes of the cell - const auto & nodes = cell.nodes(); - - // add the nodes to the map (if the mesh node is already present in the map, - // then the present discretization node is used) - auto node_0 = - node_map.insert({ nodes[0], discretization_node_type() }).first->second; - 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 }); - } - - // populate the constrained nodes - // In 1D, constraints.domain() is a set of nodes, not a mesh with cells, so we loop on - // the nodes directly - for (const auto & node : constraints.domain()) { - // get the discretization node associated with the mesh node from the map - auto it = node_map.find(node); - // add the node to the constrained nodes - constrained_nodes.insert(it->second); - } - - // all done - return; - } - }; -} - -// end of file diff --git a/lib/mito/fem/elements/seg1_embedded/IsoparametricEmbeddedSegmentP1.h b/lib/mito/fem/elements/seg1_embedded/IsoparametricEmbeddedSegmentP1.h deleted file mode 100644 index c9a6e0cb..00000000 --- a/lib/mito/fem/elements/seg1_embedded/IsoparametricEmbeddedSegmentP1.h +++ /dev/null @@ -1,143 +0,0 @@ -// -*- c++ -*- -// -// Copyright (c) 2020-2024, the MiTo Authors, all rights reserved -// - -// code guard -#pragma once - - -// DESIGN NOTES -// Class {IsoparametricEmbeddedSegmentP1} represents a first order simplex (segment) embedded in 2D -// cartesian space, equipped with linear shape functions defined in the 1D parametric space. -// The key difference from IsoparametricSegmentP1 is the non-square Jacobian (2x1 matrix) and -// the gradient computation that accounts for the embedded geometry. - - -namespace mito::fem { - - class IsoparametricEmbeddedSegmentP1 : public IsoparametricSegment<2> { - - public: - // the degree of the finite element - static constexpr int degree = 1; - // the type of shape functions (reuse from seg1, same 1D parametric shape functions) - using shape_functions_type = ShapeSegmentP1; - // the canonical element type - using canonical_element_type = typename shape_functions_type::reference_element_type; - // the parametric coordinates type - using parametric_coordinates_type = - typename canonical_element_type::parametric_coordinates_type; - // the linear 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; - - public: - // the default constructor - inline IsoparametricEmbeddedSegmentP1( - const cell_type & geometric_simplex, const coordinate_system_type & coord_system, - const connectivity_type & connectivity) : - IsoparametricSegment<2>(geometric_simplex, coord_system), - _connectivity(connectivity) - {} - - // destructor - inline ~IsoparametricEmbeddedSegmentP1() = default; - - // delete move constructor - constexpr IsoparametricEmbeddedSegmentP1(IsoparametricEmbeddedSegmentP1 &&) noexcept = - delete; - - // delete copy constructor - constexpr IsoparametricEmbeddedSegmentP1(const IsoparametricEmbeddedSegmentP1 &) = delete; - - // delete assignment operator - constexpr IsoparametricEmbeddedSegmentP1 & operator=( - const IsoparametricEmbeddedSegmentP1 &) = delete; - - // delete move assignment operator - constexpr IsoparametricEmbeddedSegmentP1 & operator=( - IsoparametricEmbeddedSegmentP1 &&) noexcept = delete; - - public: - // get the discretization nodes - constexpr auto connectivity() const noexcept -> const connectivity_type & - { - return _connectivity; - } - - // get the isoparametric mapping from parametric coordinates to physical coordinates - constexpr auto parametrization() const - { - // get the shape functions - constexpr auto phi_0 = shape_functions.shape<0>(); - 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); - } - - // get the shape function associated with local node {a} - template - requires(a >= 0 && a < n_nodes) - constexpr auto shape() const - { - // return the shape functions - return shape_functions.shape(); - } - - // get the jacobian of the isoparametric mapping from parametric to actual coordinates - 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<2, 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 2x1 matrix - return tensor::matrix_t<2, 1>{ dx_dxi[0], dx_dxi[1] }; - }); - - // and return it - return jacobian_function; - } - - // get the gradient of the a-th shape function as a function of parametric coordinates - template - 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); - // for the embedded case, the gradient is: grad(phi) = dphi/dxi * J / |J|^2 - // this follows from: dphi/dx = grad(phi)^T * J, with grad(phi) parallel to J - auto J_norm_squared = J[{ 0, 0 }] * J[{ 0, 0 }] + J[{ 1, 0 }] * J[{ 1, 0 }]; - auto J_scaled = tensor::vector_t<2>{ J[{ 0, 0 }] / J_norm_squared, - J[{ 1, 0 }] / J_norm_squared }; - // return the spatial gradients of the shape functions evaluated at {xi} - return shape_functions.dshape()(xi) * J_scaled; - }); - // and return it - return gradient_function; - } - - private: - // the discretization nodes of the simplex - const connectivity_type _connectivity; - }; - -} // namespace mito - - -// end of file diff --git a/lib/mito/fem/elements/seg1_embedded/api.h b/lib/mito/fem/elements/seg1_embedded/api.h deleted file mode 100644 index 30edb575..00000000 --- a/lib/mito/fem/elements/seg1_embedded/api.h +++ /dev/null @@ -1,21 +0,0 @@ -// -*- c++ -*- -// -// Copyright (c) 2020-2024, the MiTo Authors, all rights reserved -// - -// code guard -#pragma once - - -namespace mito::fem { - - // specialization for linear shape functions on segments embedded in 2D - template <> - struct isoparametric_simplex<1, geometry::segment_t<2>> { - using type = IsoparametricEmbeddedSegmentP1; - }; - -} - - -// end of file diff --git a/lib/mito/fem/elements/seg1_embedded/public.h b/lib/mito/fem/elements/seg1_embedded/public.h deleted file mode 100644 index 1a64d4cd..00000000 --- a/lib/mito/fem/elements/seg1_embedded/public.h +++ /dev/null @@ -1,18 +0,0 @@ -// -*- c++ -*- -// -// Copyright (c) 2020-2024, the MiTo Authors, all rights reserved -// - -// code guard -#pragma once - - -// classes implementation (reuses ShapeSegmentP1 from seg1) -#include "IsoparametricEmbeddedSegmentP1.h" -#include "DiscretizerCG.h" - -// published types and factories -#include "api.h" - - -// end of file 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/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 } From e81d32dda9ef6a202258a144af2a88b228ced47a Mon Sep 17 00:00:00 2001 From: Pawel024 Date: Tue, 27 Jan 2026 14:02:10 +0100 Subject: [PATCH 13/24] fem: expose volume form from manifold and pass it to the element classes --- lib/mito/fem/elements/IsoparametricSegment.h | 24 +++++++++++++---- lib/mito/fem/elements/IsoparametricTriangle.h | 19 ++++++++++--- lib/mito/fem/elements/seg1/DiscretizerCG.h | 21 ++++++++++----- .../elements/seg1/IsoparametricSegmentP1.h | 25 ++++++++++++----- lib/mito/fem/elements/tri1/DiscretizerCG.h | 21 ++++++++++----- .../elements/tri1/IsoparametricTriangleP1.h | 25 ++++++++++++----- lib/mito/fem/elements/tri2/DiscretizerCG.h | 26 +++++++++++------- .../elements/tri2/IsoparametricTriangleP2.h | 27 ++++++++++++++----- lib/mito/manifolds/Manifold.h | 10 +++++-- 9 files changed, 146 insertions(+), 52 deletions(-) diff --git a/lib/mito/fem/elements/IsoparametricSegment.h b/lib/mito/fem/elements/IsoparametricSegment.h index d0f8be28..f3d37b1b 100644 --- a/lib/mito/fem/elements/IsoparametricSegment.h +++ b/lib/mito/fem/elements/IsoparametricSegment.h @@ -9,13 +9,16 @@ // DESIGN NOTES // Class {IsoparametricSegment} represents a first order simplex (segment) equipped with parametric -// coordinates. The template parameter D specifies the dimension of the physical space in which -// the segment lives (D=1 for a segment in 1D, D=2 for a segment embedded in 2D, etc.). +// 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 + template class IsoparametricSegment : public utilities::Invalidatable { public: // the discretization node type @@ -29,13 +32,18 @@ namespace mito::fem { // the coordinate system type using coordinate_system_type = geometry::coordinate_system_t; // the vector type - using vector_type = tensor::vector_t; + using vector_type = tensor::vector_t; + // 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{} } {} @@ -64,6 +72,12 @@ namespace mito::fem { // 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 a216e57b..c540a1fe 100644 --- a/lib/mito/fem/elements/IsoparametricTriangle.h +++ b/lib/mito/fem/elements/IsoparametricTriangle.h @@ -8,11 +8,17 @@ // DESIGN NOTES -// Class {IsoparametricTriangle} represents a second order simplex equipped barycentric coordinates. +// Class {IsoparametricTriangle} represents a second order 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 class IsoparametricTriangle : public utilities::Invalidatable { public: // the dimension of the physical space @@ -28,14 +34,18 @@ namespace mito::fem { // the coordinate system type using coordinate_system_type = geometry::coordinate_system_t; // the vector type - using vector_type = tensor::vector_t; + using vector_type = tensor::vector_t; + // 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{} } @@ -70,6 +80,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/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 { + // agnostic of the coordinate type and volume form + template + struct Discretizer, 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; // 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 876d8cec..c444047a 100644 --- a/lib/mito/fem/elements/seg1/IsoparametricSegmentP1.h +++ b/lib/mito/fem/elements/seg1/IsoparametricSegmentP1.h @@ -8,13 +8,19 @@ // 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 first order simplex (segment) 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<1> { + template + class IsoparametricSegmentP1 : public IsoparametricSegment { + + private: + using base_type = IsoparametricSegment; public: // the degree of the finite element @@ -31,14 +37,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; + using connectivity_type = std::array; + + // 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<1>(geometric_simplex, coord_system), + const connectivity_type & connectivity, const volume_form_type & volume_form) : + base_type(geometric_simplex, coord_system, volume_form), _connectivity(connectivity) {} 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 { + // agnostic of the coordinate type and volume form + template + struct Discretizer, 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; // 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..56e7fc5e 100644 --- a/lib/mito/fem/elements/tri1/IsoparametricTriangleP1.h +++ b/lib/mito/fem/elements/tri1/IsoparametricTriangleP1.h @@ -8,13 +8,19 @@ // 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 second order 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 + class IsoparametricTriangleP1 : public IsoparametricTriangle { + + private: + using base_type = IsoparametricTriangle; public: // the degree of the finite element @@ -31,14 +37,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; + using connectivity_type = std::array; + + // 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) {} 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 { + // agnostic of the coordinate type and volume form + template + struct Discretizer, 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; + // 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; + using mesh_node_type = geometry::node_t; // 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; @@ -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..5f9b1894 100644 --- a/lib/mito/fem/elements/tri2/IsoparametricTriangleP2.h +++ b/lib/mito/fem/elements/tri2/IsoparametricTriangleP2.h @@ -8,13 +8,19 @@ // 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 second order 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 + class IsoparametricTriangleP2 : public IsoparametricTriangle { + + private: + using base_type = IsoparametricTriangle; public: // the degree of the finite element @@ -26,19 +32,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; + using connectivity_type = std::array; + + // 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) {} 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; + // 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()); } From 69a068ca791814e148639251b9e8918c0824b084 Mon Sep 17 00:00:00 2001 From: Pawel024 Date: Tue, 27 Jan 2026 14:08:14 +0100 Subject: [PATCH 14/24] fem: replace the jacobian_measure utility with element.volume_element --- lib/mito/fem/blocks/GradGradBlock.h | 4 +- lib/mito/fem/blocks/L2NormBlock.h | 4 +- lib/mito/fem/blocks/MassBlock.h | 4 +- lib/mito/fem/blocks/SourceTermBlock.h | 4 +- lib/mito/fem/blocks/public.h | 3 -- lib/mito/fem/blocks/utilities.h | 36 ----------------- .../elements/seg1/IsoparametricSegmentP1.h | 35 ++++++++++++++++- .../elements/tri1/IsoparametricTriangleP1.h | 39 +++++++++++++++++++ .../elements/tri2/IsoparametricTriangleP2.h | 38 ++++++++++++++++++ 9 files changed, 119 insertions(+), 48 deletions(-) delete mode 100644 lib/mito/fem/blocks/utilities.h diff --git a/lib/mito/fem/blocks/GradGradBlock.h b/lib/mito/fem/blocks/GradGradBlock.h index 277e2610..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 * jacobian_measure(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([&]() { diff --git a/lib/mito/fem/blocks/L2NormBlock.h b/lib/mito/fem/blocks/L2NormBlock.h index b5f86732..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 * jacobian_measure(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 8643efef..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 * jacobian_measure(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([&]() { diff --git a/lib/mito/fem/blocks/SourceTermBlock.h b/lib/mito/fem/blocks/SourceTermBlock.h index 13914a96..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 * jacobian_measure(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([&]() { diff --git a/lib/mito/fem/blocks/public.h b/lib/mito/fem/blocks/public.h index d202940c..20b5b6e9 100644 --- a/lib/mito/fem/blocks/public.h +++ b/lib/mito/fem/blocks/public.h @@ -16,9 +16,6 @@ // published types factories; this is the file you are looking for... #include "api.h" -// utilities (must be included before blocks that use jacobian_measure) -#include "utilities.h" - // classes implementation #include "AssemblyBlock.h" #include "GradGradBlock.h" diff --git a/lib/mito/fem/blocks/utilities.h b/lib/mito/fem/blocks/utilities.h deleted file mode 100644 index c46ede26..00000000 --- a/lib/mito/fem/blocks/utilities.h +++ /dev/null @@ -1,36 +0,0 @@ -// -*- c++ -*- -// -// Copyright (c) 2020-2024, the MiTo Authors, all rights reserved -// - -// code guard -#pragma once - - -namespace mito::fem::blocks { - - // Helper function to compute the Jacobian measure (integration scaling factor) - // for different Jacobian matrix shapes. - - // For square Jacobians (NxN): return |det(J)| - template - requires(tensor::square_matrix_c) - constexpr auto jacobian_measure(const matrixT & J) - { - return std::abs(tensor::determinant(J)); - } - - // For 2x1 Jacobians (embedded 1D in 2D): return sqrt(det(J^T * J)) - // This is the square root of the Gram determinant as used in the area formula - template - requires(tensor::matrix_c && matrixT::dims[0] == 2 && matrixT::dims[1] == 1) - constexpr auto jacobian_measure(const matrixT & J) - { - auto JtJ = tensor::transpose(J) * J; // 1x1 matrix - return std::sqrt(JtJ[{ 0, 0 }]); - } - -} // namespace mito - - -// end of file diff --git a/lib/mito/fem/elements/seg1/IsoparametricSegmentP1.h b/lib/mito/fem/elements/seg1/IsoparametricSegmentP1.h index c444047a..b6235803 100644 --- a/lib/mito/fem/elements/seg1/IsoparametricSegmentP1.h +++ b/lib/mito/fem/elements/seg1/IsoparametricSegmentP1.h @@ -23,6 +23,11 @@ namespace mito::fem { using base_type = IsoparametricSegment; 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 @@ -85,7 +90,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} @@ -117,6 +123,33 @@ namespace mito::fem { 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 = [&J]() { + if constexpr (D == 1) { + return tensor::vector_t<1>{ J[{ 0, 0 }] }; + } else if constexpr (D == 2) { + return tensor::vector_t<2>{ J[{ 0, 0 }], J[{ 1, 0 }] }; + } else if constexpr (D == 3) { + return tensor::vector_t<3>{ J[{ 0, 0 }], J[{ 1, 0 }], J[{ 2, 0 }] }; + } + }(); + + // 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 template requires(a >= 0 && a < n_nodes) diff --git a/lib/mito/fem/elements/tri1/IsoparametricTriangleP1.h b/lib/mito/fem/elements/tri1/IsoparametricTriangleP1.h index 56e7fc5e..82b6920c 100644 --- a/lib/mito/fem/elements/tri1/IsoparametricTriangleP1.h +++ b/lib/mito/fem/elements/tri1/IsoparametricTriangleP1.h @@ -23,6 +23,11 @@ namespace mito::fem { using base_type = IsoparametricTriangle; 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 @@ -107,6 +112,40 @@ namespace mito::fem { return jacobian_function; } + // volume element: contract the volume form with the two tangent vectors + // this follows the same pattern as Manifold::_volume + 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 = [&J]() { + if constexpr (D == 2) { + return tensor::vector_t<2>{ J[{ 0, 0 }], J[{ 1, 0 }] }; + } else if constexpr (D == 3) { + return tensor::vector_t<3>{ J[{ 0, 0 }], J[{ 1, 0 }], J[{ 2, 0 }] }; + } + }(); + + auto tangent_1 = [&J]() { + if constexpr (D == 2) { + return tensor::vector_t<2>{ J[{ 0, 1 }], J[{ 1, 1 }] }; + } else if constexpr (D == 3) { + return tensor::vector_t<3>{ J[{ 0, 1 }], J[{ 1, 1 }], J[{ 2, 1 }] }; + } + }(); + + // contract the volume form with both tangent vectors + // for N=2, include factorial (1/2!) + return (1.0 / 2.0) * this->_volume_form(x)(tangent_0, tangent_1); + }); + } + // get the gradient of the a-th shape function as a function of barycentric coordinates template requires(a >= 0 && a < n_nodes) diff --git a/lib/mito/fem/elements/tri2/IsoparametricTriangleP2.h b/lib/mito/fem/elements/tri2/IsoparametricTriangleP2.h index 5f9b1894..320f09a3 100644 --- a/lib/mito/fem/elements/tri2/IsoparametricTriangleP2.h +++ b/lib/mito/fem/elements/tri2/IsoparametricTriangleP2.h @@ -23,6 +23,11 @@ namespace mito::fem { using base_type = IsoparametricTriangle; 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 @@ -115,6 +120,39 @@ namespace mito::fem { return jacobian_function; } + // volume element: contract the volume form with the two tangent vectors + 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 = [&J]() { + if constexpr (D == 2) { + return tensor::vector_t<2>{ J[{ 0, 0 }], J[{ 1, 0 }] }; + } else if constexpr (D == 3) { + return tensor::vector_t<3>{ J[{ 0, 0 }], J[{ 1, 0 }], J[{ 2, 0 }] }; + } + }(); + + auto tangent_1 = [&J]() { + if constexpr (D == 2) { + return tensor::vector_t<2>{ J[{ 0, 1 }], J[{ 1, 1 }] }; + } else if constexpr (D == 3) { + return tensor::vector_t<3>{ J[{ 0, 1 }], J[{ 1, 1 }], J[{ 2, 1 }] }; + } + }(); + + // contract the volume form with both tangent vectors + // for N=2, include factorial (1/2!) + return (1.0 / 2.0) * this->_volume_form(x)(tangent_0, tangent_1); + }); + } + // get the gradient of the a-th shape function as a function of barycentric coordinates template requires(a >= 0 && a < n_nodes) From 9ce302ebf30a4d3d7b60cdd5cebfd7dae8e81ab6 Mon Sep 17 00:00:00 2001 From: Pawel024 Date: Tue, 27 Jan 2026 14:15:39 +0100 Subject: [PATCH 15/24] fem: update jacobian calculation in element classes to work for different physical dimensions --- lib/mito/fem/elements/IsoparametricSegment.h | 16 +++++++--- lib/mito/fem/elements/IsoparametricTriangle.h | 17 +++++++---- .../elements/seg1/IsoparametricSegmentP1.h | 18 +++++++---- .../elements/tri1/IsoparametricTriangleP1.h | 24 +++++++++++---- .../elements/tri2/IsoparametricTriangleP2.h | 30 ++++++++++++++----- 5 files changed, 77 insertions(+), 28 deletions(-) diff --git a/lib/mito/fem/elements/IsoparametricSegment.h b/lib/mito/fem/elements/IsoparametricSegment.h index f3d37b1b..9661f1b1 100644 --- a/lib/mito/fem/elements/IsoparametricSegment.h +++ b/lib/mito/fem/elements/IsoparametricSegment.h @@ -21,15 +21,23 @@ namespace mito::fem { template class IsoparametricSegment : public utilities::Invalidatable { public: + // 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; + // the discretization node type using discretization_node_type = discrete::discretization_node_t; // the underlying cell type - using cell_type = geometry::segment_t; + using cell_type = geometry::segment_t; protected: - // cartesian coordinates in dim dimensions - using coordinates_type = geometry::coordinates_t; - // the coordinate system type + // coordinate types + using coordinates_type = coordsT; using coordinate_system_type = geometry::coordinate_system_t; // the vector type using vector_type = tensor::vector_t; diff --git a/lib/mito/fem/elements/IsoparametricTriangle.h b/lib/mito/fem/elements/IsoparametricTriangle.h index c540a1fe..4c53ca4a 100644 --- a/lib/mito/fem/elements/IsoparametricTriangle.h +++ b/lib/mito/fem/elements/IsoparametricTriangle.h @@ -21,17 +21,22 @@ namespace mito::fem { template 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; + // the discretization node type using discretization_node_type = discrete::discretization_node_t; // the underlying cell type - using cell_type = geometry::triangle_t; + using cell_type = geometry::triangle_t; protected: - // cartesian coordinates in 2D - using coordinates_type = geometry::coordinates_t; - // the coordinate system type + // coordinate types + using coordinates_type = coordsT; using coordinate_system_type = geometry::coordinate_system_t; // the vector type using vector_type = tensor::vector_t; diff --git a/lib/mito/fem/elements/seg1/IsoparametricSegmentP1.h b/lib/mito/fem/elements/seg1/IsoparametricSegmentP1.h index b6235803..406e996d 100644 --- a/lib/mito/fem/elements/seg1/IsoparametricSegmentP1.h +++ b/lib/mito/fem/elements/seg1/IsoparametricSegmentP1.h @@ -104,19 +104,27 @@ 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 { // 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 + if constexpr (D == 1) { + return tensor::matrix_t<1, 1>{ dx_dxi[0] }; + } else if constexpr (D == 2) { + return tensor::matrix_t<2, 1>{ dx_dxi[0], dx_dxi[1] }; + } else if constexpr (D == 3) { + return tensor::matrix_t<3, 1>{ dx_dxi[0], dx_dxi[1], dx_dxi[2] }; + } }); // and return it diff --git a/lib/mito/fem/elements/tri1/IsoparametricTriangleP1.h b/lib/mito/fem/elements/tri1/IsoparametricTriangleP1.h index 82b6920c..4c4af44e 100644 --- a/lib/mito/fem/elements/tri1/IsoparametricTriangleP1.h +++ b/lib/mito/fem/elements/tri1/IsoparametricTriangleP1.h @@ -92,20 +92,34 @@ 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 { + // 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 - return ( - tensor::dyadic(_x0, dphi_0(xi)) + tensor::dyadic(_x1, dphi_1(xi)) - + tensor::dyadic(_x2, dphi_2(xi))); + // this is the outer product of position vectors with gradient vectors + if constexpr (D == 2) { + return ( + tensor::dyadic(this->_x0, dphi_0(xi)) + + tensor::dyadic(this->_x1, dphi_1(xi)) + + tensor::dyadic(this->_x2, dphi_2(xi))); + } else if constexpr (D == 3) { + // For embedded triangle in 3D, we need D×2 matrix + auto dxi0 = this->_x0 * dphi_0(xi)[0] + this->_x1 * dphi_1(xi)[0] + + this->_x2 * dphi_2(xi)[0]; + auto dxi1 = this->_x0 * dphi_0(xi)[1] + this->_x1 * dphi_1(xi)[1] + + this->_x2 * dphi_2(xi)[1]; + return tensor::matrix_t<3, 2>{ dxi0[0], dxi1[0], dxi0[1], + dxi1[1], dxi0[2], dxi1[2] }; + } }); // and return it diff --git a/lib/mito/fem/elements/tri2/IsoparametricTriangleP2.h b/lib/mito/fem/elements/tri2/IsoparametricTriangleP2.h index 320f09a3..46bb8685 100644 --- a/lib/mito/fem/elements/tri2/IsoparametricTriangleP2.h +++ b/lib/mito/fem/elements/tri2/IsoparametricTriangleP2.h @@ -92,14 +92,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 { + 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>(); @@ -110,10 +111,23 @@ namespace mito::fem { constexpr auto dphi_5 = shape_functions.dshape<5>(); // 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(x4, dphi_4(xi)) + tensor::dyadic(x5, dphi_5(xi))); + if constexpr (D == 2) { + return ( + 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))); + } else if constexpr (D == 3) { + // For embedded triangle in 3D + auto dxi0 = this->_x0 * dphi_0(xi)[0] + this->_x1 * dphi_1(xi)[0] + + this->_x2 * dphi_2(xi)[0] + x3 * dphi_3(xi)[0] + + x4 * dphi_4(xi)[0] + x5 * dphi_5(xi)[0]; + auto dxi1 = this->_x0 * dphi_0(xi)[1] + this->_x1 * dphi_1(xi)[1] + + this->_x2 * dphi_2(xi)[1] + x3 * dphi_3(xi)[1] + + x4 * dphi_4(xi)[1] + x5 * dphi_5(xi)[1]; + return tensor::matrix_t<3, 2>{ dxi0[0], dxi1[0], dxi0[1], + dxi1[1], dxi0[2], dxi1[2] }; + } }); // and return it From e3210e7fa70f854b5589ee9dd3575eeaa625a133 Mon Sep 17 00:00:00 2001 From: Pawel024 Date: Tue, 27 Jan 2026 15:02:07 +0100 Subject: [PATCH 16/24] tensor: add as_column_matrix factory, apply in IsoparametricSegmentP1 (mito/fem) --- lib/mito/fem/elements/seg1/IsoparametricSegmentP1.h | 8 +------- lib/mito/tensor/factories.h | 10 ++++++++++ 2 files changed, 11 insertions(+), 7 deletions(-) diff --git a/lib/mito/fem/elements/seg1/IsoparametricSegmentP1.h b/lib/mito/fem/elements/seg1/IsoparametricSegmentP1.h index 406e996d..127a6450 100644 --- a/lib/mito/fem/elements/seg1/IsoparametricSegmentP1.h +++ b/lib/mito/fem/elements/seg1/IsoparametricSegmentP1.h @@ -118,13 +118,7 @@ namespace mito::fem { auto dx_dxi = this->_x0 * dphi_0(xi) + this->_x1 * dphi_1(xi); // wrap the tangent vector as a D×1 matrix - if constexpr (D == 1) { - return tensor::matrix_t<1, 1>{ dx_dxi[0] }; - } else if constexpr (D == 2) { - return tensor::matrix_t<2, 1>{ dx_dxi[0], dx_dxi[1] }; - } else if constexpr (D == 3) { - return tensor::matrix_t<3, 1>{ dx_dxi[0], dx_dxi[1], dx_dxi[2] }; - } + return tensor::as_column_matrix(dx_dxi); }); // and return it 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 + constexpr auto as_column_matrix(const vector_t & vec) -> matrix_t + { + return [](const auto & v, std::index_sequence) { + return matrix_t{ v[Is]... }; + }(vec, std::make_index_sequence{}); + } } From 7f96cf6be45865a951727bc041c3138b14058f34 Mon Sep 17 00:00:00 2001 From: Pawel024 Date: Tue, 27 Jan 2026 15:23:07 +0100 Subject: [PATCH 17/24] fem: use pyre::tensor::col and tensor::dyadic properly to avoid repeting code --- .../elements/seg1/IsoparametricSegmentP1.h | 10 +---- .../elements/tri1/IsoparametricTriangleP1.h | 35 +++-------------- .../elements/tri2/IsoparametricTriangleP2.h | 39 ++++--------------- 3 files changed, 14 insertions(+), 70 deletions(-) diff --git a/lib/mito/fem/elements/seg1/IsoparametricSegmentP1.h b/lib/mito/fem/elements/seg1/IsoparametricSegmentP1.h index 127a6450..a791c0bf 100644 --- a/lib/mito/fem/elements/seg1/IsoparametricSegmentP1.h +++ b/lib/mito/fem/elements/seg1/IsoparametricSegmentP1.h @@ -136,15 +136,7 @@ namespace mito::fem { // get the tangent vector from the Jacobian auto J = jacobian()(xi); // extract the tangent vector from the D×1 matrix - auto tangent = [&J]() { - if constexpr (D == 1) { - return tensor::vector_t<1>{ J[{ 0, 0 }] }; - } else if constexpr (D == 2) { - return tensor::vector_t<2>{ J[{ 0, 0 }], J[{ 1, 0 }] }; - } else if constexpr (D == 3) { - return tensor::vector_t<3>{ J[{ 0, 0 }], J[{ 1, 0 }], J[{ 2, 0 }] }; - } - }(); + auto tangent = pyre::tensor::col<0>(J); // contract the volume form with the tangent vector // for N=1, no factorial needed (1/1! = 1) diff --git a/lib/mito/fem/elements/tri1/IsoparametricTriangleP1.h b/lib/mito/fem/elements/tri1/IsoparametricTriangleP1.h index 4c4af44e..e73d0554 100644 --- a/lib/mito/fem/elements/tri1/IsoparametricTriangleP1.h +++ b/lib/mito/fem/elements/tri1/IsoparametricTriangleP1.h @@ -106,20 +106,10 @@ namespace mito::fem { // compute the gradient of the isoparametric mapping // this is the outer product of position vectors with gradient vectors - if constexpr (D == 2) { - return ( - tensor::dyadic(this->_x0, dphi_0(xi)) - + tensor::dyadic(this->_x1, dphi_1(xi)) - + tensor::dyadic(this->_x2, dphi_2(xi))); - } else if constexpr (D == 3) { - // For embedded triangle in 3D, we need D×2 matrix - auto dxi0 = this->_x0 * dphi_0(xi)[0] + this->_x1 * dphi_1(xi)[0] - + this->_x2 * dphi_2(xi)[0]; - auto dxi1 = this->_x0 * dphi_0(xi)[1] + this->_x1 * dphi_1(xi)[1] - + this->_x2 * dphi_2(xi)[1]; - return tensor::matrix_t<3, 2>{ dxi0[0], dxi1[0], dxi0[1], - dxi1[1], dxi0[2], dxi1[2] }; - } + return ( + tensor::dyadic(this->_x0, dphi_0(xi)) + + tensor::dyadic(this->_x1, dphi_1(xi)) + + tensor::dyadic(this->_x2, dphi_2(xi))); }); // and return it @@ -138,21 +128,8 @@ namespace mito::fem { auto J = jacobian()(xi); // extract tangent vectors from the D×2 matrix - auto tangent_0 = [&J]() { - if constexpr (D == 2) { - return tensor::vector_t<2>{ J[{ 0, 0 }], J[{ 1, 0 }] }; - } else if constexpr (D == 3) { - return tensor::vector_t<3>{ J[{ 0, 0 }], J[{ 1, 0 }], J[{ 2, 0 }] }; - } - }(); - - auto tangent_1 = [&J]() { - if constexpr (D == 2) { - return tensor::vector_t<2>{ J[{ 0, 1 }], J[{ 1, 1 }] }; - } else if constexpr (D == 3) { - return tensor::vector_t<3>{ J[{ 0, 1 }], J[{ 1, 1 }], J[{ 2, 1 }] }; - } - }(); + auto tangent_0 = pyre::tensor::col<0>(J); + auto tangent_1 = pyre::tensor::col<1>(J); // contract the volume form with both tangent vectors // for N=2, include factorial (1/2!) diff --git a/lib/mito/fem/elements/tri2/IsoparametricTriangleP2.h b/lib/mito/fem/elements/tri2/IsoparametricTriangleP2.h index 46bb8685..a83e8018 100644 --- a/lib/mito/fem/elements/tri2/IsoparametricTriangleP2.h +++ b/lib/mito/fem/elements/tri2/IsoparametricTriangleP2.h @@ -111,23 +111,11 @@ namespace mito::fem { constexpr auto dphi_5 = shape_functions.dshape<5>(); // compute the gradient of the isoparametric mapping - if constexpr (D == 2) { - return ( - 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))); - } else if constexpr (D == 3) { - // For embedded triangle in 3D - auto dxi0 = this->_x0 * dphi_0(xi)[0] + this->_x1 * dphi_1(xi)[0] - + this->_x2 * dphi_2(xi)[0] + x3 * dphi_3(xi)[0] - + x4 * dphi_4(xi)[0] + x5 * dphi_5(xi)[0]; - auto dxi1 = this->_x0 * dphi_0(xi)[1] + this->_x1 * dphi_1(xi)[1] - + this->_x2 * dphi_2(xi)[1] + x3 * dphi_3(xi)[1] - + x4 * dphi_4(xi)[1] + x5 * dphi_5(xi)[1]; - return tensor::matrix_t<3, 2>{ dxi0[0], dxi1[0], dxi0[1], - dxi1[1], dxi0[2], dxi1[2] }; - } + return ( + 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))); }); // and return it @@ -145,21 +133,8 @@ namespace mito::fem { auto J = jacobian()(xi); // extract tangent vectors from the D×2 matrix - auto tangent_0 = [&J]() { - if constexpr (D == 2) { - return tensor::vector_t<2>{ J[{ 0, 0 }], J[{ 1, 0 }] }; - } else if constexpr (D == 3) { - return tensor::vector_t<3>{ J[{ 0, 0 }], J[{ 1, 0 }], J[{ 2, 0 }] }; - } - }(); - - auto tangent_1 = [&J]() { - if constexpr (D == 2) { - return tensor::vector_t<2>{ J[{ 0, 1 }], J[{ 1, 1 }] }; - } else if constexpr (D == 3) { - return tensor::vector_t<3>{ J[{ 0, 1 }], J[{ 1, 1 }], J[{ 2, 1 }] }; - } - }(); + auto tangent_0 = pyre::tensor::col<0>(J); + auto tangent_1 = pyre::tensor::col<1>(J); // contract the volume form with both tangent vectors // for N=2, include factorial (1/2!) From 31c7e73c68260eac9524a0ec3f8c427a9e2913f5 Mon Sep 17 00:00:00 2001 From: Pawel024 Date: Tue, 27 Jan 2026 23:54:11 +0100 Subject: [PATCH 18/24] fem: contravariant gradient computation with the riemannian formula - new elements api works --- lib/mito/fem/elements/IsoparametricSegment.h | 12 +++++ lib/mito/fem/elements/IsoparametricTriangle.h | 12 +++++ .../elements/seg1/IsoparametricSegmentP1.h | 27 +++++++---- .../elements/tri1/IsoparametricTriangleP1.h | 27 +++++++---- .../elements/tri2/IsoparametricTriangleP2.h | 27 +++++++---- lib/mito/fem/gradient.h | 48 +++++++++++++++++++ lib/mito/fem/public.h | 3 ++ 7 files changed, 132 insertions(+), 24 deletions(-) create mode 100644 lib/mito/fem/gradient.h diff --git a/lib/mito/fem/elements/IsoparametricSegment.h b/lib/mito/fem/elements/IsoparametricSegment.h index 9661f1b1..7b8a6a12 100644 --- a/lib/mito/fem/elements/IsoparametricSegment.h +++ b/lib/mito/fem/elements/IsoparametricSegment.h @@ -30,6 +30,9 @@ namespace mito::fem { // column vector) using jacobian_type = tensor::matrix_t; + // the metric space provides g, g_inv for this coordinate system + using metric_space_type = geometry::metric_space; + // the discretization node type using discretization_node_type = discrete::discretization_node_t; // the underlying cell type @@ -75,6 +78,15 @@ 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 diff --git a/lib/mito/fem/elements/IsoparametricTriangle.h b/lib/mito/fem/elements/IsoparametricTriangle.h index 4c53ca4a..94a7f085 100644 --- a/lib/mito/fem/elements/IsoparametricTriangle.h +++ b/lib/mito/fem/elements/IsoparametricTriangle.h @@ -29,6 +29,9 @@ namespace mito::fem { // the Jacobian is a D×N matrix (for a triangle: D×2) using jacobian_type = tensor::matrix_t; + // the metric space provides g, g_inv for this coordinate system + using metric_space_type = geometry::metric_space; + // the discretization node type using discretization_node_type = discrete::discretization_node_t; // the underlying cell type @@ -78,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; diff --git a/lib/mito/fem/elements/seg1/IsoparametricSegmentP1.h b/lib/mito/fem/elements/seg1/IsoparametricSegmentP1.h index a791c0bf..9291031a 100644 --- a/lib/mito/fem/elements/seg1/IsoparametricSegmentP1.h +++ b/lib/mito/fem/elements/seg1/IsoparametricSegmentP1.h @@ -21,6 +21,7 @@ namespace mito::fem { private: using base_type = IsoparametricSegment; + using metric_space_type = typename base_type::metric_space_type; public: static constexpr int D = base_type::D; // ambient dimension @@ -145,20 +146,30 @@ namespace mito::fem { } // get the gradient of the a-th shape function as a function of parametric coordinates + // returns a contravariant vector (raised index) template 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()(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(); + auto dphi_dxi_val = dphi_dxi(xi); + + // compute gradient (handles both square and embedded cases) + return fem::compute_gradient(J, g, g_inv, dphi_dxi_val); }); // and return it return gradient_function; diff --git a/lib/mito/fem/elements/tri1/IsoparametricTriangleP1.h b/lib/mito/fem/elements/tri1/IsoparametricTriangleP1.h index e73d0554..ece07580 100644 --- a/lib/mito/fem/elements/tri1/IsoparametricTriangleP1.h +++ b/lib/mito/fem/elements/tri1/IsoparametricTriangleP1.h @@ -21,6 +21,7 @@ namespace mito::fem { private: using base_type = IsoparametricTriangle; + using metric_space_type = typename base_type::metric_space_type; public: static constexpr int D = base_type::D; // ambient dimension @@ -138,20 +139,30 @@ namespace mito::fem { } // get the gradient of the a-th shape function as a function of barycentric coordinates + // returns a contravariant vector (raised index) template 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()(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(); + auto dphi_dxi_val = dphi_dxi(xi); + + // compute gradient (handles both square and embedded cases) + return fem::compute_gradient(J, g, g_inv, dphi_dxi_val); }); // and return it return gradient_function; diff --git a/lib/mito/fem/elements/tri2/IsoparametricTriangleP2.h b/lib/mito/fem/elements/tri2/IsoparametricTriangleP2.h index a83e8018..4c7a001e 100644 --- a/lib/mito/fem/elements/tri2/IsoparametricTriangleP2.h +++ b/lib/mito/fem/elements/tri2/IsoparametricTriangleP2.h @@ -21,6 +21,7 @@ namespace mito::fem { private: using base_type = IsoparametricTriangle; + using metric_space_type = typename base_type::metric_space_type; public: static constexpr int D = base_type::D; // ambient dimension @@ -143,20 +144,30 @@ namespace mito::fem { } // get the gradient of the a-th shape function as a function of barycentric coordinates + // returns a contravariant vector (raised index) template 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()(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(); + auto dphi_dxi_val = dphi_dxi(xi); + + // compute gradient (handles both square and embedded cases) + return fem::compute_gradient(J, g, g_inv, dphi_dxi_val); }); // and return it return gradient_function; 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 + 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 cc3f80db..3cccfcdd 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" From 4a08b110ba718959543c6d196d3e452cfb3bfbd6 Mon Sep 17 00:00:00 2001 From: Pawel024 Date: Wed, 28 Jan 2026 00:54:26 +0100 Subject: [PATCH 19/24] tests/fem: adapt tests to work with the new elements API --- tests/mito.lib/fem/block_grad_grad.cc | 17 ++++++++++---- .../fem/block_grad_grad_embedded_segment.cc | 16 +++++++++++-- tests/mito.lib/fem/block_grad_grad_segment.cc | 12 ++++++++-- tests/mito.lib/fem/block_mass.cc | 17 ++++++++++---- .../fem/block_mass_embedded_segment.cc | 13 +++++++++-- tests/mito.lib/fem/block_mass_segment.cc | 12 ++++++++-- tests/mito.lib/fem/domain_field.cc | 12 +++++++++- tests/mito.lib/fem/fem_field.cc | 9 ++++---- .../fem/isoparametric_embedded_segment.cc | 23 +++++++++++++++---- tests/mito.lib/fem/isoparametric_segment.cc | 12 ++++++++-- tests/mito.lib/fem/isoparametric_triangle.cc | 17 ++++++++++---- 11 files changed, 128 insertions(+), 32 deletions(-) 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(); + 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 index e7f9ec99..d523964c 100644 --- a/tests/mito.lib/fem/block_grad_grad_embedded_segment.cc +++ b/tests/mito.lib/fem/block_grad_grad_embedded_segment.cc @@ -38,9 +38,20 @@ TEST(Fem, BlockGradGradEmbeddedSegment) // 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(); + mesh.insert({ node_0, node_1 }); + + // create normal field for the submanifold (perpendicular to diagonal) + auto normal_field = mito::functions::constant( + 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, 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(); @@ -48,7 +59,8 @@ TEST(Fem, BlockGradGradEmbeddedSegment) // 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_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(); + 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(); + 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(); @@ -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(); diff --git a/tests/mito.lib/fem/block_mass_embedded_segment.cc b/tests/mito.lib/fem/block_mass_embedded_segment.cc index 03320623..e1327fc6 100644 --- a/tests/mito.lib/fem/block_mass_embedded_segment.cc +++ b/tests/mito.lib/fem/block_mass_embedded_segment.cc @@ -38,9 +38,17 @@ TEST(Fem, BlockMassEmbeddedSegment) // 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(); + mesh.insert({ node_0, node_1 }); + + // create a submanifold (1D embedded in 2D requires normal field) + auto normal_field = mito::functions::constant(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, 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(); @@ -48,7 +56,8 @@ TEST(Fem, BlockMassEmbeddedSegment) // 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(); 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(); + 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(); diff --git a/tests/mito.lib/fem/domain_field.cc b/tests/mito.lib/fem/domain_field.cc index 90dc6817..72c4f0d6 100644 --- a/tests/mito.lib/fem/domain_field.cc +++ b/tests/mito.lib/fem/domain_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; +// the type of cell +using cell_t = mito::geometry::triangle_t<2>; // the x scalar field in 2D constexpr auto x = mito::functions::component; // the y scalar field in 2D @@ -36,8 +38,16 @@ TEST(Fem, DomainField) // 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(); + 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()); // TOFIX: This syntax should also be allowed: field.localize(geometric_simplex). // However, in case, the coordinate system should be passed somehow to the function. 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; // the x scalar field in 2D constexpr auto x = mito::functions::component; // 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; + // 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 index d808e7c6..d67d224c 100644 --- a/tests/mito.lib/fem/isoparametric_embedded_segment.cc +++ b/tests/mito.lib/fem/isoparametric_embedded_segment.cc @@ -87,7 +87,7 @@ test_gradient_consistency(const auto & element) return; } -// test that the arc length computed via the Jacobian measure matches the expected value +// test that the arc length computed via the volume element matches the expected value auto test_arc_length(const auto & element, double expected_length) { @@ -97,12 +97,12 @@ test_arc_length(const auto & element, double expected_length) // 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 Jacobian measure + // compute the arc length by integrating 1 over the segment using the volume element double arc_length = 0.0; mito::tensor::constexpr_for_1([&]() { constexpr auto xi = quadrature_rule.point(q); constexpr auto w = element_t::canonical_element_type::area * quadrature_rule.weight(q); - arc_length += w * mito::fem::blocks::jacobian_measure(element.jacobian()(xi)); + arc_length += w * element.volume_element()(xi); }); // check the arc length @@ -124,8 +124,20 @@ TEST(Fem, IsoparametricEmbeddedSegment) // 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(); + 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) + auto normal_field = mito::functions::constant( + 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, 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(); @@ -133,7 +145,8 @@ TEST(Fem, IsoparametricEmbeddedSegment) // 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_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(); + 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(); + 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); From e6d974a0cfa8d6cd6866c48d8e8e3b3ac1f929c3 Mon Sep 17 00:00:00 2001 From: Pawel024 Date: Wed, 28 Jan 2026 13:54:18 +0100 Subject: [PATCH 20/24] fem: remove division by factorial in volume element computation, it's already applied by the assembly blocks --- lib/mito/fem/elements/tri1/IsoparametricTriangleP1.h | 6 +++--- lib/mito/fem/elements/tri2/IsoparametricTriangleP2.h | 5 +++-- 2 files changed, 6 insertions(+), 5 deletions(-) diff --git a/lib/mito/fem/elements/tri1/IsoparametricTriangleP1.h b/lib/mito/fem/elements/tri1/IsoparametricTriangleP1.h index ece07580..d555e30c 100644 --- a/lib/mito/fem/elements/tri1/IsoparametricTriangleP1.h +++ b/lib/mito/fem/elements/tri1/IsoparametricTriangleP1.h @@ -118,7 +118,8 @@ namespace mito::fem { } // volume element: contract the volume form with the two tangent vectors - // this follows the same pattern as Manifold::_volume + // 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) { @@ -133,8 +134,7 @@ namespace mito::fem { auto tangent_1 = pyre::tensor::col<1>(J); // contract the volume form with both tangent vectors - // for N=2, include factorial (1/2!) - return (1.0 / 2.0) * this->_volume_form(x)(tangent_0, tangent_1); + return this->_volume_form(x)(tangent_0, tangent_1); }); } diff --git a/lib/mito/fem/elements/tri2/IsoparametricTriangleP2.h b/lib/mito/fem/elements/tri2/IsoparametricTriangleP2.h index 4c7a001e..7a3c24bb 100644 --- a/lib/mito/fem/elements/tri2/IsoparametricTriangleP2.h +++ b/lib/mito/fem/elements/tri2/IsoparametricTriangleP2.h @@ -124,6 +124,8 @@ namespace mito::fem { } // 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) { @@ -138,8 +140,7 @@ namespace mito::fem { auto tangent_1 = pyre::tensor::col<1>(J); // contract the volume form with both tangent vectors - // for N=2, include factorial (1/2!) - return (1.0 / 2.0) * this->_volume_form(x)(tangent_0, tangent_1); + return this->_volume_form(x)(tangent_0, tangent_1); }); } From ff1f5c45f2872b0a32c00db1fd59c9b22f57ff98 Mon Sep 17 00:00:00 2001 From: Pawel024 Date: Wed, 28 Jan 2026 13:56:00 +0100 Subject: [PATCH 21/24] tests/fem: use a clockwise-rotated normal field to get positive volume form --- tests/mito.lib/fem/block_grad_grad_embedded_segment.cc | 4 +++- tests/mito.lib/fem/block_mass_embedded_segment.cc | 5 ++++- tests/mito.lib/fem/isoparametric_embedded_segment.cc | 7 ++++--- 3 files changed, 11 insertions(+), 5 deletions(-) diff --git a/tests/mito.lib/fem/block_grad_grad_embedded_segment.cc b/tests/mito.lib/fem/block_grad_grad_embedded_segment.cc index d523964c..d45c8d90 100644 --- a/tests/mito.lib/fem/block_grad_grad_embedded_segment.cc +++ b/tests/mito.lib/fem/block_grad_grad_embedded_segment.cc @@ -43,8 +43,10 @@ TEST(Fem, BlockGradGradEmbeddedSegment) 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( - mito::tensor::vector_t<2>{ -inv_sqrt2, inv_sqrt2 }); + mito::tensor::vector_t<2>{ inv_sqrt2, -inv_sqrt2 }); // create the submanifold auto manifold = mito::manifolds::submanifold(mesh, coord_system, normal_field); diff --git a/tests/mito.lib/fem/block_mass_embedded_segment.cc b/tests/mito.lib/fem/block_mass_embedded_segment.cc index e1327fc6..c5f4082f 100644 --- a/tests/mito.lib/fem/block_mass_embedded_segment.cc +++ b/tests/mito.lib/fem/block_mass_embedded_segment.cc @@ -43,7 +43,10 @@ TEST(Fem, BlockMassEmbeddedSegment) mesh.insert({ node_0, node_1 }); // create a submanifold (1D embedded in 2D requires normal field) - auto normal_field = mito::functions::constant(mito::tensor::vector_t<2>{ -inv_sqrt2, inv_sqrt2 }); + // 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( + mito::tensor::vector_t<2>{ inv_sqrt2, -inv_sqrt2 }); auto manifold = mito::manifolds::submanifold(mesh, coord_system, normal_field); { diff --git a/tests/mito.lib/fem/isoparametric_embedded_segment.cc b/tests/mito.lib/fem/isoparametric_embedded_segment.cc index d67d224c..ae3a22d1 100644 --- a/tests/mito.lib/fem/isoparametric_embedded_segment.cc +++ b/tests/mito.lib/fem/isoparametric_embedded_segment.cc @@ -129,9 +129,10 @@ TEST(Fem, IsoparametricEmbeddedSegment) 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) - auto normal_field = mito::functions::constant( - mito::tensor::vector_t<2>{ -0.8, 0.6 }); + // 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(mito::tensor::vector_t<2>{ 0.8, -0.6 }); // create the submanifold auto manifold = mito::manifolds::submanifold(mesh, coord_system, normal_field); From 4a1daa8d43d80b2682b8f4eb55151789d47956c4 Mon Sep 17 00:00:00 2001 From: Pawel024 Date: Wed, 28 Jan 2026 15:27:56 +0100 Subject: [PATCH 22/24] benchmarks: update the poisson benchmark to use the new elements api --- benchmarks/mito.lib/pdes/poisson.cc | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/benchmarks/mito.lib/pdes/poisson.cc b/benchmarks/mito.lib/pdes/poisson.cc index 13ceb495..6d50cac0 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; // 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; + // get the boundary mesh auto boundary_mesh = mito::mesh::boundary(mesh); From 65259dc057a5152a5648fced8758adb0388191cf Mon Sep 17 00:00:00 2001 From: Pawel024 Date: Wed, 28 Jan 2026 16:50:01 +0100 Subject: [PATCH 23/24] fem: update comments to use "n-simplex" instead of "nth order simplex" when talking about the topological dimension --- lib/mito/fem/elements/IsoparametricSegment.h | 5 ++--- lib/mito/fem/elements/IsoparametricTriangle.h | 2 +- lib/mito/fem/elements/seg1/IsoparametricSegmentP1.h | 2 +- lib/mito/fem/elements/tri1/IsoparametricTriangleP1.h | 2 +- lib/mito/fem/elements/tri2/IsoparametricTriangleP2.h | 2 +- 5 files changed, 6 insertions(+), 7 deletions(-) diff --git a/lib/mito/fem/elements/IsoparametricSegment.h b/lib/mito/fem/elements/IsoparametricSegment.h index 7b8a6a12..84e5eff4 100644 --- a/lib/mito/fem/elements/IsoparametricSegment.h +++ b/lib/mito/fem/elements/IsoparametricSegment.h @@ -8,7 +8,7 @@ // DESIGN NOTES -// Class {IsoparametricSegment} represents a first order simplex (segment) equipped with parametric +// 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) @@ -26,8 +26,7 @@ namespace mito::fem { // 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) + // the Jacobian is a D×N matrix (for a segment: D×1, a column vector) using jacobian_type = tensor::matrix_t; // the metric space provides g, g_inv for this coordinate system diff --git a/lib/mito/fem/elements/IsoparametricTriangle.h b/lib/mito/fem/elements/IsoparametricTriangle.h index 94a7f085..76cf08b0 100644 --- a/lib/mito/fem/elements/IsoparametricTriangle.h +++ b/lib/mito/fem/elements/IsoparametricTriangle.h @@ -8,7 +8,7 @@ // DESIGN NOTES -// Class {IsoparametricTriangle} represents a second order simplex equipped with barycentric +// 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) diff --git a/lib/mito/fem/elements/seg1/IsoparametricSegmentP1.h b/lib/mito/fem/elements/seg1/IsoparametricSegmentP1.h index 9291031a..ed24a93a 100644 --- a/lib/mito/fem/elements/seg1/IsoparametricSegmentP1.h +++ b/lib/mito/fem/elements/seg1/IsoparametricSegmentP1.h @@ -8,7 +8,7 @@ // DESIGN NOTES -// Class {IsoparametricSegmentP1} represents a first order simplex (segment) equipped with linear +// 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.) diff --git a/lib/mito/fem/elements/tri1/IsoparametricTriangleP1.h b/lib/mito/fem/elements/tri1/IsoparametricTriangleP1.h index d555e30c..552df9ec 100644 --- a/lib/mito/fem/elements/tri1/IsoparametricTriangleP1.h +++ b/lib/mito/fem/elements/tri1/IsoparametricTriangleP1.h @@ -8,7 +8,7 @@ // DESIGN NOTES -// Class {IsoparametricTriangleP1} represents a second order simplex equipped with linear shape +// 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.) diff --git a/lib/mito/fem/elements/tri2/IsoparametricTriangleP2.h b/lib/mito/fem/elements/tri2/IsoparametricTriangleP2.h index 7a3c24bb..8a2e72f6 100644 --- a/lib/mito/fem/elements/tri2/IsoparametricTriangleP2.h +++ b/lib/mito/fem/elements/tri2/IsoparametricTriangleP2.h @@ -8,7 +8,7 @@ // DESIGN NOTES -// Class {IsoparametricTriangleP2} represents a second order simplex equipped with quadratic shape +// 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.) From 0f27c1b70caa85666dc988519d9fdb64a437139c Mon Sep 17 00:00:00 2001 From: Pawel024 Date: Thu, 5 Feb 2026 10:45:56 +0100 Subject: [PATCH 24/24] test/fem: test integration with volume form on arbitrary segments --- .cmake/mito_tests_mito_lib.cmake | 1 + .../euclidean_submanifold_metric_2D.cc | 180 ++++++++++++++++++ 2 files changed, 181 insertions(+) create mode 100644 tests/mito.lib/geometry/euclidean_submanifold_metric_2D.cc diff --git a/.cmake/mito_tests_mito_lib.cmake b/.cmake/mito_tests_mito_lib.cmake index 215c0949..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) 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 +#include +#include +#include + + +// 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 & 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(); + + // 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 & 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(); + + // 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 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