Add support for FEM in 1D with linear isoparametric segments

.cmake/mito_tests_mito_lib.cmake

@@ -46,11 +46,15 @@ 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/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)
 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
 mito_test_driver(tests/mito.lib/io/summit_mesh_reader_2D.cc)

AUTHORS

@@ -2,3 +2,4 @@ This is a list of the MiTo copyright holders:
     Bianca Giovanardi <bianca.giovanardi@mitadel.com>
     Anwar Koshakji <anwar.koshakji@mitadel.com>
     Michael A.G. Aïvázis <michael.aivazis@mitadel.com>
+    Pawel J. Rok <pawel.rok024@gmail.com>

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<dim>;
+
+      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<coordinates_type>;
+        // 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

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"
 

lib/mito/fem/elements/public.h

@@ -17,6 +17,7 @@
 #include "api.h"
 
 // classes implementation
+#include "IsoparametricSegment.h"
 #include "IsoparametricTriangle.h"
 #include "Discretizer.h"
 

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<IsoparametricSegmentP1, discretization_t::CG> {
+        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

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 nodes
+        static constexpr int n_nodes = shape_functions_type::N;
+        // a collection of discretization nodes
+        using connectivity_type = std::array<discretization_node_type, n_nodes>;
+
+      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 <int a>
+        requires(a >= 0 && a < n_nodes)
+        constexpr auto shape() const
+        {
+            // return the shape functions
+            return shape_functions.shape<a>();
+        }
+
+        // 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 <int a>
+        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<a>()(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

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 <int a>
+        requires(a >= 0 && a < N)
+        constexpr auto shape() const
+        {
+            // return the a-th shape function
+            return std::get<a>(phi);
+        }
+
+        // get the a-th shape function's derivative as a function of parametric coordinates
+        template <int a>
+        requires(a >= 0 && a < N)
+        constexpr auto dshape() const
+        {
+            // return the a-th shape function's gradient
+            return std::get<a>(dphi);
+        }
+    };
+
+}    // namespace mito
+
+
+// end of file

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