From 4e800137da2ffa2236865a3c6092c1cf30342864 Mon Sep 17 00:00:00 2001 From: Pawel024 Date: Tue, 13 Jan 2026 17:11:59 +0100 Subject: [PATCH 1/9] fem: 1d elements - base implementation --- lib/mito/fem/elements/IsoparametricSegment.h | 74 ++++++++++ lib/mito/fem/elements/elements_library.h | 1 + lib/mito/fem/elements/public.h | 1 + lib/mito/fem/elements/seg1/DiscretizerCG.h | 65 +++++++++ .../elements/seg1/IsoparametricSegmentP1.h | 136 ++++++++++++++++++ lib/mito/fem/elements/seg1/ShapeSegmentP1.h | 58 ++++++++ lib/mito/fem/elements/seg1/api.h | 21 +++ lib/mito/fem/elements/seg1/public.h | 19 +++ 8 files changed, 375 insertions(+) create mode 100644 lib/mito/fem/elements/IsoparametricSegment.h create mode 100644 lib/mito/fem/elements/seg1/DiscretizerCG.h create mode 100644 lib/mito/fem/elements/seg1/IsoparametricSegmentP1.h create mode 100644 lib/mito/fem/elements/seg1/ShapeSegmentP1.h create mode 100644 lib/mito/fem/elements/seg1/api.h create mode 100644 lib/mito/fem/elements/seg1/public.h diff --git a/lib/mito/fem/elements/IsoparametricSegment.h b/lib/mito/fem/elements/IsoparametricSegment.h new file mode 100644 index 000000000..3c867e35d --- /dev/null +++ b/lib/mito/fem/elements/IsoparametricSegment.h @@ -0,0 +1,74 @@ +// -*- c++ -*- +// +// Copyright (c) 2020-2024, the MiTo Authors, all rights reserved +// + +// code guard +#pragma once + + +// DESIGN NOTES +// Class {IsoparametricSegment} represents a first order simplex (segment) equipped with parametric coordinates. + + +namespace mito::fem { + + 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 + using coordinates_type = geometry::coordinates_t<1, geometry::CARTESIAN>; + // the coordinate system type + using coordinate_system_type = geometry::coordinate_system_t; + // the vector type + using vector_type = tensor::vector_t<1>; + + public: + // the default constructor + constexpr IsoparametricSegment( + 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 ~IsoparametricSegment() = default; + + // delete move constructor + constexpr IsoparametricSegment(IsoparametricSegment &&) noexcept = delete; + + // delete copy constructor + constexpr IsoparametricSegment(const IsoparametricSegment &) = delete; + + // delete assignment operator + constexpr IsoparametricSegment & operator=(const IsoparametricSegment &) = delete; + + // delete move assignment operator + constexpr IsoparametricSegment & operator=(IsoparametricSegment &&) 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 d52cfc68d..ac2cb50f4 100644 --- a/lib/mito/fem/elements/elements_library.h +++ b/lib/mito/fem/elements/elements_library.h @@ -7,6 +7,7 @@ #pragma once +#include "seg1/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 40f144579..1b33bd6be 100644 --- a/lib/mito/fem/elements/public.h +++ b/lib/mito/fem/elements/public.h @@ -17,6 +17,7 @@ #include "api.h" // classes implementation +#include "IsoparametricSegment.h" #include "IsoparametricTriangle.h" #include "Discretizer.h" diff --git a/lib/mito/fem/elements/seg1/DiscretizerCG.h b/lib/mito/fem/elements/seg1/DiscretizerCG.h new file mode 100644 index 000000000..0735af4c1 --- /dev/null +++ b/lib/mito/fem/elements/seg1/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 {IsoparametricSegmentP1} 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 IsoparametricSegmentP1::discretization_node_type; + + // the connectivity type + using connectivity_type = typename IsoparametricSegmentP1::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/IsoparametricSegmentP1.h b/lib/mito/fem/elements/seg1/IsoparametricSegmentP1.h new file mode 100644 index 000000000..934a2b65f --- /dev/null +++ b/lib/mito/fem/elements/seg1/IsoparametricSegmentP1.h @@ -0,0 +1,136 @@ +// -*- c++ -*- +// +// Copyright (c) 2020-2024, the MiTo Authors, all rights reserved +// + +// code guard +#pragma once + + +// 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. + + +namespace mito::fem { + + class IsoparametricSegmentP1 : public IsoparametricSegment { + + public: + // the degree of the finite element + static constexpr int degree = 1; + // the type of 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 discretization nodes + static constexpr int n_nodes = shape_functions_type::N; + // a collection of discretization discretization nodes + using connectivity_type = std::array; + + public: + // the default constructor + inline IsoparametricSegmentP1( + const cell_type & geometric_simplex, const coordinate_system_type & coord_system, + const connectivity_type & connectivity) : + IsoparametricSegment(geometric_simplex, coord_system), + _connectivity(connectivity) + {} + + // destructor + inline ~IsoparametricSegmentP1() = default; + + // delete move constructor + constexpr IsoparametricSegmentP1(IsoparametricSegmentP1 &&) noexcept = delete; + + // delete copy constructor + constexpr IsoparametricSegmentP1(const IsoparametricSegmentP1 &) = delete; + + // delete assignment operator + constexpr IsoparametricSegmentP1 & operator=(const IsoparametricSegmentP1 &) = delete; + + // delete move assignment operator + constexpr IsoparametricSegmentP1 & operator=(IsoparametricSegmentP1 &&) 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<1> { + // get the shape functions derivatives + constexpr auto dphi_0 = shape_functions.dshape<0>(); + constexpr auto dphi_1 = shape_functions.dshape<1>(); + + // compute the jacobian of the isoparametric mapping: dx/dxi + auto dx_dxi = _x0 * dphi_0(xi) + _x1 * dphi_1(xi); + // wrap the result in a 1x1 matrix + return tensor::matrix_t<1>{ dx_dxi }; + }); + + // 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<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; + } + + private: + // the discretization nodes of the simplex + const connectivity_type _connectivity; + }; + +} // namespace mito + + +// end of file diff --git a/lib/mito/fem/elements/seg1/ShapeSegmentP1.h b/lib/mito/fem/elements/seg1/ShapeSegmentP1.h new file mode 100644 index 000000000..c4c576a54 --- /dev/null +++ b/lib/mito/fem/elements/seg1/ShapeSegmentP1.h @@ -0,0 +1,58 @@ +// -*- c++ -*- +// +// Copyright (c) 2020-2024, the MiTo Authors, all rights reserved +// + +// code guard +#pragma once + + +namespace mito::fem { + + class ShapeSegmentP1 { + + public: + // the reference element + using reference_element_type = geometry::reference_segment_t; + // the number of shape functions + static constexpr int N = 2; + + private: + // linear shape functions on the reference segment in parametric coordinates + static constexpr auto xi_0 = reference_element_type::xi<0>; + + // linear shape functions on the segment + static constexpr auto phi_0 = 1.0 - xi_0; + static constexpr auto phi_1 = xi_0; + + // the shape functions + static constexpr auto phi = std::make_tuple(phi_0, phi_1); + + // the gradients of the shape functions + static constexpr auto dphi = std::make_tuple( + fields::gradient(phi_0), fields::gradient(phi_1)); + + public: + // get the a-th shape function as a function of parametric coordinates + template + requires(a >= 0 && a < N) + constexpr auto shape() const + { + // return the a-th shape function + return std::get(phi); + } + + // get the a-th shape function's derivative as a function of parametric coordinates + template + requires(a >= 0 && a < N) + constexpr auto dshape() const + { + // return the a-th shape function's gradient + return std::get(dphi); + } + }; + +} // namespace mito + + +// end of file diff --git a/lib/mito/fem/elements/seg1/api.h b/lib/mito/fem/elements/seg1/api.h new file mode 100644 index 000000000..3fcb2b63b --- /dev/null +++ b/lib/mito/fem/elements/seg1/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 in 1D + template <> + struct isoparametric_simplex<1, geometry::segment_t<1>> { + using type = IsoparametricSegmentP1; + }; + +} + + +// end of file diff --git a/lib/mito/fem/elements/seg1/public.h b/lib/mito/fem/elements/seg1/public.h new file mode 100644 index 000000000..d3f42ddbf --- /dev/null +++ b/lib/mito/fem/elements/seg1/public.h @@ -0,0 +1,19 @@ +// -*- c++ -*- +// +// Copyright (c) 2020-2024, the MiTo Authors, all rights reserved +// + +// code guard +#pragma once + + +// classes implementation +#include "ShapeSegmentP1.h" +#include "IsoparametricSegmentP1.h" +#include "DiscretizerCG.h" + +// published types and factories +#include "api.h" + + +// end of file From cc4bedf9b9f3268c6e3595fdd81120ae4fc0e4ac Mon Sep 17 00:00:00 2001 From: Pawel024 Date: Wed, 14 Jan 2026 10:49:36 +0100 Subject: [PATCH 2/9] tests/fem: test partition of unity and gradient consistency for the P1 segment element --- .cmake/mito_tests_mito_lib.cmake | 1 + tests/mito.lib/fem/isoparametric_segment.cc | 124 ++++++++++++++++++++ 2 files changed, 125 insertions(+) create mode 100644 tests/mito.lib/fem/isoparametric_segment.cc diff --git a/.cmake/mito_tests_mito_lib.cmake b/.cmake/mito_tests_mito_lib.cmake index 88ab480e6..75d4e1575 100644 --- a/.cmake/mito_tests_mito_lib.cmake +++ b/.cmake/mito_tests_mito_lib.cmake @@ -51,6 +51,7 @@ mito_test_driver(tests/mito.lib/fem/shape_functions_triangle_construction.cc) mito_test_driver(tests/mito.lib/fem/shape_functions_triangle_p1.cc) mito_test_driver(tests/mito.lib/fem/shape_functions_triangle_p2.cc) mito_test_driver(tests/mito.lib/fem/isoparametric_triangle.cc) +mito_test_driver(tests/mito.lib/fem/isoparametric_segment.cc) # io mito_test_driver(tests/mito.lib/io/summit_mesh_reader_2D.cc) diff --git a/tests/mito.lib/fem/isoparametric_segment.cc b/tests/mito.lib/fem/isoparametric_segment.cc new file mode 100644 index 000000000..7a1e53094 --- /dev/null +++ b/tests/mito.lib/fem/isoparametric_segment.cc @@ -0,0 +1,124 @@ +// -*- c++ -*- +// +// Copyright (c) 2020-2024, the MiTo Authors, all rights reserved +// + +#include +#include + + +// the type of coordinates +using coordinates_t = mito::geometry::coordinates_t<1, 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<1>; +// 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 + 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 + EXPECT_NEAR(1.0, sum, 1.0e-15); + }); + + // 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, 3.0e-16); + }); + + // all done + return; +} + + +TEST(Fem, IsoparametricSegment) +{ + // the coordinate system + auto coord_system = coord_system_t(); + + // build nodes + auto node_0 = mito::geometry::node(coord_system, { 0.0 }); + auto node_1 = mito::geometry::node(coord_system, { 1.0 }); + + // make a geometric simplex + auto geometric_simplex = mito::geometry::segment<1>({ 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 5d7944381f765218ed89425c8ac7053fc0fd6a75 Mon Sep 17 00:00:00 2001 From: Pawel024 Date: Wed, 14 Jan 2026 10:55:50 +0100 Subject: [PATCH 3/9] tests/fem: test shape function values at nodes --- .cmake/mito_tests_mito_lib.cmake | 1 + .../fem/shape_functions_segment_p1.cc | 42 +++++++++++++++++++ 2 files changed, 43 insertions(+) create mode 100644 tests/mito.lib/fem/shape_functions_segment_p1.cc diff --git a/.cmake/mito_tests_mito_lib.cmake b/.cmake/mito_tests_mito_lib.cmake index 75d4e1575..3a0897b31 100644 --- a/.cmake/mito_tests_mito_lib.cmake +++ b/.cmake/mito_tests_mito_lib.cmake @@ -51,6 +51,7 @@ mito_test_driver(tests/mito.lib/fem/shape_functions_triangle_construction.cc) mito_test_driver(tests/mito.lib/fem/shape_functions_triangle_p1.cc) mito_test_driver(tests/mito.lib/fem/shape_functions_triangle_p2.cc) mito_test_driver(tests/mito.lib/fem/isoparametric_triangle.cc) +mito_test_driver(tests/mito.lib/fem/shape_functions_segment_p1.cc) mito_test_driver(tests/mito.lib/fem/isoparametric_segment.cc) # io diff --git a/tests/mito.lib/fem/shape_functions_segment_p1.cc b/tests/mito.lib/fem/shape_functions_segment_p1.cc new file mode 100644 index 000000000..f3388574a --- /dev/null +++ b/tests/mito.lib/fem/shape_functions_segment_p1.cc @@ -0,0 +1,42 @@ +// -*- c++ -*- +// +// Copyright (c) 2020-2024, the MiTo Authors, all rights reserved +// + +#include +#include + + +// first order shape functions type +using shape_t = mito::fem::ShapeSegmentP1; +// the parametric coordinates type +using parametric_coordinates_type = shape_t::reference_element_type::parametric_coordinates_type; + + +TEST(Fem, ShapeSegmentP1) +{ + // first order shape functions + constexpr auto element = shape_t(); + + // node 0 in parametric coordinates + constexpr auto n0 = parametric_coordinates_type{ 0.0 }; + // node 1 in parametric coordinates + constexpr auto n1 = parametric_coordinates_type{ 1.0 }; + + // the shape function associated with local node {0} + constexpr auto phi_0 = element.shape<0>(); + // check that the shape function at node 0 is 1.0 + static_assert(1.0 == phi_0(n0)); + // check that the shape function at node 1 is 0.0 + static_assert(0.0 == phi_0(n1)); + + // the shape functions at node 1 + constexpr auto phi_1 = element.shape<1>(); + // check that the shape function at node 0 is 0.0 + static_assert(0.0 == phi_1(n0)); + // check that the shape function at node 1 is 1.0 + static_assert(1.0 == phi_1(n1)); + + // all done + return; +} From 14a75d34eebc0e432afe4d313c82b31c7368adc0 Mon Sep 17 00:00:00 2001 From: Pawel024 Date: Wed, 14 Jan 2026 11:13:00 +0100 Subject: [PATCH 4/9] tests/fem: stricter static assert for the partition of unity test on a segment --- tests/mito.lib/fem/isoparametric_segment.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/mito.lib/fem/isoparametric_segment.cc b/tests/mito.lib/fem/isoparametric_segment.cc index 7a1e53094..d5f5a6ba2 100644 --- a/tests/mito.lib/fem/isoparametric_segment.cc +++ b/tests/mito.lib/fem/isoparametric_segment.cc @@ -42,14 +42,14 @@ test_partition_of_unity(const auto & element) constexpr auto xi = quadrature_rule.point(q); // compute the sum of the shape functions at {xi} for all nodes - auto sum = + 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 - EXPECT_NEAR(1.0, sum, 1.0e-15); + static_assert(1.0 == sum); }); // all done From 4dcfd3490d68b8f410dae8b93bfdb2fd4602b7ea Mon Sep 17 00:00:00 2001 From: Pawel024 Date: Wed, 14 Jan 2026 11:21:41 +0100 Subject: [PATCH 5/9] tests/fem: test stiffness matrix for the P1 segment element --- .cmake/mito_tests_mito_lib.cmake | 1 + tests/mito.lib/fem/block_grad_grad_segment.cc | 72 +++++++++++++++++++ 2 files changed, 73 insertions(+) create mode 100644 tests/mito.lib/fem/block_grad_grad_segment.cc diff --git a/.cmake/mito_tests_mito_lib.cmake b/.cmake/mito_tests_mito_lib.cmake index 3a0897b31..6f3d1c806 100644 --- a/.cmake/mito_tests_mito_lib.cmake +++ b/.cmake/mito_tests_mito_lib.cmake @@ -46,6 +46,7 @@ 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_mass.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) diff --git a/tests/mito.lib/fem/block_grad_grad_segment.cc b/tests/mito.lib/fem/block_grad_grad_segment.cc new file mode 100644 index 000000000..9e25f1c2b --- /dev/null +++ b/tests/mito.lib/fem/block_grad_grad_segment.cc @@ -0,0 +1,72 @@ +// -*- c++ -*- +// +// Copyright (c) 2020-2024, the MiTo Authors, all rights reserved +// + +#include +#include + + +// the type of coordinates +using coordinates_t = mito::geometry::coordinates_t<1, 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<1>; +// 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, BlockGradGradSegment) +{ + // the coordinate system + auto coord_system = coord_system_t(); + + // build nodes + auto node_0 = mito::geometry::node(coord_system, { 0.0 }); + auto node_1 = mito::geometry::node(coord_system, { 1.0 }); + + // make a geometric simplex + auto geometric_simplex = mito::geometry::segment<1>({ 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 }); + + // a grad-grad matrix block + auto grad_grad_block = + mito::fem::blocks::grad_grad_block(); + + // the analytical elementary stiffness matrix + 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 From 4d8f75043c9f8d1a26be7ac82ba7d7d6cf2440eb Mon Sep 17 00:00:00 2001 From: Pawel024 Date: Wed, 14 Jan 2026 11:58:00 +0100 Subject: [PATCH 6/9] tests/fem: test mass matrix for the P1 segment element --- .cmake/mito_tests_mito_lib.cmake | 1 + tests/mito.lib/fem/block_mass_segment.cc | 72 ++++++++++++++++++++++++ 2 files changed, 73 insertions(+) create mode 100644 tests/mito.lib/fem/block_mass_segment.cc diff --git a/.cmake/mito_tests_mito_lib.cmake b/.cmake/mito_tests_mito_lib.cmake index 6f3d1c806..b771cb58d 100644 --- a/.cmake/mito_tests_mito_lib.cmake +++ b/.cmake/mito_tests_mito_lib.cmake @@ -48,6 +48,7 @@ mito_test_driver(tests/mito.lib/discrete/mesh_field.cc) 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_mass.cc) +mito_test_driver(tests/mito.lib/fem/block_mass_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_mass_segment.cc b/tests/mito.lib/fem/block_mass_segment.cc new file mode 100644 index 000000000..97224cb00 --- /dev/null +++ b/tests/mito.lib/fem/block_mass_segment.cc @@ -0,0 +1,72 @@ +// -*- c++ -*- +// +// Copyright (c) 2020-2024, the MiTo Authors, all rights reserved +// + +#include +#include + + +// the type of coordinates +using coordinates_t = mito::geometry::coordinates_t<1, 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<1>; +// 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, BlockMassSegment) +{ + // the coordinate system + auto coord_system = coord_system_t(); + + // build nodes + auto node_0 = mito::geometry::node(coord_system, { 0.0 }); + auto node_1 = mito::geometry::node(coord_system, { 1.0 }); + + // make a geometric simplex + auto geometric_simplex = mito::geometry::segment<1>({ 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 }); + + // a mass matrix block + auto mass_block = mito::fem::blocks::mass_block(); + + // the analytical elementary mass matrix + 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 9858d9e7cac560a122af26e19f767eab2bb715ca Mon Sep 17 00:00:00 2001 From: Pawel024 Date: Wed, 14 Jan 2026 12:51:42 +0100 Subject: [PATCH 7/9] tests/fem: add sympy verification for the stiffness and mass matrices of a p1 segment --- tests/mito.lib/fem/verification_segment_p1.py | 38 +++++++++++++++++++ 1 file changed, 38 insertions(+) create mode 100644 tests/mito.lib/fem/verification_segment_p1.py diff --git a/tests/mito.lib/fem/verification_segment_p1.py b/tests/mito.lib/fem/verification_segment_p1.py new file mode 100644 index 000000000..de134cd19 --- /dev/null +++ b/tests/mito.lib/fem/verification_segment_p1.py @@ -0,0 +1,38 @@ +import sympy +from sympy.vector import CoordSys3D, gradient + +# define the coordinate system +C = CoordSys3D('C') + +# define the coordinate +x = C.x + +# define the shape functions +phi0 = 1 - x +phi1 = x + +# report +print("Computing the mass matrix for the 2 shape functions of a segment element:") + +# compute the mass matrix via double integration +phis = [phi0, phi1] +print("Mass matrix:") +for i, phi_i in enumerate(phis): + for j, phi_j in enumerate(phis): + integral = sympy.integrate(phi_i * phi_j, (x, 0, 1)) + print(f"M[{i},{j}] = {integral}") + +# report +print("Computing the stiffness matrix for the 2 shape functions of a segment element:") + +# compute the shape functions gradients +dphi0 = gradient(phi0) +dphi1 = gradient(phi1) + +# compute the stiffness matrix via double integration +dphis = [dphi0, dphi1] +print("Stiffness matrix:") +for i, dphi_i in enumerate(dphis): + for j, dphi_j in enumerate(dphis): + integral = sympy.integrate(dphi_i & dphi_j, (x, 0, 1)) + print(f"K[{i},{j}] = {integral}") From c7d71d75923be8b5e2ad1c39bf4f52d69380653d Mon Sep 17 00:00:00 2001 From: Pawel024 Date: Wed, 14 Jan 2026 16:35:39 +0100 Subject: [PATCH 8/9] fem: fix comment typo --- lib/mito/fem/elements/seg1/IsoparametricSegmentP1.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/lib/mito/fem/elements/seg1/IsoparametricSegmentP1.h b/lib/mito/fem/elements/seg1/IsoparametricSegmentP1.h index 934a2b65f..c4e8dde6b 100644 --- a/lib/mito/fem/elements/seg1/IsoparametricSegmentP1.h +++ b/lib/mito/fem/elements/seg1/IsoparametricSegmentP1.h @@ -28,9 +28,9 @@ namespace mito::fem { typename canonical_element_type::parametric_coordinates_type; // the linear shape functions static constexpr auto shape_functions = shape_functions_type(); - // the number of discretization discretization nodes + // the number of discretization nodes static constexpr int n_nodes = shape_functions_type::N; - // a collection of discretization discretization nodes + // a collection of discretization nodes using connectivity_type = std::array; public: From 2f2a88706269b420c5336bb18aa96d8d329847d7 Mon Sep 17 00:00:00 2001 From: Pawel024 Date: Wed, 14 Jan 2026 20:07:17 +0100 Subject: [PATCH 9/9] add my name to AUTHORS --- AUTHORS | 1 + 1 file changed, 1 insertion(+) diff --git a/AUTHORS b/AUTHORS index e1694c831..76a562312 100644 --- a/AUTHORS +++ b/AUTHORS @@ -2,3 +2,4 @@ This is a list of the MiTo copyright holders: Bianca Giovanardi Anwar Koshakji Michael A.G. Aïvázis + Pawel J. Rok