1414#include < dolfinx/fem/petsc.h>
1515#include < dolfinx/la/MatrixCSR.h>
1616#include < dolfinx/la/SparsityPattern.h>
17+ #include < dolfinx/mesh/EntityMap.h>
1718#include < map>
1819#include < memory>
1920#include < ranges>
@@ -76,37 +77,25 @@ int main(int argc, char* argv[])
7677 mesh::MeshTags<std::int32_t > cell_marker (mesh->topology (), tdim, cells,
7778 values);
7879
79- std::shared_ptr<mesh::Mesh<U>> submesh;
80- std::vector<std::int32_t > submesh_to_mesh;
80+ // We create a submesh consisting of only cells with a given tag by
81+ // calling `create_submesh`. This function also returns an
82+ // `EntityMap` object, which relates entities in the submesh to
83+ // entities in the original mesh. We will need this to assemble our
84+ // mixed-domain form.
85+ auto submesh_data = [](auto & mesh, int tdim, auto && subcells)
8186 {
82- auto [_submesh, _submesh_to_mesh, v_map, g_map]
83- = mesh::create_submesh (*mesh, tdim, cell_marker.find (2 ));
84- submesh = std::make_shared<mesh::Mesh<U>>(std::move (_submesh));
85- submesh_to_mesh = std::move (_submesh_to_mesh);
86- }
87+ auto [submesh, emap, v_map, g_map]
88+ = mesh::create_submesh (mesh, tdim, subcells);
89+ return std::pair (std::make_shared<mesh::Mesh<U>>(std::move (submesh)),
90+ std::move (emap));
91+ };
92+ auto [submesh, entity_map] = submesh_data (*mesh, tdim, cell_marker.find (2 ));
8793
8894 // We create the function space used for the trial space
8995 auto W
9096 = std::make_shared<fem::FunctionSpace<U>>(fem::create_functionspace<U>(
9197 submesh, std::make_shared<fem::FiniteElement<U>>(element)));
9298
93- // A mixed-domain form has functions defined over different meshes.
94- // The mesh associated with the measure (dx, ds, etc.) is called the
95- // integration domain. To assemble mixed-domain forms, maps must be
96- // provided taking entities in the integration domain to entities on
97- // each mesh in the form. Since one of our forms has a measure
98- // defined over `mesh` and involves a function defined over
99- // `submesh`, we must provide a map from entities in `mesh` to
100- // entities in `submesh`. This is simply the "inverse" of
101- // `submesh_to_mesh`.
102- std::vector<std::int32_t > mesh_to_submesh (num_cells_local, -1 );
103- for (std::size_t i = 0 ; i < submesh_to_mesh.size (); ++i)
104- mesh_to_submesh[submesh_to_mesh[i]] = i;
105-
106- std::map<std::shared_ptr<const mesh::Mesh<U>>,
107- std::span<const std::int32_t >>
108- entity_maps = {{submesh, mesh_to_submesh}};
109-
11099 // Next we compute the integration entities on the integration
111100 // domain `mesh`
112101 std::vector<std::int32_t > integration_entities
@@ -118,10 +107,21 @@ int main(int argc, char* argv[])
118107 subdomain_data
119108 = {{fem::IntegralType::cell, {{3 , integration_entities}}}};
120109
110+ // A mixed-domain form involves functions defined over multiple
111+ // meshes. The mesh passed to `create_form` is called the
112+ // *integration domain mesh*. To assemble a mixed-domain form, we
113+ // must supply an `EntityMap` for each additional mesh involved in
114+ // the form, relating entities in that mesh to the integration
115+ // domain mesh. In our case, `mesh` is the integration domain mesh,
116+ // and the only other mesh in our form is `submesh`. Hence, we must
117+ // provide the entity map object returned when we called
118+ // `create_submesh`, which relates entities in `submesh` to entities
119+ // in `mesh`.
120+ //
121121 // We can now create the bilinear form
122122 fem::Form<T> a_mixed
123123 = fem::create_form<T>(*form_mixed_codim0_a_mixed, {V, W}, {}, {},
124- subdomain_data, entity_maps , V->mesh ());
124+ subdomain_data, {entity_map} , V->mesh ());
125125
126126 la::SparsityPattern sp_mixed = fem::create_sparsity_pattern (a_mixed);
127127 sp_mixed.finalize ();
0 commit comments