Add a demo showing how submeshes can be used to represent boundary data

cpp/demo/mixed_poisson/main.cpp

@@ -8,6 +8,7 @@
 // * Apply boundary conditions to different fields in a mixed problem.
 // * Create integration domain data to execute finite element kernels.
 //   over subsets of the boundary.
+// * Use a submesh to represent boundary data
 //
 // The full implementation is in
 // {download}`demo_mixed_poisson/main.cpp`.
@@ -91,6 +92,7 @@
 //
 
 #include "mixed_poisson.h"
+#include <basix/cell.h>
 #include <basix/finite-element.h>
 #include <basix/mdspan.hpp>
 #include <cmath>
@@ -118,20 +120,23 @@ int main(int argc, char* argv[])
   PetscInitialize(&argc, &argv, nullptr, nullptr);
 
   {
+    mesh::CellType cell_type = mesh::CellType::triangle;
+
     // Create mesh
-    auto mesh = std::make_shared<mesh::Mesh<U>>(
-        mesh::create_rectangle<U>(MPI_COMM_WORLD, {{{0.0, 0.0}, {1.0, 1.0}}},
-                                  {32, 32}, mesh::CellType::triangle));
+    auto mesh = std::make_shared<mesh::Mesh<U>>(mesh::create_rectangle<U>(
+        MPI_COMM_WORLD, {{{0, 0}, {1, 1}}}, {32, 32}, cell_type));
 
     // Create Basix elements
-    auto RT = basix::create_element<U>(
-        basix::element::family::RT, basix::cell::type::triangle, 1,
-        basix::element::lagrange_variant::unset,
-        basix::element::dpc_variant::unset, false);
-    auto P0 = basix::create_element<U>(
-        basix::element::family::P, basix::cell::type::triangle, 0,
-        basix::element::lagrange_variant::unset,
-        basix::element::dpc_variant::unset, true);
+    basix::cell::type basix_cell_type
+        = dolfinx::mesh::cell_type_to_basix_type(cell_type);
+    auto RT
+        = basix::create_element<U>(basix::element::family::RT, basix_cell_type,
+                                   1, basix::element::lagrange_variant::unset,
+                                   basix::element::dpc_variant::unset, false);
+    auto P0
+        = basix::create_element<U>(basix::element::family::P, basix_cell_type,
+                                   0, basix::element::lagrange_variant::unset,
+                                   basix::element::dpc_variant::unset, true);
 
     // Create DOLFINx mixed element
     auto ME = std::make_shared<fem::FiniteElement<U>>(
@@ -149,7 +154,7 @@ int main(int argc, char* argv[])
     auto W0 = std::make_shared<fem::FunctionSpace<U>>(V0->collapse().first);
     auto W1 = std::make_shared<fem::FunctionSpace<U>>(V1->collapse().first);
 
-    // Create source function and interpolate sin(5x) + 1
+    // Create source function and interpolate $\sin(5x) + 1$
     auto f = std::make_shared<fem::Function<T>>(W1);
     f->interpolate(
         [](auto x) -> std::pair<std::vector<T>, std::vector<std::size_t>>
@@ -163,18 +168,7 @@ int main(int argc, char* argv[])
           return {f, {f.size()}};
         });
 
-    // Boundary condition value for u and interpolate 20y + 1
-    auto u0 = std::make_shared<fem::Function<T>>(W1);
-    u0->interpolate(
-        [](auto x) -> std::pair<std::vector<T>, std::vector<std::size_t>>
-        {
-          std::vector<T> f;
-          for (std::size_t p = 0; p < x.extent(1); ++p)
-            f.push_back(20 * x(1, p) + 1);
-          return {f, {f.size()}};
-        });
-
-    // Create boundary condition for \sigma and interpolate such that
+    // Create boundary condition for $\sigma$ and interpolate such that
     // flux = 10 (for top and bottom boundaries)
     auto g = std::make_shared<fem::Function<T>>(W0);
     g->interpolate(
@@ -200,7 +194,7 @@ int main(int argc, char* argv[])
         [](auto x)
         {
           using U = typename decltype(x)::value_type;
-          constexpr U eps = 1.0e-8;
+          constexpr U eps = 1e-8;
           std::vector<std::int8_t> marker(x.extent(1), false);
           for (std::size_t p = 0; p < x.extent(1); ++p)
           {
@@ -211,8 +205,53 @@ int main(int argc, char* argv[])
           return marker;
         });
 
-    // Compute facets with \sigma (flux) boundary condition facets, which is
-    // {all boundary facet} - {u0 boundary facets }
+    // We'd like to represent `u_0` using a function space defined only
+    // on the facets in `dfacets`. To do so, we begin by calling
+    // `create_submesh` to get a `submesh` of those facets. It also
+    // returns a map `submesh_to_mesh` whose `i`th entry is the facet in
+    // mesh corresponding to cell `i` in submesh.
+    int tdim = mesh->topology()->dim();
+    int fdim = tdim - 1;
+    std::shared_ptr<mesh::Mesh<U>> submesh;
+    std::vector<std::int32_t> submesh_to_mesh;
+    {
+      auto [_submesh, _submesh_to_mesh, v_map, g_map]
+          = mesh::create_submesh(*mesh, fdim, dfacets);
+      submesh = std::make_shared<mesh::Mesh<U>>(std::move(_submesh));
+      submesh_to_mesh = std::move(_submesh_to_mesh);
+    }
+
+    // Create an element for `u_0`
+    basix::cell::type submesh_cell_type
+        = dolfinx::mesh::cell_type_to_basix_type(
+            submesh->topology()->cell_type());
+
+    auto Qe = std::make_shared<fem::FiniteElement<U>>(
+        basix::create_element<U>(basix::element::family::P, submesh_cell_type,
+                                 1, basix::element::lagrange_variant::unset,
+                                 basix::element::dpc_variant::unset, true));
+
+    // Create a function space for `u_0` on the submesh
+    auto Q = std::make_shared<fem::FunctionSpace<U>>(
+        fem::create_functionspace<U>(submesh, Qe));
+
+    // Boundary condition value for u and interpolate $20 y + 1$
+    auto u0 = std::make_shared<fem::Function<T>>(Q);
+    u0->interpolate(
+        [](auto x) -> std::pair<std::vector<T>, std::vector<std::size_t>>
+        {
+          std::vector<T> f;
+          for (std::size_t p = 0; p < x.extent(1); ++p)
+            f.push_back(20 * x(1, p) + 1);
+          return {f, {f.size()}};
+        });
+
+    // Write u0 to file to visualise
+    io::VTKFile u0_file(MPI_COMM_WORLD, "u0.pvd", "w");
+    u0_file.write<T>({*u0}, 0);
+
+    // Compute facets with $\sigma$ (flux) boundary condition facets,
+    // which is `{all boundary facet} - {u0 boundary facets}`
     std::vector<std::int32_t> nfacets;
     std::ranges::set_difference(bfacets, dfacets, std::back_inserter(nfacets));
 
@@ -221,34 +260,51 @@ int main(int argc, char* argv[])
         = fem::locate_dofs_topological(
             *mesh->topology(), {*V0->dofmap(), *W0->dofmap()}, 1, nfacets);
 
-    // Create boundary condition for \sigma. \sigma \cdot n will be
-    // constrained to to be equal to the normal component of g. The
+    // Create boundary condition for $\sigma. $\sigma \cdot n$ will be
+    // constrained to to be equal to the normal component of $g$. The
     // boundary conditions are applied to degrees-of-freedom ndofs, and
-    // V0 is the subspace that is constrained.
+    // `V0` is the subspace that is constrained.
     fem::DirichletBC<T> bc(g, ndofs, V0);
 
-    // Create integration domain data for u0 boundary condition (applied
-    // on the ds(1) in the UFL file). First we get facet data
+    // Create integration domain data for `u0` boundary condition
+    // (applied on the `ds(1)` in the UFL file). First we get facet data
     // integration data for facets in dfacets.
     std::vector<std::int32_t> domains = fem::compute_integration_domains(
         fem::IntegralType::exterior_facet, *mesh->topology(), dfacets);
 
-    // Create data structure for the ds(1) integration domain in form
+    // Create data structure for the `ds(1)` integration domain in form
     // (see the UFL file). It is for en exterior facet integral (the key
     // in the map), and exterior facet domain marked as '1' in the UFL
-    // file, and 'domains' holds the necessary data to perform
+    // file, and `domains` holds the necessary data to perform
     // integration of selected facets.
     std::map<
         fem::IntegralType,
         std::vector<std::pair<std::int32_t, std::span<const std::int32_t>>>>
         subdomain_data{{fem::IntegralType::exterior_facet, {{1, domains}}}};
 
+    // Since we are doing a `ds(1)` integral on mesh and `u0` is defined
+    // on the `submesh`, we must provide an "entity map" relating cells
+    // in `submesh` to entities in `mesh`. This is simply the "inverse"
+    // of `submesh_to_mesh`:
+    auto facet_imap = mesh->topology()->index_map(fdim);
+    assert(facet_imap);
+    std::size_t num_facets = mesh->topology()->index_map(fdim)->size_local()
+                             + mesh->topology()->index_map(fdim)->num_ghosts();
+    std::vector<std::int32_t> mesh_to_submesh(num_facets, -1);
+    for (std::size_t i = 0; i < submesh_to_mesh.size(); ++i)
+      mesh_to_submesh[submesh_to_mesh[i]] = i;
+
+    // Create the entity map to pass to `create_form`
+    std::map<std::shared_ptr<const mesh::Mesh<U>>,
+             std::span<const std::int32_t>>
+        entity_maps = {{submesh, mesh_to_submesh}};
+
     // Define variational forms and attach he required data
     fem::Form<T> a = fem::create_form<T>(*form_mixed_poisson_a, {V, V}, {}, {},
                                          subdomain_data, {});
-    fem::Form<T> L
-        = fem::create_form<T>(*form_mixed_poisson_L, {V},
-                              {{"f", f}, {"u0", u0}}, {}, subdomain_data, {});
+    fem::Form<T> L = fem::create_form<T>(
+        *form_mixed_poisson_L, {V}, {{"f", f}, {"u0", u0}}, {}, subdomain_data,
+        entity_maps, V->mesh());
 
     // Create solution finite element Function
     auto u = std::make_shared<fem::Function<T>>(V);
@@ -272,7 +328,7 @@ int main(int argc, char* argv[])
     MatAssemblyBegin(A.mat(), MAT_FINAL_ASSEMBLY);
     MatAssemblyEnd(A.mat(), MAT_FINAL_ASSEMBLY);
 
-    // Assemble the linear form L into RHS vector
+    // Assemble the linear form `L` into RHS vector
     b.set(0);
     fem::assemble_vector(b.mutable_array(), L);
 

cpp/demo/mixed_poisson/mixed_poisson.py

@@ -3,6 +3,8 @@
 # {download}`demo_mixed_poisson/mixed_poisson.py`.  We begin by defining the
 # finite element:
 
+from basix import CellType
+from basix.cell import sub_entity_type
 from basix.ufl import element, mixed_element
 from ufl import (
     Coefficient,
@@ -16,12 +18,23 @@
     inner,
 )
 
-shape = "triangle"
-RT = element("RT", shape, 1)
-P = element("DP", shape, 0)
+# Cell type for the mesh
+msh_cell = CellType.triangle
+
+# Weakly enforced boundary data will be represented using a function space
+# defined over a submesh of the boundary. We get the submesh cell type from
+# then mesh cell type.
+submsh_cell = sub_entity_type(msh_cell, dim=1, index=0)
+
+# Define finite elements for the problem
+RT = element("RT", msh_cell, 1)
+P = element("DP", msh_cell, 0)
 ME = mixed_element([RT, P])
 
-msh = Mesh(element("Lagrange", shape, 1, shape=(2,)))
+# Define UFL mesh and submesh
+msh = Mesh(element("Lagrange", msh_cell, 1, shape=(2,)))
+submsh = Mesh(element("Lagrange", submsh_cell, 1, shape=(2,)))
+
 n = FacetNormal(msh)
 V = FunctionSpace(msh, ME)
 
@@ -30,8 +43,13 @@
 
 V0 = FunctionSpace(msh, P)
 f = Coefficient(V0)
-u0 = Coefficient(V0)
 
+# We represent boundary data using a first-degree discontinuous Lagrange space
+# defined over a submesh of the boundary
+Q = FunctionSpace(submsh, element("DP", submsh_cell, 1))
+u0 = Coefficient(Q)
+
+# Specify the weak form of the problem
 dx = Measure("dx", msh)
 ds = Measure("ds", msh)
 a = inner(sigma, tau) * dx + inner(u, div(tau)) * dx + inner(div(sigma), v) * dx