diff --git a/cpp/demo/mixed_poisson/main.cpp b/cpp/demo/mixed_poisson/main.cpp index 7e578be6785..3db8ae77bb8 100644 --- a/cpp/demo/mixed_poisson/main.cpp +++ b/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 #include #include #include @@ -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::create_rectangle(MPI_COMM_WORLD, {{{0.0, 0.0}, {1.0, 1.0}}}, - {32, 32}, mesh::CellType::triangle)); + auto mesh = std::make_shared>(mesh::create_rectangle( + MPI_COMM_WORLD, {{{0, 0}, {1, 1}}}, {32, 32}, cell_type)); // Create Basix elements - auto RT = basix::create_element( - 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( - 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(basix::element::family::RT, basix_cell_type, + 1, basix::element::lagrange_variant::unset, + basix::element::dpc_variant::unset, false); + auto P0 + = basix::create_element(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>( @@ -149,7 +154,7 @@ int main(int argc, char* argv[]) auto W0 = std::make_shared>(V0->collapse().first); auto W1 = std::make_shared>(V1->collapse().first); - // Create source function and interpolate sin(5x) + 1 + // Create source function and interpolate $\sin(5x) + 1$ auto f = std::make_shared>(W1); f->interpolate( [](auto x) -> std::pair, std::vector> @@ -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>(W1); - u0->interpolate( - [](auto x) -> std::pair, std::vector> - { - std::vector 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>(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 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> submesh; + std::vector submesh_to_mesh; + { + auto [_submesh, _submesh_to_mesh, v_map, g_map] + = mesh::create_submesh(*mesh, fdim, dfacets); + submesh = std::make_shared>(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>( + basix::create_element(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::create_functionspace(submesh, Qe)); + + // Boundary condition value for u and interpolate $20 y + 1$ + auto u0 = std::make_shared>(Q); + u0->interpolate( + [](auto x) -> std::pair, std::vector> + { + std::vector 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({*u0}, 0); + + // Compute facets with $\sigma$ (flux) boundary condition facets, + // which is `{all boundary facet} - {u0 boundary facets}` std::vector 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 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 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>>> 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 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::span> + entity_maps = {{submesh, mesh_to_submesh}}; + // Define variational forms and attach he required data fem::Form a = fem::create_form(*form_mixed_poisson_a, {V, V}, {}, {}, subdomain_data, {}); - fem::Form L - = fem::create_form(*form_mixed_poisson_L, {V}, - {{"f", f}, {"u0", u0}}, {}, subdomain_data, {}); + fem::Form L = fem::create_form( + *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>(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); diff --git a/cpp/demo/mixed_poisson/mixed_poisson.py b/cpp/demo/mixed_poisson/mixed_poisson.py index 30c5452f55d..851c82d30a6 100644 --- a/cpp/demo/mixed_poisson/mixed_poisson.py +++ b/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