FEniCS
diff --git a/‎cpp/demo/codim_0_assembly/main.cpp‎
Lines changed: 25 additions & 25 deletions b/‎cpp/demo/codim_0_assembly/main.cpp‎
Lines changed: 25 additions & 25 deletions
diff --git a/‎cpp/demo/mixed_poisson/main.cpp‎
Lines changed: 24 additions & 33 deletions b/‎cpp/demo/mixed_poisson/main.cpp‎
Lines changed: 24 additions & 33 deletions
@@ -14,6 +14,7 @@
 #include <dolfinx/fem/petsc.h>
 #include <dolfinx/la/MatrixCSR.h>
 #include <dolfinx/la/SparsityPattern.h>
+#include <dolfinx/mesh/EntityMap.h>
 #include <map>
 #include <memory>
 #include <ranges>
@@ -76,37 +77,25 @@ int main(int argc, char* argv[])
     mesh::MeshTags<std::int32_t> cell_marker(mesh->topology(), tdim, cells,
                                              values);
 
-    std::shared_ptr<mesh::Mesh<U>> submesh;
-    std::vector<std::int32_t> submesh_to_mesh;
+    // We create a submesh consisting of only cells with a given tag by
+    // calling `create_submesh`. This function also returns an
+    // `EntityMap` object, which relates entities in the submesh to
+    // entities in the original mesh. We will need this to assemble our
+    // mixed-domain form.
+    auto submesh_data = [](auto& mesh, int tdim, auto&& subcells)
     {
-      auto [_submesh, _submesh_to_mesh, v_map, g_map]
-          = mesh::create_submesh(*mesh, tdim, cell_marker.find(2));
-      submesh = std::make_shared<mesh::Mesh<U>>(std::move(_submesh));
-      submesh_to_mesh = std::move(_submesh_to_mesh);
-    }
+      auto [submesh, emap, v_map, g_map]
+          = mesh::create_submesh(mesh, tdim, subcells);
+      return std::pair(std::make_shared<mesh::Mesh<U>>(std::move(submesh)),
+                       std::move(emap));
+    };
+    auto [submesh, entity_map] = submesh_data(*mesh, tdim, cell_marker.find(2));
 
     // We create the function space used for the trial space
     auto W
         = std::make_shared<fem::FunctionSpace<U>>(fem::create_functionspace<U>(
             submesh, std::make_shared<fem::FiniteElement<U>>(element)));
 
-    // A mixed-domain form has functions defined over different meshes.
-    // The mesh associated with the measure (dx, ds, etc.) is called the
-    // integration domain. To assemble mixed-domain forms, maps must be
-    // provided taking entities in the integration domain to entities on
-    // each mesh in the form. Since one of our forms has a measure
-    // defined over `mesh` and involves a function defined over
-    // `submesh`, we must provide a map from entities in `mesh` to
-    // entities in `submesh`. This is simply the "inverse" of
-    // `submesh_to_mesh`.
-    std::vector<std::int32_t> mesh_to_submesh(num_cells_local, -1);
-    for (std::size_t i = 0; i < submesh_to_mesh.size(); ++i)
-      mesh_to_submesh[submesh_to_mesh[i]] = i;
-
-    std::map<std::shared_ptr<const mesh::Mesh<U>>,
-             std::span<const std::int32_t>>
-        entity_maps = {{submesh, mesh_to_submesh}};
-
     // Next we compute the integration entities on the integration
     // domain `mesh`
     std::vector<std::int32_t> integration_entities
@@ -118,10 +107,21 @@ int main(int argc, char* argv[])
         subdomain_data
         = {{fem::IntegralType::cell, {{3, integration_entities}}}};
 
+    // A mixed-domain form involves functions defined over multiple
+    // meshes. The mesh passed to `create_form` is called the
+    // *integration domain mesh*. To assemble a mixed-domain form, we
+    // must supply an `EntityMap` for each additional mesh involved in
+    // the form, relating entities in that mesh to the integration
+    // domain mesh. In our case, `mesh` is the integration domain mesh,
+    // and the only other mesh in our form is `submesh`. Hence, we must
+    // provide the entity map object returned when we called
+    // `create_submesh`, which relates entities in `submesh` to entities
+    // in `mesh`.
+    //
     // We can now create the bilinear form
     fem::Form<T> a_mixed
         = fem::create_form<T>(*form_mixed_codim0_a_mixed, {V, W}, {}, {},
-                              subdomain_data, entity_maps, V->mesh());
+                              subdomain_data, {entity_map}, V->mesh());
 
     la::SparsityPattern sp_mixed = fem::create_sparsity_pattern(a_mixed);
     sp_mixed.finalize();
 
@@ -207,19 +207,20 @@ int main(int argc, char* argv[])
 
     // 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.
+    // `create_submesh` to get a `submesh` of those facets.  It also returns an
+    // `EntityMap` object, which relates entities in the submesh to entities in
+    // the original mesh. We will need this to assemble our mixed-domain form.
     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_data = [](auto& mesh, int tdim, auto&& dfacets)
     {
-      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);
-    }
+      auto [submesh, e_map, v_map, g_map]
+          = mesh::create_submesh(mesh, tdim, dfacets);
+      return std::pair(std::make_shared<mesh::Mesh<U>>(std::move(submesh)),
+                       std::move(e_map));
+    };
+    auto [submesh, entity_map] = submesh_data(*mesh, fdim, dfacets);
 
     // Create an element for `u_0`
     basix::cell::type submesh_cell_type
@@ -283,35 +284,25 @@ int main(int argc, char* argv[])
         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();
-    // Since not all facets in the mesh appear in the submesh, `submesh_to_mesh`
-    // is only injective. We therefore map nonexistent facets to -1 when
-    // creating the "inverse" map mesh_to_submesh
-    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}};
+    // on the `submesh`, our form involves more than one mesh. The mesh
+    // used to define the measure and passed to `create_form` is called
+    // the integration domain mesh (here, `mesh`). To assemble our
+    // mixed-domain form, we must provide an `EntityMap` for each
+    // additional mesh in the form. In this case, the only other mesh is
+    // `submesh`. Hence, we supply the entity map returned from
+    // `create_submesh`, which relates entities in `mesh` and `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, {});
-    // Since this form involves multiple domains (i.e. both `mesh` and `submesh`
-    // for the boundary condition), we must pass the entity maps just created.
-    // We must also tell the form which domain to integrate with respect to (in
-    // this case `mesh`)
+
+    // Since this form involves multiple domains (i.e. both `mesh` and
+    // `submesh` for the boundary condition), we must pass the entity
+    // maps just created. We must also tell the form which domain to
+    // integrate with respect to (in this case `mesh`)
     fem::Form<T> L = fem::create_form<T>(
         *form_mixed_poisson_L, {V}, {{"f", f}, {"u0", u0}}, {}, subdomain_data,
-        entity_maps, V->mesh());
+        {entity_map}, V->mesh());
 
     // Create solution finite element Function
     auto u = std::make_shared<fem::Function<T>>(V);