From 4275cebb9de7bdb3a592084aa7b745b75ed16e67 Mon Sep 17 00:00:00 2001 From: Joe Dean Date: Sat, 2 Mar 2024 15:20:47 +0000 Subject: [PATCH] Add mixed-domain matrix assembly for facet integrals (codimension 0) (#3079) * Group args * Start on exterior facet integrals * Map int domain for ext facets * Pass to assembler * Use mapped cells * Sparsity for ext facet ints * Update test * Add comments * Reorder args * Add docstring * Small udpate * Fix dof trans * Improve name * Make marker into function * Format * Add helper (reduces code duplication when addin int facet integrals) * Rename * Create mixed-domain sparsity for int facet integrals * Dox fix * Fix sparsity builder * Reorder args * Group and pass different domains * Add mixed-domains for int facets * Add test * Tidy * Tidy * Docs * Tidy * Small changes * Remove stray character * Sparsity fix * Sparsity fix --------- Co-authored-by: Garth N. Wells --- cpp/dolfinx/fem/Form.h | 33 +++- cpp/dolfinx/fem/assemble_matrix_impl.h | 165 ++++++++++++++---- cpp/dolfinx/fem/sparsitybuild.cpp | 43 +++-- cpp/dolfinx/fem/sparsitybuild.h | 24 +-- cpp/dolfinx/fem/utils.h | 30 ++-- python/test/unit/fem/test_assemble_submesh.py | 55 ++++-- 6 files changed, 257 insertions(+), 93 deletions(-) diff --git a/cpp/dolfinx/fem/Form.h b/cpp/dolfinx/fem/Form.h index 1b90ede7b8c..593658a6a09 100644 --- a/cpp/dolfinx/fem/Form.h +++ b/cpp/dolfinx/fem/Form.h @@ -316,8 +316,8 @@ class Form throw std::runtime_error("No mesh entities for requested domain index."); } - /// @brief Compute the list of entity indices in `mesh` for the - /// ith integral (kernel) of a given type (i.e. cell, exterior facet, or + /// @brief Compute the list of entity indices in `mesh` for the ith + /// integral (kernel) of a given type (i.e. cell, exterior facet, or /// interior facet). /// /// @param type Integral type. @@ -339,9 +339,32 @@ class Form std::span entity_map = _entity_maps.at(msh_ptr); std::vector mapped_entities; mapped_entities.reserve(entities.size()); - std::transform(entities.begin(), entities.end(), - std::back_inserter(mapped_entities), - [&entity_map](auto e) { return entity_map[e]; }); + switch (type) + { + case IntegralType::cell: + { + std::transform(entities.begin(), entities.end(), + std::back_inserter(mapped_entities), + [&entity_map](auto e) { return entity_map[e]; }); + break; + } + case IntegralType::exterior_facet: + // Exterior and interior facets are treated the same + [[fallthrough]]; + case IntegralType::interior_facet: + { + for (std::size_t i = 0; i < entities.size(); i += 2) + { + // Add cell and the local facet index + mapped_entities.insert(mapped_entities.end(), + {entity_map[entities[i]], entities[i + 1]}); + } + break; + } + default: + throw std::runtime_error("Integral type not supported."); + } + return mapped_entities; } } diff --git a/cpp/dolfinx/fem/assemble_matrix_impl.h b/cpp/dolfinx/fem/assemble_matrix_impl.h index 5b15fef8001..4ff72297a49 100644 --- a/cpp/dolfinx/fem/assemble_matrix_impl.h +++ b/cpp/dolfinx/fem/assemble_matrix_impl.h @@ -113,8 +113,8 @@ void assemble_cells( coordinate_dofs.data(), nullptr, nullptr); // Compute A = P_0 \tilde{A} P_1^T (dof transformation) - P0(_Ae, cell_info1, c1, ndim1); // B = P0 \tilde{A} - P1T(_Ae, cell_info0, c0, ndim0); // A = B P1_T + P0(_Ae, cell_info0, c0, ndim1); // B = P0 \tilde{A} + P1T(_Ae, cell_info1, c1, ndim0); // A = B P1_T // Zero rows/columns for essential bcs auto dofs0 = std::span(dmap0.data_handle() + c0 * num_dofs0, num_dofs0); @@ -157,33 +157,74 @@ void assemble_cells( } } -/// Execute kernel over exterior facets and accumulate result in Mat +/// @brief Execute kernel over exterior facets and accumulate result in +/// a matrix. +/// @tparam T Matrix/form scalar type. +/// @param mat_set Function that accumulates computed entries into a +/// matrix. +/// @param x_dofmap Dofmap for the mesh geometry. +/// @param x Mesh geometry (coordinates). +/// @param facets Facet indices (in the integration domain mesh) to +/// execute the kernel over. +/// @param dofmap0 Test function (row) degree-of-freedom data holding +/// the (0) dofmap, (1) dofmap block size and (2) dofmap cell indices. +/// @param P0 Function that applies transformation P0.A in-place to +/// transform trial degrees-of-freedom. +/// @param dofmap1 Trial function (column) degree-of-freedom data +/// holding the (0) dofmap, (1) dofmap block size and (2) dofmap cell +/// indices. +/// @param P1T Function that applies transformation A.P1^T in-place to +/// transform trial degrees-of-freedom. +/// @param bc0 Marker for rows with Dirichlet boundary conditions applied +/// @param bc1 Marker for columns with Dirichlet boundary conditions applied +/// @param kernel Kernel function to execute over each cell. +/// @param coeffs The coefficient data array of shape (cells.size(), cstride), +/// flattened into row-major format. +/// @param cstride The coefficient stride +/// @param constants The constant data +/// @param cell_info0 The cell permutation information for the test function +/// mesh +/// @param cell_info1 The cell permutation information for the trial function +/// mesh template void assemble_exterior_facets( la::MatSet auto mat_set, mdspan2_t x_dofmap, std::span> x, - std::span facets, fem::DofTransformKernel auto P0, - mdspan2_t dofmap0, int bs0, fem::DofTransformKernel auto P1T, - mdspan2_t dofmap1, int bs1, std::span bc0, + std::span facets, + std::tuple> dofmap0, + fem::DofTransformKernel auto P0, + std::tuple> dofmap1, + fem::DofTransformKernel auto P1T, std::span bc0, std::span bc1, FEkernel auto kernel, std::span coeffs, int cstride, std::span constants, - std::span cell_info) + std::span cell_info0, + std::span cell_info1) { if (facets.empty()) return; + const auto [dmap0, bs0, facets0] = dofmap0; + const auto [dmap1, bs1, facets1] = dofmap1; + // Data structures used in assembly std::vector> coordinate_dofs(3 * x_dofmap.extent(1)); - const int num_dofs0 = dofmap0.extent(1); - const int num_dofs1 = dofmap1.extent(1); + const int num_dofs0 = dmap0.extent(1); + const int num_dofs1 = dmap1.extent(1); const int ndim0 = bs0 * num_dofs0; const int ndim1 = bs1 * num_dofs1; std::vector Ae(ndim0 * ndim1); std::span _Ae(Ae); assert(facets.size() % 2 == 0); + assert(facets0.size() == facets.size()); + assert(facets1.size() == facets.size()); for (std::size_t index = 0; index < facets.size(); index += 2) { + // Cell in the integration domain std::int32_t cell = facets[index]; + // Cell in the test function mesh + std::int32_t cell0 = facets0[index]; + // Cell in the trial function mesh + std::int32_t cell1 = facets1[index]; std::int32_t local_facet = facets[index + 1]; // Get cell coordinates/geometry @@ -201,12 +242,12 @@ void assemble_exterior_facets( kernel(Ae.data(), coeffs.data() + index / 2 * cstride, constants.data(), coordinate_dofs.data(), &local_facet, nullptr); - P0(_Ae, cell_info, cell, ndim1); - P1T(_Ae, cell_info, cell, ndim0); + P0(_Ae, cell_info0, cell0, ndim1); + P1T(_Ae, cell_info1, cell1, ndim0); // Zero rows/columns for essential bcs - auto dofs0 = std::span(dofmap0.data_handle() + cell * num_dofs0, num_dofs0); - auto dofs1 = std::span(dofmap1.data_handle() + cell * num_dofs1, num_dofs1); + auto dofs0 = std::span(dmap0.data_handle() + cell0 * num_dofs0, num_dofs0); + auto dofs1 = std::span(dmap1.data_handle() + cell1 * num_dofs1, num_dofs1); if (!bc0.empty()) { for (int i = 0; i < num_dofs0; ++i) @@ -243,22 +284,59 @@ void assemble_exterior_facets( } } -/// Execute kernel over interior facets and accumulate result in Mat +/// @brief Execute kernel over interior facets and accumulate result in a +/// matrix. +/// @tparam T Matrix/form scalar type. +/// @param mat_set Function that accumulates computed entries into a +/// matrix. +/// @param x_dofmap Dofmap for the mesh geometry. +/// @param x Mesh geometry (coordinates). +/// @param num_cell_facets Number of facets of a cell +/// @param facets Facet indices (in the integration domain mesh) to +/// execute the kernel over. +/// @param dofmap0 Test function (row) degree-of-freedom data holding +/// the (0) dofmap, (1) dofmap block size and (2) dofmap cell indices. +/// @param P0 Function that applies transformation P0.A in-place to +/// transform trial degrees-of-freedom. +/// @param dofmap1 Trial function (column) degree-of-freedom data +/// holding the (0) dofmap, (1) dofmap block size and (2) dofmap cell +/// indices. +/// @param P1T Function that applies transformation A.P1^T in-place to +/// transform trial degrees-of-freedom. +/// @param bc0 Marker for rows with Dirichlet boundary conditions applied +/// @param bc1 Marker for columns with Dirichlet boundary conditions applied +/// @param kernel Kernel function to execute over each cell. +/// @param coeffs The coefficient data array of shape (cells.size(), cstride), +/// flattened into row-major format. +/// @param cstride The coefficient stride +/// @param offsets The coefficient offsets +/// @param constants The constant data +/// @param cell_info0 The cell permutation information for the test function +/// mesh +/// @param cell_info1 The cell permutation information for the trial function +/// mesh +/// @param get_perm Function to apply facet permutations template void assemble_interior_facets( la::MatSet auto mat_set, mdspan2_t x_dofmap, std::span> x, int num_cell_facets, - std::span facets, fem::DofTransformKernel auto P0, - const DofMap& dofmap0, int bs0, fem::DofTransformKernel auto P1T, - const DofMap& dofmap1, int bs1, std::span bc0, + std::span facets, + std::tuple> dofmap0, + fem::DofTransformKernel auto P0, + std::tuple> dofmap1, + fem::DofTransformKernel auto P1T, std::span bc0, std::span bc1, FEkernel auto kernel, std::span coeffs, int cstride, std::span offsets, - std::span constants, std::span cell_info, + std::span constants, std::span cell_info0, + std::span cell_info1, const std::function& get_perm) { if (facets.empty()) return; + const auto [dmap0, bs0, facets0] = dofmap0; + const auto [dmap1, bs1, facets1] = dofmap1; + // Data structures used in assembly using X = scalar_value_type_t; std::vector coordinate_dofs(2 * x_dofmap.extent(1) * 3); @@ -273,11 +351,18 @@ void assemble_interior_facets( // Temporaries for joint dofmaps std::vector dmapjoint0, dmapjoint1; assert(facets.size() % 4 == 0); + assert(facets0.size() == facets.size()); + assert(facets1.size() == facets.size()); for (std::size_t index = 0; index < facets.size(); index += 4) { - std::array cells = {facets[index], facets[index + 2]}; - std::array local_facet - = {facets[index + 1], facets[index + 3]}; + // Cells in integration domain, test function domain and trial + // function domain + std::array cells{facets[index], facets[index + 2]}; + std::array cells0{facets0[index], facets0[index + 2]}; + std::array cells1{facets1[index], facets1[index + 2]}; + + // Local facets indices + std::array local_facet{facets[index + 1], facets[index + 3]}; // Get cell geometry auto x_dofs0 = MDSPAN_IMPL_STANDARD_NAMESPACE:: @@ -298,15 +383,15 @@ void assemble_interior_facets( } // Get dof maps for cells and pack - std::span dmap0_cell0 = dofmap0.cell_dofs(cells[0]); - std::span dmap0_cell1 = dofmap0.cell_dofs(cells[1]); + std::span dmap0_cell0 = dmap0.cell_dofs(cells0[0]); + std::span dmap0_cell1 = dmap0.cell_dofs(cells0[1]); dmapjoint0.resize(dmap0_cell0.size() + dmap0_cell1.size()); std::copy(dmap0_cell0.begin(), dmap0_cell0.end(), dmapjoint0.begin()); std::copy(dmap0_cell1.begin(), dmap0_cell1.end(), std::next(dmapjoint0.begin(), dmap0_cell0.size())); - std::span dmap1_cell0 = dofmap1.cell_dofs(cells[0]); - std::span dmap1_cell1 = dofmap1.cell_dofs(cells[1]); + std::span dmap1_cell0 = dmap1.cell_dofs(cells1[0]); + std::span dmap1_cell1 = dmap1.cell_dofs(cells1[1]); dmapjoint1.resize(dmap1_cell0.size() + dmap1_cell1.size()); std::copy(dmap1_cell0.begin(), dmap1_cell0.end(), dmapjoint1.begin()); std::copy(dmap1_cell1.begin(), dmap1_cell1.end(), @@ -336,9 +421,9 @@ void assemble_interior_facets( std::span sub_Ae0 = _Ae.subspan(bs0 * dmap0_cell0.size() * num_cols, bs0 * dmap0_cell1.size() * num_cols); - P0(_Ae, cell_info, cells[0], num_cols); - P0(sub_Ae0, cell_info, cells[1], num_cols); - P1T(_Ae, cell_info, cells[0], num_rows); + P0(_Ae, cell_info0, cells0[0], num_cols); + P0(sub_Ae0, cell_info0, cells0[1], num_cols); + P1T(_Ae, cell_info1, cells1[0], num_rows); for (int row = 0; row < num_rows; ++row) { @@ -346,7 +431,7 @@ void assemble_interior_facets( // the block matrix, so each row needs a separate span access std::span sub_Ae1 = _Ae.subspan( row * num_cols + bs1 * dmap1_cell0.size(), bs1 * dmap1_cell1.size()); - P1T(sub_Ae1, cell_info, cells[1], 1); + P1T(sub_Ae1, cell_info1, cells1[1], 1); } // Zero rows/columns for essential bcs @@ -459,10 +544,11 @@ void assemble_matrix( assert(fn); auto& [coeffs, cstride] = coefficients.at({IntegralType::exterior_facet, i}); - impl::assemble_exterior_facets(mat_set, x_dofmap, x, - a.domain(IntegralType::exterior_facet, i), - P0, dofs0, bs0, P1T, dofs1, bs1, bc0, bc1, - fn, coeffs, cstride, constants, cell_info0); + impl::assemble_exterior_facets( + mat_set, x_dofmap, x, a.domain(IntegralType::exterior_facet, i), + {dofs0, bs0, a.domain(IntegralType::exterior_facet, i, *mesh0)}, P0, + {dofs1, bs1, a.domain(IntegralType::exterior_facet, i, *mesh1)}, P1T, + bc0, bc1, fn, coeffs, cstride, constants, cell_info0, cell_info1); } if (a.num_integrals(IntegralType::interior_facet) > 0) @@ -488,11 +574,14 @@ void assemble_matrix( assert(fn); auto& [coeffs, cstride] = coefficients.at({IntegralType::interior_facet, i}); - impl::assemble_interior_facets(mat_set, x_dofmap, x, num_cell_facets, - a.domain(IntegralType::interior_facet, i), - P0, *dofmap0, bs0, P1T, *dofmap1, bs1, bc0, - bc1, fn, coeffs, cstride, c_offsets, - constants, cell_info0, get_perm); + impl::assemble_interior_facets( + mat_set, x_dofmap, x, num_cell_facets, + a.domain(IntegralType::interior_facet, i), + {*dofmap0, bs0, a.domain(IntegralType::interior_facet, i, *mesh0)}, + P0, + {*dofmap1, bs1, a.domain(IntegralType::interior_facet, i, *mesh1)}, + P1T, bc0, bc1, fn, coeffs, cstride, c_offsets, constants, cell_info0, + cell_info1, get_perm); } } } diff --git a/cpp/dolfinx/fem/sparsitybuild.cpp b/cpp/dolfinx/fem/sparsitybuild.cpp index a1d5e95441d..7a8cc1443b3 100644 --- a/cpp/dolfinx/fem/sparsitybuild.cpp +++ b/cpp/dolfinx/fem/sparsitybuild.cpp @@ -25,24 +25,37 @@ void sparsitybuild::cells( } //----------------------------------------------------------------------------- void sparsitybuild::interior_facets( - la::SparsityPattern& pattern, std::span facets, + la::SparsityPattern& pattern, + std::array, 2> cells, std::array, 2> dofmaps) { - std::array, 2> macro_dofs; - for (std::size_t index = 0; index < facets.size(); index += 2) + std::span cells0 = cells[0]; + std::span cells1 = cells[1]; + assert(cells0.size() == cells1.size()); + const DofMap& dofmap0 = dofmaps[0]; + const DofMap& dofmap1 = dofmaps[1]; + + // Iterate over facets + std::vector macro_dofs0, macro_dofs1; + for (std::size_t f = 0; f < cells[0].size(); f += 2) { - std::int32_t cell0 = facets[index]; - std::int32_t cell1 = facets[index + 1]; - for (std::size_t i = 0; i < 2; ++i) - { - auto cell_dofs0 = dofmaps[i].get().cell_dofs(cell0); - auto cell_dofs1 = dofmaps[i].get().cell_dofs(cell1); - macro_dofs[i].resize(cell_dofs0.size() + cell_dofs1.size()); - std::copy(cell_dofs0.begin(), cell_dofs0.end(), macro_dofs[i].begin()); - std::copy(cell_dofs1.begin(), cell_dofs1.end(), - std::next(macro_dofs[i].begin(), cell_dofs0.size())); - } - pattern.insert(macro_dofs[0], macro_dofs[1]); + // Test function dofs (sparsity pattern rows) + auto dofs00 = dofmap0.cell_dofs(cells0[f]); + auto dofs01 = dofmap0.cell_dofs(cells0[f + 1]); + macro_dofs0.resize(dofs00.size() + dofs01.size()); + std::copy(dofs00.begin(), dofs00.end(), macro_dofs0.begin()); + std::copy(dofs01.begin(), dofs01.end(), + std::next(macro_dofs0.begin(), dofs00.size())); + + // Trial function dofs (sparsity pattern columns) + auto dofs10 = dofmap1.cell_dofs(cells1[f]); + auto dofs11 = dofmap1.cell_dofs(cells1[f + 1]); + macro_dofs1.resize(dofs10.size() + dofs11.size()); + std::copy(dofs10.begin(), dofs10.end(), macro_dofs1.begin()); + std::copy(dofs11.begin(), dofs11.end(), + std::next(macro_dofs1.begin(), dofs10.size())); + + pattern.insert(macro_dofs0, macro_dofs1); } } //----------------------------------------------------------------------------- diff --git a/cpp/dolfinx/fem/sparsitybuild.h b/cpp/dolfinx/fem/sparsitybuild.h index 14d53b04d8e..2279cd0445f 100644 --- a/cpp/dolfinx/fem/sparsitybuild.h +++ b/cpp/dolfinx/fem/sparsitybuild.h @@ -25,10 +25,10 @@ namespace sparsitybuild { /// @brief Iterate over cells and insert entries into sparsity pattern. /// -/// Inserts the rectangular blocks of indices `dofmap[0][cells[0][i]] x -/// dofmap[1][cells[1][i]]` into the sparsity pattern, i.e. entries -/// (dofmap[0][cells[0][i]][k0], dofmap[0][cells[0][i]][k1])` will -/// appear in the sparsity pattern. +/// Inserts the rectangular blocks of indices `dofmap[0][cells[0][i]] x +/// dofmap[1][cells[1][i]]` into the sparsity pattern, i.e. entries +/// (dofmap[0][cells[0][i]][k0], dofmap[0][cells[0][i]][k1])` will +/// appear in the sparsity pattern. /// /// @param pattern Sparsity pattern to insert into. /// @param cells Lists of cells to iterate over. `cells[0]` and @@ -42,17 +42,19 @@ void cells(la::SparsityPattern& pattern, /// @brief Iterate over interior facets and insert entries into sparsity /// pattern. /// -/// Inserts the rectangular blocks of indices `[dofmap[0][cell0], -/// dofmap[0][cell1]] x [dofmap[1][cell0] + dofmap[1][cell1]]` where -/// `cell0` and `cell1` are the two cells attached to a facet. +/// Inserts the rectangular blocks of indices `[dofmap[0][cell0], +/// dofmap[0][cell1]] x [dofmap[1][cell0] + dofmap[1][cell1]]` where +/// `cell0` and `cell1` are the two cells attached to a facet. /// -/// @param[in,out] pattern Sparsity pattern to insert into -/// @param[in] facets Facets as `(cell0, cell1)` pairs (row-major) for -/// each facet. +/// @param[in,out] pattern Sparsity pattern to insert into. +/// @param[in] cells Cells to index into each dofmap. `cells[i]` is a +/// list of `(cell0, cell1)` pairs for each interior facet to index into +/// `dofmap[i]`. `cells[0]` and `cells[1]` must have the same size. /// @param[in] dofmaps Dofmaps to use in building the sparsity pattern. /// @note The sparsity pattern is not finalised. void interior_facets( - la::SparsityPattern& pattern, std::span facets, + la::SparsityPattern& pattern, + std::array, 2> cells, std::array, 2> dofmaps); } // namespace sparsitybuild diff --git a/cpp/dolfinx/fem/utils.h b/cpp/dolfinx/fem/utils.h index 1c256cef481..05c19fb3872 100644 --- a/cpp/dolfinx/fem/utils.h +++ b/cpp/dolfinx/fem/utils.h @@ -185,6 +185,16 @@ la::SparsityPattern create_sparsity_pattern(const Form& a) const std::array bs = {dofmaps[0].get().index_map_bs(), dofmaps[1].get().index_map_bs()}; + auto extract_cells = [](std::span facets) + { + assert(facets.size() % 2 == 0); + std::vector cells; + cells.reserve(facets.size() / 2); + for (std::size_t i = 0; i < facets.size(); i += 2) + cells.push_back(facets[i]); + return cells; + }; + // Create and build sparsity pattern la::SparsityPattern pattern(mesh->comm(), index_maps, bs); for (auto type : types) @@ -203,23 +213,19 @@ la::SparsityPattern create_sparsity_pattern(const Form& a) case IntegralType::interior_facet: for (int id : ids) { - std::span facets = a.domain(type, id); - std::vector f; - f.reserve(facets.size() / 2); - for (std::size_t i = 0; i < facets.size(); i += 4) - f.insert(f.end(), {facets[i], facets[i + 2]}); - sparsitybuild::interior_facets(pattern, f, {{dofmaps[0], dofmaps[1]}}); + sparsitybuild::interior_facets( + pattern, + {extract_cells(a.domain(type, id, *mesh0)), + extract_cells(a.domain(type, id, *mesh1))}, + {{dofmaps[0], dofmaps[1]}}); } break; case IntegralType::exterior_facet: for (int id : ids) { - std::span facets = a.domain(type, id); - std::vector cells; - cells.reserve(facets.size() / 2); - for (std::size_t i = 0; i < facets.size(); i += 2) - cells.push_back(facets[i]); - sparsitybuild::cells(pattern, {cells, cells}, + sparsitybuild::cells(pattern, + {extract_cells(a.domain(type, id, *mesh0)), + extract_cells(a.domain(type, id, *mesh1))}, {{dofmaps[0], dofmaps[1]}}); } break; diff --git a/python/test/unit/fem/test_assemble_submesh.py b/python/test/unit/fem/test_assemble_submesh.py index 0c80827336d..72026a9f7ea 100644 --- a/python/test/unit/fem/test_assemble_submesh.py +++ b/python/test/unit/fem/test_assemble_submesh.py @@ -118,7 +118,7 @@ def test_submesh_facet_assembly(n, k, space, ghost_mode): assert np.isclose(s_submesh, s_square_mesh) -@pytest.mark.parametrize("n", [2, 6]) +@pytest.mark.parametrize("n", [4, 6]) @pytest.mark.parametrize("k", [1, 4]) @pytest.mark.parametrize("space", ["Lagrange", "Discontinuous Lagrange"]) @pytest.mark.parametrize("ghost_mode", [GhostMode.none, GhostMode.shared_facet]) @@ -126,24 +126,54 @@ def test_mixed_dom_codim_0(n, k, space, ghost_mode): """Test assembling a form where the trial and test functions are defined on different meshes""" + def create_meshtags(msh, dim, entities, tag): + perm = np.argsort(entities) + values = np.full_like(entities, tag, dtype=np.intc) + return meshtags(msh, dim, entities[perm], values[perm]) + msh = create_rectangle( MPI.COMM_WORLD, ((0.0, 0.0), (2.0, 1.0)), (2 * n, n), ghost_mode=ghost_mode ) # Locate cells in left half of mesh and create mesh tags tdim = msh.topology.dim - cells = locate_entities(msh, tdim, lambda x: x[0] <= 1.0) - perm = np.argsort(cells) tag = 1 - values = np.full_like(cells, tag, dtype=np.intc) - ct = meshtags(msh, tdim, cells[perm], values[perm]) + cells = locate_entities(msh, tdim, lambda x: x[0] <= 1.0) + ct = create_meshtags(msh, tdim, cells, tag) + + # Locate facets on left boundary and create mesh tags + def boundary_marker(x): + return np.isclose(x[0], 0.0) - # Create integration measure on the mesh + fdim = tdim - 1 + facets = locate_entities_boundary(msh, fdim, boundary_marker) + ft = create_meshtags(msh, fdim, facets, tag) + + # Locate some interior facets and create mesh tags + def int_facets_marker(x): + dist = 1 / (2 * n) + return (x[0] > dist) & (x[0] < 1 - dist) & (x[1] > dist) & (x[1] < 1 - dist) + + int_facets = locate_entities(msh, fdim, int_facets_marker) + int_ft = create_meshtags(msh, fdim, int_facets, tag) + + # Create integration measures on the mesh dx_msh = ufl.Measure("dx", domain=msh, subdomain_data=ct) + ds_msh = ufl.Measure("ds", domain=msh, subdomain_data=ft) + dS_msh = ufl.Measure("dS", domain=msh, subdomain_data=int_ft) # Create a submesh of the left half of the mesh smsh, smsh_to_msh = create_submesh(msh, tdim, cells)[:2] + # Create some integration measures on the submesh + facets_smsh = locate_entities_boundary(smsh, fdim, boundary_marker) + ft_smsh = create_meshtags(smsh, fdim, facets_smsh, tag) + ds_smsh = ufl.Measure("ds", domain=smsh, subdomain_data=ft_smsh) + + int_facets_smsh = locate_entities(smsh, fdim, int_facets_marker) + int_ft_smsh = create_meshtags(smsh, fdim, int_facets_smsh, tag) + dS_smsh = ufl.Measure("dS", domain=smsh, subdomain_data=int_ft_smsh) + # Define function spaces over the mesh and submesh V_msh = fem.functionspace(msh, (space, k)) V_smsh = fem.functionspace(smsh, (space, k)) @@ -154,11 +184,12 @@ def test_mixed_dom_codim_0(n, k, space, ghost_mode): # Test function on the submesh w = ufl.TestFunction(V_smsh) - def ufl_form(u, v, dx): - return ufl.inner(u, v) * dx + # Define a UFL form + def ufl_form(u, v, dx, ds, dS): + return ufl.inner(u, v) * dx + ufl.inner(u, v) * ds + ufl.inner(u("+"), v("-")) * dS - # Define form to compare to and assemble - a = fem.form(ufl_form(u, v, dx_msh(tag))) + # Single-domain assembly over msh as a reference + a = fem.form(ufl_form(u, v, dx_msh(tag), ds_msh(tag), dS_msh(tag))) A = fem.assemble_matrix(a) A.scatter_reverse() @@ -166,7 +197,7 @@ def ufl_form(u, v, dx): # Entity maps must map cells in smsh (the integration domain mesh) to # cells in msh entity_maps = {msh._cpp_object: np.array(smsh_to_msh, dtype=np.int32)} - a0 = fem.form(ufl_form(u, w, ufl.dx(smsh)), entity_maps=entity_maps) + a0 = fem.form(ufl_form(u, w, ufl.dx(smsh), ds_smsh(tag), dS_smsh(tag)), entity_maps=entity_maps) A0 = fem.assemble_matrix(a0) A0.scatter_reverse() assert np.isclose(A0.squared_norm(), A.squared_norm()) @@ -180,7 +211,7 @@ def ufl_form(u, v, dx): msh_to_smsh[smsh_to_msh] = np.arange(len(smsh_to_msh)) entity_maps = {smsh._cpp_object: np.array(msh_to_smsh, dtype=np.int32)} - a1 = fem.form(ufl_form(u, w, dx_msh(tag)), entity_maps=entity_maps) + a1 = fem.form(ufl_form(u, w, dx_msh(tag), ds_msh(tag), dS_msh(tag)), entity_maps=entity_maps) A1 = fem.assemble_matrix(a1) A1.scatter_reverse() assert np.isclose(A1.squared_norm(), A.squared_norm())