diff --git a/.cmake/mito_tests_mito_lib.cmake b/.cmake/mito_tests_mito_lib.cmake index 30faf55a..3a4f1f64 100644 --- a/.cmake/mito_tests_mito_lib.cmake +++ b/.cmake/mito_tests_mito_lib.cmake @@ -35,6 +35,7 @@ mito_test_driver(tests/mito.lib/geometry/point.cc) mito_test_driver(tests/mito.lib/geometry/euclidean_metric_1D.cc) mito_test_driver(tests/mito.lib/geometry/euclidean_metric_2D.cc) mito_test_driver(tests/mito.lib/geometry/euclidean_metric_3D.cc) +mito_test_driver(tests/mito.lib/geometry/euclidean_submanifold_metric_2D.cc) mito_test_driver(tests/mito.lib/geometry/euclidean_submanifold_metric_3D.cc) mito_test_driver(tests/mito.lib/geometry/induced_metric_segment.cc) mito_test_driver(tests/mito.lib/geometry/induced_metric_triangle.cc) @@ -56,8 +57,12 @@ mito_test_driver(tests/mito.lib/discrete/mesh_field.cc) # fem mito_test_driver(tests/mito.lib/fem/block_grad_grad.cc) mito_test_driver(tests/mito.lib/fem/block_grad_grad_segment.cc) +mito_test_driver(tests/mito.lib/fem/block_grad_grad_embedded_segment.cc) +mito_test_driver(tests/mito.lib/fem/block_grad_grad_embedded_triangle.cc) mito_test_driver(tests/mito.lib/fem/block_mass.cc) mito_test_driver(tests/mito.lib/fem/block_mass_segment.cc) +mito_test_driver(tests/mito.lib/fem/block_mass_embedded_segment.cc) +mito_test_driver(tests/mito.lib/fem/block_mass_embedded_triangle.cc) mito_test_driver(tests/mito.lib/fem/block_algebra.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) @@ -66,6 +71,8 @@ mito_test_driver(tests/mito.lib/fem/isoparametric_triangle.cc) mito_test_driver(tests/mito.lib/fem/fem_field.cc) mito_test_driver(tests/mito.lib/fem/shape_functions_segment_p1.cc) mito_test_driver(tests/mito.lib/fem/isoparametric_segment.cc) +mito_test_driver(tests/mito.lib/fem/isoparametric_embedded_segment.cc) +mito_test_driver(tests/mito.lib/fem/isoparametric_embedded_triangle.cc) # io mito_test_driver(tests/mito.lib/io/summit_mesh_reader_2D.cc) @@ -189,8 +196,10 @@ mito_test_driver(tests/mito.lib/quadrature/quadrature_parametric_segment.cc) mito_test_driver(tests/mito.lib/quadrature/quadrature_parametric_triangle.cc) mito_test_driver(tests/mito.lib/quadrature/quadrature_parametric_tetrahedron.cc) mito_test_driver(tests/mito.lib/quadrature/quadrature_segment_1D.cc) +mito_test_driver(tests/mito.lib/quadrature/quadrature_segment_2D.cc) mito_test_driver(tests/mito.lib/quadrature/quadrature_segment_3D.cc) mito_test_driver(tests/mito.lib/quadrature/quadrature_triangle_2D.cc) +mito_test_driver(tests/mito.lib/quadrature/quadrature_triangle_3D.cc) mito_test_driver(tests/mito.lib/quadrature/quadrature_square_3D.cc) # integration diff --git a/.github/workflows/cmake-macos.yaml b/.github/workflows/cmake-macos.yaml index 5ea6c67b..6d7e92f2 100644 --- a/.github/workflows/cmake-macos.yaml +++ b/.github/workflows/cmake-macos.yaml @@ -82,7 +82,7 @@ jobs: repository: pyre/pyre fetch-depth: 0 path: ${{github.workspace}}/pyre - ref: tensor + ref: main - name: retrieve pyre commit SHA id: pyre-sha diff --git a/.github/workflows/cmake-ubuntu.yaml b/.github/workflows/cmake-ubuntu.yaml index 448a6e81..38edd138 100644 --- a/.github/workflows/cmake-ubuntu.yaml +++ b/.github/workflows/cmake-ubuntu.yaml @@ -57,7 +57,7 @@ jobs: repository: pyre/pyre fetch-depth: 0 path: ${{github.workspace}}/pyre - ref: tensor + ref: main - name: retrieve pyre commit SHA id: pyre-sha diff --git a/lib/mito/fem/blocks/GradientGradientBlock.h b/lib/mito/fem/blocks/GradientGradientBlock.h index a3db20a9..16b96a39 100644 --- a/lib/mito/fem/blocks/GradientGradientBlock.h +++ b/lib/mito/fem/blocks/GradientGradientBlock.h @@ -23,10 +23,6 @@ namespace mito::fem::blocks { // the type of the coefficient field using coefficient_field_type = coefficientFieldT; - public: - // instantiate the quadrature rule - static constexpr auto quadrature_rule = quadrature_rule_type(); - public: // constructor GradientGradientBlock(const coefficient_field_type & coefficient) : @@ -39,53 +35,40 @@ namespace mito::fem::blocks { requires element_of_type_c auto compute(const elementT & element) const -> elementary_shape { - // the number of nodes per element - constexpr int n_nodes = element_type::n_nodes; - - // the number of quadrature points per element - constexpr int n_quads = quadrature_rule_type::npoints; + // the parametric coordinates type + using parametric_coordinates_type = typename elementT::parametric_coordinates_type; // the elementary matrix - elementary_shape elementary_matrix{}; - - // loop on the quadrature points - tensor::constexpr_for_1([&]() { - // the parametric coordinates of the quadrature point - constexpr auto xi = quadrature_rule.point(q); - - // the coordinates of the quadrature point - auto x = element.parametrization()(xi); - - // evaluate the coefficient field at the quadrature point - auto coefficient = _coefficient(x); + return manifolds::cell_integrator(element.element()) + .integrate(mito::functions::function([&](const parametric_coordinates_type & xi) { + // the elementary contribution at quadrature point {xi} + elementary_shape elementary_matrix{}; - // the measure of the canonical simplex - constexpr auto measure = - element_type::mesh_cell_type::reference_simplex_type::measure; + // the number of nodes per element + constexpr int n_nodes = element_type::n_nodes; - // the quadrature weight at this point scaled with the area of the canonical simplex - constexpr auto w = measure * quadrature_rule.weight(q); + // the coordinates of the quadrature point + auto x = element.parametrization()(xi); - // precompute the common factor - auto factor = w * tensor::determinant(element.jacobian()(xi)); + // evaluate the coefficient at the quadrature point + auto coefficient = _coefficient(x); - // loop on the nodes of the element - tensor::constexpr_for_1([&]() { - // evaluate the spatial gradient of the element's a-th shape function at {xi} - auto dphi_a = element.template gradient()(xi); // loop on the nodes of the element - tensor::constexpr_for_1([&]() { - // evaluate the spatial gradient of the element's b-th shape function at - // {xi} - auto dphi_b = element.template gradient()(xi); - // populate the elementary contribution to the matrix - elementary_matrix[{ a, b }] += factor * dphi_a * (coefficient * dphi_b); + tensor::constexpr_for_1([&]() { + // evaluate the spatial gradient of the a-th shape function at {xi} + auto dphi_a = element.template gradient()(xi); + // loop on the nodes of the element + tensor::constexpr_for_1([&]() { + // evaluate the spatial gradient of the b-th shape function at {xi} + auto dphi_b = element.template gradient()(xi); + // populate the elementary contribution to the matrix + elementary_matrix[{ a, b }] = dphi_a * (coefficient * dphi_b); + }); }); - }); - }); - // all done - return elementary_matrix; + // all done + return elementary_matrix; + })); } private: diff --git a/lib/mito/fem/blocks/L2NormBlock.h b/lib/mito/fem/blocks/L2NormBlock.h index f7f46b49..919d78c5 100644 --- a/lib/mito/fem/blocks/L2NormBlock.h +++ b/lib/mito/fem/blocks/L2NormBlock.h @@ -16,18 +16,16 @@ namespace mito::fem::blocks { class L2NormBlock { public: - // my template parameters + // my finite element type using element_type = finiteElementT; - using elementary_shape = tensor::scalar_t; + // my quadrature rule using quadrature_rule_type = quadratureRuleT; + // my elementary shape + using elementary_shape = tensor::scalar_t; // the type of the function to compute the L2 norm of using function_type = functionT; - public: - // instantiate the quadrature rule - static constexpr auto quadrature_rule = quadrature_rule_type(); - public: // constructor L2NormBlock(const function_type & function) : _function(function) {} @@ -38,33 +36,24 @@ namespace mito::fem::blocks { requires element_of_type_c auto compute(const elementT & element) const -> elementary_shape { - // the number of quadrature points per element - constexpr int n_quads = quadrature_rule_type::npoints; - - // the elementary contribution to the L2 norm - elementary_shape norm{}; - - // loop on the quadrature points - tensor::constexpr_for_1([&]() { - // the parametric coordinates of the quadrature point - constexpr auto xi = quadrature_rule.point(q); - - // the measure of the canonical simplex - constexpr auto measure = - element_type::mesh_cell_type::reference_simplex_type::measure; + // the parametric coordinates type + using parametric_coordinates_type = typename elementT::parametric_coordinates_type; - // the quadrature weight at this point scaled with the area of the canonical simplex - constexpr auto w = measure * quadrature_rule.weight(q); + // the elementary matrix + return manifolds::cell_integrator(element.element()) + .integrate(mito::functions::function([&](const parametric_coordinates_type & xi) { + // the elementary contribution at quadrature point {xi} + elementary_shape norm{}; - // precompute the common factor - auto factor = w * tensor::determinant(element.jacobian()(xi)); + // evaluate the function at the quadrature point + auto fx = _function(xi); - // populate the elementary contribution to the matrix - norm += factor * _function(xi) * _function(xi); - }); + // assemble the elementary contribution + norm = fx * fx; - // all done - return norm; + // all done + return norm; + })); } private: diff --git a/lib/mito/fem/blocks/ValueBlock.h b/lib/mito/fem/blocks/ValueBlock.h index 678c5f05..98c0ef5c 100644 --- a/lib/mito/fem/blocks/ValueBlock.h +++ b/lib/mito/fem/blocks/ValueBlock.h @@ -25,10 +25,6 @@ namespace mito::fem::blocks { // the type of the coefficient field using coefficient_field_type = coefficientFieldT; - public: - // instantiate the quadrature rule - static constexpr auto quadrature_rule = quadrature_rule_type(); - public: // constructor ValueBlock(const coefficient_field_type & coefficient) : _coefficient(coefficient) {} @@ -39,44 +35,35 @@ namespace mito::fem::blocks { requires(element_of_type_c) auto compute(const elementT & element) const -> elementary_shape { - // the number of nodes per element - constexpr int n_nodes = element_type::n_nodes; - - // the number of quadrature points per element - constexpr int n_quads = quadrature_rule_type::npoints; + // the parametric coordinates type + using parametric_coordinates_type = typename elementT::parametric_coordinates_type; // the elementary vector - elementary_shape elementary_vector{}; - - // loop on the quadrature points - tensor::constexpr_for_1([&]() { - // the parametric coordinates of the quadrature point - constexpr auto xi = quadrature_rule.point(q); - - // the coordinates of the quadrature point - auto coord = element.parametrization()(xi); - - // the measure of the canonical simplex - constexpr auto measure = - element_type::mesh_cell_type::reference_simplex_type::measure; - - // the quadrature weight at this point scaled with the area of the canonical simplex - constexpr auto w = measure * quadrature_rule.weight(q); - - // precompute the common factor - auto factor = w * _coefficient(coord) * tensor::determinant(element.jacobian()(xi)); - - // loop on the nodes of the element - tensor::constexpr_for_1([&]() { - // evaluate the a-th shape function at {xi} - auto phi_a = element.template shape()(xi); - // populate the elementary contribution to the vector - elementary_vector[{ a }] += factor * phi_a; - }); - }); - - // all done - return elementary_vector; + return manifolds::cell_integrator(element.element()) + .integrate(mito::functions::function([&](const parametric_coordinates_type & xi) { + // the elementary contribution at quadrature point {xi} + elementary_shape elementary_vector{}; + + // the number of nodes per element + constexpr int n_nodes = element_type::n_nodes; + + // the coordinates of the quadrature point + auto x = element.parametrization()(xi); + + // evaluate the coefficient at the quadrature point + auto coefficient = _coefficient(x); + + // loop on the nodes of the element + tensor::constexpr_for_1([&]() { + // evaluate the element's a-th shape function at {xi} + const auto phi_a = element.template shape()(xi); + // populate the elementary contribution to the vector + elementary_vector[{ a }] = coefficient * phi_a; + }); + + // all done + return elementary_vector; + })); } private: diff --git a/lib/mito/fem/blocks/ValueGradientBlock.h b/lib/mito/fem/blocks/ValueGradientBlock.h index ef401123..1f24a7c6 100644 --- a/lib/mito/fem/blocks/ValueGradientBlock.h +++ b/lib/mito/fem/blocks/ValueGradientBlock.h @@ -23,10 +23,6 @@ namespace mito::fem::blocks { // the type of the coefficient field using coefficient_field_type = coefficientFieldT; - public: - // instantiate the quadrature rule - static constexpr auto quadrature_rule = quadrature_rule_type(); - public: // constructor ValueGradientBlock(const coefficient_field_type & coefficient) : _coefficient(coefficient) @@ -38,53 +34,40 @@ namespace mito::fem::blocks { requires element_of_type_c auto compute(const elementT & element) const -> elementary_shape { - // the number of nodes per element - constexpr int n_nodes = element_type::n_nodes; - - // the number of quadrature points per element - constexpr int n_quads = quadrature_rule_type::npoints; + // the parametric coordinates type + using parametric_coordinates_type = typename elementT::parametric_coordinates_type; // the elementary matrix - elementary_shape elementary_matrix{}; - - // loop on the quadrature points - tensor::constexpr_for_1([&]() { - // the parametric coordinates of the quadrature point - constexpr auto xi = quadrature_rule.point(q); - - // the coordinates of the quadrature point - auto x = element.parametrization()(xi); - - // evaluate the coefficient at the quadrature point - auto coefficient = _coefficient(x); + return manifolds::cell_integrator(element.element()) + .integrate(mito::functions::function([&](const parametric_coordinates_type & xi) { + // the elementary contribution at quadrature point {xi} + elementary_shape elementary_matrix{}; - // the measure of the canonical simplex - constexpr auto measure = - element_type::mesh_cell_type::reference_simplex_type::measure; + // the number of nodes per element + constexpr int n_nodes = element_type::n_nodes; - // the quadrature weight at this point scaled with the area of the canonical simplex - constexpr auto w = measure * quadrature_rule.weight(q); + // the coordinates of the quadrature point + auto x = element.parametrization()(xi); - // precompute the common factor - auto factor = w * tensor::determinant(element.jacobian()(xi)); + // evaluate the coefficient at the quadrature point + auto coefficient = _coefficient(x); - // loop on the nodes of the element - tensor::constexpr_for_1([&]() { - // evaluate the element's a-th shape function at {xi} - auto phi_a = element.template shape()(xi); // loop on the nodes of the element - tensor::constexpr_for_1([&]() { - // evaluate the spatial gradient of the element's b-th shape function at - // {xi} - auto dphi_b = element.template gradient()(xi); - // populate the elementary contribution to the matrix - elementary_matrix[{ a, b }] += factor * phi_a * coefficient * dphi_b; + tensor::constexpr_for_1([&]() { + // evaluate the element's a-th shape function at {xi} + const auto phi_a = element.template shape()(xi); + // loop on the nodes of the element + tensor::constexpr_for_1([&]() { + // evaluate the spatial gradient of the b-th shape function at {xi} + auto dphi_b = element.template gradient()(xi); + // populate the elementary contribution to the matrix + elementary_matrix[{ a, b }] = phi_a * coefficient * dphi_b; + }); }); - }); - }); - // all done - return elementary_matrix; + // all done + return elementary_matrix; + })); } private: diff --git a/lib/mito/fem/blocks/ValueValueBlock.h b/lib/mito/fem/blocks/ValueValueBlock.h index 4a936f96..c46813f9 100644 --- a/lib/mito/fem/blocks/ValueValueBlock.h +++ b/lib/mito/fem/blocks/ValueValueBlock.h @@ -23,10 +23,6 @@ namespace mito::fem::blocks { // the type of the coefficient field using coefficient_field_type = coefficientFieldT; - public: - // instantiate the quadrature rule - static constexpr auto quadrature_rule = quadrature_rule_type(); - public: // constructor ValueValueBlock(const coefficient_field_type & coefficient) : _coefficient(coefficient) {} @@ -37,53 +33,40 @@ namespace mito::fem::blocks { requires(element_of_type_c) auto compute(const elementT & element) const -> elementary_shape { - // the number of nodes per element - constexpr int n_nodes = element_type::n_nodes; - - // the number of quadrature points per element - constexpr int n_quads = quadrature_rule_type::npoints; + // the parametric coordinates type + using parametric_coordinates_type = typename elementT::parametric_coordinates_type; // the elementary matrix - elementary_shape elementary_matrix{}; - - // loop on the quadrature points - tensor::constexpr_for_1([&]() { - // the parametric coordinates of the quadrature point - constexpr auto xi = quadrature_rule.point(q); - - // the coordinates of the quadrature point - auto x = element.parametrization()(xi); - - // evaluate the coefficient at the quadrature point - auto coefficient = _coefficient(x); + return manifolds::cell_integrator(element.element()) + .integrate(mito::functions::function([&](const parametric_coordinates_type & xi) { + // the elementary contribution at quadrature point {xi} + elementary_shape elementary_matrix{}; - // the measure of the canonical simplex - constexpr auto measure = - element_type::mesh_cell_type::reference_simplex_type::measure; + // the number of nodes per element + constexpr int n_nodes = element_type::n_nodes; - // the quadrature weight at this point scaled with the area of the canonical simplex - constexpr auto w = measure * quadrature_rule.weight(q); + // the coordinates of the quadrature point + auto x = element.parametrization()(xi); - // precompute the common factor - auto factor = coefficient * w * tensor::determinant(element.jacobian()(xi)); + // evaluate the coefficient at the quadrature point + auto coefficient = _coefficient(x); - // loop on the nodes of the element - tensor::constexpr_for_1([&]() { - // evaluate the spatial gradient of the element's a-th shape function at {xi} - auto phi_a = element.template shape()(xi); // loop on the nodes of the element - tensor::constexpr_for_1([&]() { - // evaluate the spatial gradient of the element's b-th shape function at - // {xi} - auto phi_b = element.template shape()(xi); - // populate the elementary contribution to the matrix - elementary_matrix[{ a, b }] += factor * phi_a * phi_b; + tensor::constexpr_for_1([&]() { + // evaluate the element's a-th shape function at {xi} + const auto phi_a = element.template shape()(xi); + // loop on the nodes of the element + tensor::constexpr_for_1([&]() { + // evaluate the element's b-th shape function at {xi} + const auto phi_b = element.template shape()(xi); + // populate the elementary contribution to the matrix + elementary_matrix[{ a, b }] = coefficient * phi_a * phi_b; + }); }); - }); - }); - // all done - return elementary_matrix; + // all done + return elementary_matrix; + })); } private: diff --git a/lib/mito/fem/elements/forward.h b/lib/mito/fem/elements/forward.h index 9ecbc31a..4fe226be 100644 --- a/lib/mito/fem/elements/forward.h +++ b/lib/mito/fem/elements/forward.h @@ -13,9 +13,10 @@ namespace mito::fem { template 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; // concept of {parametrizedElementT} being compatible with finite elements of type diff --git a/lib/mito/fem/elements/seg1/DiscretizerCG.h b/lib/mito/fem/elements/seg1/DiscretizerCG.h index c77b391f..de49226e 100644 --- a/lib/mito/fem/elements/seg1/DiscretizerCG.h +++ b/lib/mito/fem/elements/seg1/DiscretizerCG.h @@ -9,9 +9,9 @@ namespace mito::fem { - // discretizer specialization for a segment in 1D with order 1 continuous Galerkin - template <> - struct Discretizer, 1>, discretization_t::CG> { + // discretizer specialization for a segment in D dimensions with order 1 continuous Galerkin + template + struct Discretizer, 1>, discretization_t::CG> { template < typename manifoldT, typename constraintsT, typename connectivity_table_type, typename map_type, typename constrained_nodes_type> @@ -21,7 +21,7 @@ namespace mito::fem { constrained_nodes_type & constrained_nodes) { // the finite element type - using finite_element_type = finite_element_family, 1>; + using finite_element_type = finite_element_family, 1>; // the discretization node type using discretization_node_type = typename finite_element_type::discretization_node_type; diff --git a/lib/mito/fem/elements/seg1/IsoparametricSegmentP1.h b/lib/mito/fem/elements/seg1/IsoparametricSegmentP1.h index 09760dcf..60bc52bb 100644 --- a/lib/mito/fem/elements/seg1/IsoparametricSegmentP1.h +++ b/lib/mito/fem/elements/seg1/IsoparametricSegmentP1.h @@ -8,8 +8,8 @@ // DESIGN NOTES -// Class {IsoparametricSegmentP1} represents a first order simplex (segment) living in 1D cartesian -// space, equipped with linear shape functions defined in the parametric space. +// Class {IsoparametricSegmentP1} represents a 1-simplex (segment) embedded in a D-dimensional +// physical space, equipped with linear shape functions defined in the parametric space. namespace mito::fem { @@ -22,6 +22,14 @@ namespace mito::fem { using parametrized_element_type = parametrizedElementT; // the underlying mesh cell type using mesh_cell_type = typename parametrized_element_type::cell_type; + // the metric volume form type + using metric_volume_form_type = typename parametrized_element_type::metric_volume_form_type; + // the coordinates type of the physical space + using coordinates_type = typename metric_volume_form_type::input_type; + // the dimension of the physical space + static constexpr int dim = coordinates_type::dim; + // the (euclidean) metric of the physical space in {coordinates_type} coordinates + using ambient_metric_type = geometry::euclidean_metric; // the degree of the finite element static constexpr int degree = 1; @@ -75,6 +83,13 @@ namespace mito::fem { return _connectivity; } + // get the parmetrized element + constexpr auto element() const noexcept + { + // all done + return _element; + } + // get the element parmetrization constexpr auto parametrization() const noexcept { @@ -94,48 +109,19 @@ namespace mito::fem { 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<1> { - // get the shape functions derivatives - constexpr auto dphi_0 = shape_functions.dshape<0>(); - constexpr auto dphi_1 = shape_functions.dshape<1>(); - - // store the coordinates of the vertices of the triangle in physical space - auto x0 = _element.parametrization()({ 0.0 }); - auto x1 = _element.parametrization()({ 1.0 }); - - // compute the jacobian of the isoparametric mapping: dx/dxi - auto dx_dxi = x0 * dphi_0(xi) + x1 * dphi_1(xi); - // wrap the result in a 1x1 matrix - return tensor::matrix_t<1>{ dx_dxi }; - }); - - // 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<1> { - // the jacobian of the mapping from the reference element to the physical - // element evaluated at {xi} - auto J = jacobian()(xi); - // the derivative of the coordinates with respect to the parametric coordinates - auto J_inv = tensor::inverse(J); - // return the spatial gradients of the shape functions evaluated at {xi} - return shape_functions.dshape()(xi) * J_inv; - }); - // and return it - return gradient_function; + // the pullback of the ambient metric along the element parametrization + auto g = + geometry::pullback_metric::field(_element.parametrization()); + // the intrinsic gradient of the a-th shape function on the parametric space + auto dphi = operators::gradient(shape_functions.shape(), g); + // push the gradient forward to the physical space with the jacobian + return _element.jacobian() * dphi; } private: diff --git a/lib/mito/fem/elements/seg1/api.h b/lib/mito/fem/elements/seg1/api.h index b61c5c27..4894e473 100644 --- a/lib/mito/fem/elements/seg1/api.h +++ b/lib/mito/fem/elements/seg1/api.h @@ -10,11 +10,11 @@ namespace mito::fem { // specialization of {finite_element_family} for first order segments - template <> - struct finite_element_family, 1> { + template + struct finite_element_family, 1> { // the dimension of the physical space - static constexpr int dim = 1; + static constexpr int dim = D; // the underlying mesh cell type using mesh_cell_type = geometry::segment_t; // the degree of the finite element diff --git a/lib/mito/fem/elements/tri1/DiscretizerCG.h b/lib/mito/fem/elements/tri1/DiscretizerCG.h index 92782591..d33a9530 100644 --- a/lib/mito/fem/elements/tri1/DiscretizerCG.h +++ b/lib/mito/fem/elements/tri1/DiscretizerCG.h @@ -10,8 +10,8 @@ namespace mito::fem { // discretizer specialization for {IsoparametricTriangleP1} with continuous Galerkin - template <> - struct Discretizer, 1>, discretization_t::CG> { + template + struct Discretizer, 1>, discretization_t::CG> { template < typename manifoldT, typename constraintsT, typename connectivity_table_type, typename map_type, typename constrained_nodes_type> @@ -21,7 +21,7 @@ namespace mito::fem { constrained_nodes_type & constrained_nodes) { // the finite element type - using finite_element_type = finite_element_family, 1>; + using finite_element_type = finite_element_family, 1>; // the discretization node type using discretization_node_type = typename finite_element_type::discretization_node_type; diff --git a/lib/mito/fem/elements/tri1/IsoparametricTriangleP1.h b/lib/mito/fem/elements/tri1/IsoparametricTriangleP1.h index c6c56e35..70bf96dc 100644 --- a/lib/mito/fem/elements/tri1/IsoparametricTriangleP1.h +++ b/lib/mito/fem/elements/tri1/IsoparametricTriangleP1.h @@ -8,8 +8,8 @@ // DESIGN NOTES -// Class {IsoparametricTriangleP1} represents a second order simplex (triangle) living in 2D -// cartesian space, equipped with linear shape functions defined in the parametric space. +// Class {IsoparametricTriangleP1} represents a 2-simplex (triangle) embedded in a D-dimensional +// physical space, equipped with linear shape functions defined in the parametric space. namespace mito::fem { @@ -22,6 +22,14 @@ namespace mito::fem { using parametrized_element_type = parametrizedElementT; // the underlying mesh cell type using mesh_cell_type = typename parametrized_element_type::cell_type; + // the metric volume form type + using metric_volume_form_type = typename parametrized_element_type::metric_volume_form_type; + // the coordinates type of the physical space + using coordinates_type = typename metric_volume_form_type::input_type; + // the dimension of the physical space + static constexpr int dim = coordinates_type::dim; + // the (euclidean) metric of the physical space in {coordinates_type} coordinates + using ambient_metric_type = geometry::euclidean_metric; // the degree of the finite element static constexpr int degree = 1; @@ -75,6 +83,13 @@ namespace mito::fem { return _connectivity; } + // get the parmetrized element + constexpr auto element() const noexcept + { + // all done + return _element; + } + // get the element parmetrization constexpr auto parametrization() const noexcept { @@ -94,51 +109,18 @@ namespace mito::fem { 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> { - // get the shape functions derivatives - constexpr auto dphi_0 = shape_functions.dshape<0>(); - constexpr auto dphi_1 = shape_functions.dshape<1>(); - constexpr auto dphi_2 = shape_functions.dshape<2>(); - - // store the coordinates of the vertices of the triangle in physical space - auto x0 = _element.parametrization()({ 0.0, 0.0 }); - auto x1 = _element.parametrization()({ 1.0, 0.0 }); - auto x2 = _element.parametrization()({ 0.0, 1.0 }); - - // compute the jacobian of the isoparametric mapping - return ( - tensor::dyadic(x0, dphi_0(xi)) + tensor::dyadic(x1, dphi_1(xi)) - + tensor::dyadic(x2, dphi_2(xi))); - }); - - // 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); - // 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; - }); - - // and return it - return gradient_function; + // the pullback of the ambient metric along the element parametrization + auto g = + geometry::pullback_metric::field(_element.parametrization()); + // the intrinsic gradient of the a-th shape function on the parametric space + auto dphi = operators::gradient(shape_functions.shape(), g); + // push the gradient forward to the physical space with the jacobian + return _element.jacobian() * dphi; } private: diff --git a/lib/mito/fem/elements/tri1/api.h b/lib/mito/fem/elements/tri1/api.h index d2efc5f9..688b638b 100644 --- a/lib/mito/fem/elements/tri1/api.h +++ b/lib/mito/fem/elements/tri1/api.h @@ -9,11 +9,12 @@ namespace mito::fem { - template <> - struct finite_element_family, 1> { + // specialization of {finite_element_family} for first order triangles + template + struct finite_element_family, 1> { // the dimension of the physical space - static constexpr int dim = 2; + static constexpr int dim = D; // the underlying mesh cell type using mesh_cell_type = geometry::triangle_t; // the degree of the finite element diff --git a/lib/mito/fem/elements/tri2/DiscretizerCG.h b/lib/mito/fem/elements/tri2/DiscretizerCG.h index b6322be3..d14e1bf8 100644 --- a/lib/mito/fem/elements/tri2/DiscretizerCG.h +++ b/lib/mito/fem/elements/tri2/DiscretizerCG.h @@ -10,8 +10,8 @@ namespace mito::fem { // discretizer specialization for {IsoparametricTriangleP2} with continuous Galerkin - template <> - struct Discretizer, 2>, discretization_t::CG> { + template + struct Discretizer, 2>, discretization_t::CG> { template < typename manifoldT, typename constraintsT, typename connectivity_table_type, typename map_type, typename constrained_nodes_type> @@ -21,7 +21,7 @@ namespace mito::fem { constrained_nodes_type & constrained_nodes) { // the finite element type - using finite_element_type = finite_element_family, 2>; + using finite_element_type = finite_element_family, 2>; // assemble the mesh node type using mesh_node_type = typename manifoldT::mesh_type::node_type; diff --git a/lib/mito/fem/elements/tri2/IsoparametricTriangleP2.h b/lib/mito/fem/elements/tri2/IsoparametricTriangleP2.h index 5f6e41aa..a3303a93 100644 --- a/lib/mito/fem/elements/tri2/IsoparametricTriangleP2.h +++ b/lib/mito/fem/elements/tri2/IsoparametricTriangleP2.h @@ -8,8 +8,8 @@ // DESIGN NOTES -// Class {IsoparametricTriangleP2} represents a second order simplex (triangle) living in 2D -// cartesian space, equipped with quadratic shape functions defined in the parametric space. +// Class {IsoparametricTriangleP2} represents a 2-simplex (triangle) embedded in a D-dimensional +// physical space, equipped with quadratic shape functions defined in the parametric space. namespace mito::fem { @@ -22,6 +22,14 @@ namespace mito::fem { using parametrized_element_type = parametrizedElementT; // the underlying mesh cell type using mesh_cell_type = typename parametrized_element_type::cell_type; + // the metric volume form type + using metric_volume_form_type = typename parametrized_element_type::metric_volume_form_type; + // the coordinates type of the physical space + using coordinates_type = typename metric_volume_form_type::input_type; + // the dimension of the physical space + static constexpr int dim = coordinates_type::dim; + // the (euclidean) metric of the physical space in {coordinates_type} coordinates + using ambient_metric_type = geometry::euclidean_metric; // the degree of the finite element static constexpr int degree = 2; @@ -82,6 +90,13 @@ namespace mito::fem { return _element.parametrization(); } + // get the parmetrized element + constexpr auto element() const noexcept + { + // all done + return _element; + } + // get the mesh cell constexpr auto cell() const noexcept -> mesh_cell_type { return _element.cell(); } @@ -94,57 +109,19 @@ namespace mito::fem { 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> { - // store the coordinates of the vertices of the triangle in physical space - auto x0 = _element.parametrization()({ 0.0, 0.0 }); - auto x1 = _element.parametrization()({ 1.0, 0.0 }); - auto x2 = _element.parametrization()({ 0.0, 1.0 }); - auto x3 = 0.5 * (x0 + x1); - auto x4 = 0.5 * (x1 + x2); - auto x5 = 0.5 * (x2 + x0); - - // get the shape functions derivatives - constexpr auto dphi_0 = shape_functions.dshape<0>(); - constexpr auto dphi_1 = shape_functions.dshape<1>(); - constexpr auto dphi_2 = shape_functions.dshape<2>(); - constexpr auto dphi_3 = shape_functions.dshape<3>(); - constexpr auto dphi_4 = shape_functions.dshape<4>(); - constexpr auto dphi_5 = shape_functions.dshape<5>(); - - // compute the jacobian 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))); - }); - - // 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); - // 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; - }); - // and return it - return gradient_function; + // the pullback of the ambient metric along the element parametrization + auto g = + geometry::pullback_metric::field(_element.parametrization()); + // the intrinsic gradient of the a-th shape function on the parametric space + auto dphi = operators::gradient(shape_functions.shape(), g); + // push the gradient forward to the physical space with the jacobian + return _element.jacobian() * dphi; } private: diff --git a/lib/mito/fem/elements/tri2/api.h b/lib/mito/fem/elements/tri2/api.h index 396b5a08..02a8e13e 100644 --- a/lib/mito/fem/elements/tri2/api.h +++ b/lib/mito/fem/elements/tri2/api.h @@ -9,11 +9,12 @@ namespace mito::fem { - template <> - struct finite_element_family, 2> { + // specialization of {finite_element_family} for second order triangles + template + struct finite_element_family, 2> { // the dimension of the physical space - static constexpr int dim = 2; + static constexpr int dim = D; // the underlying mesh cell type using mesh_cell_type = geometry::triangle_t; // the degree of the finite element diff --git a/lib/mito/manifolds/CellIntegrator.h b/lib/mito/manifolds/CellIntegrator.h new file mode 100644 index 00000000..22c72805 --- /dev/null +++ b/lib/mito/manifolds/CellIntegrator.h @@ -0,0 +1,79 @@ +// -*- c++ -*- +// +// Copyright (c) 2020-2026, the MiTo Authors, all rights reserved +// + +// code guard +#pragma once + + +namespace mito::manifolds { + + // Class {CellIntegrator} can integrate arbitrary tensors on a parametrized cell + template + class CellIntegrator { + public: + // the quadrature rule type + using quadrature_rule_type = quadratureRuleT; + // the parametrized element type + using parametrized_element_type = parametrizedElementT; + // the cell type + using cell_type = typename parametrized_element_type::cell_type; + + public: + // instantiate the quadrature rule + static constexpr auto quadrature_rule = quadrature_rule_type(); + + public: + // the constructor + constexpr CellIntegrator(const parametrized_element_type & parametrized_element) : + _parametrized_element(parametrized_element) + {} + + public: + template + auto integrate(F && f) const -> F::output_type + { + // + using result_type = F::output_type; + + // initialize the result + result_type result{}; + + // the number of quadrature points per element + constexpr int n_quads = quadrature_rule_type::npoints; + + // get cell parametrization under the manifold's coordinate system + const auto phi = _parametrized_element.parametrization(); + // compute the derivative of the cell parametrization + const auto J = functions::derivative(phi); + // get the manifold's metric volume form + const auto w = _parametrized_element.metric_volume_form(); + + // loop on the quadrature points + tensor::constexpr_for_1([&]() { + // the parametric coordinates of the quadrature point + constexpr auto xi_q = quadrature_rule.point(q); + // get the quadrature weight and scale it by the reference simplex area + constexpr auto w_q = + quadrature_rule.weight(q) * cell_type::reference_simplex_type::measure; + // construct the metric volume element at {x} by contracting the metric volume + // form with the tangent vectors at {x} + const auto dV = w(phi(xi_q))(tensor::columns(J(xi_q))); + // assemble the elementary contribution + result += f(xi_q) * w_q * dV; + }); + + // all done + return result; + } + + private: + // the parametrized element + const parametrized_element_type & _parametrized_element; + }; + +} // namespace mito + + +// end of file diff --git a/lib/mito/manifolds/ParametrizedElement.h b/lib/mito/manifolds/ParametrizedElement.h index b8b17f53..c5486078 100644 --- a/lib/mito/manifolds/ParametrizedElement.h +++ b/lib/mito/manifolds/ParametrizedElement.h @@ -54,7 +54,13 @@ namespace mito::manifolds { constexpr auto cell() const -> const cell_type & { return _cell; } // return the parametrization of this element in physical space - constexpr auto parametrization() const -> const parametrization_type & { return _parametrization; } + constexpr auto parametrization() const -> const parametrization_type & + { + return _parametrization; + } + + // return the jacobian (gradient in parametric space) of the parametrization of this element + constexpr auto jacobian() const { return operators::gradient(_parametrization); } // return the metric volume form of this element constexpr auto metric_volume_form() const -> const metric_volume_form_type & diff --git a/lib/mito/manifolds/api.h b/lib/mito/manifolds/api.h index 5740a0a6..abdb7b53 100644 --- a/lib/mito/manifolds/api.h +++ b/lib/mito/manifolds/api.h @@ -27,6 +27,14 @@ namespace mito::manifolds { const cellT & cell, const parametrizationT & parametrization, const metricVolumeFormT & metric_volume_form); + // cell integrator alias + template + using cell_integrator_t = CellIntegrator; + + // cell integrator factory + template + constexpr auto cell_integrator(const parametrizedElementT & parametrized_element); + // manifold elements view alias template using manifold_elements_view_t = ManifoldElementsView; diff --git a/lib/mito/manifolds/externals.h b/lib/mito/manifolds/externals.h index 944310ad..ded8e45d 100644 --- a/lib/mito/manifolds/externals.h +++ b/lib/mito/manifolds/externals.h @@ -12,6 +12,7 @@ // support #include "../journal.h" #include "../fields.h" +#include "../operators.h" #include "../mesh.h" #include "../geometry.h" diff --git a/lib/mito/manifolds/factories.h b/lib/mito/manifolds/factories.h index 7836f4df..5fb74677 100644 --- a/lib/mito/manifolds/factories.h +++ b/lib/mito/manifolds/factories.h @@ -28,6 +28,13 @@ namespace mito::manifolds { cell, parametrization, metric_volume_form); } + // cell integrator factory + template + constexpr auto cell_integrator(const parametrizedElementT & parametrized_element) + { + return cell_integrator_t(parametrized_element); + } + // factory manifold template constexpr auto manifold( diff --git a/lib/mito/manifolds/forward.h b/lib/mito/manifolds/forward.h index cb6154e8..3e072c55 100644 --- a/lib/mito/manifolds/forward.h +++ b/lib/mito/manifolds/forward.h @@ -18,6 +18,10 @@ namespace mito::manifolds { template class ParametrizedElement; + // cell integrator + template + class CellIntegrator; + // class manifold elements view template class ManifoldElementsView; diff --git a/lib/mito/manifolds/public.h b/lib/mito/manifolds/public.h index 78fa171a..6216bdd9 100644 --- a/lib/mito/manifolds/public.h +++ b/lib/mito/manifolds/public.h @@ -19,6 +19,7 @@ // classes implementation #include "Atlas.h" #include "ParametrizedElement.h" +#include "CellIntegrator.h" #include "ManifoldElementsView.h" #include "Manifold.h" diff --git a/lib/mito/operators/differential.h b/lib/mito/operators/differential.h index 063957bf..d89645bb 100644 --- a/lib/mito/operators/differential.h +++ b/lib/mito/operators/differential.h @@ -64,6 +64,23 @@ namespace mito::operators { return _grad(field, std::make_index_sequence{}); } + // function to compute the contravariant gradient of a scalar field with respect to a metric + template + requires( + (fields::tensor_field_c or fields::scalar_field_c) + and fields::compatible_fields_c) + constexpr auto gradient(const F & field, const G & metric_field) + { + // if the metric collapses to a scalar + if constexpr (fields::scalar_field_c) { + // divide by the metric directly + return gradient(field) / metric_field; + } else { + // multiply with the inverse of the metric + return functions::inverse(metric_field) * gradient(field); + } + } + // function to compute the divergence of a vector field template constexpr auto divergence(const F & field) diff --git a/lib/mito/quadrature/Integrator.h b/lib/mito/quadrature/Integrator.h index 6be01718..3dfabc9e 100644 --- a/lib/mito/quadrature/Integrator.h +++ b/lib/mito/quadrature/Integrator.h @@ -17,17 +17,14 @@ namespace mito::quadrature { public: // publish my template parameters using manifold_type = manifoldT; + // the cell type of the manifold using cell_type = typename manifold_type::cell_type; + // the parametric coordinates type + using parametric_coordinates_type = typename cell_type::parametric_coordinates_type; + // the reference cell type using reference_cell_type = typename manifold_type::cell_type::reference_simplex_type; - using coordinates_type = typename manifold_type::coordinates_type; - - private: - // quadrature_type, cell_type, and r identify a specific quadrature rule + // assemble to quadrature rule based on quadrature_type, cell_type, and r using quadrature_rule_type = quadrature_rule_t; - // the quadrature rule - static constexpr auto _quadratureRule = quadrature_rule_type(); - // the number of quadrature points - static constexpr int Q = quadrature_rule_type::npoints; public: // the constructor @@ -42,23 +39,10 @@ namespace mito::quadrature { for (const auto & cell : _manifold.elements()) { // get cell parametrization under the manifold's coordinate system const auto phi = cell.parametrization(); - // compute the derivative of the cell parametrization - const auto J = functions::derivative(phi); - // get the manifold's metric volume form - const auto w = cell.metric_volume_form(); - // loop on quadrature points - for (auto q = 0; q < Q; ++q) { - // get the quadrature point coordinates in physical space - const auto x_q = _quadratureRule.point(q); - // get the quadrature weight and scale it by the reference simplex area - const auto w_q = - _quadratureRule.weight(q) * cell_type::reference_simplex_type::measure; - // construct the metric volume element at {x} by contracting the metric volume - // form with the tangent vectors at {x} - const auto dV = w(phi(x_q))(tensor::columns(J(x_q))); - // assemble the elementary contribution - result += f(phi(x_q)) * w_q * dV; - } + // assemble the elementary contribution + result += manifolds::cell_integrator(cell).integrate( + mito::functions::function( + [&](const parametric_coordinates_type & xi) { return f(phi(xi)); })); } // all done 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..96b735e8 --- /dev/null +++ b/tests/mito.lib/fem/block_grad_grad_embedded_segment.cc @@ -0,0 +1,91 @@ +// -*- c++ -*- +// +// Copyright (c) 2020-2026, the MiTo Authors, all rights reserved +// + +#include +#include + + +// the type of coordinates (2D physical space) +using coordinates_t = mito::geometry::coordinates_t<2, mito::geometry::CARTESIAN>; +// the metric space type +using metric_space_t = mito::geometry::euclidean_metric_space; +// 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>; + + +TEST(Fem, BlockGradGradEmbeddedSegment) +{ + // the coordinate system + auto coord_system = mito::geometry::coordinate_system_t(); + + // an atlas under the coordinate system + auto atlas = mito::manifolds::atlas(coord_system); + + // build nodes (unit-length diagonal segment embedded in 2D) + constexpr auto inv_sqrt2 = 1.0 / std::sqrt(2.0); + auto node_0 = mito::geometry::node(coord_system, { 0.0, 0.0 }); + auto node_1 = mito::geometry::node(coord_system, { inv_sqrt2, inv_sqrt2 }); + + // make a geometric simplex + auto segment = mito::geometry::segment(node_0, node_1); + + // the normal field to the segment (perpendicular to the segment direction) + // diagonal direction is (inv_sqrt2, inv_sqrt2), so normal is (inv_sqrt2, -inv_sqrt2) + // (rotated 90° clockwise to get positive orientation with w(normal, tangent) > 0) + auto normal_field = mito::functions::constant( + mito::tensor::vector_t<2>{ inv_sqrt2, -inv_sqrt2 }); + + // strip the namespace from the placeholder for forms contractions + using mito::tensor::_; + + // the ambient metric volume form + constexpr auto w = metric_space_t::w; + + // the restriction of the metric volume form to the segment + auto wS = mito::functions::function( + [w, normal_field](const coordinates_t & x) { return w(x)(normal_field(x), _); }); + + // make a manifold element from the segment + auto element = + mito::manifolds::parametrized_element(segment, atlas.parametrization(segment), wS); + + { + // build the discretization nodes + auto discretization_node_0 = discretization_node_t(); + auto discretization_node_1 = discretization_node_t(); + + // the degree of the finite element + constexpr int degree = 1; + // assemble the finite element type + using finite_element_t = mito::fem::finite_element_family; + + // a finite element on the embedded segment + auto element_p1 = mito::fem::finite_element( + element, { discretization_node_0, discretization_node_1 }); + + // a grad-grad matrix block + auto coefficient = mito::functions::identity(); + constexpr int doe_diffusion = 2 * (finite_element_t::degree - 1); + auto grad_grad_block = + mito::fem::blocks::grad_grad_block(coefficient); + + // the analytical elementary stiffness matrix (same as 1D for unit-length segment) + auto analytical_block = mito::tensor::matrix_t<2>{ 1.0, -1.0, -1.0, 1.0 }; + + // compute the elementary contribution of the block + auto computed_block = grad_grad_block.compute(element_p1); + + // compute the error + auto error = mito::tensor::norm(computed_block - analytical_block); + + // check the error is zero to machine precision + EXPECT_DOUBLE_EQ(0.0, error); + } + + // all done + return; +} diff --git a/tests/mito.lib/fem/block_grad_grad_embedded_triangle.cc b/tests/mito.lib/fem/block_grad_grad_embedded_triangle.cc new file mode 100644 index 00000000..ebddffc3 --- /dev/null +++ b/tests/mito.lib/fem/block_grad_grad_embedded_triangle.cc @@ -0,0 +1,94 @@ +// -*- c++ -*- +// +// Copyright (c) 2020-2026, the MiTo Authors, all rights reserved +// + +#include +#include + + +// the type of coordinates (3D physical space) +using coordinates_t = mito::geometry::coordinates_t<3, mito::geometry::CARTESIAN>; +// the metric space type +using metric_space_t = mito::geometry::euclidean_metric_space; +// the type of discretization node +using discretization_node_t = mito::discrete::discretization_node_t; +// the type of cell (triangle embedded in 3D) +using cell_t = mito::geometry::triangle_t<3>; + + +TEST(Fem, BlockGradGradEmbeddedTriangle) +{ + // the coordinate system + auto coord_system = mito::geometry::coordinate_system_t(); + + // an atlas under the coordinate system + auto atlas = mito::manifolds::atlas(coord_system); + + // build nodes for a unit right triangle rigidly rotated out of the xy-plane + // (the legs {t_1} = (0.6, 0.0, 0.8) and {t_2} = (0.0, 1.0, 0.0) are unit and orthogonal) + auto node_0 = mito::geometry::node(coord_system, { 0.0, 0.0, 0.0 }); + auto node_1 = mito::geometry::node(coord_system, { 0.6, 0.0, 0.8 }); + auto node_2 = mito::geometry::node(coord_system, { 0.0, 1.0, 0.0 }); + + // make a geometric simplex + auto triangle = mito::geometry::triangle(node_0, node_1, node_2); + + // the normal field to the triangle: {n} = {t_1} x {t_2} normalized, + // oriented such that w(n, t_1, t_2) > 0 + auto normal_field = + mito::functions::constant(mito::tensor::vector_t<3>{ -0.8, 0.0, 0.6 }); + + // strip the namespace from the placeholder for forms contractions + using mito::tensor::_; + + // the ambient metric volume form + constexpr auto w = metric_space_t::w; + + // the restriction of the metric volume form to the triangle + auto wS = mito::functions::function( + [w, normal_field](const coordinates_t & x) { return w(x)(normal_field(x), _, _); }); + + // make a manifold element from the triangle + auto element = + mito::manifolds::parametrized_element(triangle, atlas.parametrization(triangle), wS); + + { + // build the discretization nodes + auto discretization_node_0 = discretization_node_t(); + auto discretization_node_1 = discretization_node_t(); + auto discretization_node_2 = discretization_node_t(); + + // the degree of the finite element + constexpr int degree = 1; + // assemble the finite element type + using finite_element_t = mito::fem::finite_element_family; + + // a finite element on the embedded triangle + auto element_p1 = mito::fem::finite_element( + element, { discretization_node_0, discretization_node_1, discretization_node_2 }); + + // a grad-grad matrix block + auto coefficient = mito::functions::identity(); + constexpr int doe_diffusion = 2 * (finite_element_t::degree - 1); + auto grad_grad_block = + mito::fem::blocks::grad_grad_block(coefficient); + + // the analytical elementary stiffness matrix (the stiffness matrix is invariant under + // rigid rotations, so it is the same as for the unit right triangle in 2D) + auto analytical_block = 1.0 / 2.0 * mito::tensor::matrix_t<3>{ 2.0, -1.0, -1.0, -1.0, 1.0, + 0.0, -1.0, 0.0, 1.0 }; + + // compute the elementary contribution of the block + auto computed_block = grad_grad_block.compute(element_p1); + + // compute the error + auto error = mito::tensor::norm(computed_block - analytical_block); + + // check the error is zero to machine precision + EXPECT_NEAR(0.0, error, 1.0e-15); + } + + // all done + return; +} diff --git a/tests/mito.lib/fem/block_mass_embedded_segment.cc b/tests/mito.lib/fem/block_mass_embedded_segment.cc new file mode 100644 index 00000000..e43a92e7 --- /dev/null +++ b/tests/mito.lib/fem/block_mass_embedded_segment.cc @@ -0,0 +1,90 @@ +// -*- c++ -*- +// +// Copyright (c) 2020-2026, the MiTo Authors, all rights reserved +// + +#include +#include + + +// the type of coordinates (2D physical space) +using coordinates_t = mito::geometry::coordinates_t<2, mito::geometry::CARTESIAN>; +// the metric space type +using metric_space_t = mito::geometry::euclidean_metric_space; +// 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>; + + +TEST(Fem, BlockMassEmbeddedSegment) +{ + // the coordinate system + auto coord_system = mito::geometry::coordinate_system_t(); + + // an atlas under the coordinate system + auto atlas = mito::manifolds::atlas(coord_system); + + // build nodes (unit-length diagonal segment embedded in 2D) + constexpr auto inv_sqrt2 = 1.0 / std::sqrt(2.0); + auto node_0 = mito::geometry::node(coord_system, { 0.0, 0.0 }); + auto node_1 = mito::geometry::node(coord_system, { inv_sqrt2, inv_sqrt2 }); + + // make a geometric simplex + auto segment = mito::geometry::segment(node_0, node_1); + + // the normal field to the segment (perpendicular to the segment direction) + // diagonal direction is (inv_sqrt2, inv_sqrt2), so normal is (inv_sqrt2, -inv_sqrt2) + // (rotated 90° clockwise to get positive orientation with w(normal, tangent) > 0) + auto normal_field = mito::functions::constant( + mito::tensor::vector_t<2>{ inv_sqrt2, -inv_sqrt2 }); + + // strip the namespace from the placeholder for forms contractions + using mito::tensor::_; + + // the ambient metric volume form + constexpr auto w = metric_space_t::w; + + // the restriction of the metric volume form to the segment + auto wS = mito::functions::function( + [w, normal_field](const coordinates_t & x) { return w(x)(normal_field(x), _); }); + + // make a manifold element from the segment + auto element = + mito::manifolds::parametrized_element(segment, atlas.parametrization(segment), wS); + + { + // build the discretization nodes + auto discretization_node_0 = discretization_node_t(); + auto discretization_node_1 = discretization_node_t(); + + // the degree of the finite element + constexpr int degree = 1; + // assemble the finite element type + using finite_element_t = mito::fem::finite_element_family; + + // a finite element on the embedded segment + auto element_p1 = mito::fem::finite_element( + element, { discretization_node_0, discretization_node_1 }); + + // a mass matrix block + auto density = mito::functions::one; + constexpr int doe_mass = 2 * finite_element_t::degree; + auto mass_block = mito::fem::blocks::value_value_block(density); + + // the analytical elementary mass matrix (same as 1D for unit-length segment) + auto analytical_block = 1.0 / 6.0 * mito::tensor::matrix_t<2>{ 2.0, 1.0, 1.0, 2.0 }; + + // compute the elementary contribution of the block + auto computed_block = mass_block.compute(element_p1); + + // compute the error + auto error = mito::tensor::norm(computed_block - analytical_block); + + // check the error is reasonable + EXPECT_NEAR(0.0, error, 1.5e-16); + } + + // all done + return; +} diff --git a/tests/mito.lib/fem/block_mass_embedded_triangle.cc b/tests/mito.lib/fem/block_mass_embedded_triangle.cc new file mode 100644 index 00000000..22991d6d --- /dev/null +++ b/tests/mito.lib/fem/block_mass_embedded_triangle.cc @@ -0,0 +1,93 @@ +// -*- c++ -*- +// +// Copyright (c) 2020-2026, the MiTo Authors, all rights reserved +// + +#include +#include + + +// the type of coordinates (3D physical space) +using coordinates_t = mito::geometry::coordinates_t<3, mito::geometry::CARTESIAN>; +// the metric space type +using metric_space_t = mito::geometry::euclidean_metric_space; +// the type of discretization node +using discretization_node_t = mito::discrete::discretization_node_t; +// the type of cell (triangle embedded in 3D) +using cell_t = mito::geometry::triangle_t<3>; + + +TEST(Fem, BlockMassEmbeddedTriangle) +{ + // the coordinate system + auto coord_system = mito::geometry::coordinate_system_t(); + + // an atlas under the coordinate system + auto atlas = mito::manifolds::atlas(coord_system); + + // build nodes for a unit right triangle rigidly rotated out of the xy-plane + // (the legs {t_1} = (0.6, 0.0, 0.8) and {t_2} = (0.0, 1.0, 0.0) are unit and orthogonal) + auto node_0 = mito::geometry::node(coord_system, { 0.0, 0.0, 0.0 }); + auto node_1 = mito::geometry::node(coord_system, { 0.6, 0.0, 0.8 }); + auto node_2 = mito::geometry::node(coord_system, { 0.0, 1.0, 0.0 }); + + // make a geometric simplex + auto triangle = mito::geometry::triangle(node_0, node_1, node_2); + + // the normal field to the triangle: {n} = {t_1} x {t_2} normalized, + // oriented such that w(n, t_1, t_2) > 0 + auto normal_field = + mito::functions::constant(mito::tensor::vector_t<3>{ -0.8, 0.0, 0.6 }); + + // strip the namespace from the placeholder for forms contractions + using mito::tensor::_; + + // the ambient metric volume form + constexpr auto w = metric_space_t::w; + + // the restriction of the metric volume form to the triangle + auto wS = mito::functions::function( + [w, normal_field](const coordinates_t & x) { return w(x)(normal_field(x), _, _); }); + + // make a manifold element from the triangle + auto element = + mito::manifolds::parametrized_element(triangle, atlas.parametrization(triangle), wS); + + { + // build the discretization nodes + auto discretization_node_0 = discretization_node_t(); + auto discretization_node_1 = discretization_node_t(); + auto discretization_node_2 = discretization_node_t(); + + // the degree of the finite element + constexpr int degree = 1; + // assemble the finite element type + using finite_element_t = mito::fem::finite_element_family; + + // a finite element on the embedded triangle + auto element_p1 = mito::fem::finite_element( + element, { discretization_node_0, discretization_node_1, discretization_node_2 }); + + // a mass matrix block + auto density = mito::functions::one; + constexpr int doe_mass = 2 * finite_element_t::degree; + auto mass_block = mito::fem::blocks::value_value_block(density); + + // the analytical elementary mass matrix (the mass matrix is invariant under rigid + // rotations, so it is the same as for the unit right triangle in 2D) + auto analytical_block = + 1.0 / 24.0 * mito::tensor::matrix_t<3>{ 2.0, 1.0, 1.0, 1.0, 2.0, 1.0, 1.0, 1.0, 2.0 }; + + // compute the elementary contribution of the block + auto computed_block = mass_block.compute(element_p1); + + // compute the error + auto error = mito::tensor::norm(computed_block - analytical_block); + + // check the error is zero to machine precision + EXPECT_NEAR(0.0, error, 1.0e-15); + } + + // all done + return; +} diff --git a/tests/mito.lib/fem/isoparametric_embedded_segment.cc b/tests/mito.lib/fem/isoparametric_embedded_segment.cc new file mode 100644 index 00000000..6c0c6692 --- /dev/null +++ b/tests/mito.lib/fem/isoparametric_embedded_segment.cc @@ -0,0 +1,163 @@ +// -*- c++ -*- +// +// Copyright (c) 2020-2026, the MiTo Authors, all rights reserved +// + +#include +#include + + +// the type of coordinates (2D physical space) +using coordinates_t = mito::geometry::coordinates_t<2, mito::geometry::CARTESIAN>; +// the metric space type +using metric_space_t = mito::geometry::euclidean_metric_space; +// the type of discretization node +using discretization_node_t = mito::discrete::discretization_node_t; +// the type of cell (segment embedded in 2D) +using cell_t = mito::geometry::segment_t<2>; +// the reference simplex +using reference_simplex_t = cell_t::reference_simplex_type; +// Gauss quadrature on segments with degree of exactness 2 +using quadrature_rule_t = + mito::quadrature::quadrature_rule_t; + + +// 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 that the sum of the shape functions gradients is the zero vector + EXPECT_NEAR(0.0, mito::tensor::norm(sum), 3.0e-16); + }); + + // all done + return; +} + +TEST(Fem, IsoparametricEmbeddedSegment) +{ + // the coordinate system + auto coord_system = mito::geometry::coordinate_system_t(); + + // an atlas under the coordinate system + auto atlas = mito::manifolds::atlas(coord_system); + + // build nodes for a diagonal segment from (0,0) to (3,4) - length 5 + auto node_0 = mito::geometry::node(coord_system, { 0.0, 0.0 }); + auto node_1 = mito::geometry::node(coord_system, { 3.0, 4.0 }); + + // make a geometric simplex + auto segment = mito::geometry::segment(node_0, node_1); + + // the normal field to the segment (perpendicular to the segment direction) + // segment direction is (3,4)/5 = (0.6, 0.8), so normal is (0.8, -0.6) + // (rotated 90° clockwise to get positive orientation with w(normal, tangent) > 0) + auto normal_field = + mito::functions::constant(mito::tensor::vector_t<2>{ 0.8, -0.6 }); + + // strip the namespace from the placeholder for forms contractions + using mito::tensor::_; + + // the ambient metric volume form + constexpr auto w = metric_space_t::w; + + // the restriction of the metric volume form to the segment + auto wS = mito::functions::function( + [w, normal_field](const coordinates_t & x) { return w(x)(normal_field(x), _); }); + + // make a manifold element from the segment + auto element = + mito::manifolds::parametrized_element(segment, atlas.parametrization(segment), wS); + + { + // build the discretization nodes + auto discretization_node_0 = discretization_node_t(); + auto discretization_node_1 = discretization_node_t(); + + // the degree of the finite element + constexpr int degree = 1; + // assemble the finite element type + using finite_element_t = mito::fem::finite_element_family; + + // first order isoparametric finite element on the embedded segment + auto element_p1 = mito::fem::finite_element( + element, { discretization_node_0, discretization_node_1 }); + + // check that first order shape functions are a partition of unity + test_partition_of_unity(element_p1); + + // check that the gradients of first order shape functions sum to 0.0 + test_gradient_consistency(element_p1); + + // check the gradient values at the midpoint for the (0,0) to (3,4) segment + { + auto xi = reference_simplex_t::parametric_coordinates_type{ 0.5 }; + auto grad_0 = element_p1.gradient<0>()(xi); + auto grad_1 = element_p1.gradient<1>()(xi); + + // for linear shape functions, the gradients should be constant and opposite + // the gradient should be in the direction of the segment: (3/5, 4/5) / 5 = (3/25, 4/25) + // phi_0 decreases from 1 to 0, so grad_0 = (-3/25, -4/25) + // phi_1 increases from 0 to 1, so grad_1 = (3/25, 4/25) + EXPECT_NEAR(-3.0 / 25.0, grad_0[0], 1.0e-15); + EXPECT_NEAR(-4.0 / 25.0, grad_0[1], 1.0e-15); + EXPECT_NEAR(3.0 / 25.0, grad_1[0], 1.0e-15); + EXPECT_NEAR(4.0 / 25.0, grad_1[1], 1.0e-15); + } + } + + // all done + return; +} diff --git a/tests/mito.lib/fem/isoparametric_embedded_triangle.cc b/tests/mito.lib/fem/isoparametric_embedded_triangle.cc new file mode 100644 index 00000000..9495e14f --- /dev/null +++ b/tests/mito.lib/fem/isoparametric_embedded_triangle.cc @@ -0,0 +1,175 @@ +// -*- c++ -*- +// +// Copyright (c) 2020-2026, the MiTo Authors, all rights reserved +// + +#include +#include + + +// the type of coordinates (3D physical space) +using coordinates_t = mito::geometry::coordinates_t<3, mito::geometry::CARTESIAN>; +// the metric space type +using metric_space_t = mito::geometry::euclidean_metric_space; +// the type of discretization node +using discretization_node_t = mito::discrete::discretization_node_t; +// the type of cell (triangle embedded in 3D) +using cell_t = mito::geometry::triangle_t<3>; +// the reference simplex +using reference_simplex_t = cell_t::reference_simplex_type; +// Gauss quadrature on triangles with degree of exactness 2 +using quadrature_rule_t = + mito::quadrature::quadrature_rule_t; + + +// 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 that the sum of the shape functions gradients is the zero vector + EXPECT_NEAR(0.0, mito::tensor::norm(sum), 3.0e-16); + }); + + // all done + return; +} + +TEST(Fem, IsoparametricEmbeddedTriangle) +{ + // the coordinate system + auto coord_system = mito::geometry::coordinate_system_t(); + + // an atlas under the coordinate system + auto atlas = mito::manifolds::atlas(coord_system); + + // build nodes for a unit right triangle rigidly rotated out of the xy-plane + // (the legs {t_1} = (0.6, 0.0, 0.8) and {t_2} = (0.0, 1.0, 0.0) are unit and orthogonal) + auto node_0 = mito::geometry::node(coord_system, { 0.0, 0.0, 0.0 }); + auto node_1 = mito::geometry::node(coord_system, { 0.6, 0.0, 0.8 }); + auto node_2 = mito::geometry::node(coord_system, { 0.0, 1.0, 0.0 }); + + // make a geometric simplex + auto triangle = mito::geometry::triangle(node_0, node_1, node_2); + + // the normal field to the triangle: {n} = {t_1} x {t_2} normalized, + // oriented such that w(n, t_1, t_2) > 0 + auto normal_field = + mito::functions::constant(mito::tensor::vector_t<3>{ -0.8, 0.0, 0.6 }); + + // strip the namespace from the placeholder for forms contractions + using mito::tensor::_; + + // the ambient metric volume form + constexpr auto w = metric_space_t::w; + + // the restriction of the metric volume form to the triangle + auto wS = mito::functions::function( + [w, normal_field](const coordinates_t & x) { return w(x)(normal_field(x), _, _); }); + + // make a manifold element from the triangle + auto element = + mito::manifolds::parametrized_element(triangle, atlas.parametrization(triangle), wS); + + { + // build the discretization nodes + auto discretization_node_0 = discretization_node_t(); + auto discretization_node_1 = discretization_node_t(); + auto discretization_node_2 = discretization_node_t(); + + // the degree of the finite element + constexpr int degree = 1; + // assemble the finite element type + using finite_element_t = mito::fem::finite_element_family; + + // first order isoparametric finite element on the embedded triangle + auto element_p1 = mito::fem::finite_element( + element, { discretization_node_0, discretization_node_1, discretization_node_2 }); + + // check that first order shape functions are a partition of unity + test_partition_of_unity(element_p1); + + // check that the gradients of first order shape functions sum to 0.0 + test_gradient_consistency(element_p1); + + // check the gradient values at the barycenter + { + auto xi = reference_simplex_t::parametric_coordinates_type{ 1.0 / 3.0, 1.0 / 3.0 }; + auto grad_0 = element_p1.gradient<0>()(xi); + auto grad_1 = element_p1.gradient<1>()(xi); + auto grad_2 = element_p1.gradient<2>()(xi); + + // the tangent vectors (the triangle legs) and the normal vector + constexpr auto t_1 = mito::tensor::vector_t<3>{ 0.6, 0.0, 0.8 }; + constexpr auto t_2 = mito::tensor::vector_t<3>{ 0.0, 1.0, 0.0 }; + constexpr auto n = mito::tensor::vector_t<3>{ -0.8, 0.0, 0.6 }; + + // the gradients lie in the plane of the triangle + EXPECT_NEAR(0.0, grad_0 * n, 1.0e-15); + EXPECT_NEAR(0.0, grad_1 * n, 1.0e-15); + EXPECT_NEAR(0.0, grad_2 * n, 1.0e-15); + + // the directional derivatives along the legs recover the parametric derivatives + EXPECT_NEAR(-1.0, grad_0 * t_1, 1.0e-15); + EXPECT_NEAR(-1.0, grad_0 * t_2, 1.0e-15); + EXPECT_NEAR(1.0, grad_1 * t_1, 1.0e-15); + EXPECT_NEAR(0.0, grad_1 * t_2, 1.0e-15); + EXPECT_NEAR(0.0, grad_2 * t_1, 1.0e-15); + EXPECT_NEAR(1.0, grad_2 * t_2, 1.0e-15); + } + } + + // all done + return; +} diff --git a/tests/mito.lib/geometry/euclidean_submanifold_metric_2D.cc b/tests/mito.lib/geometry/euclidean_submanifold_metric_2D.cc new file mode 100644 index 00000000..52d8257a --- /dev/null +++ b/tests/mito.lib/geometry/euclidean_submanifold_metric_2D.cc @@ -0,0 +1,180 @@ +// -*- c++ -*- +// +// Copyright (c) 2020-2026, the MiTo Authors, all rights reserved +// + +#include +#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 vector + auto director0 = x1 - x0; + + // compute signed length of segment + auto length = w(director0); + + // all done + return length; +} + + +TEST(Tensor, EuclideanSubmanifoldMetric2D) +{ + // the basis one-forms + constexpr auto dx = mito::tensor::one_form(e_x); + constexpr auto dy = mito::tensor::one_form(e_y); + + // the coordinate system + auto coord_system = mito::geometry::coordinate_system(); + + // pick two sets of coordinates + constexpr auto x_0 = mito::geometry::cartesian::coordinates({ 0.0, 0.0 }); + constexpr auto x_1 = mito::geometry::cartesian::coordinates({ 3.0, 4.0 }); + + // assemble the normal vector to the submanifold from tangent + constexpr auto tangent_vector = x_1 - x_0; + constexpr auto normal_vector = + mito::tensor::vector_t<2>{ tangent_vector[1], -tangent_vector[0] } + / mito::tensor::norm(tangent_vector); + + // the 2D metric volume element + constexpr auto w = mito::tensor::wedge(dx, dy); + + // the 1D restriction of the 2D metric volume element + constexpr auto wS = w(normal_vector, _); + + // build nodes of a line segment + auto node_0 = mito::geometry::node(coord_system, x_0); + auto node_1 = mito::geometry::node(coord_system, x_1); + + // check that the length of the line segment is correct + EXPECT_DOUBLE_EQ(5.0, length(wS, coord_system, node_0, node_1)); + + // check that the odd permutation of the vertices gives a negative length + EXPECT_DOUBLE_EQ(-5.0, length(wS, coord_system, node_1, node_0)); + + // the normal form to the submanifold + constexpr auto normal_form = mito::tensor::one_form(normal_vector); + + // rebuild the volume form as wV = wedge(normal_form, wS) + constexpr auto wV = mito::tensor::wedge(normal_form, wS); + + // check that wV coincides with w + EXPECT_DOUBLE_EQ(w(e_x, e_y), wV(e_x, e_y)); +} + + +// helper to test a single segment: computes the restricted volume form and verifies length +auto +test_segment( + mito::geometry::coordinate_system_t & 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 diff --git a/tests/mito.lib/quadrature/quadrature_segment_2D.cc b/tests/mito.lib/quadrature/quadrature_segment_2D.cc new file mode 100644 index 00000000..4ffa944c --- /dev/null +++ b/tests/mito.lib/quadrature/quadrature_segment_2D.cc @@ -0,0 +1,72 @@ +// -*- c++ -*- +// +// Copyright (c) 2020-2026, the MiTo Authors, all rights reserved +// + +#include +#include + + +// strip the namespace +using mito::tensor::_; +using mito::quadrature::GAUSS; + +// alias for a set of cartesian coordinates in 2D +using coordinates_t = mito::geometry::coordinates_t<2, mito::geometry::CARTESIAN>; +// the euclidean metric space type +using metric_space_t = mito::geometry::euclidean_metric_space; + + +TEST(Quadrature, Segment2D) +{ + // make a channel + journal::info_t channel("tests.quadrature_segment_2D"); + + // the coordinate system + auto coord_system = mito::geometry::coordinate_system(); + + // build nodes for a diagonal segment from (0,0) to (3,4) - length 5 + auto node_0 = mito::geometry::node(coord_system, { 0.0, 0.0 }); + auto node_1 = mito::geometry::node(coord_system, { 3.0, 4.0 }); + + // make a geometric simplex + auto segment = mito::geometry::segment(node_0, node_1); + + // a mesh with {segment} + auto mesh = mito::mesh::mesh>(); + mesh.insert(segment); + + // the normal field to the segment (perpendicular to the segment direction) + // segment direction is (3,4)/5 = (0.6, 0.8), so normal is (0.8, -0.6) + // (rotated 90° clockwise to get positive orientation with w(normal, tangent) > 0) + auto normal_field = + mito::functions::constant(mito::tensor::vector_t<2>{ 0.8, -0.6 }); + + // the ambient metric volume form + constexpr auto w = metric_space_t::w; + + // the restriction of the metric volume form to the segment + auto wS = mito::functions::function( + [w, normal_field](const coordinates_t & x) { return w(x)(normal_field(x), _); }); + + // This instantiates a quad rule on the cells (pairing cell type and degree of exactness) + auto manifold = mito::manifolds::manifold(mesh, coord_system, wS); + + // a scalar field + auto f_one = mito::functions::constant(1.0); + + // a zero degree-of-exactness integrator on the submanifold + auto integrator = mito::quadrature::integrator(manifold); + + // integrate {f_one} on the submanifold + auto result = integrator.integrate(f_one); + // the exact result + auto exact = 5.0; + // report + channel << "result: " << result << journal::endl; + channel << "exact: " << exact << journal::endl; + // expect a decent match with the exact solution + EXPECT_DOUBLE_EQ(result, exact); +} + +// end of file diff --git a/tests/mito.lib/quadrature/quadrature_triangle_3D.cc b/tests/mito.lib/quadrature/quadrature_triangle_3D.cc new file mode 100644 index 00000000..11956a4d --- /dev/null +++ b/tests/mito.lib/quadrature/quadrature_triangle_3D.cc @@ -0,0 +1,74 @@ +// -*- c++ -*- +// +// Copyright (c) 2020-2026, the MiTo Authors, all rights reserved +// + +#include +#include + + +// strip the namespace +using mito::tensor::_; +using mito::quadrature::GAUSS; + +// alias for a set of cartesian coordinates in 3D +using coordinates_t = mito::geometry::coordinates_t<3, mito::geometry::CARTESIAN>; +// the euclidean metric space type +using metric_space_t = mito::geometry::euclidean_metric_space; + + +TEST(Quadrature, Triangle3D) +{ + // make a channel + journal::info_t channel("tests.quadrature_triangle_3D"); + + // the coordinate system + auto coord_system = mito::geometry::coordinate_system(); + + // create nodes + // build nodes for a unit right triangle rigidly rotated out of the xy-plane + // (the legs {t_1} = (0.6, 0.0, 0.8) and {t_2} = (0.0, 1.0, 0.0) are unit and orthogonal) + auto node_0 = mito::geometry::node(coord_system, { 0.0, 0.0, 0.0 }); + auto node_1 = mito::geometry::node(coord_system, { 0.6, 0.0, 0.8 }); + auto node_2 = mito::geometry::node(coord_system, { 0.0, 1.0, 0.0 }); + + // make a geometric simplex + auto triangle = mito::geometry::triangle(node_0, node_1, node_2); + + // a mesh with {triangle} + auto mesh = mito::mesh::mesh>(); + mesh.insert(triangle); + + // the normal field to the triangle: {n} = {t_1} x {t_2} normalized, + // oriented such that w(n, t_1, t_2) > 0 + auto normal_field = + mito::functions::constant(mito::tensor::vector_t<3>{ -0.8, 0.0, 0.6 }); + + // the ambient metric volume form + constexpr auto w = metric_space_t::w; + + // the restriction of the metric volume form to the triangle + auto wS = mito::functions::function( + [w, normal_field](const coordinates_t & x) { return w(x)(normal_field(x), _, _); }); + + // This instantiates a quad rule on the cells (pairing cell type and degree of exactness) + auto manifold = mito::manifolds::manifold(mesh, coord_system, wS); + + // a scalar field + auto f_one = mito::functions::constant(1.0); + + // a zero degree-of-exactness integrator on the submanifold + auto integrator = mito::quadrature::integrator(manifold); + + // integrate {f_one} on the submanifold + auto result = integrator.integrate(f_one); + // the exact result + auto exact = 0.5; + // report + channel << "result: " << result << journal::endl; + channel << "exact: " << exact << journal::endl; + // expect a decent match with the exact solution + EXPECT_DOUBLE_EQ(result, exact); +} + +// end of file