From a50f9782e1cdd8fbfde3166b3f9258a09c6ef5e5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= Date: Wed, 3 Sep 2025 09:16:41 +0000 Subject: [PATCH 01/37] Unify assemblers of vertex, exterior facet and ridge integrals (new) --- cpp/dolfinx/fem/Form.h | 3 +- cpp/dolfinx/fem/assemble_vector_impl.h | 215 ++++++++++--------------- cpp/dolfinx/fem/utils.cpp | 16 +- cpp/dolfinx/fem/utils.h | 120 ++++++++++++-- python/dolfinx/fem/forms.py | 5 + python/dolfinx/wrappers/fem.cpp | 3 +- python/test/unit/fem/test_assembler.py | 66 ++++++++ 7 files changed, 283 insertions(+), 145 deletions(-) diff --git a/cpp/dolfinx/fem/Form.h b/cpp/dolfinx/fem/Form.h index 8a361dcf0c..7787d75211 100644 --- a/cpp/dolfinx/fem/Form.h +++ b/cpp/dolfinx/fem/Form.h @@ -40,7 +40,8 @@ enum class IntegralType : std::int8_t cell = 0, ///< Cell exterior_facet = 1, ///< Exterior facet interior_facet = 2, ///< Interior facet - vertex = 3 ///< Vertex + vertex = 3, ///< Vertex + ridge = 4 ///< Ridge }; /// @brief Represents integral data, containing the kernel, and a list diff --git a/cpp/dolfinx/fem/assemble_vector_impl.h b/cpp/dolfinx/fem/assemble_vector_impl.h index 3e8c1c2a13..46b449102b 100644 --- a/cpp/dolfinx/fem/assemble_vector_impl.h +++ b/cpp/dolfinx/fem/assemble_vector_impl.h @@ -740,19 +740,19 @@ void assemble_cells( /// `(cells.size(), coeffs_per_cell)`. /// @param[in] cell_info0 The cell permutation information for the test /// function mesh. -/// @param[in] perms Facet permutation integer. Empty if facet +/// @param[in] perms Entity permutation integer. Empty if entity /// permutations are not required. template ::value_type> requires std::is_same_v::value_type, T> -void assemble_exterior_facets( +void assemble_exterior_entities( fem::DofTransformKernel auto P0, V&& b, mdspan2_t x_dofmap, md::mdspan, md::extents> x, md::mdspan> - facets, + entities, std::tuple>> @@ -762,24 +762,24 @@ void assemble_exterior_facets( std::span cell_info0, md::mdspan> perms) { - if (facets.empty()) + if (entities.empty()) return; - const auto [dmap, bs, facets0] = dofmap; + const auto [dmap, bs, entities0] = dofmap; assert(_bs < 0 or _bs == bs); // Create data structures used in assembly const int num_dofs = dmap.extent(1); std::vector> cdofs(3 * x_dofmap.extent(1)); std::vector be(bs * num_dofs); - assert(facets0.size() == facets.size()); - for (std::size_t f = 0; f < facets.extent(0); ++f) + assert(entities0.size() == entities.size()); + for (std::size_t f = 0; f < entities.extent(0); ++f) { // Cell in the integration domain, local facet index relative to the // integration domain cell, and cell in the test function mesh - std::int32_t cell = facets(f, 0); - std::int32_t local_facet = facets(f, 1); - std::int32_t cell0 = facets0(f, 0); + std::int32_t cell = entities(f, 0); + std::int32_t local_entity = entities(f, 1); + std::int32_t cell0 = entities0(f, 0); // Get cell coordinates/geometry auto x_dofs = md::submdspan(x_dofmap, cell, md::full_extent); @@ -787,12 +787,12 @@ void assemble_exterior_facets( std::copy_n(&x(x_dofs[i], 0), 3, std::next(cdofs.begin(), 3 * i)); // Permutations - std::uint8_t perm = perms.empty() ? 0 : perms(cell, local_facet); + std::uint8_t perm = perms.empty() ? 0 : perms(cell, local_entity); // Tabulate element vector std::ranges::fill(be, 0); kernel(be.data(), &coeffs(f, 0), constants.data(), cdofs.data(), - &local_facet, &perm, nullptr); + &local_entity, &perm, nullptr); P0(be, cell_info0, cell0, 1); // Add element vector to global vector @@ -941,94 +941,6 @@ void assemble_interior_facets( } } -/// @brief Execute kernel over a set of vertices and accumulate result in -/// vector. -/// -/// @tparam T Scalar type -/// @tparam _bs Block size of the form test function dof map. If less -/// than zero the block size is determined at runtime. If `_bs` is -/// positive the block size is used as a compile-time constant, which -/// has performance benefits. -/// @param[in] P0 Function that applies transformation `P0.b` in-place -/// to `b` to transform test degrees-of-freedom. -/// @param[in,out] b Array to accumulate into. -/// @param[in] x_dofmap Dofmap for the mesh geometry. -/// @param[in] x Mesh geometry (coordinates). -/// @param[in] vertices Vertex indices `(vertices.size(), 2)` - first entry -/// holds the index of the cell adjacent to the vertex, and the second -/// stores the local index of the vertex within the cell. -/// @param[in] dofmap Test function (row) degree-of-freedom data holding -/// the (0) dofmap, (1) dofmap block size and (2) dofmap cell indices. -/// @param[in] kernel Kernel function to execute over each cell. -/// @param[in] constants Constant coefficient data in the kernel. -/// @param[in] coeffs Coefficient data in the kernel. It has shape -/// `(vertices.size(), num_cell_coeffs)`. `coeffs(i, j)` is the `j`th -/// coefficient for cell `i`. -/// @param[in] cell_info0 Cell permutation information for the test -/// function mesh. -template -void assemble_vertices( - fem::DofTransformKernel auto P0, std::span b, mdspan2_t x_dofmap, - md::mdspan, - md::extents> - x, - md::mdspan> - vertices, - std::tuple>> - dofmap, - FEkernel auto kernel, std::span constants, - md::mdspan> coeffs, - std::span cell_info0) -{ - if (vertices.empty()) - return; - - const auto [dmap, bs, vertices0] = dofmap; - assert(_bs < 0 or _bs == bs); - - // Create data structures used in assembly - std::vector> cdofs(3 * x_dofmap.extent(1)); - std::vector be(bs * dmap.extent(1)); - - // Iterate over active vertices - for (std::size_t index = 0; index < vertices.extent(0); ++index) - { - // Integration domain cell, local index, and test function cell - std::int32_t cell = vertices(index, 0); - std::int32_t local_index = vertices(index, 1); - std::int32_t c0 = vertices0(index, 0); - - // Get cell coordinates/geometry - auto x_dofs = md::submdspan(x_dofmap, cell, md::full_extent); - for (std::size_t i = 0; i < x_dofs.size(); ++i) - std::copy_n(&x(x_dofs[i], 0), 3, std::next(cdofs.begin(), 3 * i)); - - // Tabulate vector for vertex - std::ranges::fill(be, 0); - kernel(be.data(), &coeffs(index, 0), constants.data(), cdofs.data(), - &local_index, nullptr, nullptr); - P0(be, cell_info0, c0, 1); - - // Scatter vertex vector to 'global' vector array - auto dofs = md::submdspan(dmap, c0, md::full_extent); - if constexpr (_bs > 0) - { - for (std::size_t i = 0; i < dofs.size(); ++i) - for (int k = 0; k < _bs; ++k) - b[_bs * dofs[i] + k] += be[_bs * i + k]; - } - else - { - for (std::size_t i = 0; i < dofs.size(); ++i) - for (int k = 0; k < bs; ++k) - b[bs * dofs[i] + k] += be[bs * i + k]; - } - } -} - /// Modify RHS vector to account for boundary condition such that: /// /// b <- b - alpha * A.(x_bc - x0) @@ -1368,7 +1280,7 @@ void assemble_vector( } } - md::mdspan> perms; + md::mdspan> facet_perms; if (L.needs_facet_permutations()) { mesh::CellType cell_type = mesh->topology()->cell_types()[cell_type_idx]; @@ -1377,8 +1289,8 @@ void assemble_vector( mesh->topology_mutable()->create_entity_permutations(); const std::vector& p = mesh->topology()->get_facet_permutations(); - perms = md::mdspan(p.data(), p.size() / num_facets_per_cell, - num_facets_per_cell); + facet_perms = md::mdspan(p.data(), p.size() / num_facets_per_cell, + num_facets_per_cell); } using mdspanx2_t @@ -1398,24 +1310,24 @@ void assemble_vector( assert((facets.size() / 2) * cstride == coeffs.size()); if (bs == 1) { - impl::assemble_exterior_facets<1>( + impl::assemble_exterior_entities<1>( P0, b, x_dofmap, x, facets, {dofs, bs, facets1}, fn, constants, md::mdspan(coeffs.data(), facets.extent(0), cstride), cell_info0, - perms); + facet_perms); } else if (bs == 3) { - impl::assemble_exterior_facets<3>( + impl::assemble_exterior_entities<3>( P0, b, x_dofmap, x, facets, {dofs, bs, facets1}, fn, constants, md::mdspan(coeffs.data(), facets.size() / 2, cstride), cell_info0, - perms); + facet_perms); } else { - impl::assemble_exterior_facets( + impl::assemble_exterior_entities( P0, b, x_dofmap, x, facets, {dofs, bs, facets1}, fn, constants, md::mdspan(coeffs.data(), facets.size() / 2, cstride), cell_info0, - perms); + facet_perms); } } @@ -1444,7 +1356,7 @@ void assemble_vector( mdspanx22_t(facets1.data(), facets1.size() / 4, 2, 2)}, fn, constants, mdspanx2x_t(coeffs.data(), facets.size() / 4, 2, cstride), - cell_info0, perms); + cell_info0, facet_perms); } else if (bs == 3) { @@ -1455,7 +1367,7 @@ void assemble_vector( mdspanx22_t(facets1.data(), facets1.size() / 4, 2, 2)}, fn, constants, mdspanx2x_t(coeffs.data(), facets.size() / 4, 2, cstride), - cell_info0, perms); + cell_info0, facet_perms); } else { @@ -1466,34 +1378,85 @@ void assemble_vector( mdspanx22_t(facets1.data(), facets1.size() / 4, 2, 2)}, fn, constants, mdspanx2x_t(coeffs.data(), facets.size() / 4, 2, cstride), - cell_info0, perms); + cell_info0, facet_perms); } } + md::mdspan> vertex_perms; for (int i = 0; i < L.num_integrals(IntegralType::vertex, 0); ++i) { auto fn = L.kernel(IntegralType::vertex, i, 0); assert(fn); - std::span vertices = L.domain(IntegralType::vertex, i, cell_type_idx); - std::span vertices0 - = L.domain_arg(IntegralType::vertex, 0, i, cell_type_idx); + std::span v = L.domain(IntegralType::vertex, i, cell_type_idx); + mdspanx2_t vertices(v.data(), v.size() / 2, 2); + std::span v1 = L.domain_arg(IntegralType::vertex, 0, i, cell_type_idx); + mdspanx2_t vertices1(v1.data(), v1.size() / 2, 2); auto& [coeffs, cstride] = coefficients.at({IntegralType::vertex, i}); assert(vertices.size() * cstride == coeffs.size()); + if (bs == 1) + { + impl::assemble_exterior_entities<1>( + P0, b, x_dofmap, x, vertices, {dofs, bs, vertices1}, fn, constants, + md::mdspan(coeffs.data(), vertices.extent(0), cstride), cell_info0, + facet_perms); + } + else if (bs == 3) + { + impl::assemble_exterior_entities<3>( + P0, b, x_dofmap, x, vertices, {dofs, bs, vertices1}, fn, constants, + md::mdspan(coeffs.data(), vertices.size() / 2, cstride), cell_info0, + facet_perms); + } + else + { + impl::assemble_exterior_entities( + P0, b, x_dofmap, x, vertices, {dofs, bs, vertices1}, fn, constants, + md::mdspan(coeffs.data(), vertices.size() / 2, cstride), cell_info0, + facet_perms); + } + } + + // NOTE: Empty for now, but I guess we require this for mixed-dimensional + // assembly + md::mdspan> ridge_perms; + for (int i = 0; i < L.num_integrals(IntegralType::ridge, 0); ++i) + { + auto fn = L.kernel(IntegralType::ridge, i, 0); + assert(fn); - impl::assemble_vertices( - P0, b, x_dofmap, x, - md::mdspan>( - vertices.data(), vertices.size() / 2, 2), - {dofs, bs, - md::mdspan>( - vertices0.data(), vertices0.size() / 2, 2)}, - fn, constants, - md::mdspan(coeffs.data(), vertices.size() / 2, cstride), cell_info0); + std::span r = L.domain(IntegralType::ridge, i, cell_type_idx); + mdspanx2_t ridges(r.data(), r.size() / 2, 2); + std::span r1 = L.domain_arg(IntegralType::ridge, 0, i, cell_type_idx); + mdspanx2_t ridges1(r1.data(), r1.size() / 2, 2); + + auto& [coeffs, cstride] = coefficients.at({IntegralType::ridge, i}); + + assert(ridges.size() * cstride == coeffs.size()); + + if (bs == 1) + { + impl::assemble_exterior_entities<1>( + P0, b, x_dofmap, x, ridges, {dofs, bs, ridges1}, fn, constants, + md::mdspan(coeffs.data(), ridges.extent(0), cstride), cell_info0, + facet_perms); + } + else if (bs == 3) + { + impl::assemble_exterior_entities<3>( + P0, b, x_dofmap, x, ridges, {dofs, bs, ridges1}, fn, constants, + md::mdspan(coeffs.data(), ridges.size() / 2, cstride), cell_info0, + facet_perms); + } + else + { + impl::assemble_exterior_entities( + P0, b, x_dofmap, x, ridges, {dofs, bs, ridges1}, fn, constants, + md::mdspan(coeffs.data(), ridges.size() / 2, cstride), cell_info0, + facet_perms); + } } } } diff --git a/cpp/dolfinx/fem/utils.cpp b/cpp/dolfinx/fem/utils.cpp index ec7a8501ae..c11baebff7 100644 --- a/cpp/dolfinx/fem/utils.cpp +++ b/cpp/dolfinx/fem/utils.cpp @@ -154,6 +154,9 @@ fem::compute_integration_domains(fem::IntegralType integral_type, case IntegralType::vertex: dim = 0; break; + case IntegralType::ridge: + dim = tdim - 2; + break; default: throw std::runtime_error( "Cannot compute integration domains. Integral type not supported."); @@ -253,11 +256,22 @@ fem::compute_integration_domains(fem::IntegralType integral_type, auto [v_to_c, c_to_v] = get_connectivities(0); for (auto vertex : entities) { - std::array pair = impl::get_cell_vertex_pairs<1>( + std::array pair = impl::get_cell_entity_pairs<1>( vertex, v_to_c->links(vertex), *c_to_v); entity_data.insert(entity_data.end(), pair.begin(), pair.end()); } + break; + } + case IntegralType::ridge: + { + auto [r_to_c, c_to_r] = get_connectivities(tdim - 2); + for (auto entity : entities) + { + std::array pair = impl::get_cell_entity_pairs<1>( + entity, r_to_c->links(entity), *c_to_r); + entity_data.insert(entity_data.end(), pair.begin(), pair.end()); + } } } diff --git a/cpp/dolfinx/fem/utils.h b/cpp/dolfinx/fem/utils.h index 22436c8e50..d9860b821a 100644 --- a/cpp/dolfinx/fem/utils.h +++ b/cpp/dolfinx/fem/utils.h @@ -92,17 +92,17 @@ get_cell_facet_pairs(std::int32_t f, std::span cells, return cell_local_facet_pairs; } -/// Helper function to get an array of of (cell, local_vertex) pairs -/// corresponding to a given vertex index. -/// @note If the vertex is connected to multiple cells, the first one is picked. -/// @param[in] v vertex index -/// @param[in] cells List of cells incident to the vertex -/// @param[in] c_to_v Cell to vertex connectivity -/// @return Vector of (cell, local_vertex) pairs +/// Helper function to get an array of of (cell, local_entity) pairs +/// corresponding to a given entity index. +/// @note If the entity is connected to multiple cells, the first one is picked. +/// @param[in] e entity index +/// @param[in] cells List of cells incident to the entity +/// @param[in] c_to_e Cell to entity connectivity +/// @return Vector of (cell, local_entity) pairs template std::array -get_cell_vertex_pairs(std::int32_t v, std::span cells, - const graph::AdjacencyList& c_to_v) +get_cell_entity_pairs(std::int32_t e, std::span cells, + const graph::AdjacencyList& c_to_e) { static_assert(num_cells == 1); // Patch assembly not supported. @@ -111,11 +111,11 @@ get_cell_vertex_pairs(std::int32_t v, std::span cells, // Use first cell for assembly over by default std::int32_t cell = cells[0]; - // Find local index of vertex within cell - auto cell_vertices = c_to_v.links(cell); - auto it = std::ranges::find(cell_vertices, v); - assert(it != cell_vertices.end()); - std::int32_t local_index = std::distance(cell_vertices.begin(), it); + // Find local index of entity within cell + auto cell_entities = c_to_e.links(cell); + auto it = std::ranges::find(cell_entities, e); + assert(it != cell_entities.end()); + std::int32_t local_index = std::distance(cell_entities.begin(), it); return {cell, local_index}; } @@ -490,7 +490,7 @@ Form create_form_factory( // integral offsets. Since the UFL forms for each type of cell should be // the same, I think this assumption is OK. const int* integral_offsets = ufcx_forms[0].get().form_integral_offsets; - std::array num_integrals_type; + std::array num_integrals_type; for (std::size_t i = 0; i < num_integrals_type.size(); ++i) num_integrals_type[i] = integral_offsets[i + 1] - integral_offsets[i]; @@ -510,6 +510,14 @@ Form create_form_factory( mesh->topology_mutable()->create_connectivity(tdim, tdim - 1); } + // Create ridges, if required + if (num_integrals_type[ridge] > 0) + { + mesh->topology_mutable()->create_entities(tdim - 2); + mesh->topology_mutable()->create_connectivity(tdim - 2, tdim); + mesh->topology_mutable()->create_connectivity(tdim, tdim - 2); + } + // Get list of integral IDs, and load tabulate tensor into memory for // each std::map, integral_data> integrals; @@ -800,7 +808,7 @@ Form create_form_factory( cell_and_vertex.reserve(2 * vertices_range.size()); for (std::int32_t vertex : vertices_range) { - std::array pair = impl::get_cell_vertex_pairs<1>( + std::array pair = impl::get_cell_entity_pairs<1>( vertex, v_to_c->links(vertex), *c_to_v); cell_and_vertex.insert(cell_and_vertex.end(), pair.begin(), @@ -838,6 +846,86 @@ Form create_form_factory( } } + // Attach ridge kernels + { + for (std::size_t form_idx = 0; form_idx < ufcx_forms.size(); ++form_idx) + { + const ufcx_form& form = ufcx_forms[form_idx]; + + std::span ids(form.form_integral_ids + integral_offsets[ridge], + num_integrals_type[ridge]); + auto sd = subdomains.find(IntegralType::ridge); + for (int i = 0; i < num_integrals_type[ridge]; ++i) + { + const int id = ids[i]; + ufcx_integral* integral + = form.form_integrals[integral_offsets[ridge] + i]; + assert(integral); + check_geometry_hash(*integral, form_idx); + + std::vector active_coeffs; + for (int j = 0; j < form.num_coefficients; ++j) + { + if (integral->enabled_coefficients[j]) + active_coeffs.push_back(j); + } + + impl::kernel_t k = impl::extract_kernel(integral); + assert(k); + + // Build list of entities to assembler over + auto r_to_c = topology->connectivity(tdim - 2, tdim); + assert(r_to_c); + auto c_to_r = topology->connectivity(tdim, tdim - 2); + assert(c_to_r); + + // pack for a range of vertices a flattened list of cell index c_i and + // local ridge index l_i: + // [c_0, l_0, ..., c_n, l_n] + auto get_cells_and_ridges = [r_to_c, c_to_r](auto ridge_range) + { + std::vector cell_and_ridge; + cell_and_ridge.reserve(2 * ridge_range.size()); + for (std::int32_t ridge : ridge_range) + { + std::array pair = impl::get_cell_entity_pairs<1>( + ridge, r_to_c->links(ridge), *c_to_r); + + cell_and_ridge.insert(cell_and_ridge.end(), pair.begin(), + pair.end()); + } + assert(cell_and_ridge.size() == 2 * ridge_range.size()); + return cell_and_ridge; + }; + + if (id == -1) + { + // Default ridge kernel operates on all (owned) ridges + std::int32_t num_ridges = topology->index_map(tdim - 2)->size_local(); + std::vector cells_and_ridges + = get_cells_and_ridges(std::ranges::views::iota(0, num_ridges)); + + integrals.insert({{IntegralType::ridge, i, form_idx}, + {k, cells_and_ridges, active_coeffs}}); + } + else if (sd != subdomains.end()) + { + // NOTE: This requires that pairs are sorted + auto it = std::ranges::lower_bound(sd->second, id, std::less<>{}, + [](auto& a) { return a.first; }); + if (it != sd->second.end() and it->first == id) + { + integrals.insert({{IntegralType::ridge, i, form_idx}, + {k, + std::vector(it->second.begin(), + it->second.end()), + active_coeffs}}); + } + } + } + } + } + return Form(spaces, std::move(integrals), mesh, coefficients, constants, needs_facet_permutations, entity_maps); } diff --git a/python/dolfinx/fem/forms.py b/python/dolfinx/fem/forms.py index 9967e3eadf..f54ce4617c 100644 --- a/python/dolfinx/fem/forms.py +++ b/python/dolfinx/fem/forms.py @@ -169,6 +169,10 @@ def get_integration_domains( subdomain._cpp_object.topology.create_connectivity(0, tdim) subdomain._cpp_object.topology.create_connectivity(tdim, 0) + if integral_type is IntegralType.ridge: + tdim = subdomain.topology.dim + subdomain._cpp_object.topology.create_connectivity(tdim - 2, tdim) + subdomain._cpp_object.topology.create_connectivity(tdim, tdim - 2) # Compute integration domains only for each subdomain id in # the integrals. If a process has no integral entities, # insert an empty array. @@ -219,6 +223,7 @@ def form_cpp_class( "exterior_facet": IntegralType.exterior_facet, "interior_facet": IntegralType.interior_facet, "vertex": IntegralType.vertex, + "ridge": IntegralType.ridge, } diff --git a/python/dolfinx/wrappers/fem.cpp b/python/dolfinx/wrappers/fem.cpp index 1acc71e53a..ee55f9f368 100644 --- a/python/dolfinx/wrappers/fem.cpp +++ b/python/dolfinx/wrappers/fem.cpp @@ -1273,7 +1273,8 @@ void fem(nb::module_& m) "exterior facet integral") .value("interior_facet", dolfinx::fem::IntegralType::interior_facet, "exterior facet integral") - .value("vertex", dolfinx::fem::IntegralType::vertex, "vertex integral"); + .value("vertex", dolfinx::fem::IntegralType::vertex, "vertex integral") + .value("ridge", dolfinx::fem::IntegralType::ridge, "ridge integral"); declare_function_space(m, "float32"); declare_function_space(m, "float64"); diff --git a/python/test/unit/fem/test_assembler.py b/python/test/unit/fem/test_assembler.py index 39fe8077fb..83f72d25ab 100644 --- a/python/test/unit/fem/test_assembler.py +++ b/python/test/unit/fem/test_assembler.py @@ -1866,3 +1866,69 @@ def check_vertex_integral_against_sum(form, vertices, weighted=False): assert expected_value_l == pytest.approx( value_l.array[: expected_value_l.size], abs=5e4 * np.finfo(rdtype).eps ) + + +@pytest.mark.parametrize( + "cell_type", + [ + mesh.CellType.triangle, + mesh.CellType.quadrilateral, + #mesh.CellType.tetrahedron, + # # mesh.CellType.pyramid, + # mesh.CellType.prism, + # mesh.CellType.hexahedron, + ], +) +@pytest.mark.parametrize("ghost_mode", [mesh.GhostMode.none, mesh.GhostMode.shared_facet]) +@pytest.mark.parametrize( + "dtype", + [ + np.float32, + np.float64, + pytest.param( + np.complex64, + marks=pytest.mark.skipif( + os.name == "nt", reason="win32 platform does not support C99 _Complex numbers" + ), + ), + pytest.param( + np.complex128, + marks=pytest.mark.skipif( + os.name == "nt", reason="win32 platform does not support C99 _Complex numbers" + ), + ), + ], +) +def test_ridge_integrals(cell_type, ghost_mode, dtype): + comm = MPI.COMM_WORLD + rdtype = np.real(dtype(0)).dtype + + msh = None + cell_dim = mesh.cell_dim(cell_type) + if cell_dim == 2: + msh = mesh.create_unit_square( + comm, 4, 4, cell_type=cell_type, ghost_mode=ghost_mode, dtype=rdtype + ) + elif cell_dim == 3: + msh = mesh.create_unit_cube( + comm, 4, 4, 4, cell_type=cell_type, ghost_mode=ghost_mode, dtype=rdtype + ) + else: + raise RuntimeError("Bad dimension") + gdim = msh.geometry.dim + V = dolfinx.fem.functionspace(msh, ("Lagrange", 3)) + + x = ufl.SpatialCoordinate(msh) + u = ufl.TestFunction(V) + + integrand = ufl.conj(u) * x[gdim -1] + dr = ufl.Measure("dr", domain=msh) + F = dolfinx.fem.form(integrand * dr, dtype=dtype) + b = dolfinx.fem.assemble_vector(F) + + if cell_dim == 2: + dP = ufl.Measure("dP", domain=msh) + Fp = dolfinx.fem.form(integrand * dP, dtype=dtype) + b_p = dolfinx.fem.assemble_vector(Fp) + tol = np.finfo(rdtype).eps + np.testing.assert_allclose(b.array, b_p.array, atol=tol, rtol=tol) \ No newline at end of file From 5167f239e68a5322cf88d412eb8fe7b732c15eb9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= Date: Wed, 3 Sep 2025 09:50:29 +0000 Subject: [PATCH 02/37] Assemble scalar added --- cpp/dolfinx/fem/assemble_scalar_impl.h | 98 ++++++++++++-------------- cpp/dolfinx/fem/assemble_vector_impl.h | 3 +- python/dolfinx/fem/forms.py | 3 +- python/test/unit/fem/test_assembler.py | 85 ++++++++++++++++++++-- 4 files changed, 128 insertions(+), 61 deletions(-) diff --git a/cpp/dolfinx/fem/assemble_scalar_impl.h b/cpp/dolfinx/fem/assemble_scalar_impl.h index 3f060e7853..0887179430 100644 --- a/cpp/dolfinx/fem/assemble_scalar_impl.h +++ b/cpp/dolfinx/fem/assemble_scalar_impl.h @@ -57,32 +57,32 @@ T assemble_cells(mdspan2_t x_dofmap, return value; } -/// Execute kernel over exterior facets and accumulate result +/// Execute kernel over exterior entities and accumulate result template -T assemble_exterior_facets( +T assemble_exterior_entities( mdspan2_t x_dofmap, md::mdspan, md::extents> x, md::mdspan> - facets, + entities, FEkernel auto fn, std::span constants, md::mdspan> coeffs, md::mdspan> perms) { T value(0); - if (facets.empty()) + if (entities.empty()) return value; // Create data structures used in assembly std::vector> cdofs(3 * x_dofmap.extent(1)); // Iterate over all facets - for (std::size_t f = 0; f < facets.extent(0); ++f) + for (std::size_t f = 0; f < entities.extent(0); ++f) { - std::int32_t cell = facets(f, 0); - std::int32_t local_facet = facets(f, 1); + std::int32_t cell = entities(f, 0); + std::int32_t local_entity = entities(f, 1); // Get cell coordinates/geometry auto x_dofs = md::submdspan(x_dofmap, cell, md::full_extent); @@ -90,8 +90,8 @@ T assemble_exterior_facets( std::copy_n(&x(x_dofs[i], 0), 3, std::next(cdofs.begin(), 3 * i)); // Permutations - std::uint8_t perm = perms.empty() ? 0 : perms(cell, local_facet); - fn(&value, &coeffs(f, 0), constants.data(), cdofs.data(), &local_facet, + std::uint8_t perm = perms.empty() ? 0 : perms(cell, local_entity); + fn(&value, &coeffs(f, 0), constants.data(), cdofs.data(), &local_entity, &perm, nullptr); } @@ -150,43 +150,6 @@ T assemble_interior_facets( return value; } -/// Assemble functional over vertices -template -T assemble_vertices(mdspan2_t x_dofmap, - md::mdspan, - md::extents> - x, - md::mdspan> - vertices, - FEkernel auto fn, std::span constants, - md::mdspan> coeffs) -{ - T value(0); - if (vertices.empty()) - return value; - - // Create data structures used in assembly - std::vector> cdofs(3 * x_dofmap.extent(1)); - - // Iterate over all cells - for (std::size_t index = 0; index < vertices.extent(0); ++index) - { - std::int32_t cell = vertices(index, 0); - std::int32_t local_vertex_index = vertices(index, 1); - - // Get cell coordinates/geometry - auto x_dofs = md::submdspan(x_dofmap, cell, md::full_extent); - for (std::size_t i = 0; i < x_dofs.size(); ++i) - std::copy_n(&x(x_dofs[i], 0), 3, std::next(cdofs.begin(), 3 * i)); - - fn(&value, &coeffs(index, 0), constants.data(), cdofs.data(), - &local_vertex_index, nullptr, nullptr); - } - - return value; -} - /// Assemble functional into an scalar with provided mesh geometry. template T assemble_scalar( @@ -217,14 +180,14 @@ T assemble_scalar( mesh::CellType cell_type = mesh->topology()->cell_type(); int num_facets_per_cell = mesh::cell_num_entities(cell_type, mesh->topology()->dim() - 1); - md::mdspan> perms; + md::mdspan> facet_perms; if (M.needs_facet_permutations()) { mesh->topology_mutable()->create_entity_permutations(); const std::vector& p = mesh->topology()->get_facet_permutations(); - perms = md::mdspan(p.data(), p.size() / num_facets_per_cell, - num_facets_per_cell); + facet_perms = md::mdspan(p.data(), p.size() / num_facets_per_cell, + num_facets_per_cell); } for (int i = 0; i < M.num_integrals(IntegralType::exterior_facet, 0); ++i) @@ -241,13 +204,14 @@ T assemble_scalar( constexpr std::size_t shape1 = 2 * num_adjacent_cells; assert((facets.size() / 2) * cstride == coeffs.size()); - value += impl::assemble_exterior_facets( + value += impl::assemble_exterior_entities( x_dofmap, x, md::mdspan>( facets.data(), facets.size() / shape1, 2), fn, constants, - md::mdspan(coeffs.data(), facets.size() / shape1, cstride), perms); + md::mdspan(coeffs.data(), facets.size() / shape1, cstride), + facet_perms); } for (int i = 0; i < M.num_integrals(IntegralType::interior_facet, 0); ++i) @@ -272,9 +236,10 @@ T assemble_scalar( md::mdspan>( coeffs.data(), facets.size() / shape1, 2, cstride), - perms); + facet_perms); } + md::mdspan> vertex_perms; for (int i = 0; i < M.num_integrals(IntegralType::vertex, 0); ++i) { auto fn = M.kernel(IntegralType::vertex, i, 0); @@ -290,13 +255,38 @@ T assemble_scalar( // Two values per adj. cell (cell index and local vertex index). constexpr std::size_t shape1 = 2 * num_adjacent_cells; - value += impl::assemble_vertices( + value += impl::assemble_exterior_entities( x_dofmap, x, md::mdspan>( vertices.data(), vertices.size() / shape1, shape1), fn, constants, - md::mdspan(coeffs.data(), vertices.size() / shape1, cstride)); + md::mdspan(coeffs.data(), vertices.size() / shape1, cstride), + vertex_perms); + } + + md::mdspan> ridge_perms; + for (int i = 0; i < M.num_integrals(IntegralType::ridge, 0); ++i) + { + auto fn = M.kernel(IntegralType::ridge, i, 0); + assert(fn); + + auto& [coeffs, cstride] = coefficients.at({IntegralType::ridge, i}); + + std::span ridges = M.domain(IntegralType::ridge, i, 0); + assert(ridges.size() * cstride == coeffs.size()); + + constexpr std::size_t num_adjacent_cells = 1; + // Two values per adj. cell (cell index and local ridge index). + constexpr std::size_t shape1 = 2 * num_adjacent_cells; + value += impl::assemble_exterior_entities( + x_dofmap, x, + md::mdspan>( + ridges.data(), ridges.size() / shape1, shape1), + fn, constants, + md::mdspan(coeffs.data(), ridges.size() / shape1, cstride), + ridge_perms); } return value; diff --git a/cpp/dolfinx/fem/assemble_vector_impl.h b/cpp/dolfinx/fem/assemble_vector_impl.h index 46b449102b..2d1ccbe811 100644 --- a/cpp/dolfinx/fem/assemble_vector_impl.h +++ b/cpp/dolfinx/fem/assemble_vector_impl.h @@ -718,7 +718,8 @@ void assemble_cells( } } -/// @brief Execute kernel over cells and accumulate result in vector. +/// @brief Execute kernel over exterior entities and accumulate result in +/// vector. /// /// @tparam T Scalar type. /// @tparam _bs The block size of the form test function dof map. If diff --git a/python/dolfinx/fem/forms.py b/python/dolfinx/fem/forms.py index f54ce4617c..1dc913e5f4 100644 --- a/python/dolfinx/fem/forms.py +++ b/python/dolfinx/fem/forms.py @@ -173,6 +173,7 @@ def get_integration_domains( tdim = subdomain.topology.dim subdomain._cpp_object.topology.create_connectivity(tdim - 2, tdim) subdomain._cpp_object.topology.create_connectivity(tdim, tdim - 2) + # Compute integration domains only for each subdomain id in # the integrals. If a process has no integral entities, # insert an empty array. @@ -374,7 +375,7 @@ def _form(form): # Extract subdomain ids from ufcx_form subdomain_ids = {type: [] for type in sd.get(domain).keys()} - integral_offsets = [ufcx_form.form_integral_offsets[i] for i in range(5)] + integral_offsets = [ufcx_form.form_integral_offsets[i] for i in range(6)] for i in range(len(integral_offsets) - 1): integral_type = IntegralType(i) for j in range(integral_offsets[i], integral_offsets[i + 1]): diff --git a/python/test/unit/fem/test_assembler.py b/python/test/unit/fem/test_assembler.py index 83f72d25ab..72c1447d5e 100644 --- a/python/test/unit/fem/test_assembler.py +++ b/python/test/unit/fem/test_assembler.py @@ -1873,8 +1873,8 @@ def check_vertex_integral_against_sum(form, vertices, weighted=False): [ mesh.CellType.triangle, mesh.CellType.quadrilateral, - #mesh.CellType.tetrahedron, - # # mesh.CellType.pyramid, + # mesh.CellType.tetrahedron, + # mesh.CellType.pyramid, # mesh.CellType.prism, # mesh.CellType.hexahedron, ], @@ -1899,7 +1899,7 @@ def check_vertex_integral_against_sum(form, vertices, weighted=False): ), ], ) -def test_ridge_integrals(cell_type, ghost_mode, dtype): +def test_ridge_integrals_rank1(cell_type, ghost_mode, dtype): comm = MPI.COMM_WORLD rdtype = np.real(dtype(0)).dtype @@ -1921,7 +1921,7 @@ def test_ridge_integrals(cell_type, ghost_mode, dtype): x = ufl.SpatialCoordinate(msh) u = ufl.TestFunction(V) - integrand = ufl.conj(u) * x[gdim -1] + integrand = ufl.conj(u) * x[gdim - 1] dr = ufl.Measure("dr", domain=msh) F = dolfinx.fem.form(integrand * dr, dtype=dtype) b = dolfinx.fem.assemble_vector(F) @@ -1931,4 +1931,79 @@ def test_ridge_integrals(cell_type, ghost_mode, dtype): Fp = dolfinx.fem.form(integrand * dP, dtype=dtype) b_p = dolfinx.fem.assemble_vector(Fp) tol = np.finfo(rdtype).eps - np.testing.assert_allclose(b.array, b_p.array, atol=tol, rtol=tol) \ No newline at end of file + np.testing.assert_allclose(b.array, b_p.array, atol=tol, rtol=tol) + + +@pytest.mark.parametrize( + "cell_type", + [ + mesh.CellType.triangle, + mesh.CellType.quadrilateral, + mesh.CellType.tetrahedron, + # mesh.CellType.pyramid, + # mesh.CellType.prism, + mesh.CellType.hexahedron, + ], +) +@pytest.mark.parametrize("ghost_mode", [mesh.GhostMode.none, mesh.GhostMode.shared_facet]) +@pytest.mark.parametrize( + "dtype", + [ + np.float32, + np.float64, + pytest.param( + np.complex64, + marks=pytest.mark.skipif( + os.name == "nt", reason="win32 platform does not support C99 _Complex numbers" + ), + ), + pytest.param( + np.complex128, + marks=pytest.mark.skipif( + os.name == "nt", reason="win32 platform does not support C99 _Complex numbers" + ), + ), + ], +) +def test_ridge_integrals_rank0(cell_type, ghost_mode, dtype): + comm = MPI.COMM_WORLD + rdtype = np.real(dtype(0)).dtype + + msh = None + cell_dim = mesh.cell_dim(cell_type) + if cell_dim == 2: + msh = mesh.create_unit_square( + comm, 4, 4, cell_type=cell_type, ghost_mode=ghost_mode, dtype=rdtype + ) + elif cell_dim == 3: + msh = mesh.create_unit_cube( + comm, 4, 4, 4, cell_type=cell_type, ghost_mode=ghost_mode, dtype=rdtype + ) + else: + raise RuntimeError("Bad dimension") + gdim = msh.geometry.dim + + x = ufl.SpatialCoordinate(msh) + + def marked_ridges(x): + return np.isclose(x[0], 1) & np.isclose(x[1], 0.5) + + exterior_ridges = dolfinx.mesh.locate_entities_boundary( + msh, msh.topology.dim - 2, marked_ridges + ) + et = dolfinx.mesh.meshtags( + msh, msh.topology.dim - 2, exterior_ridges, np.full_like(exterior_ridges, 33) + ) + dr = ufl.Measure("dr", domain=msh, subdomain_data=et, subdomain_id=33) + integrand = x[0] + x[1] + x[gdim - 1] ** 2 + J_compiled = dolfinx.fem.form(integrand * dr, dtype=dtype) + J = dolfinx.fem.assemble_scalar(J_compiled) + J_sum = msh.comm.allreduce(J, op=MPI.SUM) + if cell_dim == 2: + ref_sol = 1 + 1 / 2 + (1 / 2) ** 2 + elif cell_dim == 3: + ref_sol = 1 + 1 / 2 + 1 / 3 + else: + raise RuntimeError("Bad dimension") + tol = 50 * np.finfo(rdtype).eps + assert np.isclose(J_sum, ref_sol, atol=tol, rtol=tol) From 804523bcf5f29b14778e141371efe19da621e4df Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= Date: Wed, 3 Sep 2025 09:51:19 +0000 Subject: [PATCH 03/37] Use vertex_perms --- cpp/dolfinx/fem/assemble_vector_impl.h | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/cpp/dolfinx/fem/assemble_vector_impl.h b/cpp/dolfinx/fem/assemble_vector_impl.h index 2d1ccbe811..d5eac078e1 100644 --- a/cpp/dolfinx/fem/assemble_vector_impl.h +++ b/cpp/dolfinx/fem/assemble_vector_impl.h @@ -1402,21 +1402,21 @@ void assemble_vector( impl::assemble_exterior_entities<1>( P0, b, x_dofmap, x, vertices, {dofs, bs, vertices1}, fn, constants, md::mdspan(coeffs.data(), vertices.extent(0), cstride), cell_info0, - facet_perms); + vertex_perms); } else if (bs == 3) { impl::assemble_exterior_entities<3>( P0, b, x_dofmap, x, vertices, {dofs, bs, vertices1}, fn, constants, md::mdspan(coeffs.data(), vertices.size() / 2, cstride), cell_info0, - facet_perms); + vertex_perms); } else { impl::assemble_exterior_entities( P0, b, x_dofmap, x, vertices, {dofs, bs, vertices1}, fn, constants, md::mdspan(coeffs.data(), vertices.size() / 2, cstride), cell_info0, - facet_perms); + vertex_perms); } } From a8e0738985c981a67eabd9bfa1f7ed05ba6d0520 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= Date: Wed, 3 Sep 2025 10:01:32 +0000 Subject: [PATCH 04/37] Fix docs --- cpp/dolfinx/fem/assemble_vector_impl.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cpp/dolfinx/fem/assemble_vector_impl.h b/cpp/dolfinx/fem/assemble_vector_impl.h index d5eac078e1..2b2ce04c1d 100644 --- a/cpp/dolfinx/fem/assemble_vector_impl.h +++ b/cpp/dolfinx/fem/assemble_vector_impl.h @@ -731,7 +731,7 @@ void assemble_cells( /// @param[in,out] b The vector to accumulate into. /// @param[in] x_dofmap Dofmap for the mesh geometry. /// @param[in] x Mesh geometry (coordinates). -/// @param[in] facets Facets (in the integration domain mesh) to execute +/// @param[in] entities Entities (in the integration domain mesh) to execute /// the kernel over. /// @param[in] dofmap Test function (row) degree-of-freedom data holding /// the (0) dofmap, (1) dofmap block size and (2) dofmap cell indices. From 992b7ddbcd580d78c7fa15f3bb521f7cb3b34621 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= Date: Wed, 3 Sep 2025 10:29:51 +0000 Subject: [PATCH 05/37] Fix correct usage of perms --- cpp/dolfinx/fem/assemble_vector_impl.h | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/cpp/dolfinx/fem/assemble_vector_impl.h b/cpp/dolfinx/fem/assemble_vector_impl.h index 2b2ce04c1d..22d498a56c 100644 --- a/cpp/dolfinx/fem/assemble_vector_impl.h +++ b/cpp/dolfinx/fem/assemble_vector_impl.h @@ -1442,21 +1442,21 @@ void assemble_vector( impl::assemble_exterior_entities<1>( P0, b, x_dofmap, x, ridges, {dofs, bs, ridges1}, fn, constants, md::mdspan(coeffs.data(), ridges.extent(0), cstride), cell_info0, - facet_perms); + ridge_perms); } else if (bs == 3) { impl::assemble_exterior_entities<3>( P0, b, x_dofmap, x, ridges, {dofs, bs, ridges1}, fn, constants, md::mdspan(coeffs.data(), ridges.size() / 2, cstride), cell_info0, - facet_perms); + ridge_perms); } else { impl::assemble_exterior_entities( P0, b, x_dofmap, x, ridges, {dofs, bs, ridges1}, fn, constants, md::mdspan(coeffs.data(), ridges.size() / 2, cstride), cell_info0, - facet_perms); + ridge_perms); } } } From 0d24551424f2f05d40c9c67f7a4cc0066616ddc6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= Date: Wed, 3 Sep 2025 10:51:36 +0000 Subject: [PATCH 06/37] add 3D vector test --- python/test/unit/fem/test_assembler.py | 140 +++++++++++++++++++++---- 1 file changed, 118 insertions(+), 22 deletions(-) diff --git a/python/test/unit/fem/test_assembler.py b/python/test/unit/fem/test_assembler.py index 72c1447d5e..5d7e5cf959 100644 --- a/python/test/unit/fem/test_assembler.py +++ b/python/test/unit/fem/test_assembler.py @@ -1873,10 +1873,6 @@ def check_vertex_integral_against_sum(form, vertices, weighted=False): [ mesh.CellType.triangle, mesh.CellType.quadrilateral, - # mesh.CellType.tetrahedron, - # mesh.CellType.pyramid, - # mesh.CellType.prism, - # mesh.CellType.hexahedron, ], ) @pytest.mark.parametrize("ghost_mode", [mesh.GhostMode.none, mesh.GhostMode.shared_facet]) @@ -1899,22 +1895,14 @@ def check_vertex_integral_against_sum(form, vertices, weighted=False): ), ], ) -def test_ridge_integrals_rank1(cell_type, ghost_mode, dtype): +def test_ridge_integrals_rank1_2D(cell_type, ghost_mode, dtype): comm = MPI.COMM_WORLD rdtype = np.real(dtype(0)).dtype msh = None - cell_dim = mesh.cell_dim(cell_type) - if cell_dim == 2: - msh = mesh.create_unit_square( - comm, 4, 4, cell_type=cell_type, ghost_mode=ghost_mode, dtype=rdtype - ) - elif cell_dim == 3: - msh = mesh.create_unit_cube( - comm, 4, 4, 4, cell_type=cell_type, ghost_mode=ghost_mode, dtype=rdtype - ) - else: - raise RuntimeError("Bad dimension") + msh = mesh.create_unit_square( + comm, 4, 4, cell_type=cell_type, ghost_mode=ghost_mode, dtype=rdtype + ) gdim = msh.geometry.dim V = dolfinx.fem.functionspace(msh, ("Lagrange", 3)) @@ -1926,12 +1914,11 @@ def test_ridge_integrals_rank1(cell_type, ghost_mode, dtype): F = dolfinx.fem.form(integrand * dr, dtype=dtype) b = dolfinx.fem.assemble_vector(F) - if cell_dim == 2: - dP = ufl.Measure("dP", domain=msh) - Fp = dolfinx.fem.form(integrand * dP, dtype=dtype) - b_p = dolfinx.fem.assemble_vector(Fp) - tol = np.finfo(rdtype).eps - np.testing.assert_allclose(b.array, b_p.array, atol=tol, rtol=tol) + dP = ufl.Measure("dP", domain=msh) + Fp = dolfinx.fem.form(integrand * dP, dtype=dtype) + b_p = dolfinx.fem.assemble_vector(Fp) + tol = np.finfo(rdtype).eps + np.testing.assert_allclose(b.array, b_p.array, atol=tol, rtol=tol) @pytest.mark.parametrize( @@ -2007,3 +1994,112 @@ def marked_ridges(x): raise RuntimeError("Bad dimension") tol = 50 * np.finfo(rdtype).eps assert np.isclose(J_sum, ref_sol, atol=tol, rtol=tol) + + +@pytest.mark.parametrize( + "cell_type", + [ + mesh.CellType.tetrahedron, + # mesh.CellType.pyramid, + # mesh.CellType.prism, + mesh.CellType.hexahedron, + ], +) +@pytest.mark.parametrize("ghost_mode", [mesh.GhostMode.none, mesh.GhostMode.shared_facet]) +@pytest.mark.parametrize( + "dtype", + [ + np.float32, + np.float64, + pytest.param( + np.complex64, + marks=pytest.mark.skipif( + os.name == "nt", reason="win32 platform does not support C99 _Complex numbers" + ), + ), + pytest.param( + np.complex128, + marks=pytest.mark.skipif( + os.name == "nt", reason="win32 platform does not support C99 _Complex numbers" + ), + ), + ], +) +def test_ridge_integrals_rank1_3D(cell_type, ghost_mode, dtype): + comm = MPI.COMM_WORLD + rdtype = np.real(dtype(0)).dtype + + msh = None + N = 4 + msh = mesh.create_unit_cube( + comm, N, N, N, cell_type=cell_type, ghost_mode=ghost_mode, dtype=rdtype + ) + gdim = msh.geometry.dim + + el = ("Lagrange", 3, (gdim,)) + V = dolfinx.fem.functionspace(msh, el) + + x = ufl.SpatialCoordinate(msh) + v = ufl.TestFunction(V) + + def marked_ridges(x): + return np.isclose(x[0], 1) & np.isclose(x[1], 0.5) + + exterior_ridges = dolfinx.mesh.locate_entities_boundary( + msh, msh.topology.dim - 2, marked_ridges + ) + et = dolfinx.mesh.meshtags( + msh, msh.topology.dim - 2, exterior_ridges, np.full_like(exterior_ridges, 33) + ) + + def f(x): + return x[0] + ufl.sin(x[1]), x[2] + x[1], x[0] ** 2 - 3 * x[1] + + def integrand(x, v): + expr = ufl.as_vector(f(x)) + return ufl.inner(expr, v) + + metadata = {"quadrature_degree": 10} + dr = ufl.Measure("dr", domain=msh, subdomain_data=et, subdomain_id=33, metadata=metadata) + F = dolfinx.fem.form(integrand(x, v) * dr, dtype=dtype) + b = dolfinx.fem.assemble_vector(F) + b.scatter_reverse(la.InsertMode.add) + b.scatter_forward() + + # Create reference solution on unit interval + assert gdim == 3 + if comm.rank == 0: + nodes = np.zeros((N + 1, gdim), dtype=rdtype) + nodes[:, 0] = 1 + nodes[:, 1] = 0.5 + nodes[:, 2] = np.linspace(0, 1, N + 1) + connectivity = ( + np.repeat(np.arange(nodes.shape[0]), 2)[1:-1] + .reshape(nodes.shape[0] - 1, 2) + .astype(np.int64) + ) + else: + nodes = np.zeros((0, gdim), dtype=rdtype) + connectivity = np.zeros((0, 2), dtype=np.int64) + + c_el = ufl.Mesh( + basix.ufl.element("Lagrange", basix.CellType.interval, 1, shape=(nodes.shape[1],)) + ) + line_mesh = dolfinx.mesh.create_mesh( + MPI.COMM_WORLD, + x=nodes, + cells=connectivity, + e=c_el, + partitioner=dolfinx.mesh.create_cell_partitioner(ghost_mode), + ) + + Q = dolfinx.fem.functionspace(line_mesh, el) + q = ufl.TestFunction(Q) + + x_e = ufl.SpatialCoordinate(line_mesh) + F_ref = dolfinx.fem.form(integrand(x_e, q) * ufl.dx(domain=line_mesh, metadata=metadata)) + b_ref = dolfinx.fem.assemble_vector(F_ref) + b_ref.scatter_reverse(la.InsertMode.add) + b_ref.scatter_forward() + tol = 10 * np.finfo(rdtype).eps + assert np.isclose(la.norm(b), la.norm(b_ref), rtol=tol, atol=tol) From 848ec5be9113091a059f97e7be25c9d63eecf684 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= Date: Wed, 3 Sep 2025 12:37:26 +0000 Subject: [PATCH 07/37] Remove code duplication --- cpp/dolfinx/fem/assemble_scalar_impl.h | 99 ++++++++------------------ 1 file changed, 29 insertions(+), 70 deletions(-) diff --git a/cpp/dolfinx/fem/assemble_scalar_impl.h b/cpp/dolfinx/fem/assemble_scalar_impl.h index 0887179430..eee96b3618 100644 --- a/cpp/dolfinx/fem/assemble_scalar_impl.h +++ b/cpp/dolfinx/fem/assemble_scalar_impl.h @@ -190,30 +190,6 @@ T assemble_scalar( num_facets_per_cell); } - for (int i = 0; i < M.num_integrals(IntegralType::exterior_facet, 0); ++i) - { - auto fn = M.kernel(IntegralType::exterior_facet, i, 0); - assert(fn); - auto& [coeffs, cstride] - = coefficients.at({IntegralType::exterior_facet, i}); - - std::span facets = M.domain(IntegralType::exterior_facet, i, 0); - - constexpr std::size_t num_adjacent_cells = 1; - // Two values per each adj. cell (cell index and local facet index). - constexpr std::size_t shape1 = 2 * num_adjacent_cells; - - assert((facets.size() / 2) * cstride == coeffs.size()); - value += impl::assemble_exterior_entities( - x_dofmap, x, - md::mdspan>( - facets.data(), facets.size() / shape1, 2), - fn, constants, - md::mdspan(coeffs.data(), facets.size() / shape1, cstride), - facet_perms); - } - for (int i = 0; i < M.num_integrals(IntegralType::interior_facet, 0); ++i) { auto fn = M.kernel(IntegralType::interior_facet, i, 0); @@ -239,54 +215,37 @@ T assemble_scalar( facet_perms); } - md::mdspan> vertex_perms; - for (int i = 0; i < M.num_integrals(IntegralType::vertex, 0); ++i) + for (fem::IntegralType itg_type : + {fem::IntegralType::exterior_facet, fem::IntegralType::vertex, + fem::IntegralType::ridge}) { - auto fn = M.kernel(IntegralType::vertex, i, 0); - assert(fn); - - auto& [coeffs, cstride] = coefficients.at({IntegralType::vertex, i}); + md::mdspan> perms; + if (itg_type == fem::IntegralType::exterior_facet && !facet_perms.empty()) + perms = facet_perms; + else + perms = md::mdspan>(); - std::span vertices - = M.domain(IntegralType::vertex, i, 0); - assert(vertices.size() * cstride == coeffs.size()); - - constexpr std::size_t num_adjacent_cells = 1; - // Two values per adj. cell (cell index and local vertex index). - constexpr std::size_t shape1 = 2 * num_adjacent_cells; - - value += impl::assemble_exterior_entities( - x_dofmap, x, - md::mdspan>( - vertices.data(), vertices.size() / shape1, shape1), - fn, constants, - md::mdspan(coeffs.data(), vertices.size() / shape1, cstride), - vertex_perms); - } - - md::mdspan> ridge_perms; - for (int i = 0; i < M.num_integrals(IntegralType::ridge, 0); ++i) - { - auto fn = M.kernel(IntegralType::ridge, i, 0); - assert(fn); - - auto& [coeffs, cstride] = coefficients.at({IntegralType::ridge, i}); - - std::span ridges = M.domain(IntegralType::ridge, i, 0); - assert(ridges.size() * cstride == coeffs.size()); - - constexpr std::size_t num_adjacent_cells = 1; - // Two values per adj. cell (cell index and local ridge index). - constexpr std::size_t shape1 = 2 * num_adjacent_cells; - value += impl::assemble_exterior_entities( - x_dofmap, x, - md::mdspan>( - ridges.data(), ridges.size() / shape1, shape1), - fn, constants, - md::mdspan(coeffs.data(), ridges.size() / shape1, cstride), - ridge_perms); + for (int i = 0; i < M.num_integrals(itg_type, 0); ++i) + { + auto fn = M.kernel(itg_type, i, 0); + assert(fn); + auto& [coeffs, cstride] = coefficients.at({itg_type, i}); + + std::span entities = M.domain(itg_type, i, 0); + + constexpr std::size_t num_adjacent_cells = 1; + // Two values per each adj. cell (cell index and local entity index). + constexpr std::size_t shape1 = 2 * num_adjacent_cells; + + assert((entities.size() / 2) * cstride == coeffs.size()); + value += impl::assemble_exterior_entities( + x_dofmap, x, + md::mdspan>( + entities.data(), entities.size() / shape1, 2), + fn, constants, + md::mdspan(coeffs.data(), entities.size() / shape1, cstride), perms); + } } return value; From 433b89e8deef60eca589cce3bda4953ea472b978 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= Date: Wed, 3 Sep 2025 13:32:38 +0000 Subject: [PATCH 08/37] Huge duplication removal --- cpp/dolfinx/fem/assemble_vector_impl.h | 144 +++++++------------------ cpp/dolfinx/fem/pack.h | 39 ++----- python/test/unit/fem/test_assembler.py | 34 ++++-- 3 files changed, 71 insertions(+), 146 deletions(-) diff --git a/cpp/dolfinx/fem/assemble_vector_impl.h b/cpp/dolfinx/fem/assemble_vector_impl.h index 22d498a56c..52d1830893 100644 --- a/cpp/dolfinx/fem/assemble_vector_impl.h +++ b/cpp/dolfinx/fem/assemble_vector_impl.h @@ -1298,40 +1298,6 @@ void assemble_vector( = md::mdspan>; - for (int i = 0; i < L.num_integrals(IntegralType::exterior_facet, 0); ++i) - { - auto fn = L.kernel(IntegralType::exterior_facet, i, 0); - assert(fn); - auto& [coeffs, cstride] - = coefficients.at({IntegralType::exterior_facet, i}); - std::span f = L.domain(IntegralType::exterior_facet, i, 0); - mdspanx2_t facets(f.data(), f.size() / 2, 2); - std::span f1 = L.domain_arg(IntegralType::exterior_facet, 0, i, 0); - mdspanx2_t facets1(f1.data(), f1.size() / 2, 2); - assert((facets.size() / 2) * cstride == coeffs.size()); - if (bs == 1) - { - impl::assemble_exterior_entities<1>( - P0, b, x_dofmap, x, facets, {dofs, bs, facets1}, fn, constants, - md::mdspan(coeffs.data(), facets.extent(0), cstride), cell_info0, - facet_perms); - } - else if (bs == 3) - { - impl::assemble_exterior_entities<3>( - P0, b, x_dofmap, x, facets, {dofs, bs, facets1}, fn, constants, - md::mdspan(coeffs.data(), facets.size() / 2, cstride), cell_info0, - facet_perms); - } - else - { - impl::assemble_exterior_entities( - P0, b, x_dofmap, x, facets, {dofs, bs, facets1}, fn, constants, - md::mdspan(coeffs.data(), facets.size() / 2, cstride), cell_info0, - facet_perms); - } - } - for (int i = 0; i < L.num_integrals(IntegralType::interior_facet, 0); ++i) { using mdspanx22_t @@ -1383,80 +1349,48 @@ void assemble_vector( } } - md::mdspan> vertex_perms; - for (int i = 0; i < L.num_integrals(IntegralType::vertex, 0); ++i) + for (fem::IntegralType itg_type : + {fem::IntegralType::exterior_facet, fem::IntegralType::vertex, + fem::IntegralType::ridge}) { - auto fn = L.kernel(IntegralType::vertex, i, 0); - assert(fn); - - std::span v = L.domain(IntegralType::vertex, i, cell_type_idx); - mdspanx2_t vertices(v.data(), v.size() / 2, 2); - std::span v1 = L.domain_arg(IntegralType::vertex, 0, i, cell_type_idx); - mdspanx2_t vertices1(v1.data(), v1.size() / 2, 2); - - auto& [coeffs, cstride] = coefficients.at({IntegralType::vertex, i}); - - assert(vertices.size() * cstride == coeffs.size()); - if (bs == 1) - { - impl::assemble_exterior_entities<1>( - P0, b, x_dofmap, x, vertices, {dofs, bs, vertices1}, fn, constants, - md::mdspan(coeffs.data(), vertices.extent(0), cstride), cell_info0, - vertex_perms); - } - else if (bs == 3) - { - impl::assemble_exterior_entities<3>( - P0, b, x_dofmap, x, vertices, {dofs, bs, vertices1}, fn, constants, - md::mdspan(coeffs.data(), vertices.size() / 2, cstride), cell_info0, - vertex_perms); - } + md::mdspan> perms; + if (itg_type == fem::IntegralType::exterior_facet && !facet_perms.empty()) + perms = facet_perms; else + perms = md::mdspan>(); + for (int i = 0; i < L.num_integrals(itg_type, 0); ++i) { - impl::assemble_exterior_entities( - P0, b, x_dofmap, x, vertices, {dofs, bs, vertices1}, fn, constants, - md::mdspan(coeffs.data(), vertices.size() / 2, cstride), cell_info0, - vertex_perms); - } - } - - // NOTE: Empty for now, but I guess we require this for mixed-dimensional - // assembly - md::mdspan> ridge_perms; - for (int i = 0; i < L.num_integrals(IntegralType::ridge, 0); ++i) - { - auto fn = L.kernel(IntegralType::ridge, i, 0); - assert(fn); - - std::span r = L.domain(IntegralType::ridge, i, cell_type_idx); - mdspanx2_t ridges(r.data(), r.size() / 2, 2); - std::span r1 = L.domain_arg(IntegralType::ridge, 0, i, cell_type_idx); - mdspanx2_t ridges1(r1.data(), r1.size() / 2, 2); - - auto& [coeffs, cstride] = coefficients.at({IntegralType::ridge, i}); - - assert(ridges.size() * cstride == coeffs.size()); - - if (bs == 1) - { - impl::assemble_exterior_entities<1>( - P0, b, x_dofmap, x, ridges, {dofs, bs, ridges1}, fn, constants, - md::mdspan(coeffs.data(), ridges.extent(0), cstride), cell_info0, - ridge_perms); - } - else if (bs == 3) - { - impl::assemble_exterior_entities<3>( - P0, b, x_dofmap, x, ridges, {dofs, bs, ridges1}, fn, constants, - md::mdspan(coeffs.data(), ridges.size() / 2, cstride), cell_info0, - ridge_perms); - } - else - { - impl::assemble_exterior_entities( - P0, b, x_dofmap, x, ridges, {dofs, bs, ridges1}, fn, constants, - md::mdspan(coeffs.data(), ridges.size() / 2, cstride), cell_info0, - ridge_perms); + auto fn = L.kernel(itg_type, i, 0); + assert(fn); + auto& [coeffs, cstride] = coefficients.at({itg_type, i}); + std::span e = L.domain(itg_type, i, 0); + mdspanx2_t entities(e.data(), e.size() / 2, 2); + std::span e1 = L.domain_arg(itg_type, 0, i, 0); + mdspanx2_t entities1(e1.data(), e1.size() / 2, 2); + assert((entities.size() / 2) * cstride == coeffs.size()); + if (bs == 1) + { + impl::assemble_exterior_entities<1>( + P0, b, x_dofmap, x, entities, {dofs, bs, entities1}, fn, + constants, md::mdspan(coeffs.data(), entities.extent(0), cstride), + cell_info0, perms); + } + else if (bs == 3) + { + impl::assemble_exterior_entities<3>( + P0, b, x_dofmap, x, entities, {dofs, bs, entities1}, fn, + constants, + md::mdspan(coeffs.data(), entities.size() / 2, cstride), + cell_info0, perms); + } + else + { + impl::assemble_exterior_entities( + P0, b, x_dofmap, x, entities, {dofs, bs, entities1}, fn, + constants, + md::mdspan(coeffs.data(), entities.size() / 2, cstride), + cell_info0, perms); + } } } } diff --git a/cpp/dolfinx/fem/pack.h b/cpp/dolfinx/fem/pack.h index 9b08e845d8..65777820ee 100644 --- a/cpp/dolfinx/fem/pack.h +++ b/cpp/dolfinx/fem/pack.h @@ -179,11 +179,8 @@ allocate_coefficient_storage(const Form& form, IntegralType integral_type, const std::vector offsets = form.coefficient_offsets(); cstride = offsets.back(); num_entities = form.domain(integral_type, id, 0).size(); - if (integral_type == IntegralType::exterior_facet - or integral_type == IntegralType::interior_facet) - { + if (integral_type != IntegralType::cell) num_entities /= 2; - } } return {std::vector(num_entities * cstride), cstride}; @@ -270,28 +267,6 @@ void pack_coefficients(const Form& form, } break; } - case IntegralType::exterior_facet: - { - // Iterate over coefficients coefficients that are active in - // exterior facet integrals - for (int coeff : form.active_coeffs(IntegralType::exterior_facet, id)) - { - auto mesh = coefficients[coeff]->function_space()->mesh(); - std::span facets_b - = form.domain_coeff(IntegralType::exterior_facet, id, coeff); - md::mdspan> - facets(facets_b.data(), facets_b.size() / 2, 2); - auto cells = md::submdspan(facets, md::full_extent, 0); - - std::span cell_info - = impl::get_cell_orientation_info(*coefficients[coeff]); - impl::pack_coefficient_entity(std::span(c), cstride, - *coefficients[coeff], cell_info, cells, - offsets[coeff]); - } - break; - } case IntegralType::interior_facet: { // Iterate over coefficients that are active in interior @@ -322,25 +297,27 @@ void pack_coefficients(const Form& form, } break; } + case IntegralType::exterior_facet: case IntegralType::vertex: + case IntegralType::ridge: { // Iterate over coefficients that are active in vertex integrals - for (int coeff : form.active_coeffs(IntegralType::vertex, id)) + for (int coeff : form.active_coeffs(integral_type, id)) { // Get coefficient mesh auto mesh = coefficients[coeff]->function_space()->mesh(); assert(mesh); - std::span vertices_b - = form.domain_coeff(IntegralType::vertex, id, coeff); + std::span entitites_b + = form.domain_coeff(integral_type, id, coeff); md::mdspan> - vertices(vertices_b.data(), vertices_b.size() / 2, 2); + entities(entitites_b.data(), entitites_b.size() / 2, 2); std::span cell_info = impl::get_cell_orientation_info(*coefficients[coeff]); impl::pack_coefficient_entity( std::span(c), cstride, *coefficients[coeff], cell_info, - md::submdspan(vertices, md::full_extent, 0), offsets[coeff]); + md::submdspan(entities, md::full_extent, 0), offsets[coeff]); } break; } diff --git a/python/test/unit/fem/test_assembler.py b/python/test/unit/fem/test_assembler.py index 5d7e5cf959..ec47f17a30 100644 --- a/python/test/unit/fem/test_assembler.py +++ b/python/test/unit/fem/test_assembler.py @@ -1996,6 +1996,7 @@ def marked_ridges(x): assert np.isclose(J_sum, ref_sol, atol=tol, rtol=tol) +@pytest.mark.parametrize("coefficient", [True, False]) @pytest.mark.parametrize( "cell_type", [ @@ -2025,11 +2026,10 @@ def marked_ridges(x): ), ], ) -def test_ridge_integrals_rank1_3D(cell_type, ghost_mode, dtype): +def test_ridge_integrals_rank1_3D(cell_type, ghost_mode, dtype, coefficient): comm = MPI.COMM_WORLD rdtype = np.real(dtype(0)).dtype - msh = None N = 4 msh = mesh.create_unit_cube( comm, N, N, N, cell_type=cell_type, ghost_mode=ghost_mode, dtype=rdtype @@ -2037,7 +2037,9 @@ def test_ridge_integrals_rank1_3D(cell_type, ghost_mode, dtype): gdim = msh.geometry.dim el = ("Lagrange", 3, (gdim,)) - V = dolfinx.fem.functionspace(msh, el) + + volume_element = basix.ufl.element(el[0], msh.basix_cell(), el[1], shape=el[2], dtype=rdtype) + V = dolfinx.fem.functionspace(msh, volume_element) x = ufl.SpatialCoordinate(msh) v = ufl.TestFunction(V) @@ -2052,12 +2054,17 @@ def marked_ridges(x): msh, msh.topology.dim - 2, exterior_ridges, np.full_like(exterior_ridges, 33) ) - def f(x): - return x[0] + ufl.sin(x[1]), x[2] + x[1], x[0] ** 2 - 3 * x[1] + def f(mod, x): + return x[0] + mod.sin(x[1]), x[2] + x[1], x[0] ** 2 - 3 * x[1] def integrand(x, v): - expr = ufl.as_vector(f(x)) - return ufl.inner(expr, v) + if coefficient: + Z = dolfinx.fem.functionspace(v.ufl_function_space().mesh, ("Lagrange", 1, (gdim,))) + z = dolfinx.fem.Function(Z, dtype=dtype) + z.interpolate(lambda x: f(np, x)) + else: + z = ufl.as_vector(f(ufl, x)) + return ufl.inner(z, v) metadata = {"quadrature_degree": 10} dr = ufl.Measure("dr", domain=msh, subdomain_data=et, subdomain_id=33, metadata=metadata) @@ -2083,7 +2090,9 @@ def integrand(x, v): connectivity = np.zeros((0, 2), dtype=np.int64) c_el = ufl.Mesh( - basix.ufl.element("Lagrange", basix.CellType.interval, 1, shape=(nodes.shape[1],)) + basix.ufl.element( + "Lagrange", basix.CellType.interval, 1, shape=(nodes.shape[1],), dtype=dtype + ), ) line_mesh = dolfinx.mesh.create_mesh( MPI.COMM_WORLD, @@ -2093,11 +2102,16 @@ def integrand(x, v): partitioner=dolfinx.mesh.create_cell_partitioner(ghost_mode), ) - Q = dolfinx.fem.functionspace(line_mesh, el) + line_element = basix.ufl.element( + el[0], line_mesh.basix_cell(), el[1], shape=el[2], dtype=rdtype + ) + Q = dolfinx.fem.functionspace(line_mesh, line_element) q = ufl.TestFunction(Q) x_e = ufl.SpatialCoordinate(line_mesh) - F_ref = dolfinx.fem.form(integrand(x_e, q) * ufl.dx(domain=line_mesh, metadata=metadata)) + F_ref = dolfinx.fem.form( + integrand(x_e, q) * ufl.dx(domain=line_mesh, metadata=metadata), dtype=dtype + ) b_ref = dolfinx.fem.assemble_vector(F_ref) b_ref.scatter_reverse(la.InsertMode.add) b_ref.scatter_forward() From 6016f7b5b4d7762fa39e9c4cc3bbc55cdb9c6a42 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= Date: Wed, 3 Sep 2025 13:43:27 +0000 Subject: [PATCH 09/37] Remove more duplication --- cpp/dolfinx/fem/utils.cpp | 36 +++++++----------------------------- 1 file changed, 7 insertions(+), 29 deletions(-) diff --git a/cpp/dolfinx/fem/utils.cpp b/cpp/dolfinx/fem/utils.cpp index c11baebff7..0aadd5fec6 100644 --- a/cpp/dolfinx/fem/utils.cpp +++ b/cpp/dolfinx/fem/utils.cpp @@ -202,22 +202,6 @@ fem::compute_integration_domains(fem::IntegralType integral_type, entity_data.insert(entity_data.begin(), entities.begin(), entities.end()); break; } - case IntegralType::exterior_facet: - { - auto [f_to_c, c_to_f] = get_connectivities(tdim - 1); - // Create list of tagged boundary facets - const std::vector bfacets = mesh::exterior_facet_indices(topology); - std::vector facets; - std::ranges::set_intersection(entities, bfacets, - std::back_inserter(facets)); - for (auto f : facets) - { - // Get the facet as a pair of (cell, local facet) - auto facet = impl::get_cell_facet_pairs<1>(f, f_to_c->links(f), *c_to_f); - entity_data.insert(entity_data.end(), facet.begin(), facet.end()); - } - break; - } case IntegralType::interior_facet: { auto [f_to_c, c_to_f] = get_connectivities(tdim - 1); @@ -251,28 +235,22 @@ fem::compute_integration_domains(fem::IntegralType integral_type, } break; } + case IntegralType::exterior_facet: case IntegralType::vertex: - { - auto [v_to_c, c_to_v] = get_connectivities(0); - for (auto vertex : entities) - { - std::array pair = impl::get_cell_entity_pairs<1>( - vertex, v_to_c->links(vertex), *c_to_v); - - entity_data.insert(entity_data.end(), pair.begin(), pair.end()); - } - break; - } case IntegralType::ridge: { - auto [r_to_c, c_to_r] = get_connectivities(tdim - 2); + auto [e_to_c, c_to_e] = get_connectivities(dim); for (auto entity : entities) { std::array pair = impl::get_cell_entity_pairs<1>( - entity, r_to_c->links(entity), *c_to_r); + entity, e_to_c->links(entity), *c_to_e); entity_data.insert(entity_data.end(), pair.begin(), pair.end()); } + break; } + default: + throw std::runtime_error( + "Cannot compute integration domains. Integral type not supported."); } return entity_data; From fca8d496059b7c0069151f6899c0d0fd82edc214 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= Date: Wed, 3 Sep 2025 14:14:06 +0000 Subject: [PATCH 10/37] Even more duplication --- cpp/dolfinx/fem/utils.h | 289 ++++++++++++---------------------------- 1 file changed, 82 insertions(+), 207 deletions(-) diff --git a/cpp/dolfinx/fem/utils.h b/cpp/dolfinx/fem/utils.h index d9860b821a..637e85f2c8 100644 --- a/cpp/dolfinx/fem/utils.h +++ b/cpp/dolfinx/fem/utils.h @@ -601,78 +601,9 @@ Form create_form_factory( } } - // Attach exterior facet kernels - std::vector default_facets_ext; - { - std::span ids(ufcx_forms[0].get().form_integral_ids - + integral_offsets[exterior_facet], - num_integrals_type[exterior_facet]); - auto sd = subdomains.find(IntegralType::exterior_facet); - for (std::size_t form_idx = 0; form_idx < ufcx_forms.size(); ++form_idx) - { - const ufcx_form& ufcx_form = ufcx_forms[form_idx]; - for (int i = 0; i < num_integrals_type[exterior_facet]; ++i) - { - const int id = ids[i]; - ufcx_integral* integral - = ufcx_form.form_integrals[integral_offsets[exterior_facet] + i]; - assert(integral); - check_geometry_hash(*integral, form_idx); - - std::vector active_coeffs; - for (int j = 0; j < ufcx_form.num_coefficients; ++j) - { - if (integral->enabled_coefficients[j]) - active_coeffs.push_back(j); - } - - impl::kernel_t k = impl::extract_kernel(integral); - - // Build list of entities to assembler over - const std::vector bfacets = mesh::exterior_facet_indices(*topology); - auto f_to_c = topology->connectivity(tdim - 1, tdim); - assert(f_to_c); - auto c_to_f = topology->connectivity(tdim, tdim - 1); - assert(c_to_f); - if (id == -1) - { - // Default kernel, operates on all (owned) exterior facets - default_facets_ext.reserve(2 * bfacets.size()); - for (std::int32_t f : bfacets) - { - // There will only be one pair for an exterior facet integral - std::array pair - = impl::get_cell_facet_pairs<1>(f, f_to_c->links(f), *c_to_f); - default_facets_ext.insert(default_facets_ext.end(), pair.begin(), - pair.end()); - } - integrals.insert({{IntegralType::exterior_facet, i, form_idx}, - {k, default_facets_ext, active_coeffs}}); - } - else if (sd != subdomains.end()) - { - // NOTE: This requires that pairs are sorted - auto it = std::ranges::lower_bound(sd->second, id, std::less<>{}, - [](auto& a) { return a.first; }); - if (it != sd->second.end() and it->first == id) - { - integrals.insert({{IntegralType::exterior_facet, i, form_idx}, - {k, - std::vector(it->second.begin(), - it->second.end()), - active_coeffs}}); - } - } - - if (integral->needs_facet_permutations) - needs_facet_permutations = true; - } - } - } - // Attach interior facet kernels - std::vector default_facets_int; { + std::vector default_facets_int; std::span ids(ufcx_forms[0].get().form_integral_ids + integral_offsets[interior_facet], num_integrals_type[interior_facet]); @@ -765,162 +696,106 @@ Form create_form_factory( } } - // Attach vertex kernels + // Attach exterior entity integrals { - for (std::size_t form_idx = 0; form_idx < ufcx_forms.size(); ++form_idx) + for (IntegralType itg_type : {IntegralType::exterior_facet, + IntegralType::vertex, IntegralType::ridge}) { - const ufcx_form& form = ufcx_forms[form_idx]; + std::size_t dim; + if (itg_type == IntegralType::exterior_facet) + dim = tdim - 1; + else if (itg_type == IntegralType::ridge) + dim = tdim - 2; + else if (itg_type == IntegralType::vertex) + dim = 0; + else + throw std::runtime_error("Unsupported integral type"); - std::span ids(form.form_integral_ids - + integral_offsets[vertex], - num_integrals_type[vertex]); - auto sd = subdomains.find(IntegralType::vertex); - for (int i = 0; i < num_integrals_type[vertex]; ++i) + auto get_default_integration_entities + = [dim](const mesh::Topology& topology, IntegralType itg_type) { - const int id = ids[i]; - ufcx_integral* integral - = form.form_integrals[integral_offsets[vertex] + i]; - assert(integral); - check_geometry_hash(*integral, form_idx); - - std::vector active_coeffs; - for (int j = 0; j < form.num_coefficients; ++j) + if (itg_type == IntegralType::exterior_facet) { - if (integral->enabled_coefficients[j]) - active_coeffs.push_back(j); + // Integrate over all owned exterior facets + return mesh::exterior_facet_indices(topology); } - - impl::kernel_t k = impl::extract_kernel(integral); - assert(k); - - // Build list of entities to assembler over - auto v_to_c = topology->connectivity(0, tdim); - assert(v_to_c); - auto c_to_v = topology->connectivity(tdim, 0); - assert(c_to_v); - - // pack for a range of vertices a flattened list of cell index c_i and - // local vertex index l_i: - // [c_0, l_0, ..., c_n, l_n] - auto get_cells_and_vertices = [v_to_c, c_to_v](auto vertices_range) - { - std::vector cell_and_vertex; - cell_and_vertex.reserve(2 * vertices_range.size()); - for (std::int32_t vertex : vertices_range) - { - std::array pair = impl::get_cell_entity_pairs<1>( - vertex, v_to_c->links(vertex), *c_to_v); - - cell_and_vertex.insert(cell_and_vertex.end(), pair.begin(), - pair.end()); - } - assert(cell_and_vertex.size() == 2 * vertices_range.size()); - return cell_and_vertex; - }; - - if (id == -1) + else { + // Integrate over all owned entities // Default vertex kernel operates on all (owned) vertices - std::int32_t num_vertices = topology->index_map(0)->size_local(); - std::vector cells_and_vertices = get_cells_and_vertices( - std::ranges::views::iota(0, num_vertices)); - - integrals.insert({{IntegralType::vertex, i, form_idx}, - {k, cells_and_vertices, active_coeffs}}); - } - else if (sd != subdomains.end()) - { - // NOTE: This requires that pairs are sorted - auto it = std::ranges::lower_bound(sd->second, id, std::less<>{}, - [](auto& a) { return a.first; }); - if (it != sd->second.end() and it->first == id) - { - integrals.insert({{IntegralType::vertex, i, form_idx}, - {k, - std::vector(it->second.begin(), - it->second.end()), - active_coeffs}}); - } + std::int32_t num_entities = topology.index_map(dim)->size_local(); + std::vector entities(num_entities); + std::iota(entities.begin(), entities.end(), 0); + return entities; } - } - } - } + }; - // Attach ridge kernels - { - for (std::size_t form_idx = 0; form_idx < ufcx_forms.size(); ++form_idx) - { - const ufcx_form& form = ufcx_forms[form_idx]; + std::vector default_entities_ext; - std::span ids(form.form_integral_ids + integral_offsets[ridge], - num_integrals_type[ridge]); - auto sd = subdomains.find(IntegralType::ridge); - for (int i = 0; i < num_integrals_type[ridge]; ++i) + std::span ids(ufcx_forms[0].get().form_integral_ids + + integral_offsets[(std::int8_t)itg_type], + num_integrals_type[(std::int8_t)itg_type]); + auto sd = subdomains.find(itg_type); + for (std::size_t form_idx = 0; form_idx < ufcx_forms.size(); ++form_idx) { - const int id = ids[i]; - ufcx_integral* integral - = form.form_integrals[integral_offsets[ridge] + i]; - assert(integral); - check_geometry_hash(*integral, form_idx); - - std::vector active_coeffs; - for (int j = 0; j < form.num_coefficients; ++j) - { - if (integral->enabled_coefficients[j]) - active_coeffs.push_back(j); - } - - impl::kernel_t k = impl::extract_kernel(integral); - assert(k); - - // Build list of entities to assembler over - auto r_to_c = topology->connectivity(tdim - 2, tdim); - assert(r_to_c); - auto c_to_r = topology->connectivity(tdim, tdim - 2); - assert(c_to_r); - - // pack for a range of vertices a flattened list of cell index c_i and - // local ridge index l_i: - // [c_0, l_0, ..., c_n, l_n] - auto get_cells_and_ridges = [r_to_c, c_to_r](auto ridge_range) + const ufcx_form& ufcx_form = ufcx_forms[form_idx]; + for (int i = 0; i < num_integrals_type[(std::int8_t)itg_type]; ++i) { - std::vector cell_and_ridge; - cell_and_ridge.reserve(2 * ridge_range.size()); - for (std::int32_t ridge : ridge_range) + const int id = ids[i]; + ufcx_integral* integral + = ufcx_form.form_integrals[integral_offsets[(std::int8_t)itg_type] + + i]; + assert(integral); + check_geometry_hash(*integral, form_idx); + + std::vector active_coeffs; + for (int j = 0; j < ufcx_form.num_coefficients; ++j) { - std::array pair = impl::get_cell_entity_pairs<1>( - ridge, r_to_c->links(ridge), *c_to_r); - - cell_and_ridge.insert(cell_and_ridge.end(), pair.begin(), - pair.end()); + if (integral->enabled_coefficients[j]) + active_coeffs.push_back(j); } - assert(cell_and_ridge.size() == 2 * ridge_range.size()); - return cell_and_ridge; - }; - if (id == -1) - { - // Default ridge kernel operates on all (owned) ridges - std::int32_t num_ridges = topology->index_map(tdim - 2)->size_local(); - std::vector cells_and_ridges - = get_cells_and_ridges(std::ranges::views::iota(0, num_ridges)); + impl::kernel_t k = impl::extract_kernel(integral); - integrals.insert({{IntegralType::ridge, i, form_idx}, - {k, cells_and_ridges, active_coeffs}}); - } - else if (sd != subdomains.end()) - { - // NOTE: This requires that pairs are sorted - auto it = std::ranges::lower_bound(sd->second, id, std::less<>{}, - [](auto& a) { return a.first; }); - if (it != sd->second.end() and it->first == id) + // Build list of entities to assembler over + auto e_to_c = topology->connectivity(dim, tdim); + assert(e_to_c); + auto c_to_e = topology->connectivity(tdim, dim); + assert(c_to_e); + if (id == -1) { - integrals.insert({{IntegralType::ridge, i, form_idx}, - {k, - std::vector(it->second.begin(), - it->second.end()), - active_coeffs}}); + std::vector default_entities + = get_default_integration_entities(*topology, itg_type); + // Default kernel + default_entities_ext.reserve(2 * default_entities.size()); + for (std::int32_t e : default_entities) + { + // There will only be one pair for an exterior facet integral + std::array pair = impl::get_cell_entity_pairs<1>( + e, e_to_c->links(e), *c_to_e); + default_entities_ext.insert(default_entities_ext.end(), + pair.begin(), pair.end()); + } + integrals.insert({{itg_type, i, form_idx}, + {k, default_entities_ext, active_coeffs}}); } + else if (sd != subdomains.end()) + { + // NOTE: This requires that pairs are sorted + auto it = std::ranges::lower_bound(sd->second, id, std::less<>{}, + [](auto& a) { return a.first; }); + if (it != sd->second.end() and it->first == id) + { + integrals.insert({{itg_type, i, form_idx}, + {k, + std::vector(it->second.begin(), + it->second.end()), + active_coeffs}}); + } + } + + if (integral->needs_facet_permutations) + needs_facet_permutations = true; } } } From 34d2b997b44fe4c33bb5f7adacd51c4da1272da8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= Date: Wed, 3 Sep 2025 14:14:39 +0000 Subject: [PATCH 11/37] Use real dtype for basix element --- python/test/unit/fem/test_assembler.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/test/unit/fem/test_assembler.py b/python/test/unit/fem/test_assembler.py index ec47f17a30..dc39463c7d 100644 --- a/python/test/unit/fem/test_assembler.py +++ b/python/test/unit/fem/test_assembler.py @@ -2091,7 +2091,7 @@ def integrand(x, v): c_el = ufl.Mesh( basix.ufl.element( - "Lagrange", basix.CellType.interval, 1, shape=(nodes.shape[1],), dtype=dtype + "Lagrange", basix.CellType.interval, 1, shape=(nodes.shape[1],), dtype=rdtype ), ) line_mesh = dolfinx.mesh.create_mesh( From 18b8907a2feae0c658f49ece3bc0f02756ccb7e7 Mon Sep 17 00:00:00 2001 From: jorgensd Date: Wed, 3 Sep 2025 16:02:55 +0000 Subject: [PATCH 12/37] Use switch --- cpp/dolfinx/fem/utils.h | 20 ++++++++++++++++---- 1 file changed, 16 insertions(+), 4 deletions(-) diff --git a/cpp/dolfinx/fem/utils.h b/cpp/dolfinx/fem/utils.h index 637e85f2c8..daf2f892bf 100644 --- a/cpp/dolfinx/fem/utils.h +++ b/cpp/dolfinx/fem/utils.h @@ -702,14 +702,26 @@ Form create_form_factory( IntegralType::vertex, IntegralType::ridge}) { std::size_t dim; - if (itg_type == IntegralType::exterior_facet) + switch (itg_type) + { + case IntegralType::exterior_facet: + { dim = tdim - 1; - else if (itg_type == IntegralType::ridge) + break; + } + case IntegralType::ridge: + { dim = tdim - 2; - else if (itg_type == IntegralType::vertex) + break; + } + case IntegralType::vertex: + { dim = 0; - else + break; + } + default: throw std::runtime_error("Unsupported integral type"); + } auto get_default_integration_entities = [dim](const mesh::Topology& topology, IntegralType itg_type) From bc2e17c611b9748bfe8c5f768dec8f8827969d27 Mon Sep 17 00:00:00 2001 From: jorgensd Date: Wed, 3 Sep 2025 16:51:15 +0000 Subject: [PATCH 13/37] Use more switches --- cpp/dolfinx/fem/assemble_scalar_impl.h | 10 +++++++--- cpp/dolfinx/fem/assemble_vector_impl.h | 10 +++++++--- 2 files changed, 14 insertions(+), 6 deletions(-) diff --git a/cpp/dolfinx/fem/assemble_scalar_impl.h b/cpp/dolfinx/fem/assemble_scalar_impl.h index eee96b3618..b6f58606f1 100644 --- a/cpp/dolfinx/fem/assemble_scalar_impl.h +++ b/cpp/dolfinx/fem/assemble_scalar_impl.h @@ -220,10 +220,14 @@ T assemble_scalar( fem::IntegralType::ridge}) { md::mdspan> perms; - if (itg_type == fem::IntegralType::exterior_facet && !facet_perms.empty()) + switch (itg_type) + { + case fem::IntegralType::exterior_facet: + { perms = facet_perms; - else - perms = md::mdspan>(); + break; + } + } for (int i = 0; i < M.num_integrals(itg_type, 0); ++i) { diff --git a/cpp/dolfinx/fem/assemble_vector_impl.h b/cpp/dolfinx/fem/assemble_vector_impl.h index 52d1830893..27379f0670 100644 --- a/cpp/dolfinx/fem/assemble_vector_impl.h +++ b/cpp/dolfinx/fem/assemble_vector_impl.h @@ -1354,10 +1354,14 @@ void assemble_vector( fem::IntegralType::ridge}) { md::mdspan> perms; - if (itg_type == fem::IntegralType::exterior_facet && !facet_perms.empty()) + switch (itg_type) + { + case fem::IntegralType::exterior_facet: + { perms = facet_perms; - else - perms = md::mdspan>(); + break; + } + } for (int i = 0; i < L.num_integrals(itg_type, 0); ++i) { auto fn = L.kernel(itg_type, i, 0); From 4097d9fc4a5515d2dfc619989ee7f6938e3cd600 Mon Sep 17 00:00:00 2001 From: jorgensd Date: Wed, 3 Sep 2025 17:12:23 +0000 Subject: [PATCH 14/37] Add default --- cpp/dolfinx/fem/assemble_scalar_impl.h | 2 ++ cpp/dolfinx/fem/assemble_vector_impl.h | 2 ++ 2 files changed, 4 insertions(+) diff --git a/cpp/dolfinx/fem/assemble_scalar_impl.h b/cpp/dolfinx/fem/assemble_scalar_impl.h index b6f58606f1..e535815683 100644 --- a/cpp/dolfinx/fem/assemble_scalar_impl.h +++ b/cpp/dolfinx/fem/assemble_scalar_impl.h @@ -227,6 +227,8 @@ T assemble_scalar( perms = facet_perms; break; } + default: + break; } for (int i = 0; i < M.num_integrals(itg_type, 0); ++i) diff --git a/cpp/dolfinx/fem/assemble_vector_impl.h b/cpp/dolfinx/fem/assemble_vector_impl.h index 27379f0670..3161354a3a 100644 --- a/cpp/dolfinx/fem/assemble_vector_impl.h +++ b/cpp/dolfinx/fem/assemble_vector_impl.h @@ -1361,6 +1361,8 @@ void assemble_vector( perms = facet_perms; break; } + default: + break; } for (int i = 0; i < L.num_integrals(itg_type, 0); ++i) { From cd71b1c32bee40f5aecb6c056bf67708a6cd535a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= Date: Thu, 4 Sep 2025 07:36:49 +0000 Subject: [PATCH 15/37] add vertex and ridge matrix integrals + lifting --- cpp/dolfinx/fem/assemble_matrix_impl.h | 121 ++++++++++++++----------- cpp/dolfinx/fem/assemble_vector_impl.h | 108 ++++++++++++---------- 2 files changed, 130 insertions(+), 99 deletions(-) diff --git a/cpp/dolfinx/fem/assemble_matrix_impl.h b/cpp/dolfinx/fem/assemble_matrix_impl.h index 90351a05b7..f2404b9e61 100644 --- a/cpp/dolfinx/fem/assemble_matrix_impl.h +++ b/cpp/dolfinx/fem/assemble_matrix_impl.h @@ -156,7 +156,7 @@ void assemble_cells_matrix( } } -/// @brief Execute kernel over exterior facets and accumulate result in +/// @brief Execute kernel over exterior entities and accumulate result in /// a matrix. /// /// @tparam T Matrix/form scalar type. @@ -164,8 +164,8 @@ void assemble_cells_matrix( /// matrix. /// @param[in] x_dofmap Dofmap for the mesh geometry. /// @param[in] x Mesh geometry (coordinates). -/// @param[in] facets Facet indices (in the integration domain mesh) to -/// execute the kernel over. +/// @param[in] entities Integration entities (in the integration domain mesh) to +/// execute the kernel over. These are pairs (cell, local entity index) /// @param[in] dofmap0 Test function (row) degree-of-freedom data /// holding the (0) dofmap, (1) dofmap block size and (2) dofmap cell /// indices. @@ -191,14 +191,14 @@ void assemble_cells_matrix( /// @param[in] perms Facet permutation integer. Empty if facet /// permutations are not required. template -void assemble_exterior_facets( +void assemble_exterior_entities( la::MatSet auto mat_set, mdspan2_t x_dofmap, md::mdspan, md::extents> x, md::mdspan> - facets, + entities, std::tuple>> @@ -215,11 +215,11 @@ void assemble_exterior_facets( std::span cell_info1, md::mdspan> perms) { - if (facets.empty()) + if (entities.empty()) return; - const auto [dmap0, bs0, facets0] = dofmap0; - const auto [dmap1, bs1, facets1] = dofmap1; + const auto [dmap0, bs0, entities0] = dofmap0; + const auto [dmap1, bs1, entities1] = dofmap1; // Data structures used in assembly std::vector> cdofs(3 * x_dofmap.extent(1)); @@ -228,17 +228,17 @@ void assemble_exterior_facets( const int ndim0 = bs0 * num_dofs0; const int ndim1 = bs1 * num_dofs1; std::vector Ae(ndim0 * ndim1); - assert(facets0.size() == facets.size()); - assert(facets1.size() == facets.size()); - for (std::size_t f = 0; f < facets.extent(0); ++f) + assert(entities0.size() == entities.size()); + assert(entities1.size() == entities.size()); + for (std::size_t f = 0; f < entities.extent(0); ++f) { // Cell in the integration domain, local facet index relative to the // integration domain cell, and cells in the test and trial function // meshes - std::int32_t cell = facets(f, 0); - std::int32_t local_facet = facets(f, 1); - std::int32_t cell0 = facets0(f, 0); - std::int32_t cell1 = facets1(f, 0); + std::int32_t cell = entities(f, 0); + std::int32_t local_entity = entities(f, 1); + std::int32_t cell0 = entities0(f, 0); + std::int32_t cell1 = entities1(f, 0); // Get cell coordinates/geometry auto x_dofs = md::submdspan(x_dofmap, cell, md::full_extent); @@ -246,12 +246,12 @@ void assemble_exterior_facets( std::copy_n(&x(x_dofs[i], 0), 3, std::next(cdofs.begin(), 3 * i)); // Permutations - std::uint8_t perm = perms.empty() ? 0 : perms(cell, local_facet); + std::uint8_t perm = perms.empty() ? 0 : perms(cell, local_entity); // Tabulate tensor std::ranges::fill(Ae, 0); kernel(Ae.data(), &coeffs(f, 0), constants.data(), cdofs.data(), - &local_facet, &perm, nullptr); + &local_entity, &perm, nullptr); P0(Ae, cell_info0, cell0, ndim1); P1T(Ae, cell_info1, cell1, ndim0); @@ -605,7 +605,7 @@ void assemble_matrix( cell_info0, cell_info1); } - md::mdspan> perms; + md::mdspan> facet_perms; if (a.needs_facet_permutations()) { mesh::CellType cell_type = mesh->topology()->cell_types()[cell_type_idx]; @@ -614,40 +614,8 @@ void assemble_matrix( mesh->topology_mutable()->create_entity_permutations(); const std::vector& p = mesh->topology()->get_facet_permutations(); - perms = md::mdspan(p.data(), p.size() / num_facets_per_cell, - num_facets_per_cell); - } - - for (int i = 0; - i < a.num_integrals(IntegralType::exterior_facet, cell_type_idx); ++i) - { - if (num_cell_types > 1) - { - throw std::runtime_error("Exterior facet integrals with mixed " - "topology aren't supported yet"); - } - - using mdspanx2_t - = md::mdspan>; - - auto fn = a.kernel(IntegralType::exterior_facet, i, 0); - assert(fn); - auto& [coeffs, cstride] - = coefficients.at({IntegralType::exterior_facet, i}); - - std::span f = a.domain(IntegralType::exterior_facet, i, 0); - mdspanx2_t facets(f.data(), f.size() / 2, 2); - std::span f0 = a.domain_arg(IntegralType::exterior_facet, 0, i, 0); - mdspanx2_t facets0(f0.data(), f0.size() / 2, 2); - std::span f1 = a.domain_arg(IntegralType::exterior_facet, 1, i, 0); - mdspanx2_t facets1(f1.data(), f1.size() / 2, 2); - assert((facets.size() / 2) * cstride == coeffs.size()); - impl::assemble_exterior_facets( - mat_set, x_dofmap, x, facets, {dofs0, bs0, facets0}, P0, - {dofs1, bs1, facets1}, P1T, bc0, bc1, fn, - md::mdspan(coeffs.data(), facets.extent(0), cstride), constants, - cell_info0, cell_info1, perms); + facet_perms = md::mdspan(p.data(), p.size() / num_facets_per_cell, + num_facets_per_cell); } for (int i = 0; @@ -685,7 +653,54 @@ void assemble_matrix( mdspanx22_t(facets1.data(), facets1.size() / 4, 2, 2)}, P1T, bc0, bc1, fn, mdspanx2x_t(coeffs.data(), facets.size() / 4, 2, cstride), constants, - cell_info0, cell_info1, perms); + cell_info0, cell_info1, facet_perms); + } + + for (fem::IntegralType itg_type : + {fem::IntegralType::exterior_facet, fem::IntegralType::vertex, + fem::IntegralType::ridge}) + { + md::mdspan> perms; + switch (itg_type) + { + case fem::IntegralType::exterior_facet: + { + perms = facet_perms; + break; + } + default: + break; + } + + for (int i = 0; i < a.num_integrals(itg_type, cell_type_idx); ++i) + { + if (num_cell_types > 1) + { + throw std::runtime_error("Exterior facet integrals with mixed " + "topology aren't supported yet"); + } + + using mdspanx2_t + = md::mdspan>; + + auto fn = a.kernel(itg_type, i, 0); + assert(fn); + auto& [coeffs, cstride] = coefficients.at({itg_type, i}); + + std::span e = a.domain(itg_type, i, 0); + mdspanx2_t entities(e.data(), e.size() / 2, 2); + std::span e0 = a.domain_arg(itg_type, 0, i, 0); + mdspanx2_t entities0(e0.data(), e0.size() / 2, 2); + std::span e1 = a.domain_arg(itg_type, 1, i, 0); + mdspanx2_t entities1(e1.data(), e1.size() / 2, 2); + assert((entities.size() / 2) * cstride == coeffs.size()); + impl::assemble_exterior_entities( + mat_set, x_dofmap, x, entities, {dofs0, bs0, entities0}, P0, + {dofs1, bs1, entities1}, P1T, bc0, bc1, fn, + md::mdspan(coeffs.data(), entities.extent(0), cstride), constants, + cell_info0, cell_info1, perms); + } } } } diff --git a/cpp/dolfinx/fem/assemble_vector_impl.h b/cpp/dolfinx/fem/assemble_vector_impl.h index 3161354a3a..9ab29f87bf 100644 --- a/cpp/dolfinx/fem/assemble_vector_impl.h +++ b/cpp/dolfinx/fem/assemble_vector_impl.h @@ -229,9 +229,9 @@ void _lift_bc_cells( /// @param[in,out] b Vector to modify. /// @param[in] x_dofmap Degree-of-freedom map for the mesh geometry. /// @param[in] x Mesh geometry (coordinates). -/// @param[in] kernel Kernel function to execute over each facet. -/// @param[in] facets Facets to execute the kernel over, where for the -/// `i`th facet `facets(i, 0)` is the attached cell and `facets(i, 1)` +/// @param[in] kernel Kernel function to execute over each entity. +/// @param[in] entities Entities to execute the kernel over, where for the +/// `i`th entity `enities(i, 0)` is the attached cell and `entities(i, 1)` /// is the local index of the facet relative to the cell. /// @param[in] dofmap0 Test function (row) degree-of-freedom data /// holding the (0) dofmap, (1) dofmap block size and (2) dofmap @@ -263,7 +263,7 @@ void _lift_bc_cells( template ::value_type> requires std::is_same_v::value_type, T> -void _lift_bc_exterior_facets( +void _lift_bc_exterior_entities( V&& b, mdspan2_t x_dofmap, md::mdspan, md::extents> @@ -271,7 +271,7 @@ void _lift_bc_exterior_facets( FEkernel auto kernel, md::mdspan> - facets, + entities, std::tuple>> @@ -288,11 +288,11 @@ void _lift_bc_exterior_facets( std::span bc_markers1, std::span x0, T alpha, md::mdspan> perms) { - if (facets.empty()) + if (entities.empty()) return; - const auto [dmap0, bs0, facets0] = dofmap0; - const auto [dmap1, bs1, facets1] = dofmap1; + const auto [dmap0, bs0, entities0] = dofmap0; + const auto [dmap1, bs1, entities1] = dofmap1; const int num_rows = bs0 * dmap0.extent(1); const int num_cols = bs1 * dmap1.extent(1); @@ -300,18 +300,18 @@ void _lift_bc_exterior_facets( // Data structures used in bc application std::vector> cdofs(3 * x_dofmap.extent(1)); std::vector Ae(num_rows * num_cols), be(num_rows); - assert(facets0.size() == facets.size()); - assert(facets1.size() == facets.size()); - for (std::size_t index = 0; index < facets.extent(0); ++index) + assert(entities0.size() == entities.size()); + assert(entities1.size() == entities.size()); + for (std::size_t index = 0; index < entities.extent(0); ++index) { // Cell in integration domain, test function and trial function // meshes - std::int32_t cell = facets(index, 0); - std::int32_t cell0 = facets0(index, 0); - std::int32_t cell1 = facets1(index, 0); + std::int32_t cell = entities(index, 0); + std::int32_t cell0 = entities0(index, 0); + std::int32_t cell1 = entities1(index, 0); - // Local facet index - std::int32_t local_facet = facets(index, 1); + // Local entity index + std::int32_t local_entity = entities(index, 1); // Get dof maps for cell auto dofs1 = md::submdspan(dmap1, cell1, md::full_extent); @@ -342,10 +342,10 @@ void _lift_bc_exterior_facets( auto dofs0 = md::submdspan(dmap0, cell0, md::full_extent); // Permutations - std::uint8_t perm = perms.empty() ? 0 : perms(cell, local_facet); + std::uint8_t perm = perms.empty() ? 0 : perms(cell, local_entity); std::ranges::fill(Ae, 0); kernel(Ae.data(), &coeffs(index, 0), constants.data(), cdofs.data(), - &local_facet, &perm, nullptr); + &local_entity, &perm, nullptr); P0(Ae, cell_info0, cell0, num_cols); P1T(Ae, cell_info1, cell1, num_rows); @@ -1045,7 +1045,7 @@ void lift_bc(V&& b, const Form& a, mdspan2_t x_dofmap, } } - md::mdspan> perms; + md::mdspan> facet_perms; if (a.needs_facet_permutations()) { mesh::CellType cell_type = mesh->topology()->cell_type(); @@ -1054,32 +1054,8 @@ void lift_bc(V&& b, const Form& a, mdspan2_t x_dofmap, mesh->topology_mutable()->create_entity_permutations(); const std::vector& p = mesh->topology()->get_facet_permutations(); - perms = md::mdspan(p.data(), p.size() / num_facets_per_cell, - num_facets_per_cell); - } - - for (int i = 0; i < a.num_integrals(IntegralType::exterior_facet, 0); ++i) - { - auto kernel = a.kernel(IntegralType::exterior_facet, i, 0); - assert(kernel); - auto& [coeffs, cstride] - = coefficients.at({IntegralType::exterior_facet, i}); - - using mdspanx2_t - = md::mdspan>; - std::span f = a.domain(IntegralType::exterior_facet, i, 0); - mdspanx2_t facets(f.data(), f.size() / 2, 2); - std::span f0 = a.domain_arg(IntegralType::exterior_facet, 0, i, 0); - mdspanx2_t facets0(f0.data(), f0.size() / 2, 2); - std::span f1 = a.domain_arg(IntegralType::exterior_facet, 1, i, 0); - mdspanx2_t facets1(f1.data(), f1.size() / 2, 2); - assert(coeffs.size() == facets.extent(0) * cstride); - _lift_bc_exterior_facets( - b, x_dofmap, x, kernel, facets, {dofmap0, bs0, facets0}, P0, - {dofmap1, bs1, facets1}, P1T, constants, - md::mdspan(coeffs.data(), facets.extent(0), cstride), cell_info0, - cell_info1, bc_values1, bc_markers1, x0, alpha, perms); + facet_perms = md::mdspan(p.data(), p.size() / num_facets_per_cell, + num_facets_per_cell); } for (int i = 0; i < a.num_integrals(IntegralType::interior_facet, 0); ++i) @@ -1105,7 +1081,47 @@ void lift_bc(V&& b, const Form& a, mdspan2_t x_dofmap, b, x_dofmap, x, kernel, facets, {dofmap0, bs0, facets0}, P0, {dofmap1, bs1, facets1}, P1T, constants, mdspanx2x_t(coeffs.data(), facets.extent(0), 2, cstride), cell_info0, - cell_info1, bc_values1, bc_markers1, x0, alpha, perms); + cell_info1, bc_values1, bc_markers1, x0, alpha, facet_perms); + } + + for (fem::IntegralType itg_type : + {fem::IntegralType::exterior_facet, fem::IntegralType::vertex, + fem::IntegralType::ridge}) + { + md::mdspan> perms; + switch (itg_type) + { + case fem::IntegralType::exterior_facet: + { + perms = facet_perms; + break; + } + default: + break; + } + + for (int i = 0; i < a.num_integrals(itg_type, 0); ++i) + { + auto kernel = a.kernel(itg_type, i, 0); + assert(kernel); + auto& [coeffs, cstride] = coefficients.at({itg_type, i}); + + using mdspanx2_t + = md::mdspan>; + std::span e = a.domain(itg_type, i, 0); + mdspanx2_t entities(e.data(), e.size() / 2, 2); + std::span e0 = a.domain_arg(itg_type, 0, i, 0); + mdspanx2_t entities0(e0.data(), e0.size() / 2, 2); + std::span e1 = a.domain_arg(itg_type, 1, i, 0); + mdspanx2_t entities1(e1.data(), e1.size() / 2, 2); + assert(coeffs.size() == entities.extent(0) * cstride); + _lift_bc_exterior_entities( + b, x_dofmap, x, kernel, entities, {dofmap0, bs0, entities0}, P0, + {dofmap1, bs1, entities1}, P1T, constants, + md::mdspan(coeffs.data(), entities.extent(0), cstride), cell_info0, + cell_info1, bc_values1, bc_markers1, x0, alpha, perms); + } } } From 12342723a857e717476fe628fa1f2a33c0581cda Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= Date: Thu, 4 Sep 2025 07:52:11 +0000 Subject: [PATCH 16/37] Adapt lifting test to check matrix and lifting of ridge and vertex integrals --- cpp/dolfinx/fem/utils.h | 2 ++ python/test/unit/fem/test_assembler.py | 29 ++++++++++++++++++++++++-- 2 files changed, 29 insertions(+), 2 deletions(-) diff --git a/cpp/dolfinx/fem/utils.h b/cpp/dolfinx/fem/utils.h index daf2f892bf..e2e20c70ad 100644 --- a/cpp/dolfinx/fem/utils.h +++ b/cpp/dolfinx/fem/utils.h @@ -291,6 +291,8 @@ void build_sparsity_pattern(la::SparsityPattern& pattern, const Form& a) } break; case IntegralType::exterior_facet: + case IntegralType::ridge: + case IntegralType::vertex: for (int i = 0; i < a.num_integrals(type, cell_type_idx); ++i) { sparsitybuild::cells(pattern, diff --git a/python/test/unit/fem/test_assembler.py b/python/test/unit/fem/test_assembler.py index dc39463c7d..82b86324c4 100644 --- a/python/test/unit/fem/test_assembler.py +++ b/python/test/unit/fem/test_assembler.py @@ -48,7 +48,7 @@ locate_entities_boundary, meshtags, ) -from ufl import derivative, dS, ds, dx, inner +from ufl import derivative, dP, dr, dS, ds, dx, inner from ufl.geometry import SpatialCoordinate dtype_parametrize = pytest.mark.parametrize( @@ -224,7 +224,32 @@ def test_assembly_bcs(self, mode): mesh = create_unit_square(MPI.COMM_WORLD, 12, 12, ghost_mode=mode) V = functionspace(mesh, ("Lagrange", 1)) u, v = ufl.TrialFunction(V), ufl.TestFunction(V) - a = form(inner(u, v) * dx + inner(u, v) * ds) + + exterior_ridges = locate_entities_boundary( + mesh, mesh.topology.dim - 2, lambda x: np.full(x.shape[1], True, dtype=bool) + ) + ridge_tags = dolfinx.mesh.meshtags( + mesh, + mesh.topology.dim - 2, + exterior_ridges, + np.ones_like(exterior_ridges, dtype=np.int32), + ) + dr_exterior = dr(subdomain_data=ridge_tags, subdomain_id=1) + exterior_vertices = locate_entities_boundary( + mesh, 0, lambda x: np.isclose(x[0], 0.0) | np.isclose(x[1], 1.0) + ) + vertex_tags = dolfinx.mesh.meshtags( + mesh, 0, exterior_vertices, np.ones_like(exterior_vertices, dtype=np.int32) + ) + dP_exterior = dP(subdomain_data=vertex_tags, subdomain_id=1) + x = ufl.SpatialCoordinate(mesh) + f = ufl.exp(x[0] + x[1]) + a = form( + inner(u, v) * dx + + inner(u, v) * ds + - inner(u, v) * dr_exterior + + f * inner(u, v) * dP_exterior + ) L = form(inner(1.0, v) * dx) bdofsV = locate_dofs_geometrical(V, lambda x: np.isclose(x[0], 0.0) | np.isclose(x[0], 1.0)) From 1d32d8acc69d470fa3ab4cd71b87686b9bf87752 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= Date: Tue, 9 Sep 2025 12:44:01 +0000 Subject: [PATCH 17/37] Merge main into branch --- .github/workflows/build-wheels.yml | 2 +- .github/workflows/ccpp.yml | 4 +-- .github/workflows/macos.yml | 2 +- .github/workflows/windows.yml | 6 ++-- cpp/dolfinx/fem/assemble_scalar_impl.h | 44 +++++++++++++------------- 5 files changed, 29 insertions(+), 29 deletions(-) diff --git a/.github/workflows/build-wheels.yml b/.github/workflows/build-wheels.yml index 5a582900a5..e5ddabf038 100644 --- a/.github/workflows/build-wheels.yml +++ b/.github/workflows/build-wheels.yml @@ -72,7 +72,7 @@ jobs: CIBW_ENVIRONMENT: PIP_EXTRA_INDEX_URL=file:///project/simple PETSC_DIR=/usr/local MAKEFLAGS=-j2 steps: - - uses: actions/setup-python@v5 + - uses: actions/setup-python@v6 - name: Install Python dependencies run: python -m pip install cibuildwheel simple503 wheel diff --git a/.github/workflows/ccpp.yml b/.github/workflows/ccpp.yml index 087519a8de..6e9069c2a3 100644 --- a/.github/workflows/ccpp.yml +++ b/.github/workflows/ccpp.yml @@ -24,7 +24,7 @@ jobs: uses: actions/checkout@v5 - name: Set up Python - uses: actions/setup-python@v5 + uses: actions/setup-python@v6 with: python-version: '3.12' @@ -87,7 +87,7 @@ jobs: libhdf5-mpi-dev liblapack-dev libparmetis-dev libpugixml-dev \ libspdlog-dev mpi-default-dev ninja-build pkg-config - name: Set up Python - uses: actions/setup-python@v5 + uses: actions/setup-python@v6 with: python-version: "3.12" diff --git a/.github/workflows/macos.yml b/.github/workflows/macos.yml index f073088068..6310af71a8 100644 --- a/.github/workflows/macos.yml +++ b/.github/workflows/macos.yml @@ -34,7 +34,7 @@ jobs: echo DYLD_LIBRARY_PATH=/usr/local/lib >> $GITHUB_ENV - name: Set up Python - uses: actions/setup-python@v5 + uses: actions/setup-python@v6 id: cp3 with: python-version: "3.12" diff --git a/.github/workflows/windows.yml b/.github/workflows/windows.yml index f8340e8ad8..5bc75252f4 100644 --- a/.github/workflows/windows.yml +++ b/.github/workflows/windows.yml @@ -43,7 +43,7 @@ jobs: -Source "${{ env.FEED_URL }}" - name: Checkout DOLFINx - uses: actions/checkout@v4 + uses: actions/checkout@v5 with: path: dolfinx-src @@ -52,12 +52,12 @@ jobs: cat dolfinx-src/.github/workflows/fenicsx-refs.env >> $GITHUB_ENV - name: Set up Python - uses: actions/setup-python@v5 + uses: actions/setup-python@v6 with: python-version: "3.12" - name: Checkout UFL - uses: actions/checkout@v4 + uses: actions/checkout@v5 with: repository: "fenics/ufl" path: ufl diff --git a/cpp/dolfinx/fem/assemble_scalar_impl.h b/cpp/dolfinx/fem/assemble_scalar_impl.h index e535815683..418f72b691 100644 --- a/cpp/dolfinx/fem/assemble_scalar_impl.h +++ b/cpp/dolfinx/fem/assemble_scalar_impl.h @@ -29,14 +29,14 @@ T assemble_cells(mdspan2_t x_dofmap, x, std::span cells, FEkernel auto fn, std::span constants, - md::mdspan> coeffs) + md::mdspan> coeffs, + std::span> cdofs_b) { T value(0); if (cells.empty()) return value; - // Create data structures used in assembly - std::vector> cdofs(3 * x_dofmap.extent(1)); + assert(cdofs_b.size() >= 3 * x_dofmap.extent(1)); // Iterate over all cells for (std::size_t index = 0; index < cells.size(); ++index) @@ -46,11 +46,9 @@ T assemble_cells(mdspan2_t x_dofmap, // Get cell coordinates/geometry auto x_dofs = md::submdspan(x_dofmap, c, md::full_extent); for (std::size_t i = 0; i < x_dofs.size(); ++i) - { - std::copy_n(&x(x_dofs[i], 0), 3, std::next(cdofs.begin(), 3 * i)); - } + std::copy_n(&x(x_dofs[i], 0), 3, std::next(cdofs_b.begin(), 3 * i)); - fn(&value, &coeffs(index, 0), constants.data(), cdofs.data(), nullptr, + fn(&value, &coeffs(index, 0), constants.data(), cdofs_b.data(), nullptr, nullptr, nullptr); } @@ -69,14 +67,14 @@ T assemble_exterior_entities( entities, FEkernel auto fn, std::span constants, md::mdspan> coeffs, - md::mdspan> perms) + md::mdspan> perms, + std::span> cdofs_b) { T value(0); if (entities.empty()) return value; - // Create data structures used in assembly - std::vector> cdofs(3 * x_dofmap.extent(1)); + assert(cdofs_b.size() >= 3 * x_dofmap.extent(1)); // Iterate over all facets for (std::size_t f = 0; f < entities.extent(0); ++f) @@ -87,11 +85,11 @@ T assemble_exterior_entities( // Get cell coordinates/geometry auto x_dofs = md::submdspan(x_dofmap, cell, md::full_extent); for (std::size_t i = 0; i < x_dofs.size(); ++i) - std::copy_n(&x(x_dofs[i], 0), 3, std::next(cdofs.begin(), 3 * i)); + std::copy_n(&x(x_dofs[i], 0), 3, std::next(cdofs_b.begin(), 3 * i)); // Permutations std::uint8_t perm = perms.empty() ? 0 : perms(cell, local_entity); - fn(&value, &coeffs(f, 0), constants.data(), cdofs.data(), &local_entity, + fn(&value, &coeffs(f, 0), constants.data(), cdofs_b.data(), &local_entity, &perm, nullptr); } @@ -112,18 +110,17 @@ T assemble_interior_facets( md::mdspan> coeffs, - md::mdspan> perms) + md::mdspan> perms, + std::span> cdofs_b) { T value(0); if (facets.empty()) return value; // Create data structures used in assembly - using X = scalar_value_t; - std::vector cdofs(2 * x_dofmap.extent(1) * 3); - std::span cdofs0(cdofs.data(), x_dofmap.extent(1) * 3); - std::span cdofs1(cdofs.data() + x_dofmap.extent(1) * 3, - x_dofmap.extent(1) * 3); + assert(cdofs_b.size() >= 2 * 3 * x_dofmap.extent(1)); + auto cdofs0 = cdofs_b.first(3 * x_dofmap.extent(1)); + auto cdofs1 = cdofs_b.last(3 * x_dofmap.extent(1)); // Iterate over all facets for (std::size_t f = 0; f < facets.extent(0); ++f) @@ -143,7 +140,7 @@ T assemble_interior_facets( ? std::array{0, 0} : std::array{perms(cells[0], local_facet[0]), perms(cells[1], local_facet[1])}; - fn(&value, &coeffs(f, 0, 0), constants.data(), cdofs.data(), + fn(&value, &coeffs(f, 0, 0), constants.data(), cdofs_b.data(), local_facet.data(), perm.data(), nullptr); } @@ -164,6 +161,8 @@ T assemble_scalar( std::shared_ptr> mesh = M.mesh(); assert(mesh); + std::vector> cdofs_b(2 * 3 * x_dofmap.extent(1)); + T value = 0; for (int i = 0; i < M.num_integrals(IntegralType::cell, 0); ++i) { @@ -174,7 +173,7 @@ T assemble_scalar( assert(cells.size() * cstride == coeffs.size()); value += impl::assemble_cells( x_dofmap, x, cells, fn, constants, - md::mdspan(coeffs.data(), cells.size(), cstride)); + md::mdspan(coeffs.data(), cells.size(), cstride), cdofs_b); } mesh::CellType cell_type = mesh->topology()->cell_type(); @@ -212,7 +211,7 @@ T assemble_scalar( md::mdspan>( coeffs.data(), facets.size() / shape1, 2, cstride), - facet_perms); + facet_perms, cdofs_b); } for (fem::IntegralType itg_type : @@ -250,7 +249,8 @@ T assemble_scalar( md::extents>( entities.data(), entities.size() / shape1, 2), fn, constants, - md::mdspan(coeffs.data(), entities.size() / shape1, cstride), perms); + md::mdspan(coeffs.data(), entities.size() / shape1, cstride), perms, + cdofs_b); } } From 28c64f9400041c90c6b4ed19472055ad10955831 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= Date: Tue, 9 Sep 2025 14:23:27 +0000 Subject: [PATCH 18/37] Improve name of inner-most assembly function --- cpp/dolfinx/fem/assemble_matrix_impl.h | 6 +++--- cpp/dolfinx/fem/assemble_scalar_impl.h | 6 +++--- cpp/dolfinx/fem/assemble_vector_impl.h | 10 +++++----- 3 files changed, 11 insertions(+), 11 deletions(-) diff --git a/cpp/dolfinx/fem/assemble_matrix_impl.h b/cpp/dolfinx/fem/assemble_matrix_impl.h index f2404b9e61..6cbf865a70 100644 --- a/cpp/dolfinx/fem/assemble_matrix_impl.h +++ b/cpp/dolfinx/fem/assemble_matrix_impl.h @@ -156,7 +156,7 @@ void assemble_cells_matrix( } } -/// @brief Execute kernel over exterior entities and accumulate result in +/// @brief Execute kernel over a set of entities and accumulate result in /// a matrix. /// /// @tparam T Matrix/form scalar type. @@ -191,7 +191,7 @@ void assemble_cells_matrix( /// @param[in] perms Facet permutation integer. Empty if facet /// permutations are not required. template -void assemble_exterior_entities( +void assemble_entities_over_cells( la::MatSet auto mat_set, mdspan2_t x_dofmap, md::mdspan, md::extents> @@ -695,7 +695,7 @@ void assemble_matrix( std::span e1 = a.domain_arg(itg_type, 1, i, 0); mdspanx2_t entities1(e1.data(), e1.size() / 2, 2); assert((entities.size() / 2) * cstride == coeffs.size()); - impl::assemble_exterior_entities( + impl::assemble_entities_over_cells( mat_set, x_dofmap, x, entities, {dofs0, bs0, entities0}, P0, {dofs1, bs1, entities1}, P1T, bc0, bc1, fn, md::mdspan(coeffs.data(), entities.extent(0), cstride), constants, diff --git a/cpp/dolfinx/fem/assemble_scalar_impl.h b/cpp/dolfinx/fem/assemble_scalar_impl.h index 418f72b691..aafaeb34fa 100644 --- a/cpp/dolfinx/fem/assemble_scalar_impl.h +++ b/cpp/dolfinx/fem/assemble_scalar_impl.h @@ -55,9 +55,9 @@ T assemble_cells(mdspan2_t x_dofmap, return value; } -/// Execute kernel over exterior entities and accumulate result +/// Execute kernel over a set of entities and accumulate result template -T assemble_exterior_entities( +T assemble_entities_over_cells( mdspan2_t x_dofmap, md::mdspan, md::extents> @@ -243,7 +243,7 @@ T assemble_scalar( constexpr std::size_t shape1 = 2 * num_adjacent_cells; assert((entities.size() / 2) * cstride == coeffs.size()); - value += impl::assemble_exterior_entities( + value += impl::assemble_entities_over_cells( x_dofmap, x, md::mdspan>( diff --git a/cpp/dolfinx/fem/assemble_vector_impl.h b/cpp/dolfinx/fem/assemble_vector_impl.h index 9ab29f87bf..55e6024980 100644 --- a/cpp/dolfinx/fem/assemble_vector_impl.h +++ b/cpp/dolfinx/fem/assemble_vector_impl.h @@ -718,7 +718,7 @@ void assemble_cells( } } -/// @brief Execute kernel over exterior entities and accumulate result in +/// @brief Execute kernel over a set of entities and accumulate result in /// vector. /// /// @tparam T Scalar type. @@ -746,7 +746,7 @@ void assemble_cells( template ::value_type> requires std::is_same_v::value_type, T> -void assemble_exterior_entities( +void assemble_entities_over_cells( fem::DofTransformKernel auto P0, V&& b, mdspan2_t x_dofmap, md::mdspan, md::extents> @@ -1392,14 +1392,14 @@ void assemble_vector( assert((entities.size() / 2) * cstride == coeffs.size()); if (bs == 1) { - impl::assemble_exterior_entities<1>( + impl::assemble_entities_over_cells<1>( P0, b, x_dofmap, x, entities, {dofs, bs, entities1}, fn, constants, md::mdspan(coeffs.data(), entities.extent(0), cstride), cell_info0, perms); } else if (bs == 3) { - impl::assemble_exterior_entities<3>( + impl::assemble_entities_over_cells<3>( P0, b, x_dofmap, x, entities, {dofs, bs, entities1}, fn, constants, md::mdspan(coeffs.data(), entities.size() / 2, cstride), @@ -1407,7 +1407,7 @@ void assemble_vector( } else { - impl::assemble_exterior_entities( + impl::assemble_entities_over_cells( P0, b, x_dofmap, x, entities, {dofs, bs, entities1}, fn, constants, md::mdspan(coeffs.data(), entities.size() / 2, cstride), From d02191524099f2d209c2d4ec6db575cfa1547954 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20Schartum=20Dokken?= Date: Tue, 9 Sep 2025 19:42:16 +0200 Subject: [PATCH 19/37] Apply suggestions from code review Co-authored-by: Garth N. Wells --- cpp/dolfinx/fem/assemble_matrix_impl.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cpp/dolfinx/fem/assemble_matrix_impl.h b/cpp/dolfinx/fem/assemble_matrix_impl.h index 6cbf865a70..f6b22b3ca6 100644 --- a/cpp/dolfinx/fem/assemble_matrix_impl.h +++ b/cpp/dolfinx/fem/assemble_matrix_impl.h @@ -156,7 +156,7 @@ void assemble_cells_matrix( } } -/// @brief Execute kernel over a set of entities and accumulate result in +/// @brief Execute kernel over entities of codimension > 1 and accumulate result in /// a matrix. /// /// @tparam T Matrix/form scalar type. From 5e986917d31156ba250f7c91b7dadfb612945574 Mon Sep 17 00:00:00 2001 From: jorgensd Date: Tue, 9 Sep 2025 17:56:10 +0000 Subject: [PATCH 20/37] Rename "IntegralType::exterior_facet" to "IntegralType::facet" + address other comments by Garth --- cpp/demo/mixed_poisson/main.cpp | 4 +- cpp/dolfinx/fem/Form.h | 8 +-- cpp/dolfinx/fem/assemble_matrix_impl.h | 37 +++++------ cpp/dolfinx/fem/assemble_scalar_impl.h | 35 +++++------ cpp/dolfinx/fem/assemble_vector_impl.h | 61 ++++++++----------- cpp/dolfinx/fem/pack.h | 2 +- cpp/dolfinx/fem/utils.cpp | 4 +- cpp/dolfinx/fem/utils.h | 19 +++--- python/dolfinx/fem/forms.py | 4 +- python/dolfinx/wrappers/fem.cpp | 4 +- .../test_assemble_mesh_independent_form.py | 4 +- python/test/unit/fem/test_assemble_submesh.py | 6 +- python/test/unit/fem/test_assembler.py | 8 +-- 13 files changed, 95 insertions(+), 101 deletions(-) diff --git a/cpp/demo/mixed_poisson/main.cpp b/cpp/demo/mixed_poisson/main.cpp index 65eeb8fa39..ed8e638bc2 100644 --- a/cpp/demo/mixed_poisson/main.cpp +++ b/cpp/demo/mixed_poisson/main.cpp @@ -271,7 +271,7 @@ int main(int argc, char* argv[]) // (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); + fem::IntegralType::facet, *mesh->topology(), dfacets); // Create data structure for the `ds(1)` integration domain in form // (see the UFL file). It is for en exterior facet integral (the key @@ -281,7 +281,7 @@ int main(int argc, char* argv[]) std::map< fem::IntegralType, std::vector>>> - subdomain_data{{fem::IntegralType::exterior_facet, {{1, domains}}}}; + subdomain_data{{fem::IntegralType::facet, {{1, domains}}}}; // Since we are doing a `ds(1)` integral on mesh and `u0` is defined // on the `submesh`, our form involves more than one mesh. The mesh diff --git a/cpp/dolfinx/fem/Form.h b/cpp/dolfinx/fem/Form.h index 7787d75211..b23598d178 100644 --- a/cpp/dolfinx/fem/Form.h +++ b/cpp/dolfinx/fem/Form.h @@ -38,7 +38,7 @@ class Function; enum class IntegralType : std::int8_t { cell = 0, ///< Cell - exterior_facet = 1, ///< Exterior facet + facet = 1, ///< Facet interior_facet = 2, ///< Interior facet vertex = 3, ///< Vertex ridge = 4 ///< Ridge @@ -277,7 +277,7 @@ class Form std::vector e; if (type == IntegralType::cell) e = emap.sub_topology_to_topology(itg.entities, inverse); - else if (type == IntegralType::exterior_facet + else if (type == IntegralType::facet or type == IntegralType::interior_facet) { const mesh::Topology topology = *_mesh->topology(); @@ -320,7 +320,7 @@ class Form std::vector e; if (type == IntegralType::cell) e = emap.sub_topology_to_topology(integral.entities, inverse); - else if (type == IntegralType::exterior_facet + else if (type == IntegralType::facet or type == IntegralType::interior_facet) { const mesh::Topology topology = *_mesh->topology(); @@ -456,7 +456,7 @@ class Form /// integrated over by a given integral (kernel). /// /// - For IntegralType::cell, returns a list of cell indices. - /// - For IntegralType::exterior_facet, returns a list with shape + /// - For IntegralType::facet, returns a list with shape /// `(num_facets, 2)`, where `[cell_index, 0]` is the cell index and /// `[cell_index, 1]` is the local facet index relative to the cell. /// - For IntegralType::interior_facet the shape is `(num_facets, 4)`, diff --git a/cpp/dolfinx/fem/assemble_matrix_impl.h b/cpp/dolfinx/fem/assemble_matrix_impl.h index f6b22b3ca6..66c0c12bb2 100644 --- a/cpp/dolfinx/fem/assemble_matrix_impl.h +++ b/cpp/dolfinx/fem/assemble_matrix_impl.h @@ -156,8 +156,16 @@ void assemble_cells_matrix( } } -/// @brief Execute kernel over entities of codimension > 1 and accumulate result in -/// a matrix. +/// @brief Execute kernel over entities of codimension > 1 and accumulate result +/// in a matrix. +/// +/// Each entity is represented by (i) a cell that the entity is attached to +/// and (ii) the local index of the entity with respect to the cell. The +/// kernel is executed for each entity. The kernel can access data +/// (e.g., coefficients, basis functions) associated with the attached cell. +/// However, entities may be attached to more than one cell. This function +/// therefore computes 'one-sided' integrals, i.e. evaluates integrals as seen +/// from cell used to define the entity. /// /// @tparam T Matrix/form scalar type. /// @param[in] mat_set Function that accumulates computed entries into a @@ -191,7 +199,7 @@ void assemble_cells_matrix( /// @param[in] perms Facet permutation integer. Empty if facet /// permutations are not required. template -void assemble_entities_over_cells( +void assemble_entities( la::MatSet auto mat_set, mdspan2_t x_dofmap, md::mdspan, md::extents> @@ -656,21 +664,14 @@ void assemble_matrix( cell_info0, cell_info1, facet_perms); } - for (fem::IntegralType itg_type : - {fem::IntegralType::exterior_facet, fem::IntegralType::vertex, - fem::IntegralType::ridge}) + for (auto itg_type : {fem::IntegralType::facet, fem::IntegralType::vertex, + fem::IntegralType::ridge}) { - md::mdspan> perms; - switch (itg_type) - { - case fem::IntegralType::exterior_facet: - { - perms = facet_perms; - break; - } - default: - break; - } + md::mdspan> perms + = itg_type == fem::IntegralType::facet + ? facet_perms + : md::mdspan>{}; for (int i = 0; i < a.num_integrals(itg_type, cell_type_idx); ++i) { @@ -695,7 +696,7 @@ void assemble_matrix( std::span e1 = a.domain_arg(itg_type, 1, i, 0); mdspanx2_t entities1(e1.data(), e1.size() / 2, 2); assert((entities.size() / 2) * cstride == coeffs.size()); - impl::assemble_entities_over_cells( + impl::assemble_entities( mat_set, x_dofmap, x, entities, {dofs0, bs0, entities0}, P0, {dofs1, bs1, entities1}, P1T, bc0, bc1, fn, md::mdspan(coeffs.data(), entities.extent(0), cstride), constants, diff --git a/cpp/dolfinx/fem/assemble_scalar_impl.h b/cpp/dolfinx/fem/assemble_scalar_impl.h index aafaeb34fa..892ec46b6c 100644 --- a/cpp/dolfinx/fem/assemble_scalar_impl.h +++ b/cpp/dolfinx/fem/assemble_scalar_impl.h @@ -55,9 +55,18 @@ T assemble_cells(mdspan2_t x_dofmap, return value; } -/// Execute kernel over a set of entities and accumulate result +/// @brief Execute kernel over entities of codimension > 1 and accumulate result +/// in a scalar. +/// +/// Each entity is represented by (i) a cell that the entity is attached to +/// and (ii) the local index of the entity with respect to the cell. The +/// kernel is executed for each entity. The kernel can access data +/// (e.g., coefficients, basis functions) associated with the attached cell. +/// However, entities may be attached to more than one cell. This function +/// therefore computes 'one-sided' integrals, i.e. evaluates integrals as seen +/// from cell used to define the entity. template -T assemble_entities_over_cells( +T assemble_entities( mdspan2_t x_dofmap, md::mdspan, md::extents> @@ -214,21 +223,13 @@ T assemble_scalar( facet_perms, cdofs_b); } - for (fem::IntegralType itg_type : - {fem::IntegralType::exterior_facet, fem::IntegralType::vertex, - fem::IntegralType::ridge}) + for (auto itg_type : {fem::IntegralType::facet, fem::IntegralType::vertex, + fem::IntegralType::ridge}) { - md::mdspan> perms; - switch (itg_type) - { - case fem::IntegralType::exterior_facet: - { - perms = facet_perms; - break; - } - default: - break; - } + md::mdspan> perms + = (itg_type == fem::IntegralType::facet) + ? facet_perms + : md::mdspan>{}; for (int i = 0; i < M.num_integrals(itg_type, 0); ++i) { @@ -243,7 +244,7 @@ T assemble_scalar( constexpr std::size_t shape1 = 2 * num_adjacent_cells; assert((entities.size() / 2) * cstride == coeffs.size()); - value += impl::assemble_entities_over_cells( + value += impl::assemble_entities( x_dofmap, x, md::mdspan>( diff --git a/cpp/dolfinx/fem/assemble_vector_impl.h b/cpp/dolfinx/fem/assemble_vector_impl.h index 55e6024980..d6d52a2ccf 100644 --- a/cpp/dolfinx/fem/assemble_vector_impl.h +++ b/cpp/dolfinx/fem/assemble_vector_impl.h @@ -718,8 +718,16 @@ void assemble_cells( } } -/// @brief Execute kernel over a set of entities and accumulate result in -/// vector. +/// @brief Execute kernel over entities of codimension > 1 and accumulate result +/// in a matrix. +/// +/// Each entity is represented by (i) a cell that the entity is attached to +/// and (ii) the local index of the entity with respect to the cell. The +/// kernel is executed for each entity. The kernel can access data +/// (e.g., coefficients, basis functions) associated with the attached cell. +/// However, entities may be attached to more than one cell. This function +/// therefore computes 'one-sided' integrals, i.e. evaluates integrals as seen +/// from cell used to define the entity. /// /// @tparam T Scalar type. /// @tparam _bs The block size of the form test function dof map. If @@ -746,7 +754,7 @@ void assemble_cells( template ::value_type> requires std::is_same_v::value_type, T> -void assemble_entities_over_cells( +void assemble_entities( fem::DofTransformKernel auto P0, V&& b, mdspan2_t x_dofmap, md::mdspan, md::extents> @@ -1084,21 +1092,13 @@ void lift_bc(V&& b, const Form& a, mdspan2_t x_dofmap, cell_info1, bc_values1, bc_markers1, x0, alpha, facet_perms); } - for (fem::IntegralType itg_type : - {fem::IntegralType::exterior_facet, fem::IntegralType::vertex, - fem::IntegralType::ridge}) + for (auto itg_type : {fem::IntegralType::facet, fem::IntegralType::vertex, + fem::IntegralType::ridge}) { - md::mdspan> perms; - switch (itg_type) - { - case fem::IntegralType::exterior_facet: - { - perms = facet_perms; - break; - } - default: - break; - } + md::mdspan> perms + = (itg_type == fem::IntegralType::facet) + ? facet_perms + : md::mdspan>{}; for (int i = 0; i < a.num_integrals(itg_type, 0); ++i) { @@ -1365,21 +1365,14 @@ void assemble_vector( } } - for (fem::IntegralType itg_type : - {fem::IntegralType::exterior_facet, fem::IntegralType::vertex, - fem::IntegralType::ridge}) + for (auto itg_type : {fem::IntegralType::facet, fem::IntegralType::vertex, + fem::IntegralType::ridge}) { - md::mdspan> perms; - switch (itg_type) - { - case fem::IntegralType::exterior_facet: - { - perms = facet_perms; - break; - } - default: - break; - } + md::mdspan> perms + = (itg_type == fem::IntegralType::facet) + ? facet_perms + : md::mdspan>{}; for (int i = 0; i < L.num_integrals(itg_type, 0); ++i) { auto fn = L.kernel(itg_type, i, 0); @@ -1392,14 +1385,14 @@ void assemble_vector( assert((entities.size() / 2) * cstride == coeffs.size()); if (bs == 1) { - impl::assemble_entities_over_cells<1>( + impl::assemble_entities<1>( P0, b, x_dofmap, x, entities, {dofs, bs, entities1}, fn, constants, md::mdspan(coeffs.data(), entities.extent(0), cstride), cell_info0, perms); } else if (bs == 3) { - impl::assemble_entities_over_cells<3>( + impl::assemble_entities<3>( P0, b, x_dofmap, x, entities, {dofs, bs, entities1}, fn, constants, md::mdspan(coeffs.data(), entities.size() / 2, cstride), @@ -1407,7 +1400,7 @@ void assemble_vector( } else { - impl::assemble_entities_over_cells( + impl::assemble_entities( P0, b, x_dofmap, x, entities, {dofs, bs, entities1}, fn, constants, md::mdspan(coeffs.data(), entities.size() / 2, cstride), diff --git a/cpp/dolfinx/fem/pack.h b/cpp/dolfinx/fem/pack.h index 65777820ee..bae72ee375 100644 --- a/cpp/dolfinx/fem/pack.h +++ b/cpp/dolfinx/fem/pack.h @@ -297,7 +297,7 @@ void pack_coefficients(const Form& form, } break; } - case IntegralType::exterior_facet: + case IntegralType::facet: case IntegralType::vertex: case IntegralType::ridge: { diff --git a/cpp/dolfinx/fem/utils.cpp b/cpp/dolfinx/fem/utils.cpp index 0aadd5fec6..ff47af04c8 100644 --- a/cpp/dolfinx/fem/utils.cpp +++ b/cpp/dolfinx/fem/utils.cpp @@ -145,7 +145,7 @@ fem::compute_integration_domains(fem::IntegralType integral_type, case IntegralType::cell: dim = tdim; break; - case IntegralType::exterior_facet: + case IntegralType::facet: dim = tdim - 1; break; case IntegralType::interior_facet: @@ -235,7 +235,7 @@ fem::compute_integration_domains(fem::IntegralType integral_type, } break; } - case IntegralType::exterior_facet: + case IntegralType::facet: case IntegralType::vertex: case IntegralType::ridge: { diff --git a/cpp/dolfinx/fem/utils.h b/cpp/dolfinx/fem/utils.h index e2e20c70ad..00fc135aa7 100644 --- a/cpp/dolfinx/fem/utils.h +++ b/cpp/dolfinx/fem/utils.h @@ -146,13 +146,13 @@ get_cell_entity_pairs(std::int32_t e, std::span cells, /// @param[in] entities List of mesh entities. Depending on the `IntegralType` /// these are associated with different entities: /// `IntegralType::cell`: cells -/// `IntegralType::exterior_facet`: facets +/// `IntegralType::facet`: facets /// `IntegralType::interior_facet`: facets /// `IntegralType::vertex`: vertices /// @return List of integration entity data, depending on the `IntegralType` the /// data per entity has different layouts /// `IntegralType::cell`: cell -/// `IntegralType::exterior_facet`: (cell, local_facet) +/// `IntegralType::facet`: (cell, local_facet) /// `IntegralType::interior_facet`: (cell, local_facet) /// `IntegralType::vertex`: (cell, local_vertex) std::vector @@ -239,7 +239,7 @@ void build_sparsity_pattern(la::SparsityPattern& pattern, const Form& a) const std::set types = a.integral_types(); if (types.find(IntegralType::interior_facet) != types.end() - or types.find(IntegralType::exterior_facet) != types.end()) + or types.find(IntegralType::facet) != types.end()) { // FIXME: cleanup these calls? Some of the happen internally again. int tdim = mesh->topology()->dim(); @@ -290,7 +290,7 @@ void build_sparsity_pattern(la::SparsityPattern& pattern, const Form& a) {{dofmaps[0], dofmaps[1]}}); } break; - case IntegralType::exterior_facet: + case IntegralType::facet: case IntegralType::ridge: case IntegralType::vertex: for (int i = 0; i < a.num_integrals(type, cell_type_idx); ++i) @@ -504,8 +504,7 @@ Form create_form_factory( } // Create facets, if required - if (num_integrals_type[exterior_facet] > 0 - or num_integrals_type[interior_facet] > 0) + if (num_integrals_type[facet] > 0 or num_integrals_type[interior_facet] > 0) { mesh->topology_mutable()->create_entities(tdim - 1); mesh->topology_mutable()->create_connectivity(tdim - 1, tdim); @@ -700,13 +699,13 @@ Form create_form_factory( // Attach exterior entity integrals { - for (IntegralType itg_type : {IntegralType::exterior_facet, - IntegralType::vertex, IntegralType::ridge}) + for (IntegralType itg_type : + {IntegralType::facet, IntegralType::vertex, IntegralType::ridge}) { std::size_t dim; switch (itg_type) { - case IntegralType::exterior_facet: + case IntegralType::facet: { dim = tdim - 1; break; @@ -728,7 +727,7 @@ Form create_form_factory( auto get_default_integration_entities = [dim](const mesh::Topology& topology, IntegralType itg_type) { - if (itg_type == IntegralType::exterior_facet) + if (itg_type == IntegralType::facet) { // Integrate over all owned exterior facets return mesh::exterior_facet_indices(topology); diff --git a/python/dolfinx/fem/forms.py b/python/dolfinx/fem/forms.py index bcc1cf9908..325876b756 100644 --- a/python/dolfinx/fem/forms.py +++ b/python/dolfinx/fem/forms.py @@ -157,7 +157,7 @@ def get_integration_domains( else: domains = [] if not isinstance(subdomain, list): - if integral_type in (IntegralType.exterior_facet, IntegralType.interior_facet): + if integral_type in (IntegralType.facet, IntegralType.interior_facet): tdim = subdomain.topology.dim subdomain._cpp_object.topology.create_connectivity(tdim - 1, tdim) subdomain._cpp_object.topology.create_connectivity(tdim, tdim - 1) @@ -222,7 +222,7 @@ def form_cpp_class( _ufl_to_dolfinx_domain = { "cell": IntegralType.cell, - "exterior_facet": IntegralType.exterior_facet, + "facet": IntegralType.facet, "interior_facet": IntegralType.interior_facet, "vertex": IntegralType.vertex, "ridge": IntegralType.ridge, diff --git a/python/dolfinx/wrappers/fem.cpp b/python/dolfinx/wrappers/fem.cpp index ee55f9f368..bc18235163 100644 --- a/python/dolfinx/wrappers/fem.cpp +++ b/python/dolfinx/wrappers/fem.cpp @@ -772,7 +772,7 @@ void declare_form(nb::module_& m, std::string type) case dolfinx::fem::IntegralType::cell: return nb::ndarray(_d.data(), {_d.size()}); - case dolfinx::fem::IntegralType::exterior_facet: + case dolfinx::fem::IntegralType::facet: { return nb::ndarray( _d.data(), {_d.size() / 2, 2}); @@ -1269,7 +1269,7 @@ void fem(nb::module_& m) nb::enum_(m, "_IntegralType") .value("cell", dolfinx::fem::IntegralType::cell, "cell integral") - .value("exterior_facet", dolfinx::fem::IntegralType::exterior_facet, + .value("facet", dolfinx::fem::IntegralType::facet, "exterior facet integral") .value("interior_facet", dolfinx::fem::IntegralType::interior_facet, "exterior facet integral") diff --git a/python/test/unit/fem/test_assemble_mesh_independent_form.py b/python/test/unit/fem/test_assemble_mesh_independent_form.py index bfff112d93..5d20745bd3 100644 --- a/python/test/unit/fem/test_assemble_mesh_independent_form.py +++ b/python/test/unit/fem/test_assemble_mesh_independent_form.py @@ -122,9 +122,9 @@ def g(x): wh.interpolate(g) facet_entities = dolfinx.fem.compute_integration_domains( - dolfinx.fem.IntegralType.exterior_facet, mesh.topology, facets + dolfinx.fem.IntegralType.facet, mesh.topology, facets ) - subdomains = {dolfinx.fem.IntegralType.exterior_facet: [(subdomain_id, facet_entities)]} + subdomains = {dolfinx.fem.IntegralType.facet: [(subdomain_id, facet_entities)]} form = dolfinx.fem.create_form( compiled_form, [Vh], mesh, subdomains, {w: wh}, {}, [entity_map] diff --git a/python/test/unit/fem/test_assemble_submesh.py b/python/test/unit/fem/test_assemble_submesh.py index 0334eb7087..e56eb7b5cb 100644 --- a/python/test/unit/fem/test_assemble_submesh.py +++ b/python/test/unit/fem/test_assemble_submesh.py @@ -158,7 +158,7 @@ def interior_marker(x): def a_ufl(u, v, f, g, measure): "Helper function to create a UFL bilinear form. The form depends on the integral type" - if measure.integral_type() == "cell" or measure.integral_type() == "exterior_facet": + if measure.integral_type() == "cell" or measure.integral_type() == "facet": return ufl.inner(f * g * u, v) * measure else: assert measure.integral_type() == "interior_facet" @@ -167,7 +167,7 @@ def a_ufl(u, v, f, g, measure): def L_ufl(v, f, g, measure): "Helper function to create a UFL linear form. The form depends on the integral type" - if measure.integral_type() == "cell" or measure.integral_type() == "exterior_facet": + if measure.integral_type() == "cell" or measure.integral_type() == "facet": return ufl.inner(f * g, v) * measure else: assert measure.integral_type() == "interior_facet" @@ -175,7 +175,7 @@ def L_ufl(v, f, g, measure): def M_ufl(f, g, measure): - if measure.integral_type() == "cell" or measure.integral_type() == "exterior_facet": + if measure.integral_type() == "cell" or measure.integral_type() == "facet": return f * g * measure else: assert measure.integral_type() == "interior_facet" diff --git a/python/test/unit/fem/test_assembler.py b/python/test/unit/fem/test_assembler.py index 82b86324c4..f1824d0654 100644 --- a/python/test/unit/fem/test_assembler.py +++ b/python/test/unit/fem/test_assembler.py @@ -1737,8 +1737,8 @@ def check_vertex_integral_against_sum(form, vertices, weighted=False): # b) With create_form vertices = np.arange(num_vertices) - fem.compute_integration_domains(fem.IntegralType.exterior_facet, msh.topology, vertices) - subdomains = {fem.IntegralType.exterior_facet: [(0, cell_vertex_pairs)]} + fem.compute_integration_domains(fem.IntegralType.facet, msh.topology, vertices) + subdomains = {fem.IntegralType.facet: [(0, cell_vertex_pairs)]} compiled_form = fem.compile_form( comm, x[0] * ufl.dP, form_compiler_options={"scalar_type": dtype} @@ -1877,8 +1877,8 @@ def check_vertex_integral_against_sum(form, vertices, weighted=False): # b) With create_form vertices = np.arange(num_vertices) - fem.compute_integration_domains(fem.IntegralType.exterior_facet, msh.topology, vertices) - subdomains = {fem.IntegralType.exterior_facet: [(0, cell_vertex_pairs)]} + fem.compute_integration_domains(fem.IntegralType.facet, msh.topology, vertices) + subdomains = {fem.IntegralType.facet: [(0, cell_vertex_pairs)]} compiled_form = fem.compile_form( comm, x[0] * v * ufl.dP, form_compiler_options={"scalar_type": dtype} From 796da8fe503767ffa7b9961a7a061aead711e3cb Mon Sep 17 00:00:00 2001 From: jorgensd Date: Tue, 9 Sep 2025 18:05:21 +0000 Subject: [PATCH 21/37] Fix spacing --- cpp/dolfinx/fem/utils.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/cpp/dolfinx/fem/utils.h b/cpp/dolfinx/fem/utils.h index 00fc135aa7..389259e4c1 100644 --- a/cpp/dolfinx/fem/utils.h +++ b/cpp/dolfinx/fem/utils.h @@ -146,13 +146,13 @@ get_cell_entity_pairs(std::int32_t e, std::span cells, /// @param[in] entities List of mesh entities. Depending on the `IntegralType` /// these are associated with different entities: /// `IntegralType::cell`: cells -/// `IntegralType::facet`: facets +/// `IntegralType::facet`: facets /// `IntegralType::interior_facet`: facets /// `IntegralType::vertex`: vertices /// @return List of integration entity data, depending on the `IntegralType` the /// data per entity has different layouts /// `IntegralType::cell`: cell -/// `IntegralType::facet`: (cell, local_facet) +/// `IntegralType::facet`: (cell, local_facet) /// `IntegralType::interior_facet`: (cell, local_facet) /// `IntegralType::vertex`: (cell, local_vertex) std::vector From 82e2c89eb5dc75769289c37fc70f28dc58312aa1 Mon Sep 17 00:00:00 2001 From: jorgensd Date: Tue, 9 Sep 2025 18:10:09 +0000 Subject: [PATCH 22/37] Use renaming --- cpp/dolfinx/fem/utils.h | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/cpp/dolfinx/fem/utils.h b/cpp/dolfinx/fem/utils.h index 389259e4c1..2a036f4639 100644 --- a/cpp/dolfinx/fem/utils.h +++ b/cpp/dolfinx/fem/utils.h @@ -504,7 +504,9 @@ Form create_form_factory( } // Create facets, if required - if (num_integrals_type[facet] > 0 or num_integrals_type[interior_facet] > 0) + // NOTE: exterior_facet and interior_facet is declared in ufcx.h + if (num_integrals_type[exterior_facet] > 0 + or num_integrals_type[interior_facet] > 0) { mesh->topology_mutable()->create_entities(tdim - 1); mesh->topology_mutable()->create_connectivity(tdim - 1, tdim); From c9298ea7c7b65214434e7c59f259c31f9807ffb8 Mon Sep 17 00:00:00 2001 From: jorgensd Date: Tue, 9 Sep 2025 18:21:11 +0000 Subject: [PATCH 23/37] Simplify --- cpp/dolfinx/fem/assemble_scalar_impl.h | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/cpp/dolfinx/fem/assemble_scalar_impl.h b/cpp/dolfinx/fem/assemble_scalar_impl.h index 892ec46b6c..92555c4caa 100644 --- a/cpp/dolfinx/fem/assemble_scalar_impl.h +++ b/cpp/dolfinx/fem/assemble_scalar_impl.h @@ -239,18 +239,15 @@ T assemble_scalar( std::span entities = M.domain(itg_type, i, 0); - constexpr std::size_t num_adjacent_cells = 1; // Two values per each adj. cell (cell index and local entity index). - constexpr std::size_t shape1 = 2 * num_adjacent_cells; - assert((entities.size() / 2) * cstride == coeffs.size()); value += impl::assemble_entities( x_dofmap, x, md::mdspan>( - entities.data(), entities.size() / shape1, 2), + entities.data(), entities.size() / 2, 2), fn, constants, - md::mdspan(coeffs.data(), entities.size() / shape1, cstride), perms, + md::mdspan(coeffs.data(), entities.size() / 2, cstride), perms, cdofs_b); } } From 6a389dafc1e0a7a5ecbbf1e794457d2e16666d97 Mon Sep 17 00:00:00 2001 From: jorgensd Date: Tue, 9 Sep 2025 18:22:23 +0000 Subject: [PATCH 24/37] Rename lifting --- cpp/dolfinx/fem/assemble_vector_impl.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/cpp/dolfinx/fem/assemble_vector_impl.h b/cpp/dolfinx/fem/assemble_vector_impl.h index d6d52a2ccf..53a751ae1b 100644 --- a/cpp/dolfinx/fem/assemble_vector_impl.h +++ b/cpp/dolfinx/fem/assemble_vector_impl.h @@ -263,7 +263,7 @@ void _lift_bc_cells( template ::value_type> requires std::is_same_v::value_type, T> -void _lift_bc_exterior_entities( +void _lift_bc_entities( V&& b, mdspan2_t x_dofmap, md::mdspan, md::extents> @@ -1116,7 +1116,7 @@ void lift_bc(V&& b, const Form& a, mdspan2_t x_dofmap, std::span e1 = a.domain_arg(itg_type, 1, i, 0); mdspanx2_t entities1(e1.data(), e1.size() / 2, 2); assert(coeffs.size() == entities.extent(0) * cstride); - _lift_bc_exterior_entities( + _lift_bc_entities( b, x_dofmap, x, kernel, entities, {dofmap0, bs0, entities0}, P0, {dofmap1, bs1, entities1}, P1T, constants, md::mdspan(coeffs.data(), entities.extent(0), cstride), cell_info0, From f7cc45c21165cae740f359444ce07427a29ff6aa Mon Sep 17 00:00:00 2001 From: jorgensd Date: Tue, 9 Sep 2025 18:24:47 +0000 Subject: [PATCH 25/37] Remove outdated comment --- cpp/dolfinx/fem/utils.h | 1 - 1 file changed, 1 deletion(-) diff --git a/cpp/dolfinx/fem/utils.h b/cpp/dolfinx/fem/utils.h index 2a036f4639..9f0ec71f4b 100644 --- a/cpp/dolfinx/fem/utils.h +++ b/cpp/dolfinx/fem/utils.h @@ -737,7 +737,6 @@ Form create_form_factory( else { // Integrate over all owned entities - // Default vertex kernel operates on all (owned) vertices std::int32_t num_entities = topology.index_map(dim)->size_local(); std::vector entities(num_entities); std::iota(entities.begin(), entities.end(), 0); From a1ff6ffa6bb3dee901a1014de51d79398aebf246 Mon Sep 17 00:00:00 2001 From: jorgensd Date: Tue, 9 Sep 2025 18:28:50 +0000 Subject: [PATCH 26/37] Fix typo --- cpp/dolfinx/common/sort.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cpp/dolfinx/common/sort.h b/cpp/dolfinx/common/sort.h index 4ee49fd936..6aad81ebc8 100644 --- a/cpp/dolfinx/common/sort.h +++ b/cpp/dolfinx/common/sort.h @@ -67,7 +67,7 @@ struct __radix_sort /// @code /// std::array a{2, 3, 1}; /// std::array i{0, 1, 2}; - /// dolfixn::radix_sort(i, [&](auto i){ return a[i]; }); // yields i = {2, 0, + /// dolfinx::radix_sort(i, [&](auto i){ return a[i]; }); // yields i = {2, 0, /// 1} and a[i] = {1, 2, 3}; /// @endcode /// @tparam R Type of range to be sorted. From 11418439630ad69fb7abf9e2078dbd0c0e3f89c8 Mon Sep 17 00:00:00 2001 From: jorgensd Date: Tue, 9 Sep 2025 18:29:06 +0000 Subject: [PATCH 27/37] Explicit lambda function declaration --- cpp/dolfinx/fem/utils.h | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/cpp/dolfinx/fem/utils.h b/cpp/dolfinx/fem/utils.h index 9f0ec71f4b..e0590814bb 100644 --- a/cpp/dolfinx/fem/utils.h +++ b/cpp/dolfinx/fem/utils.h @@ -726,7 +726,9 @@ Form create_form_factory( throw std::runtime_error("Unsupported integral type"); } - auto get_default_integration_entities + const std::function(const mesh::Topology&, + IntegralType)> + get_default_integration_entities = [dim](const mesh::Topology& topology, IntegralType itg_type) { if (itg_type == IntegralType::facet) From 5965b12764101cb57254c7ff9d9459d9b15239c3 Mon Sep 17 00:00:00 2001 From: jorgensd Date: Tue, 9 Sep 2025 19:01:01 +0000 Subject: [PATCH 28/37] Revert IntegralType::facet backt to IntegralType::exterior_facet --- cpp/demo/mixed_poisson/main.cpp | 4 ++-- cpp/dolfinx/fem/Form.h | 8 ++++---- cpp/dolfinx/fem/assemble_matrix_impl.h | 6 +++--- cpp/dolfinx/fem/assemble_scalar_impl.h | 6 +++--- cpp/dolfinx/fem/assemble_vector_impl.h | 12 ++++++------ cpp/dolfinx/fem/pack.h | 2 +- cpp/dolfinx/fem/utils.cpp | 4 ++-- cpp/dolfinx/fem/utils.h | 16 ++++++++-------- python/dolfinx/wrappers/fem.cpp | 4 ++-- 9 files changed, 31 insertions(+), 31 deletions(-) diff --git a/cpp/demo/mixed_poisson/main.cpp b/cpp/demo/mixed_poisson/main.cpp index ed8e638bc2..65eeb8fa39 100644 --- a/cpp/demo/mixed_poisson/main.cpp +++ b/cpp/demo/mixed_poisson/main.cpp @@ -271,7 +271,7 @@ int main(int argc, char* argv[]) // (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::facet, *mesh->topology(), dfacets); + fem::IntegralType::exterior_facet, *mesh->topology(), dfacets); // Create data structure for the `ds(1)` integration domain in form // (see the UFL file). It is for en exterior facet integral (the key @@ -281,7 +281,7 @@ int main(int argc, char* argv[]) std::map< fem::IntegralType, std::vector>>> - subdomain_data{{fem::IntegralType::facet, {{1, domains}}}}; + 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`, our form involves more than one mesh. The mesh diff --git a/cpp/dolfinx/fem/Form.h b/cpp/dolfinx/fem/Form.h index b23598d178..6bbb426e5c 100644 --- a/cpp/dolfinx/fem/Form.h +++ b/cpp/dolfinx/fem/Form.h @@ -38,7 +38,7 @@ class Function; enum class IntegralType : std::int8_t { cell = 0, ///< Cell - facet = 1, ///< Facet + exterior_facet = 1, ///< Facet interior_facet = 2, ///< Interior facet vertex = 3, ///< Vertex ridge = 4 ///< Ridge @@ -277,7 +277,7 @@ class Form std::vector e; if (type == IntegralType::cell) e = emap.sub_topology_to_topology(itg.entities, inverse); - else if (type == IntegralType::facet + else if (type == IntegralType::exterior_facet or type == IntegralType::interior_facet) { const mesh::Topology topology = *_mesh->topology(); @@ -320,7 +320,7 @@ class Form std::vector e; if (type == IntegralType::cell) e = emap.sub_topology_to_topology(integral.entities, inverse); - else if (type == IntegralType::facet + else if (type == IntegralType::exterior_facet or type == IntegralType::interior_facet) { const mesh::Topology topology = *_mesh->topology(); @@ -456,7 +456,7 @@ class Form /// integrated over by a given integral (kernel). /// /// - For IntegralType::cell, returns a list of cell indices. - /// - For IntegralType::facet, returns a list with shape + /// - For IntegralType::exterior_facet, returns a list with shape /// `(num_facets, 2)`, where `[cell_index, 0]` is the cell index and /// `[cell_index, 1]` is the local facet index relative to the cell. /// - For IntegralType::interior_facet the shape is `(num_facets, 4)`, diff --git a/cpp/dolfinx/fem/assemble_matrix_impl.h b/cpp/dolfinx/fem/assemble_matrix_impl.h index 66c0c12bb2..9d98225898 100644 --- a/cpp/dolfinx/fem/assemble_matrix_impl.h +++ b/cpp/dolfinx/fem/assemble_matrix_impl.h @@ -664,11 +664,11 @@ void assemble_matrix( cell_info0, cell_info1, facet_perms); } - for (auto itg_type : {fem::IntegralType::facet, fem::IntegralType::vertex, - fem::IntegralType::ridge}) + for (auto itg_type : {fem::IntegralType::exterior_facet, + fem::IntegralType::vertex, fem::IntegralType::ridge}) { md::mdspan> perms - = itg_type == fem::IntegralType::facet + = itg_type == fem::IntegralType::exterior_facet ? facet_perms : md::mdspan>{}; diff --git a/cpp/dolfinx/fem/assemble_scalar_impl.h b/cpp/dolfinx/fem/assemble_scalar_impl.h index 92555c4caa..f810f92797 100644 --- a/cpp/dolfinx/fem/assemble_scalar_impl.h +++ b/cpp/dolfinx/fem/assemble_scalar_impl.h @@ -223,11 +223,11 @@ T assemble_scalar( facet_perms, cdofs_b); } - for (auto itg_type : {fem::IntegralType::facet, fem::IntegralType::vertex, - fem::IntegralType::ridge}) + for (auto itg_type : {fem::IntegralType::exterior_facet, + fem::IntegralType::vertex, fem::IntegralType::ridge}) { md::mdspan> perms - = (itg_type == fem::IntegralType::facet) + = (itg_type == fem::IntegralType::exterior_facet) ? facet_perms : md::mdspan>{}; diff --git a/cpp/dolfinx/fem/assemble_vector_impl.h b/cpp/dolfinx/fem/assemble_vector_impl.h index 53a751ae1b..1a73307f63 100644 --- a/cpp/dolfinx/fem/assemble_vector_impl.h +++ b/cpp/dolfinx/fem/assemble_vector_impl.h @@ -1092,11 +1092,11 @@ void lift_bc(V&& b, const Form& a, mdspan2_t x_dofmap, cell_info1, bc_values1, bc_markers1, x0, alpha, facet_perms); } - for (auto itg_type : {fem::IntegralType::facet, fem::IntegralType::vertex, - fem::IntegralType::ridge}) + for (auto itg_type : {fem::IntegralType::exterior_facet, + fem::IntegralType::vertex, fem::IntegralType::ridge}) { md::mdspan> perms - = (itg_type == fem::IntegralType::facet) + = (itg_type == fem::IntegralType::exterior_facet) ? facet_perms : md::mdspan>{}; @@ -1365,11 +1365,11 @@ void assemble_vector( } } - for (auto itg_type : {fem::IntegralType::facet, fem::IntegralType::vertex, - fem::IntegralType::ridge}) + for (auto itg_type : {fem::IntegralType::exterior_facet, + fem::IntegralType::vertex, fem::IntegralType::ridge}) { md::mdspan> perms - = (itg_type == fem::IntegralType::facet) + = (itg_type == fem::IntegralType::exterior_facet) ? facet_perms : md::mdspan>{}; diff --git a/cpp/dolfinx/fem/pack.h b/cpp/dolfinx/fem/pack.h index bae72ee375..65777820ee 100644 --- a/cpp/dolfinx/fem/pack.h +++ b/cpp/dolfinx/fem/pack.h @@ -297,7 +297,7 @@ void pack_coefficients(const Form& form, } break; } - case IntegralType::facet: + case IntegralType::exterior_facet: case IntegralType::vertex: case IntegralType::ridge: { diff --git a/cpp/dolfinx/fem/utils.cpp b/cpp/dolfinx/fem/utils.cpp index ff47af04c8..0aadd5fec6 100644 --- a/cpp/dolfinx/fem/utils.cpp +++ b/cpp/dolfinx/fem/utils.cpp @@ -145,7 +145,7 @@ fem::compute_integration_domains(fem::IntegralType integral_type, case IntegralType::cell: dim = tdim; break; - case IntegralType::facet: + case IntegralType::exterior_facet: dim = tdim - 1; break; case IntegralType::interior_facet: @@ -235,7 +235,7 @@ fem::compute_integration_domains(fem::IntegralType integral_type, } break; } - case IntegralType::facet: + case IntegralType::exterior_facet: case IntegralType::vertex: case IntegralType::ridge: { diff --git a/cpp/dolfinx/fem/utils.h b/cpp/dolfinx/fem/utils.h index e0590814bb..00fdc6211b 100644 --- a/cpp/dolfinx/fem/utils.h +++ b/cpp/dolfinx/fem/utils.h @@ -146,13 +146,13 @@ get_cell_entity_pairs(std::int32_t e, std::span cells, /// @param[in] entities List of mesh entities. Depending on the `IntegralType` /// these are associated with different entities: /// `IntegralType::cell`: cells -/// `IntegralType::facet`: facets +/// `IntegralType::exterior_facet`: facets /// `IntegralType::interior_facet`: facets /// `IntegralType::vertex`: vertices /// @return List of integration entity data, depending on the `IntegralType` the /// data per entity has different layouts /// `IntegralType::cell`: cell -/// `IntegralType::facet`: (cell, local_facet) +/// `IntegralType::exterior_facet`: (cell, local_facet) /// `IntegralType::interior_facet`: (cell, local_facet) /// `IntegralType::vertex`: (cell, local_vertex) std::vector @@ -239,7 +239,7 @@ void build_sparsity_pattern(la::SparsityPattern& pattern, const Form& a) const std::set types = a.integral_types(); if (types.find(IntegralType::interior_facet) != types.end() - or types.find(IntegralType::facet) != types.end()) + or types.find(IntegralType::exterior_facet) != types.end()) { // FIXME: cleanup these calls? Some of the happen internally again. int tdim = mesh->topology()->dim(); @@ -290,7 +290,7 @@ void build_sparsity_pattern(la::SparsityPattern& pattern, const Form& a) {{dofmaps[0], dofmaps[1]}}); } break; - case IntegralType::facet: + case IntegralType::exterior_facet: case IntegralType::ridge: case IntegralType::vertex: for (int i = 0; i < a.num_integrals(type, cell_type_idx); ++i) @@ -701,13 +701,13 @@ Form create_form_factory( // Attach exterior entity integrals { - for (IntegralType itg_type : - {IntegralType::facet, IntegralType::vertex, IntegralType::ridge}) + for (IntegralType itg_type : {IntegralType::exterior_facet, + IntegralType::vertex, IntegralType::ridge}) { std::size_t dim; switch (itg_type) { - case IntegralType::facet: + case IntegralType::exterior_facet: { dim = tdim - 1; break; @@ -731,7 +731,7 @@ Form create_form_factory( get_default_integration_entities = [dim](const mesh::Topology& topology, IntegralType itg_type) { - if (itg_type == IntegralType::facet) + if (itg_type == IntegralType::exterior_facet) { // Integrate over all owned exterior facets return mesh::exterior_facet_indices(topology); diff --git a/python/dolfinx/wrappers/fem.cpp b/python/dolfinx/wrappers/fem.cpp index bc18235163..5903efa15b 100644 --- a/python/dolfinx/wrappers/fem.cpp +++ b/python/dolfinx/wrappers/fem.cpp @@ -772,7 +772,7 @@ void declare_form(nb::module_& m, std::string type) case dolfinx::fem::IntegralType::cell: return nb::ndarray(_d.data(), {_d.size()}); - case dolfinx::fem::IntegralType::facet: + case dolfinx::fem::IntegralType::exterior_facet: { return nb::ndarray( _d.data(), {_d.size() / 2, 2}); @@ -1269,7 +1269,7 @@ void fem(nb::module_& m) nb::enum_(m, "_IntegralType") .value("cell", dolfinx::fem::IntegralType::cell, "cell integral") - .value("facet", dolfinx::fem::IntegralType::facet, + .value("facet", dolfinx::fem::IntegralType::exterior_facet, "exterior facet integral") .value("interior_facet", dolfinx::fem::IntegralType::interior_facet, "exterior facet integral") From daadfa10c3a3565f3e8a6022e3567fcae0d05c2b Mon Sep 17 00:00:00 2001 From: jorgensd Date: Tue, 9 Sep 2025 19:14:36 +0000 Subject: [PATCH 29/37] More reverting --- python/dolfinx/fem/forms.py | 4 ++-- python/dolfinx/wrappers/fem.cpp | 2 +- .../test/unit/fem/test_assemble_mesh_independent_form.py | 4 ++-- python/test/unit/fem/test_assembler.py | 8 ++++---- 4 files changed, 9 insertions(+), 9 deletions(-) diff --git a/python/dolfinx/fem/forms.py b/python/dolfinx/fem/forms.py index 325876b756..a9a03254ab 100644 --- a/python/dolfinx/fem/forms.py +++ b/python/dolfinx/fem/forms.py @@ -157,7 +157,7 @@ def get_integration_domains( else: domains = [] if not isinstance(subdomain, list): - if integral_type in (IntegralType.facet, IntegralType.interior_facet): + if integral_type in (IntegralType.exterior_facet, IntegralType.interior_facet): tdim = subdomain.topology.dim subdomain._cpp_object.topology.create_connectivity(tdim - 1, tdim) subdomain._cpp_object.topology.create_connectivity(tdim, tdim - 1) @@ -222,7 +222,7 @@ def form_cpp_class( _ufl_to_dolfinx_domain = { "cell": IntegralType.cell, - "facet": IntegralType.facet, + "facet": IntegralType.exterior_facet, "interior_facet": IntegralType.interior_facet, "vertex": IntegralType.vertex, "ridge": IntegralType.ridge, diff --git a/python/dolfinx/wrappers/fem.cpp b/python/dolfinx/wrappers/fem.cpp index 5903efa15b..ee55f9f368 100644 --- a/python/dolfinx/wrappers/fem.cpp +++ b/python/dolfinx/wrappers/fem.cpp @@ -1269,7 +1269,7 @@ void fem(nb::module_& m) nb::enum_(m, "_IntegralType") .value("cell", dolfinx::fem::IntegralType::cell, "cell integral") - .value("facet", dolfinx::fem::IntegralType::exterior_facet, + .value("exterior_facet", dolfinx::fem::IntegralType::exterior_facet, "exterior facet integral") .value("interior_facet", dolfinx::fem::IntegralType::interior_facet, "exterior facet integral") diff --git a/python/test/unit/fem/test_assemble_mesh_independent_form.py b/python/test/unit/fem/test_assemble_mesh_independent_form.py index 5d20745bd3..bfff112d93 100644 --- a/python/test/unit/fem/test_assemble_mesh_independent_form.py +++ b/python/test/unit/fem/test_assemble_mesh_independent_form.py @@ -122,9 +122,9 @@ def g(x): wh.interpolate(g) facet_entities = dolfinx.fem.compute_integration_domains( - dolfinx.fem.IntegralType.facet, mesh.topology, facets + dolfinx.fem.IntegralType.exterior_facet, mesh.topology, facets ) - subdomains = {dolfinx.fem.IntegralType.facet: [(subdomain_id, facet_entities)]} + subdomains = {dolfinx.fem.IntegralType.exterior_facet: [(subdomain_id, facet_entities)]} form = dolfinx.fem.create_form( compiled_form, [Vh], mesh, subdomains, {w: wh}, {}, [entity_map] diff --git a/python/test/unit/fem/test_assembler.py b/python/test/unit/fem/test_assembler.py index f1824d0654..82b86324c4 100644 --- a/python/test/unit/fem/test_assembler.py +++ b/python/test/unit/fem/test_assembler.py @@ -1737,8 +1737,8 @@ def check_vertex_integral_against_sum(form, vertices, weighted=False): # b) With create_form vertices = np.arange(num_vertices) - fem.compute_integration_domains(fem.IntegralType.facet, msh.topology, vertices) - subdomains = {fem.IntegralType.facet: [(0, cell_vertex_pairs)]} + fem.compute_integration_domains(fem.IntegralType.exterior_facet, msh.topology, vertices) + subdomains = {fem.IntegralType.exterior_facet: [(0, cell_vertex_pairs)]} compiled_form = fem.compile_form( comm, x[0] * ufl.dP, form_compiler_options={"scalar_type": dtype} @@ -1877,8 +1877,8 @@ def check_vertex_integral_against_sum(form, vertices, weighted=False): # b) With create_form vertices = np.arange(num_vertices) - fem.compute_integration_domains(fem.IntegralType.facet, msh.topology, vertices) - subdomains = {fem.IntegralType.facet: [(0, cell_vertex_pairs)]} + fem.compute_integration_domains(fem.IntegralType.exterior_facet, msh.topology, vertices) + subdomains = {fem.IntegralType.exterior_facet: [(0, cell_vertex_pairs)]} compiled_form = fem.compile_form( comm, x[0] * v * ufl.dP, form_compiler_options={"scalar_type": dtype} From 6662a079c79a6d4c333e6d1049a0b1b08ad27a36 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20Schartum=20Dokken?= Date: Tue, 9 Sep 2025 21:16:26 +0200 Subject: [PATCH 30/37] More reverting --- python/dolfinx/fem/forms.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/dolfinx/fem/forms.py b/python/dolfinx/fem/forms.py index a9a03254ab..bcc1cf9908 100644 --- a/python/dolfinx/fem/forms.py +++ b/python/dolfinx/fem/forms.py @@ -222,7 +222,7 @@ def form_cpp_class( _ufl_to_dolfinx_domain = { "cell": IntegralType.cell, - "facet": IntegralType.exterior_facet, + "exterior_facet": IntegralType.exterior_facet, "interior_facet": IntegralType.interior_facet, "vertex": IntegralType.vertex, "ridge": IntegralType.ridge, From 3b01e9322788052b4f7b7a104a87de46911f5cc7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20Schartum=20Dokken?= Date: Tue, 9 Sep 2025 21:19:24 +0200 Subject: [PATCH 31/37] Revert spacing --- cpp/dolfinx/fem/utils.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/cpp/dolfinx/fem/utils.h b/cpp/dolfinx/fem/utils.h index 00fdc6211b..9b100a9c08 100644 --- a/cpp/dolfinx/fem/utils.h +++ b/cpp/dolfinx/fem/utils.h @@ -146,13 +146,13 @@ get_cell_entity_pairs(std::int32_t e, std::span cells, /// @param[in] entities List of mesh entities. Depending on the `IntegralType` /// these are associated with different entities: /// `IntegralType::cell`: cells -/// `IntegralType::exterior_facet`: facets +/// `IntegralType::exterior_facet`: facets /// `IntegralType::interior_facet`: facets /// `IntegralType::vertex`: vertices /// @return List of integration entity data, depending on the `IntegralType` the /// data per entity has different layouts /// `IntegralType::cell`: cell -/// `IntegralType::exterior_facet`: (cell, local_facet) +/// `IntegralType::exterior_facet`: (cell, local_facet) /// `IntegralType::interior_facet`: (cell, local_facet) /// `IntegralType::vertex`: (cell, local_vertex) std::vector From 6f53f75c53b1756dfe6af2fc9b6035085bb0604f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20Schartum=20Dokken?= Date: Tue, 9 Sep 2025 21:19:52 +0200 Subject: [PATCH 32/37] Revert spacing --- cpp/dolfinx/fem/utils.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cpp/dolfinx/fem/utils.h b/cpp/dolfinx/fem/utils.h index 9b100a9c08..094cf3e55b 100644 --- a/cpp/dolfinx/fem/utils.h +++ b/cpp/dolfinx/fem/utils.h @@ -152,7 +152,7 @@ get_cell_entity_pairs(std::int32_t e, std::span cells, /// @return List of integration entity data, depending on the `IntegralType` the /// data per entity has different layouts /// `IntegralType::cell`: cell -/// `IntegralType::exterior_facet`: (cell, local_facet) +/// `IntegralType::exterior_facet`: (cell, local_facet) /// `IntegralType::interior_facet`: (cell, local_facet) /// `IntegralType::vertex`: (cell, local_vertex) std::vector From 93059d2287d3167c816c8adb69d85d5c2b25928a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20Schartum=20Dokken?= Date: Tue, 9 Sep 2025 21:21:46 +0200 Subject: [PATCH 33/37] Further reverting --- python/test/unit/fem/test_assemble_submesh.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/python/test/unit/fem/test_assemble_submesh.py b/python/test/unit/fem/test_assemble_submesh.py index e56eb7b5cb..61e649b9f9 100644 --- a/python/test/unit/fem/test_assemble_submesh.py +++ b/python/test/unit/fem/test_assemble_submesh.py @@ -158,7 +158,7 @@ def interior_marker(x): def a_ufl(u, v, f, g, measure): "Helper function to create a UFL bilinear form. The form depends on the integral type" - if measure.integral_type() == "cell" or measure.integral_type() == "facet": + if measure.integral_type() == "cell" or measure.integral_type() == "exterior_facet": return ufl.inner(f * g * u, v) * measure else: assert measure.integral_type() == "interior_facet" @@ -167,7 +167,7 @@ def a_ufl(u, v, f, g, measure): def L_ufl(v, f, g, measure): "Helper function to create a UFL linear form. The form depends on the integral type" - if measure.integral_type() == "cell" or measure.integral_type() == "facet": + if measure.integral_type() == "cell" or measure.integral_type() == "exterior_facet": return ufl.inner(f * g, v) * measure else: assert measure.integral_type() == "interior_facet" From 34a990ca4a8feb01c34bae937fb962e0e38026f8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20Schartum=20Dokken?= Date: Tue, 9 Sep 2025 21:47:59 +0200 Subject: [PATCH 34/37] Yet another revert --- python/test/unit/fem/test_assemble_submesh.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/test/unit/fem/test_assemble_submesh.py b/python/test/unit/fem/test_assemble_submesh.py index 61e649b9f9..0334eb7087 100644 --- a/python/test/unit/fem/test_assemble_submesh.py +++ b/python/test/unit/fem/test_assemble_submesh.py @@ -175,7 +175,7 @@ def L_ufl(v, f, g, measure): def M_ufl(f, g, measure): - if measure.integral_type() == "cell" or measure.integral_type() == "facet": + if measure.integral_type() == "cell" or measure.integral_type() == "exterior_facet": return f * g * measure else: assert measure.integral_type() == "interior_facet" From 347af6242656acb0a9c43246d26af3e0a69d4eaa Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20Schartum=20Dokken?= Date: Fri, 12 Sep 2025 17:30:44 +0200 Subject: [PATCH 35/37] Apply suggestions from code review Co-authored-by: Michal Habera --- cpp/dolfinx/fem/assemble_matrix_impl.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cpp/dolfinx/fem/assemble_matrix_impl.h b/cpp/dolfinx/fem/assemble_matrix_impl.h index 9d98225898..0afe4150f0 100644 --- a/cpp/dolfinx/fem/assemble_matrix_impl.h +++ b/cpp/dolfinx/fem/assemble_matrix_impl.h @@ -240,7 +240,7 @@ void assemble_entities( assert(entities1.size() == entities.size()); for (std::size_t f = 0; f < entities.extent(0); ++f) { - // Cell in the integration domain, local facet index relative to the + // Cell in the integration domain, local entity index relative to the // integration domain cell, and cells in the test and trial function // meshes std::int32_t cell = entities(f, 0); From deeaaf364f6a8c3ad2ccc863f07bab41c5b03fd9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20Schartum=20Dokken?= Date: Fri, 26 Sep 2025 09:39:46 +0200 Subject: [PATCH 36/37] Apply suggestions from code review MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Paul T. Kühner <56360279+schnellerhase@users.noreply.github.com> --- cpp/dolfinx/fem/assemble_matrix_impl.h | 4 ++-- cpp/dolfinx/fem/assemble_scalar_impl.h | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/cpp/dolfinx/fem/assemble_matrix_impl.h b/cpp/dolfinx/fem/assemble_matrix_impl.h index 0afe4150f0..3b6148f9e9 100644 --- a/cpp/dolfinx/fem/assemble_matrix_impl.h +++ b/cpp/dolfinx/fem/assemble_matrix_impl.h @@ -156,7 +156,7 @@ void assemble_cells_matrix( } } -/// @brief Execute kernel over entities of codimension > 1 and accumulate result +/// @brief Execute kernel over entities of codimension ≥ 1 and accumulate result /// in a matrix. /// /// Each entity is represented by (i) a cell that the entity is attached to @@ -196,7 +196,7 @@ void assemble_cells_matrix( /// function mesh. /// @param[in] cell_info1 Cell permutation information for the trial /// function mesh. -/// @param[in] perms Facet permutation integer. Empty if facet +/// @param[in] perms Entity permutation integer. Empty if entity /// permutations are not required. template void assemble_entities( diff --git a/cpp/dolfinx/fem/assemble_scalar_impl.h b/cpp/dolfinx/fem/assemble_scalar_impl.h index f810f92797..37166b46ea 100644 --- a/cpp/dolfinx/fem/assemble_scalar_impl.h +++ b/cpp/dolfinx/fem/assemble_scalar_impl.h @@ -55,7 +55,7 @@ T assemble_cells(mdspan2_t x_dofmap, return value; } -/// @brief Execute kernel over entities of codimension > 1 and accumulate result +/// @brief Execute kernel over entities of codimension ≥ 1 and accumulate result /// in a scalar. /// /// Each entity is represented by (i) a cell that the entity is attached to From 678463a3e7d262f30665dd1caf85522e7d47c82c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20Schartum=20Dokken?= Date: Fri, 26 Sep 2025 09:40:48 +0200 Subject: [PATCH 37/37] Apply suggestions from code review MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Paul T. Kühner <56360279+schnellerhase@users.noreply.github.com> --- cpp/dolfinx/fem/assemble_vector_impl.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cpp/dolfinx/fem/assemble_vector_impl.h b/cpp/dolfinx/fem/assemble_vector_impl.h index 1a73307f63..f57e646dd9 100644 --- a/cpp/dolfinx/fem/assemble_vector_impl.h +++ b/cpp/dolfinx/fem/assemble_vector_impl.h @@ -718,7 +718,7 @@ void assemble_cells( } } -/// @brief Execute kernel over entities of codimension > 1 and accumulate result +/// @brief Execute kernel over entities of codimension ≥ 1 and accumulate result /// in a matrix. /// /// Each entity is represented by (i) a cell that the entity is attached to