Embedded Interfaces

.cmake/mito_tests_mito_lib.cmake

@@ -100,6 +100,7 @@ endif()
 if(WITH_PETSC)
     mito_test_driver(tests/mito.lib/matrix_solvers/petsc_initialize_finalize.cc)
     mito_test_driver(tests/mito.lib/matrix_solvers/petsc_ksp.cc)
+    mito_test_driver(tests/mito.lib/fem/hybrid_dimensional_diffusion_continuous.cc)
 endif()
 
 # tensor

lib/mito/constraints/Dirichlet.h

@@ -18,6 +18,8 @@ namespace mito::constraints {
       public:
         using domain_type = meshT;
         using function_type = fieldT;
+        // the type of the values prescribed by the constraint
+        using value_type = typename function_type::output_type;
         using node_type = typename domain_type::node_type;
         using nodes_type = std::set<node_type>;
 

lib/mito/fem/FunctionSpace.h

@@ -38,8 +38,11 @@ namespace mito::fem {
 
         // the discretization node type
         using discretization_node_type = typename finite_element_type::discretization_node_type;
-        // the constrained nodes type
-        using constrained_nodes_type = std::set<discretization_node_type>;
+        // the type of the values prescribed by the constraints
+        using constraint_value_type = typename constraints_type::value_type;
+        // the constrained values type (map from constrained node to prescribed value)
+        using constrained_values_type =
+            std::map<discretization_node_type, constraint_value_type>;
         // the type of a map between the mesh nodes and discretization nodes
         using map_type = std::unordered_map<
             mesh_node_type, discretization_node_type, utilities::hash_function<mesh_node_type>>;
@@ -68,7 +71,24 @@ namespace mito::fem {
             // TODO: merge the discretization type with the finite element type
             // discretize the manifold subject to the constraints
             discretize<finite_element_type, discretizationT>(
-                manifold, constraints, _connectivity_table, _node_map, _constrained_nodes);
+                manifold, constraints, _connectivity_table, _node_map, _constrained_values);
+        }
+
+        // the constructor with a pre-populated node map (for coupled problems that share
+        // discretization nodes with another function space)
+        template <discretization_t discretizationT = discretization_t::CG>
+        constexpr FunctionSpace(
+            const manifold_type & manifold, const constraints_type & constraints,
+            const map_type & shared_node_map) :
+            _manifold(manifold),
+            _constraints(constraints),
+            _connectivity_table(),
+            _node_map(shared_node_map)
+        {
+            // discretize the manifold subject to the constraints; mesh nodes already present in
+            // the node map reuse their discretization nodes
+            discretize<finite_element_type, discretizationT>(
+                manifold, constraints, _connectivity_table, _node_map, _constrained_values);
         }
 
         // destructor
@@ -113,10 +133,10 @@ namespace mito::fem {
             return _constraints;
         }
 
-        // get the constrained nodes
-        constexpr auto constrained_nodes() const noexcept -> const constrained_nodes_type &
+        // get the constrained nodes and their prescribed values
+        constexpr auto constrained_values() const noexcept -> const constrained_values_type &
         {
-            return _constrained_nodes;
+            return _constrained_values;
         }
 
         // accessor for the node map
@@ -156,8 +176,8 @@ namespace mito::fem {
         // the constraints
         const constraints_type & _constraints;
 
-        // the constrained nodes
-        constrained_nodes_type _constrained_nodes;
+        // the constrained nodes and their prescribed values
+        constrained_values_type _constrained_values;
 
         // the connectivity table of the finite elements
         connectivity_table_type _connectivity_table;

lib/mito/fem/api.h

@@ -41,14 +41,19 @@ namespace mito::fem {
     constexpr auto weakform(const lhsBlockT & lhs_block, const rhsBlockT & rhs_block);
 
     // discrete system alias
-    template <class functionSpaceT, class weakformT, class linearSystemT>
-    using discrete_system_t = DiscreteSystem<functionSpaceT, weakformT, linearSystemT>;
+    template <class linearSystemT, contribution_c... contributionTs>
+    using discrete_system_t = DiscreteSystem<linearSystemT, contributionTs...>;
 
-    // discrete system factory
-    template <class linearSystemT, class functionSpaceT, class weakformT>
+    // discrete system factory (one contribution per function space)
+    template <class linearSystemT, contribution_c... contributionTs>
     constexpr auto discrete_system(
-        const functionSpaceT & function_space, const weakformT & weakform,
-        const std::string & label);
+        const std::string & label, const contributionTs &... contributions);
+
+    // discrete system factory (single function space)
+    template <class linearSystemT, function_space_c functionSpaceT, class weakformT>
+    constexpr auto discrete_system(
+        const std::string & label, const functionSpaceT & function_space,
+        const weakformT & weakform);
 }
 
 

lib/mito/fem/elements/Discretizer.h

@@ -13,19 +13,19 @@ namespace mito::fem {
     struct Discretizer {
         template <
             typename manifoldT, typename constraintsT, typename connectivity_table_type,
-            typename map_type, typename constrained_nodes_type>
+            typename map_type, typename constrained_values_type>
         static void apply(
             const manifoldT &, const constraintsT &, connectivity_table_type &, map_type &,
-            constrained_nodes_type &);
+            constrained_values_type &);
     };
 
     template <class elementT, discretization_t discretizationT>
     auto discretize(
         const auto & manifold, const auto & constraints, auto & connectivity, auto & node_map,
-        auto & constrained_nodes)
+        auto & constrained_values)
     {
         Discretizer<elementT, discretizationT>::apply(
-            manifold, constraints, connectivity, node_map, constrained_nodes);
+            manifold, constraints, connectivity, node_map, constrained_values);
     }
 
 }

lib/mito/fem/elements/seg1/DiscretizerCG.h

@@ -14,11 +14,11 @@ namespace mito::fem {
     struct Discretizer<finite_element_family<geometry::segment_t<D>, 1>, discretization_t::CG> {
         template <
             typename manifoldT, typename constraintsT, typename connectivity_table_type,
-            typename map_type, typename constrained_nodes_type>
+            typename map_type, typename constrained_values_type>
         static void apply(
             const manifoldT & manifold, const constraintsT & constraints,
             connectivity_table_type & connectivity, map_type & node_map,
-            constrained_nodes_type & constrained_nodes)
+            constrained_values_type & constrained_values)
         {
             // the finite element type
             using finite_element_type = finite_element_family<geometry::segment_t<D>, 1>;
@@ -49,14 +49,20 @@ namespace mito::fem {
                 connectivity.emplace(cell.simplex().id(), connectivity_type{ node_0, node_1 });
             }
 
-            // populate the constrained nodes
+            // get the coordinate system and the constraint function
+            const auto & coord_system = manifold.coordinate_system();
+            const auto & function = constraints.function();
+
+            // populate the constrained nodes with the values of the constraint function
             // 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);
+                // add the node to the constrained nodes with the value of the constraint function
+                // at the node coordinates
+                constrained_values.insert(
+                    { it->second, function(coord_system.coordinates(node->point())) });
             }
 
             // all done

lib/mito/fem/elements/tri1/DiscretizerCG.h

@@ -14,11 +14,11 @@ namespace mito::fem {
     struct Discretizer<finite_element_family<geometry::triangle_t<D>, 1>, discretization_t::CG> {
         template <
             typename manifoldT, typename constraintsT, typename connectivity_table_type,
-            typename map_type, typename constrained_nodes_type>
+            typename map_type, typename constrained_values_type>
         static void apply(
             const manifoldT & manifold, const constraintsT & constraints,
             connectivity_table_type & connectivity, map_type & node_map,
-            constrained_nodes_type & constrained_nodes)
+            constrained_values_type & constrained_values)
         {
             // the finite element type
             using finite_element_type = finite_element_family<geometry::triangle_t<D>, 1>;
@@ -52,13 +52,19 @@ namespace mito::fem {
                     cell.simplex().id(), connectivity_type{ node_0, node_1, node_2 });
             }
 
-            // populate the constrained nodes
+            // get the coordinate system and the constraint function
+            const auto & coord_system = manifold.coordinate_system();
+            const auto & function = constraints.function();
+
+            // populate the constrained nodes with the values of the constraint function
             for (const auto & cell : constraints.domain().cells()) {
                 for (const auto & node : cell.nodes()) {
                     // 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);
+                    // add the node to the constrained nodes with the value of the constraint
+                    // function at the node coordinates
+                    constrained_values.insert(
+                        { it->second, function(coord_system.coordinates(node->point())) });
                 }
             }
 

lib/mito/fem/elements/tri2/DiscretizerCG.h

@@ -14,11 +14,11 @@ namespace mito::fem {
     struct Discretizer<finite_element_family<geometry::triangle_t<D>, 2>, discretization_t::CG> {
         template <
             typename manifoldT, typename constraintsT, typename connectivity_table_type,
-            typename map_type, typename constrained_nodes_type>
+            typename map_type, typename constrained_values_type>
         static void apply(
             const manifoldT & manifold, const constraintsT & constraints,
             connectivity_table_type & connectivity, map_type & node_map,
-            constrained_nodes_type & constrained_nodes)
+            constrained_values_type & constrained_values)
         {
             // the finite element type
             using finite_element_type = finite_element_family<geometry::triangle_t<D>, 2>;
@@ -82,22 +82,33 @@ namespace mito::fem {
                     connectivity_type{ node_0, node_1, node_2, node_3, node_4, node_5 });
             }
 
-            // populate the constrained nodes
+            // get the coordinate system and the constraint function
+            const auto & coord_system = manifold.coordinate_system();
+            const auto & function = constraints.function();
+
+            // populate the constrained nodes with the values of the constraint function
             for (const auto & cell : constraints.domain().cells()) {
                 for (const auto & node : cell.nodes()) {
                     // 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);
+                    // add the node to the constrained nodes with the value of the constraint
+                    // function at the node coordinates
+                    constrained_values.insert(
+                        { it->second, function(coord_system.coordinates(node->point())) });
                 }
                 auto node_0 = node_map.at(cell.nodes()[0]);
                 auto node_1 = node_map.at(cell.nodes()[1]);
                 auto ordered_nodes = (node_0.id() < node_1.id()) ?
                                          std::array{ node_0.id(), node_1.id() } :
                                          std::array{ node_1.id(), node_0.id() };
                 auto node = mid_nodes_map.at(ordered_nodes);
-                // add the node to the constrained nodes
-                constrained_nodes.insert(node);
+                // the coordinates of the two vertices of the boundary cell
+                auto coord_0 = coord_system.coordinates(cell.nodes()[0]->point());
+                auto coord_1 = coord_system.coordinates(cell.nodes()[1]->point());
+                // add the mid node to the constrained nodes with the value of the constraint
+                // function at the midpoint of the boundary cell
+                constrained_values.insert(
+                    { node, function(coord_0 + 0.5 * (coord_1 - coord_0)) });
             }
 
             // all done

lib/mito/fem/factories.h

@@ -25,21 +25,47 @@ namespace mito::fem {
         return function_space_t<elementT, manifoldT, constraintsT>(manifold, constraints);
     }
 
+    // function space factory with a pre-populated node map (for coupled problems that share
+    // discretization nodes with another function space)
+    template <
+        class elementT, manifolds::manifold_c manifoldT, constraints::constraint_c constraintsT>
+    // require compatibility between the manifold cell and the finite element cell
+    requires(
+        std::is_same_v<typename manifoldT::mesh_type::cell_type, typename elementT::mesh_cell_type>)
+    constexpr auto function_space(
+        const manifoldT & manifold, const constraintsT & constraints,
+        const typename function_space_t<elementT, manifoldT, constraintsT>::map_type &
+            shared_node_map)
+    {
+        // build a function space on the manifold, reusing the discretization nodes of the mesh
+        // nodes already present in {shared_node_map}
+        return function_space_t<elementT, manifoldT, constraintsT>(
+            manifold, constraints, shared_node_map);
+    }
+
     // weakform factory
     template <class lhsBlockT, class rhsBlockT>
     constexpr auto weakform(const lhsBlockT & lhs_block, const rhsBlockT & rhs_block)
     {
         return weakform_t<lhsBlockT, rhsBlockT>(lhs_block, rhs_block);
     }
 
-    // discrete system factory
-    template <class linearSystemT, class functionSpaceT, class weakformT>
+    // discrete system factory (one contribution per function space)
+    template <class linearSystemT, contribution_c... contributionTs>
+    constexpr auto discrete_system(
+        const std::string & label, const contributionTs &... contributions)
+    {
+        return discrete_system_t<linearSystemT, contributionTs...>(label, contributions...);
+    }
+
+    // discrete system factory (single function space)
+    template <class linearSystemT, function_space_c functionSpaceT, class weakformT>
     constexpr auto discrete_system(
         const std::string & label, const functionSpaceT & function_space,
         const weakformT & weakform)
     {
-        return discrete_system_t<functionSpaceT, weakformT, linearSystemT>(
-            label, function_space, weakform);
+        return discrete_system_t<linearSystemT, Contribution<functionSpaceT, weakformT>>(
+            label, Contribution<functionSpaceT, weakformT>{ function_space, weakform });
     }
 }
 

lib/mito/fem/forward.h

@@ -35,8 +35,22 @@ namespace mito::fem {
     requires compatible_assembly_blocks_c<lhsBlockT, rhsBlockT>
     class Weakform;
 
+    // a contribution to a discrete system: a weakform to be assembled on a function space,
+    // scaled by a constant coefficient
+    template <function_space_c functionSpaceT, class weakformT>
+    struct Contribution;
+
+    // concept of a contribution
+    template <class C>
+    concept contribution_c = requires(C c) {
+        // require that C only binds to {Contribution} specializations
+        []<function_space_c functionSpaceT, class weakformT>(
+            const Contribution<functionSpaceT, weakformT> &) {
+        }(c);
+    };
+
     // class discrete system
-    template <function_space_c functionSpaceT, class weakformT, class linearSystemT>
+    template <class linearSystemT, contribution_c... contributionTs>
     class DiscreteSystem;
 
     // class domain field

lib/mito/manifolds/Manifold.h

@@ -73,6 +73,12 @@ namespace mito::manifolds {
         // accessor for the mesh
         constexpr auto mesh() const noexcept -> const mesh_type & { return _mesh; }
 
+        // accessor for the coordinate system
+        constexpr auto coordinate_system() const noexcept -> const coordinate_system_type &
+        {
+            return _atlas.coordinate_system();
+        }
+
         // return an iterable view of the manifold elements
         constexpr auto elements() const noexcept { return manifold_elements_view_type{ *this }; }