From 2d5a4c27e48b5fbe0ec013a9147d5aa69ae93122 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= Date: Wed, 10 Sep 2025 12:56:03 +0200 Subject: [PATCH 01/30] Change facet type --- cpp/dolfinx/fem/utils.h | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/cpp/dolfinx/fem/utils.h b/cpp/dolfinx/fem/utils.h index 22436c8e50d..48cab5be495 100644 --- a/cpp/dolfinx/fem/utils.h +++ b/cpp/dolfinx/fem/utils.h @@ -502,7 +502,7 @@ Form create_form_factory( } // Create facets, if required - if (num_integrals_type[exterior_facet] > 0 + if (num_integrals_type[facet] > 0 or num_integrals_type[interior_facet] > 0) { mesh->topology_mutable()->create_entities(tdim - 1); @@ -597,17 +597,17 @@ Form create_form_factory( std::vector default_facets_ext; { std::span ids(ufcx_forms[0].get().form_integral_ids - + integral_offsets[exterior_facet], - num_integrals_type[exterior_facet]); + + integral_offsets[facet], + num_integrals_type[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) + for (int i = 0; i < num_integrals_type[facet]; ++i) { const int id = ids[i]; ufcx_integral* integral - = ufcx_form.form_integrals[integral_offsets[exterior_facet] + i]; + = ufcx_form.form_integrals[integral_offsets[facet] + i]; assert(integral); check_geometry_hash(*integral, form_idx); From 0c572f4ca409dcb239e1f58c57e04d994a8a59f3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= Date: Wed, 10 Sep 2025 11:41:59 +0000 Subject: [PATCH 02/30] Change exeterior_facet to facet --- .github/workflows/fenicsx-refs.env | 3 ++- cpp/demo/mixed_poisson/main.cpp | 4 ++-- cpp/dolfinx/fem/Form.h | 8 +++---- cpp/dolfinx/fem/assemble_matrix_impl.h | 15 ++++++------ cpp/dolfinx/fem/assemble_scalar_impl.h | 9 ++++--- cpp/dolfinx/fem/assemble_vector_impl.h | 24 +++++++++---------- cpp/dolfinx/fem/pack.h | 8 +++---- cpp/dolfinx/fem/utils.cpp | 4 ++-- cpp/dolfinx/fem/utils.h | 17 +++++++------ python/dolfinx/fem/forms.py | 12 ++++++---- python/dolfinx/wrappers/fem.cpp | 5 ++-- .../test_assemble_mesh_independent_form.py | 4 ++-- python/test/unit/fem/test_assembler.py | 8 +++---- 13 files changed, 59 insertions(+), 62 deletions(-) diff --git a/.github/workflows/fenicsx-refs.env b/.github/workflows/fenicsx-refs.env index 22c7e44a7e7..5a671128ff2 100644 --- a/.github/workflows/fenicsx-refs.env +++ b/.github/workflows/fenicsx-refs.env @@ -3,4 +3,5 @@ basix_ref=main ufl_repository=FEniCS/ufl ufl_ref=main ffcx_repository=FEniCS/ffcx -ffcx_ref=main +ffcx_ref=dokken/integral-type + diff --git a/cpp/demo/mixed_poisson/main.cpp b/cpp/demo/mixed_poisson/main.cpp index 65eeb8fa395..ed8e638bc25 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 8a361dcf0c8..f00fb2a8b7d 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 }; @@ -276,7 +276,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(); @@ -319,7 +319,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(); @@ -455,7 +455,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 90351a05b7a..66b5a58883c 100644 --- a/cpp/dolfinx/fem/assemble_matrix_impl.h +++ b/cpp/dolfinx/fem/assemble_matrix_impl.h @@ -618,8 +618,8 @@ void assemble_matrix( num_facets_per_cell); } - for (int i = 0; - i < a.num_integrals(IntegralType::exterior_facet, cell_type_idx); ++i) + for (int i = 0; i < a.num_integrals(IntegralType::facet, cell_type_idx); + ++i) { if (num_cell_types > 1) { @@ -631,16 +631,15 @@ void assemble_matrix( = md::mdspan>; - auto fn = a.kernel(IntegralType::exterior_facet, i, 0); + auto fn = a.kernel(IntegralType::facet, i, 0); assert(fn); - auto& [coeffs, cstride] - = coefficients.at({IntegralType::exterior_facet, i}); + auto& [coeffs, cstride] = coefficients.at({IntegralType::facet, i}); - std::span f = a.domain(IntegralType::exterior_facet, i, 0); + std::span f = a.domain(IntegralType::facet, i, 0); mdspanx2_t facets(f.data(), f.size() / 2, 2); - std::span f0 = a.domain_arg(IntegralType::exterior_facet, 0, i, 0); + std::span f0 = a.domain_arg(IntegralType::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); + std::span f1 = a.domain_arg(IntegralType::facet, 1, i, 0); mdspanx2_t facets1(f1.data(), f1.size() / 2, 2); assert((facets.size() / 2) * cstride == coeffs.size()); impl::assemble_exterior_facets( diff --git a/cpp/dolfinx/fem/assemble_scalar_impl.h b/cpp/dolfinx/fem/assemble_scalar_impl.h index 1ade0a9792e..3ee1aa2ea94 100644 --- a/cpp/dolfinx/fem/assemble_scalar_impl.h +++ b/cpp/dolfinx/fem/assemble_scalar_impl.h @@ -226,14 +226,13 @@ T assemble_scalar( num_facets_per_cell); } - for (int i = 0; i < M.num_integrals(IntegralType::exterior_facet, 0); ++i) + for (int i = 0; i < M.num_integrals(IntegralType::facet, 0); ++i) { - auto fn = M.kernel(IntegralType::exterior_facet, i, 0); + auto fn = M.kernel(IntegralType::facet, i, 0); assert(fn); - auto& [coeffs, cstride] - = coefficients.at({IntegralType::exterior_facet, i}); + auto& [coeffs, cstride] = coefficients.at({IntegralType::facet, i}); - std::span facets = M.domain(IntegralType::exterior_facet, i, 0); + std::span facets = M.domain(IntegralType::facet, i, 0); // Two values per each adjacent cell (cell index and local facet // index) diff --git a/cpp/dolfinx/fem/assemble_vector_impl.h b/cpp/dolfinx/fem/assemble_vector_impl.h index 3e8c1c2a13d..87353e641ef 100644 --- a/cpp/dolfinx/fem/assemble_vector_impl.h +++ b/cpp/dolfinx/fem/assemble_vector_impl.h @@ -1145,21 +1145,20 @@ void lift_bc(V&& b, const Form& a, mdspan2_t x_dofmap, num_facets_per_cell); } - for (int i = 0; i < a.num_integrals(IntegralType::exterior_facet, 0); ++i) + for (int i = 0; i < a.num_integrals(IntegralType::facet, 0); ++i) { - auto kernel = a.kernel(IntegralType::exterior_facet, i, 0); + auto kernel = a.kernel(IntegralType::facet, i, 0); assert(kernel); - auto& [coeffs, cstride] - = coefficients.at({IntegralType::exterior_facet, i}); + auto& [coeffs, cstride] = coefficients.at({IntegralType::facet, i}); using mdspanx2_t = md::mdspan>; - std::span f = a.domain(IntegralType::exterior_facet, i, 0); + std::span f = a.domain(IntegralType::facet, i, 0); mdspanx2_t facets(f.data(), f.size() / 2, 2); - std::span f0 = a.domain_arg(IntegralType::exterior_facet, 0, i, 0); + std::span f0 = a.domain_arg(IntegralType::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); + std::span f1 = a.domain_arg(IntegralType::facet, 1, i, 0); mdspanx2_t facets1(f1.data(), f1.size() / 2, 2); assert(coeffs.size() == facets.extent(0) * cstride); _lift_bc_exterior_facets( @@ -1385,15 +1384,14 @@ void assemble_vector( = md::mdspan>; - for (int i = 0; i < L.num_integrals(IntegralType::exterior_facet, 0); ++i) + for (int i = 0; i < L.num_integrals(IntegralType::facet, 0); ++i) { - auto fn = L.kernel(IntegralType::exterior_facet, i, 0); + auto fn = L.kernel(IntegralType::facet, i, 0); assert(fn); - auto& [coeffs, cstride] - = coefficients.at({IntegralType::exterior_facet, i}); - std::span f = L.domain(IntegralType::exterior_facet, i, 0); + auto& [coeffs, cstride] = coefficients.at({IntegralType::facet, i}); + std::span f = L.domain(IntegralType::facet, i, 0); mdspanx2_t facets(f.data(), f.size() / 2, 2); - std::span f1 = L.domain_arg(IntegralType::exterior_facet, 0, i, 0); + std::span f1 = L.domain_arg(IntegralType::facet, 0, i, 0); mdspanx2_t facets1(f1.data(), f1.size() / 2, 2); assert((facets.size() / 2) * cstride == coeffs.size()); if (bs == 1) diff --git a/cpp/dolfinx/fem/pack.h b/cpp/dolfinx/fem/pack.h index 9b08e845d81..a943d872c90 100644 --- a/cpp/dolfinx/fem/pack.h +++ b/cpp/dolfinx/fem/pack.h @@ -179,7 +179,7 @@ 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 + if (integral_type == IntegralType::facet or integral_type == IntegralType::interior_facet) { num_entities /= 2; @@ -270,15 +270,15 @@ void pack_coefficients(const Form& form, } break; } - case IntegralType::exterior_facet: + case IntegralType::facet: { // Iterate over coefficients coefficients that are active in // exterior facet integrals - for (int coeff : form.active_coeffs(IntegralType::exterior_facet, id)) + for (int coeff : form.active_coeffs(IntegralType::facet, id)) { auto mesh = coefficients[coeff]->function_space()->mesh(); std::span facets_b - = form.domain_coeff(IntegralType::exterior_facet, id, coeff); + = form.domain_coeff(IntegralType::facet, id, coeff); md::mdspan> facets(facets_b.data(), facets_b.size() / 2, 2); diff --git a/cpp/dolfinx/fem/utils.cpp b/cpp/dolfinx/fem/utils.cpp index ec7a8501ae3..4b36d5dd18e 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: @@ -199,7 +199,7 @@ fem::compute_integration_domains(fem::IntegralType integral_type, entity_data.insert(entity_data.begin(), entities.begin(), entities.end()); break; } - case IntegralType::exterior_facet: + case IntegralType::facet: { auto [f_to_c, c_to_f] = get_connectivities(tdim - 1); // Create list of tagged boundary facets diff --git a/cpp/dolfinx/fem/utils.h b/cpp/dolfinx/fem/utils.h index 48cab5be495..6ee160191aa 100644 --- a/cpp/dolfinx/fem/utils.h +++ b/cpp/dolfinx/fem/utils.h @@ -146,13 +146,13 @@ get_cell_vertex_pairs(std::int32_t v, 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: for (int i = 0; i < a.num_integrals(type, cell_type_idx); ++i) { sparsitybuild::cells(pattern, @@ -502,8 +502,7 @@ Form create_form_factory( } // Create facets, if required - if (num_integrals_type[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); @@ -599,7 +598,7 @@ Form create_form_factory( std::span ids(ufcx_forms[0].get().form_integral_ids + integral_offsets[facet], num_integrals_type[facet]); - auto sd = subdomains.find(IntegralType::exterior_facet); + auto sd = subdomains.find(IntegralType::facet); for (std::size_t form_idx = 0; form_idx < ufcx_forms.size(); ++form_idx) { const ufcx_form& ufcx_form = ufcx_forms[form_idx]; @@ -638,7 +637,7 @@ Form create_form_factory( default_facets_ext.insert(default_facets_ext.end(), pair.begin(), pair.end()); } - integrals.insert({{IntegralType::exterior_facet, i, form_idx}, + integrals.insert({{IntegralType::facet, i, form_idx}, {k, default_facets_ext, active_coeffs}}); } else if (sd != subdomains.end()) @@ -648,7 +647,7 @@ Form create_form_factory( [](auto& a) { return a.first; }); if (it != sd->second.end() and it->first == id) { - integrals.insert({{IntegralType::exterior_facet, i, form_idx}, + integrals.insert({{IntegralType::facet, i, form_idx}, {k, std::vector(it->second.begin(), it->second.end()), diff --git a/python/dolfinx/fem/forms.py b/python/dolfinx/fem/forms.py index 3c0e246ed11..4e5da546039 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) @@ -217,7 +217,7 @@ def form_cpp_class( _ufl_to_dolfinx_domain = { "cell": IntegralType.cell, - "exterior_facet": IntegralType.exterior_facet, + "exterior_facet": IntegralType.facet, "interior_facet": IntegralType.interior_facet, "vertex": IntegralType.vertex, } @@ -369,17 +369,19 @@ def _form(form): constants = [c._cpp_object for c in form.constants()] # Extract subdomain ids from ufcx_form - subdomain_ids = {type: [] for type in sd.get(domain).keys()} + subdomain_ids = {_ufl_to_dolfinx_domain[type]: [] for type in sd.get(domain).keys()} integral_offsets = [ufcx_form.form_integral_offsets[i] for i in range(5)] for i in range(len(integral_offsets) - 1): integral_type = IntegralType(i) for j in range(integral_offsets[i], integral_offsets[i + 1]): - subdomain_ids[integral_type.name].append(ufcx_form.form_integral_ids[j]) + subdomain_ids[integral_type].append(ufcx_form.form_integral_ids[j]) # Subdomain markers (possibly empty list for some integral types) subdomains = { _ufl_to_dolfinx_domain[key]: get_integration_domains( - _ufl_to_dolfinx_domain[key], subdomain_data[0], subdomain_ids[key] + _ufl_to_dolfinx_domain[key], + subdomain_data[0], + subdomain_ids[_ufl_to_dolfinx_domain[key]], ) for (key, subdomain_data) in sd.get(domain).items() } diff --git a/python/dolfinx/wrappers/fem.cpp b/python/dolfinx/wrappers/fem.cpp index 1acc71e53a2..595c3b69341 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,8 +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, - "exterior facet integral") + .value("facet", dolfinx::fem::IntegralType::facet, "facet integral") .value("interior_facet", dolfinx::fem::IntegralType::interior_facet, "exterior facet integral") .value("vertex", dolfinx::fem::IntegralType::vertex, "vertex 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 bfff112d936..5d20745bd31 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_assembler.py b/python/test/unit/fem/test_assembler.py index 39fe8077fb0..72de8437822 100644 --- a/python/test/unit/fem/test_assembler.py +++ b/python/test/unit/fem/test_assembler.py @@ -1712,8 +1712,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} @@ -1852,8 +1852,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 75ddd17a323b776c27d2217fabaf279722825d2f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= Date: Wed, 10 Sep 2025 12:08:13 +0000 Subject: [PATCH 03/30] Add debug to understand why ffcx is not picked up on windows --- .github/workflows/windows.yml | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/.github/workflows/windows.yml b/.github/workflows/windows.yml index 5bc75252f4c..a9e23cb35d7 100644 --- a/.github/workflows/windows.yml +++ b/.github/workflows/windows.yml @@ -96,6 +96,11 @@ jobs: run: | python -m pip -v install --check-build-dependencies --no-build-isolation --no-cache-dir .[ci] --config-settings=cmake.args=-DBasix_DIR="D:/a/dolfinx/basix-install/lib/cmake/basix" --config-settings=cmake.build-type="Release" + - name: Debug env variables + run: | + cat dolfinx-src/.github/workflows/fenicsx-refs.env + cat ${{ env.ffcx_ref}} + - name: Checkout FFCx uses: actions/checkout@v5 with: From b9a6670fc6719d96a7f935654a338ca02eb26972 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= Date: Wed, 10 Sep 2025 12:21:04 +0000 Subject: [PATCH 04/30] use echo --- .github/workflows/windows.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/windows.yml b/.github/workflows/windows.yml index a9e23cb35d7..1784292a8d4 100644 --- a/.github/workflows/windows.yml +++ b/.github/workflows/windows.yml @@ -99,7 +99,7 @@ jobs: - name: Debug env variables run: | cat dolfinx-src/.github/workflows/fenicsx-refs.env - cat ${{ env.ffcx_ref}} + echo ${{ env.ffcx_ref}} - name: Checkout FFCx uses: actions/checkout@v5 From 9a7e92e14a3bdd40ab9e1e7b5fcd02b11202279b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= Date: Wed, 10 Sep 2025 12:27:50 +0000 Subject: [PATCH 05/30] Fix space --- .github/workflows/windows.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/windows.yml b/.github/workflows/windows.yml index 1784292a8d4..84d9d62ecf7 100644 --- a/.github/workflows/windows.yml +++ b/.github/workflows/windows.yml @@ -99,7 +99,7 @@ jobs: - name: Debug env variables run: | cat dolfinx-src/.github/workflows/fenicsx-refs.env - echo ${{ env.ffcx_ref}} + echo ${{ env.ffcx_ref }} - name: Checkout FFCx uses: actions/checkout@v5 From 1dfaef4875c1018e9b6306d8a44955db49541d8d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= Date: Wed, 10 Sep 2025 12:49:28 +0000 Subject: [PATCH 06/30] Fix windows CI --- .github/workflows/windows.yml | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/.github/workflows/windows.yml b/.github/workflows/windows.yml index 84d9d62ecf7..0214ded3d17 100644 --- a/.github/workflows/windows.yml +++ b/.github/workflows/windows.yml @@ -49,7 +49,7 @@ jobs: - name: Load environment variables run: | - cat dolfinx-src/.github/workflows/fenicsx-refs.env >> $GITHUB_ENV + cat dolfinx-src/.github/workflows/fenicsx-refs.env >> $env:GITHUB_ENV - name: Set up Python uses: actions/setup-python@v6 @@ -96,10 +96,6 @@ jobs: run: | python -m pip -v install --check-build-dependencies --no-build-isolation --no-cache-dir .[ci] --config-settings=cmake.args=-DBasix_DIR="D:/a/dolfinx/basix-install/lib/cmake/basix" --config-settings=cmake.build-type="Release" - - name: Debug env variables - run: | - cat dolfinx-src/.github/workflows/fenicsx-refs.env - echo ${{ env.ffcx_ref }} - name: Checkout FFCx uses: actions/checkout@v5 From c3e4f8110a687ee0f54ec250e0064b5137ca4c34 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= Date: Wed, 10 Sep 2025 15:40:00 +0000 Subject: [PATCH 07/30] Replace IntegralType with structure storing the necessary information --- cpp/demo/codim_0_assembly/main.cpp | 4 +- cpp/demo/custom_kernel/main.cpp | 4 +- cpp/demo/mixed_poisson/main.cpp | 6 +- cpp/dolfinx/fem/CMakeLists.txt | 1 + cpp/dolfinx/fem/Form.cpp | 24 ++++++++ cpp/dolfinx/fem/Form.h | 46 +++++++------- cpp/dolfinx/fem/assemble_matrix_impl.h | 43 +++++++------ cpp/dolfinx/fem/assemble_scalar_impl.h | 38 ++++++------ cpp/dolfinx/fem/assemble_vector_impl.h | 83 ++++++++++++++------------ cpp/dolfinx/fem/pack.h | 58 +++++------------- cpp/dolfinx/fem/utils.cpp | 43 ++++--------- cpp/dolfinx/fem/utils.h | 66 ++++++++++---------- python/dolfinx/fem/forms.py | 19 +++--- python/dolfinx/wrappers/fem.cpp | 29 +++++---- 14 files changed, 226 insertions(+), 238 deletions(-) create mode 100644 cpp/dolfinx/fem/Form.cpp diff --git a/cpp/demo/codim_0_assembly/main.cpp b/cpp/demo/codim_0_assembly/main.cpp index 5be7086e5e2..d227c01f24f 100644 --- a/cpp/demo/codim_0_assembly/main.cpp +++ b/cpp/demo/codim_0_assembly/main.cpp @@ -100,12 +100,12 @@ int main(int argc, char* argv[]) // domain `mesh` std::vector integration_entities = fem::compute_integration_domains( - fem::IntegralType::cell, *mesh->topology(), cell_marker.find(2)); + fem::IntegralType(0, 1), *mesh->topology(), cell_marker.find(2)); std::map< fem::IntegralType, std::vector>>> subdomain_data - = {{fem::IntegralType::cell, {{3, integration_entities}}}}; + = {{fem::IntegralType(0, 1), {{3, integration_entities}}}}; // A mixed-domain form involves functions defined over multiple // meshes. The mesh passed to `create_form` is called the diff --git a/cpp/demo/custom_kernel/main.cpp b/cpp/demo/custom_kernel/main.cpp index 6ddcfacabfb..8bdb1836e21 100644 --- a/cpp/demo/custom_kernel/main.cpp +++ b/cpp/demo/custom_kernel/main.cpp @@ -75,7 +75,7 @@ double assemble_matrix0(std::shared_ptr> V, { // Kernel data (ID, kernel function, cell indices to execute over) std::map integrals{ - std::pair{std::tuple{fem::IntegralType::cell, 0, 0}, + std::pair{std::tuple{fem::IntegralType(0, 1), 0, 0}, fem::integral_data(kernel, cells, std::vector{})}}; fem::Form a({V, V}, integrals, V->mesh(), {}, {}, false, {}); @@ -105,7 +105,7 @@ double assemble_vector0(std::shared_ptr> V, { auto mesh = V->mesh(); std::map integrals{ - std::pair{std::tuple{fem::IntegralType::cell, 0, 0}, + std::pair{std::tuple{fem::IntegralType(0, 1), 0, 0}, fem::integral_data(kernel, cells, std::vector{})}}; fem::Form L({V}, integrals, mesh, {}, {}, false, {}); auto dofmap = V->dofmap(); diff --git a/cpp/demo/mixed_poisson/main.cpp b/cpp/demo/mixed_poisson/main.cpp index ed8e638bc25..a6f330affe7 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); + IntegralType(1, 1), *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 @@ -279,9 +279,9 @@ int main(int argc, char* argv[]) // file, and `domains` holds the necessary data to perform // integration of selected facets. std::map< - fem::IntegralType, + std::int32_t, std::vector>>> - subdomain_data{{fem::IntegralType::facet, {{1, domains}}}}; + subdomain_data{{IntegralType(1, 1), {{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/CMakeLists.txt b/cpp/dolfinx/fem/CMakeLists.txt index 548df96a011..d72c37a5390 100644 --- a/cpp/dolfinx/fem/CMakeLists.txt +++ b/cpp/dolfinx/fem/CMakeLists.txt @@ -34,6 +34,7 @@ target_sources( ${CMAKE_CURRENT_SOURCE_DIR}/DofMap.cpp ${CMAKE_CURRENT_SOURCE_DIR}/ElementDofLayout.cpp ${CMAKE_CURRENT_SOURCE_DIR}/FiniteElement.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/Form.cpp ${CMAKE_CURRENT_SOURCE_DIR}/dofmapbuilder.cpp ${CMAKE_CURRENT_SOURCE_DIR}/petsc.cpp ${CMAKE_CURRENT_SOURCE_DIR}/sparsitybuild.cpp diff --git a/cpp/dolfinx/fem/Form.cpp b/cpp/dolfinx/fem/Form.cpp new file mode 100644 index 00000000000..c8aa0854b83 --- /dev/null +++ b/cpp/dolfinx/fem/Form.cpp @@ -0,0 +1,24 @@ +// Copyright (C) 2025 Jørgen S. Dokken +// +// This file is part of DOLFINx (https://www.fenicsproject.org) +// +// SPDX-License-Identifier: LGPL-3.0-or-later + +#include "Form.h" + +bool dolfinx::fem::operator<(const dolfinx::fem::IntegralType& self, + const dolfinx::fem::IntegralType& other) +{ + // First, compare codim + if (self.codim < other.codim) + { + return true; + } + // If codims are equal, compare num_cells + if (self.codim == other.codim) + { + return self.num_cells < other.num_cells; + } + // Otherwise, this object is not less than 'other' + return false; +} diff --git a/cpp/dolfinx/fem/Form.h b/cpp/dolfinx/fem/Form.h index f00fb2a8b7d..8416a231e7b 100644 --- a/cpp/dolfinx/fem/Form.h +++ b/cpp/dolfinx/fem/Form.h @@ -35,14 +35,24 @@ template class Function; /// @brief Type of integral -enum class IntegralType : std::int8_t + +struct IntegralType { - cell = 0, ///< Cell - facet = 1, ///< Facet - interior_facet = 2, ///< Interior facet - vertex = 3 ///< Vertex + std::int32_t codim; ///< Codimension of the integral + std::int32_t num_cells; ///< Number of cells in the integral + + /// @brief Equality operator for integral type + /// @param other The other IntegralType to compare with + /// @return True if the integral types are equal, false otherwise + bool operator==(const IntegralType& other) const + { + return codim == other.codim && num_cells == other.num_cells; + } }; +// Overload the less-than operator +bool operator<(const IntegralType& self, const IntegralType& other); + /// @brief Represents integral data, containing the kernel, and a list /// of entities to integrate over and the indicies of the coefficient /// functions (relative to the Form) active for this integral. @@ -274,10 +284,9 @@ class Form { auto [type, idx, kernel_idx] = key; std::vector e; - if (type == IntegralType::cell) + if (type.codim == 0) e = emap.sub_topology_to_topology(itg.entities, inverse); - else if (type == IntegralType::facet - or type == IntegralType::interior_facet) + else if (type.codim = 1) { const mesh::Topology topology = *_mesh->topology(); int tdim = topology.dim(); @@ -317,10 +326,9 @@ class Form bool inverse = emap.sub_topology() == mesh0->topology(); std::vector e; - if (type == IntegralType::cell) + if (type.codim == 0) e = emap.sub_topology_to_topology(integral.entities, inverse); - else if (type == IntegralType::facet - or type == IntegralType::interior_facet) + else if (type.codim == 1) { const mesh::Topology topology = *_mesh->topology(); int tdim = topology.dim(); @@ -454,16 +462,12 @@ class Form /// These are the entities in the mesh returned by ::mesh that are /// integrated over by a given integral (kernel). /// - /// - For IntegralType::cell, returns a list of cell indices. - /// - 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)`, - /// where `[cell_index, 0]` is one attached cell and `[cell_index, 1]` - /// is the is the local facet index relative to the cell, and - /// `[cell_index, 2]` is the other one attached cell and `[cell_index, 1]` - /// is the is the local facet index relative to this cell. Storage - /// is row-major. + /// - For IntegralType of codim 0 returns a list of cell indices. + /// - For any other codim return a set of tuples `[cell_index, + /// local_entity_index]`. + /// - If the number of cells in the integral type, return as many tuples per + /// entity as there are number of cells in the integral type. Storage is + /// row-major. /// /// @param[in] type Integral type. /// @param[in] idx Integral index in flattened list of integral kernels. diff --git a/cpp/dolfinx/fem/assemble_matrix_impl.h b/cpp/dolfinx/fem/assemble_matrix_impl.h index 66b5a58883c..9ee0c1cbf79 100644 --- a/cpp/dolfinx/fem/assemble_matrix_impl.h +++ b/cpp/dolfinx/fem/assemble_matrix_impl.h @@ -588,15 +588,15 @@ void assemble_matrix( cell_info0 = std::span(mesh0->topology()->get_cell_permutation_info()); cell_info1 = std::span(mesh1->topology()->get_cell_permutation_info()); } - - for (int i = 0; i < a.num_integrals(IntegralType::cell, cell_type_idx); ++i) + IntegralType cell_integral = IntegralType(0, 1); + for (int i = 0; i < a.num_integrals(cell_integral, cell_type_idx); ++i) { - auto fn = a.kernel(IntegralType::cell, i, cell_type_idx); + auto fn = a.kernel(cell_integral, i, cell_type_idx); assert(fn); - std::span cells = a.domain(IntegralType::cell, i, cell_type_idx); - std::span cells0 = a.domain_arg(IntegralType::cell, 0, i, cell_type_idx); - std::span cells1 = a.domain_arg(IntegralType::cell, 1, i, cell_type_idx); - auto& [coeffs, cstride] = coefficients.at({IntegralType::cell, i}); + std::span cells = a.domain(cell_integral, i, cell_type_idx); + std::span cells0 = a.domain_arg(cell_integral, 0, i, cell_type_idx); + std::span cells1 = a.domain_arg(cell_integral, 1, i, cell_type_idx); + auto& [coeffs, cstride] = coefficients.at({cell_integral, i}); assert(cells.size() * cstride == coeffs.size()); impl::assemble_cells_matrix( mat_set, x_dofmap, x, cells, {dofs0, bs0, cells0}, P0, @@ -618,8 +618,8 @@ void assemble_matrix( num_facets_per_cell); } - for (int i = 0; i < a.num_integrals(IntegralType::facet, cell_type_idx); - ++i) + IntegralType facet = IntegralType(1, 1); + for (int i = 0; i < a.num_integrals(facet, cell_type_idx); ++i) { if (num_cell_types > 1) { @@ -631,15 +631,15 @@ void assemble_matrix( = md::mdspan>; - auto fn = a.kernel(IntegralType::facet, i, 0); + auto fn = a.kernel(facet, i, 0); assert(fn); - auto& [coeffs, cstride] = coefficients.at({IntegralType::facet, i}); + auto& [coeffs, cstride] = coefficients.at({facet, i}); - std::span f = a.domain(IntegralType::facet, i, 0); + std::span f = a.domain(facet, i, 0); mdspanx2_t facets(f.data(), f.size() / 2, 2); - std::span f0 = a.domain_arg(IntegralType::facet, 0, i, 0); + std::span f0 = a.domain_arg(facet, 0, i, 0); mdspanx2_t facets0(f0.data(), f0.size() / 2, 2); - std::span f1 = a.domain_arg(IntegralType::facet, 1, i, 0); + std::span f1 = a.domain_arg(facet, 1, i, 0); mdspanx2_t facets1(f1.data(), f1.size() / 2, 2); assert((facets.size() / 2) * cstride == coeffs.size()); impl::assemble_exterior_facets( @@ -649,8 +649,8 @@ void assemble_matrix( cell_info0, cell_info1, perms); } - for (int i = 0; - i < a.num_integrals(IntegralType::interior_facet, cell_type_idx); ++i) + IntegralType interior_facet = IntegralType(1, 2); + for (int i = 0; i < a.num_integrals(interior_facet, cell_type_idx); ++i) { if (num_cell_types > 1) { @@ -665,14 +665,13 @@ void assemble_matrix( = md::mdspan>; - auto fn = a.kernel(IntegralType::interior_facet, i, 0); + auto fn = a.kernel(interior_facet, i, 0); assert(fn); - auto& [coeffs, cstride] - = coefficients.at({IntegralType::interior_facet, i}); + auto& [coeffs, cstride] = coefficients.at({interior_facet, i}); - std::span facets = a.domain(IntegralType::interior_facet, i, 0); - std::span facets0 = a.domain_arg(IntegralType::interior_facet, 0, i, 0); - std::span facets1 = a.domain_arg(IntegralType::interior_facet, 1, i, 0); + std::span facets = a.domain(interior_facet, i, 0); + std::span facets0 = a.domain_arg(interior_facet, 0, i, 0); + std::span facets1 = a.domain_arg(interior_facet, 1, i, 0); assert((facets.size() / 4) * 2 * cstride == coeffs.size()); impl::assemble_interior_facets( mat_set, x_dofmap, x, diff --git a/cpp/dolfinx/fem/assemble_scalar_impl.h b/cpp/dolfinx/fem/assemble_scalar_impl.h index 3ee1aa2ea94..0a6c0c1e982 100644 --- a/cpp/dolfinx/fem/assemble_scalar_impl.h +++ b/cpp/dolfinx/fem/assemble_scalar_impl.h @@ -201,12 +201,13 @@ T assemble_scalar( 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) + IntegralType cell_integral = IntegralType(0, 1); + for (int i = 0; i < M.num_integrals(cell_integral, 0); ++i) { - auto fn = M.kernel(IntegralType::cell, i, 0); + auto fn = M.kernel(cell_integral, i, 0); assert(fn); - auto& [coeffs, cstride] = coefficients.at({IntegralType::cell, i}); - std::span cells = M.domain(IntegralType::cell, i, 0); + auto& [coeffs, cstride] = coefficients.at({cell_integral, i}); + std::span cells = M.domain(cell_integral, i, 0); assert(cells.size() * cstride == coeffs.size()); value += impl::assemble_cells( x_dofmap, x, cells, fn, constants, @@ -226,13 +227,14 @@ T assemble_scalar( num_facets_per_cell); } - for (int i = 0; i < M.num_integrals(IntegralType::facet, 0); ++i) + IntegralType facet_type = IntegralType(1, 1); + for (int i = 0; i < M.num_integrals(facet_type, 0); ++i) { - auto fn = M.kernel(IntegralType::facet, i, 0); + auto fn = M.kernel(facet_type, i, 0); assert(fn); - auto& [coeffs, cstride] = coefficients.at({IntegralType::facet, i}); + auto& [coeffs, cstride] = coefficients.at({facet_type, i}); - std::span facets = M.domain(IntegralType::facet, i, 0); + std::span facets = M.domain(facet_type, i, 0); // Two values per each adjacent cell (cell index and local facet // index) @@ -250,13 +252,13 @@ T assemble_scalar( cdofs_b); } - for (int i = 0; i < M.num_integrals(IntegralType::interior_facet, 0); ++i) + IntegralType interior_facet_type = IntegralType(1, 2); + for (int i = 0; i < M.num_integrals(interior_facet_type, 0); ++i) { - auto fn = M.kernel(IntegralType::interior_facet, i, 0); + auto fn = M.kernel(interior_facet_type, i, 0); assert(fn); - auto& [coeffs, cstride] - = coefficients.at({IntegralType::interior_facet, i}); - std::span facets = M.domain(IntegralType::interior_facet, i, 0); + auto& [coeffs, cstride] = coefficients.at({interior_facet_type, i}); + std::span facets = M.domain(interior_facet_type, i, 0); constexpr std::size_t num_adjacent_cells = 2; // Two values per each adj. cell (cell index and local facet index). @@ -275,15 +277,15 @@ T assemble_scalar( perms, cdofs_b); } - for (int i = 0; i < M.num_integrals(IntegralType::vertex, 0); ++i) + IntegralType vertex_type = IntegralType(1, 2); + for (int i = 0; i < M.num_integrals(vertex_type, 0); ++i) { - auto fn = M.kernel(IntegralType::vertex, i, 0); + auto fn = M.kernel(vertex_type, i, 0); assert(fn); - auto& [coeffs, cstride] = coefficients.at({IntegralType::vertex, i}); + auto& [coeffs, cstride] = coefficients.at({vertex_type, i}); - std::span vertices - = M.domain(IntegralType::vertex, i, 0); + std::span vertices = M.domain(vertex_type, i, 0); assert(vertices.size() * cstride == coeffs.size()); constexpr std::size_t num_adjacent_cells = 1; diff --git a/cpp/dolfinx/fem/assemble_vector_impl.h b/cpp/dolfinx/fem/assemble_vector_impl.h index 87353e641ef..77a86d9ca8d 100644 --- a/cpp/dolfinx/fem/assemble_vector_impl.h +++ b/cpp/dolfinx/fem/assemble_vector_impl.h @@ -1100,14 +1100,15 @@ void lift_bc(V&& b, const Form& a, mdspan2_t x_dofmap, = element1->template dof_transformation_right_fn( doftransform::transpose); - for (int i = 0; i < a.num_integrals(IntegralType::cell, 0); ++i) + IntegralType cell_integral_type = IntegralType(0, 1); + for (int i = 0; i < a.num_integrals(cell_integral_type, 0); ++i) { - auto kernel = a.kernel(IntegralType::cell, i, 0); + auto kernel = a.kernel(cell_integral_type, i, 0); assert(kernel); - auto& [_coeffs, cstride] = coefficients.at({IntegralType::cell, i}); - std::span cells = a.domain(IntegralType::cell, i, 0); - std::span cells0 = a.domain_arg(IntegralType::cell, 0, i, 0); - std::span cells1 = a.domain_arg(IntegralType::cell, 1, i, 0); + auto& [_coeffs, cstride] = coefficients.at({cell_integral_type, i}); + std::span cells = a.domain(cell_integral_type, i, 0); + std::span cells0 = a.domain_arg(cell_integral_type, 0, i, 0); + std::span cells1 = a.domain_arg(cell_integral_type, 1, i, 0); assert(_coeffs.size() == cells.size() * cstride); auto coeffs = md::mdspan(_coeffs.data(), cells.size(), cstride); if (bs0 == 1 and bs1 == 1) @@ -1145,20 +1146,21 @@ void lift_bc(V&& b, const Form& a, mdspan2_t x_dofmap, num_facets_per_cell); } - for (int i = 0; i < a.num_integrals(IntegralType::facet, 0); ++i) + IntegralType facet_integral_type = IntegralType(1, 1); + for (int i = 0; i < a.num_integrals(facet_integral_type, 0); ++i) { - auto kernel = a.kernel(IntegralType::facet, i, 0); + auto kernel = a.kernel(facet_integral_type, i, 0); assert(kernel); - auto& [coeffs, cstride] = coefficients.at({IntegralType::facet, i}); + auto& [coeffs, cstride] = coefficients.at({facet_integral_type, i}); using mdspanx2_t = md::mdspan>; - std::span f = a.domain(IntegralType::facet, i, 0); + std::span f = a.domain(facet_integral_type, i, 0); mdspanx2_t facets(f.data(), f.size() / 2, 2); - std::span f0 = a.domain_arg(IntegralType::facet, 0, i, 0); + std::span f0 = a.domain_arg(facet_integral_type, 0, i, 0); mdspanx2_t facets0(f0.data(), f0.size() / 2, 2); - std::span f1 = a.domain_arg(IntegralType::facet, 1, i, 0); + std::span f1 = a.domain_arg(facet_integral_type, 1, i, 0); mdspanx2_t facets1(f1.data(), f1.size() / 2, 2); assert(coeffs.size() == facets.extent(0) * cstride); _lift_bc_exterior_facets( @@ -1168,12 +1170,13 @@ void lift_bc(V&& b, const Form& a, mdspan2_t x_dofmap, cell_info1, bc_values1, bc_markers1, x0, alpha, perms); } - for (int i = 0; i < a.num_integrals(IntegralType::interior_facet, 0); ++i) + IntegralType interior_facet_integral_type = IntegralType(1, 2); + for (int i = 0; i < a.num_integrals(interior_facet_integral_type, 0); ++i) { - auto kernel = a.kernel(IntegralType::interior_facet, i, 0); + auto kernel = a.kernel(interior_facet_integral_type, i, 0); assert(kernel); auto& [coeffs, cstride] - = coefficients.at({IntegralType::interior_facet, i}); + = coefficients.at({interior_facet_integral_type, i}); using mdspanx22_t = md::mdspan& a, mdspan2_t x_dofmap, using mdspanx2x_t = md::mdspan>; - std::span f = a.domain(IntegralType::interior_facet, i, 0); + std::span f = a.domain(interior_facet_integral_type, i, 0); mdspanx22_t facets(f.data(), f.size() / 4, 2, 2); - std::span f0 = a.domain_arg(IntegralType::interior_facet, 0, i, 0); + std::span f0 = a.domain_arg(interior_facet_integral_type, 0, i, 0); mdspanx22_t facets0(f0.data(), f0.size() / 4, 2, 2); - std::span f1 = a.domain_arg(IntegralType::interior_facet, 1, i, 0); + std::span f1 = a.domain_arg(interior_facet_integral_type, 1, i, 0); mdspanx22_t facets1(f1.data(), f1.size() / 4, 2, 2); _lift_bc_interior_facets( b, x_dofmap, x, kernel, facets, {dofmap0, bs0, facets0}, P0, @@ -1339,13 +1342,14 @@ void assemble_vector( cell_info0 = std::span(mesh0->topology()->get_cell_permutation_info()); } - for (int i = 0; i < L.num_integrals(IntegralType::cell, 0); ++i) + IntegralType cell_integral_type = IntegralType(0, 1); + for (int i = 0; i < L.num_integrals(cell_integral_type, 0); ++i) { - auto fn = L.kernel(IntegralType::cell, i, cell_type_idx); + auto fn = L.kernel(cell_integral_type, i, cell_type_idx); assert(fn); - std::span cells = L.domain(IntegralType::cell, i, cell_type_idx); - std::span cells0 = L.domain_arg(IntegralType::cell, 0, i, cell_type_idx); - auto& [coeffs, cstride] = coefficients.at({IntegralType::cell, i}); + std::span cells = L.domain(cell_integral_type, i, cell_type_idx); + std::span cells0 = L.domain_arg(cell_integral_type, 0, i, cell_type_idx); + auto& [coeffs, cstride] = coefficients.at({cell_integral_type, i}); assert(cells.size() * cstride == coeffs.size()); if (bs == 1) { @@ -1384,14 +1388,15 @@ void assemble_vector( = md::mdspan>; - for (int i = 0; i < L.num_integrals(IntegralType::facet, 0); ++i) + IntegralType facet_integral_type = IntegralType(1, 1); + for (int i = 0; i < L.num_integrals(facet_integral_type, 0); ++i) { - auto fn = L.kernel(IntegralType::facet, i, 0); + auto fn = L.kernel(facet_integral_type, i, 0); assert(fn); - auto& [coeffs, cstride] = coefficients.at({IntegralType::facet, i}); - std::span f = L.domain(IntegralType::facet, i, 0); + auto& [coeffs, cstride] = coefficients.at({facet_integral_type, i}); + std::span f = L.domain(facet_integral_type, i, 0); mdspanx2_t facets(f.data(), f.size() / 2, 2); - std::span f1 = L.domain_arg(IntegralType::facet, 0, i, 0); + std::span f1 = L.domain_arg(facet_integral_type, 0, i, 0); mdspanx2_t facets1(f1.data(), f1.size() / 2, 2); assert((facets.size() / 2) * cstride == coeffs.size()); if (bs == 1) @@ -1417,7 +1422,8 @@ void assemble_vector( } } - for (int i = 0; i < L.num_integrals(IntegralType::interior_facet, 0); ++i) + IntegralType interior_facet_integral_type = IntegralType(1, 2); + for (int i = 0; i < L.num_integrals(interior_facet_integral_type, 0); ++i) { using mdspanx22_t = md::mdspan>; - auto fn = L.kernel(IntegralType::interior_facet, i, 0); + auto fn = L.kernel(interior_facet_integral_type, i, 0); assert(fn); auto& [coeffs, cstride] - = coefficients.at({IntegralType::interior_facet, i}); - std::span facets = L.domain(IntegralType::interior_facet, i, 0); - std::span facets1 = L.domain_arg(IntegralType::interior_facet, 0, i, 0); + = coefficients.at({interior_facet_integral_type, i}); + std::span facets = L.domain(interior_facet_integral_type, i, 0); + std::span facets1 = L.domain_arg(interior_facet_integral_type, 0, i, 0); assert((facets.size() / 4) * 2 * cstride == coeffs.size()); if (bs == 1) { @@ -1468,16 +1474,17 @@ void assemble_vector( } } - for (int i = 0; i < L.num_integrals(IntegralType::vertex, 0); ++i) + IntegralType vertex_integral_type = IntegralType(2, 1); + for (int i = 0; i < L.num_integrals(vertex_integral_type, 0); ++i) { - auto fn = L.kernel(IntegralType::vertex, i, 0); + auto fn = L.kernel(vertex_integral_type, i, 0); assert(fn); - std::span vertices = L.domain(IntegralType::vertex, i, cell_type_idx); + std::span vertices = L.domain(vertex_integral_type, i, cell_type_idx); std::span vertices0 - = L.domain_arg(IntegralType::vertex, 0, i, cell_type_idx); + = L.domain_arg(vertex_integral_type, 0, i, cell_type_idx); - auto& [coeffs, cstride] = coefficients.at({IntegralType::vertex, i}); + auto& [coeffs, cstride] = coefficients.at({vertex_integral_type, i}); assert(vertices.size() * cstride == coeffs.size()); diff --git a/cpp/dolfinx/fem/pack.h b/cpp/dolfinx/fem/pack.h index a943d872c90..07b747cf373 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::facet - or integral_type == IntegralType::interior_facet) - { + if (integral_type.codim != 0) num_entities /= 2; - } } return {std::vector(num_entities * cstride), cstride}; @@ -237,12 +234,10 @@ void pack_coefficients(const Form& form, int cstride = coeff_data.second; if (!coefficients.empty()) { - switch (integral_type) - { - case IntegralType::cell: + if (integral_type.codim == 0) { // Iterate over coefficients that are active in cell integrals - for (int coeff : form.active_coeffs(IntegralType::cell, id)) + for (int coeff : form.active_coeffs(IntegralType(0, 1), id)) { // Get coefficient mesh auto mesh = coefficients[coeff]->function_space()->mesh(); @@ -260,7 +255,7 @@ void pack_coefficients(const Form& form, } std::span cells_b - = form.domain_coeff(IntegralType::cell, id, coeff); + = form.domain_coeff(IntegralType(0, 1), id, coeff); md::mdspan cells(cells_b.data(), cells_b.size()); std::span cell_info = impl::get_cell_orientation_info(*coefficients[coeff]); @@ -268,21 +263,21 @@ void pack_coefficients(const Form& form, *coefficients[coeff], cell_info, cells, offsets[coeff]); } - break; } - case IntegralType::facet: + else if (integral_type.num_cells == 1) { // Iterate over coefficients coefficients that are active in // exterior facet integrals - for (int coeff : form.active_coeffs(IntegralType::facet, id)) + for (int coeff : + form.active_coeffs(IntegralType(integral_type.codim, 1), id)) { auto mesh = coefficients[coeff]->function_space()->mesh(); - std::span facets_b - = form.domain_coeff(IntegralType::facet, id, coeff); + std::span entities_b = form.domain_coeff( + IntegralType(integral_type.codim, 1), id, coeff); md::mdspan> - facets(facets_b.data(), facets_b.size() / 2, 2); - auto cells = md::submdspan(facets, md::full_extent, 0); + entities(entities_b.data(), entities_b.size() / 2, 2); + auto cells = md::submdspan(entities, md::full_extent, 0); std::span cell_info = impl::get_cell_orientation_info(*coefficients[coeff]); @@ -290,17 +285,16 @@ void pack_coefficients(const Form& form, *coefficients[coeff], cell_info, cells, offsets[coeff]); } - break; } - case IntegralType::interior_facet: + else if ((integral_type.codim == 1) and (integral_type.num_cells == 2)) { // Iterate over coefficients that are active in interior // facet integrals - for (int coeff : form.active_coeffs(IntegralType::interior_facet, id)) + for (int coeff : form.active_coeffs(IntegralType(1, 2), id)) { auto mesh = coefficients[coeff]->function_space()->mesh(); std::span facets_b - = form.domain_coeff(IntegralType::interior_facet, id, coeff); + = form.domain_coeff(IntegralType(1, 2), id, coeff); md::mdspan> facets(facets_b.data(), facets_b.size() / 4, 4); @@ -320,31 +314,9 @@ void pack_coefficients(const Form& form, *coefficients[coeff], cell_info, cells1, offsets[coeff] + offsets[coeff + 1]); } - break; } - case IntegralType::vertex: + else { - // Iterate over coefficients that are active in vertex integrals - for (int coeff : form.active_coeffs(IntegralType::vertex, id)) - { - // Get coefficient mesh - auto mesh = coefficients[coeff]->function_space()->mesh(); - assert(mesh); - - std::span vertices_b - = form.domain_coeff(IntegralType::vertex, id, coeff); - md::mdspan> - vertices(vertices_b.data(), vertices_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]); - } - break; - } - default: throw std::runtime_error( "Could not pack coefficient. Integral type not supported."); } diff --git a/cpp/dolfinx/fem/utils.cpp b/cpp/dolfinx/fem/utils.cpp index 4b36d5dd18e..7ef7a8f9bb8 100644 --- a/cpp/dolfinx/fem/utils.cpp +++ b/cpp/dolfinx/fem/utils.cpp @@ -133,31 +133,14 @@ std::vector fem::get_constant_names(const ufcx_form& ufcx_form) } //----------------------------------------------------------------------------- std::vector -fem::compute_integration_domains(fem::IntegralType integral_type, +fem::compute_integration_domains(const fem::IntegralType& integral_type, const mesh::Topology& topology, std::span entities) { const int tdim = topology.dim(); - int dim = -1; - switch (integral_type) - { - case IntegralType::cell: - dim = tdim; - break; - case IntegralType::facet: - dim = tdim - 1; - break; - case IntegralType::interior_facet: - dim = tdim - 1; - break; - case IntegralType::vertex: - dim = 0; - break; - default: - throw std::runtime_error( - "Cannot compute integration domains. Integral type not supported."); - } + int dim = tdim - integral_type.codim; + assert(dim >= 0); { // Create span of the owned entities (leaves off any ghosts) @@ -192,14 +175,11 @@ fem::compute_integration_domains(fem::IntegralType integral_type, }; std::vector entity_data; - switch (integral_type) - { - case IntegralType::cell: + if (integral_type.codim == 0) { entity_data.insert(entity_data.begin(), entities.begin(), entities.end()); - break; } - case IntegralType::facet: + else if (integral_type.codim == 1 and integral_type.num_cells == 1) { auto [f_to_c, c_to_f] = get_connectivities(tdim - 1); // Create list of tagged boundary facets @@ -213,9 +193,8 @@ fem::compute_integration_domains(fem::IntegralType integral_type, 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: + else if (integral_type.codim == 1 and integral_type.num_cells == 2) { auto [f_to_c, c_to_f] = get_connectivities(tdim - 1); @@ -246,9 +225,8 @@ fem::compute_integration_domains(fem::IntegralType integral_type, "Use \"shared facet\" ghost mode when creating the mesh."); } } - break; } - case IntegralType::vertex: + else if (integral_type.codim == tdim and integral_type.num_cells == 1) { auto [v_to_c, c_to_v] = get_connectivities(0); for (auto vertex : entities) @@ -259,8 +237,13 @@ fem::compute_integration_domains(fem::IntegralType integral_type, entity_data.insert(entity_data.end(), pair.begin(), pair.end()); } } + else + { + throw std::runtime_error("Unsupported integral type with codim " + + std::to_string(integral_type.codim) + + " and num_cells " + + std::to_string(integral_type.num_cells)); } - return entity_data; } //----------------------------------------------------------------------------- diff --git a/cpp/dolfinx/fem/utils.h b/cpp/dolfinx/fem/utils.h index 6ee160191aa..bedbe0208bb 100644 --- a/cpp/dolfinx/fem/utils.h +++ b/cpp/dolfinx/fem/utils.h @@ -144,19 +144,11 @@ get_cell_vertex_pairs(std::int32_t v, std::span cells, /// @param[in] integral_type Integral type. /// @param[in] topology Mesh topology. /// @param[in] entities List of mesh entities. Depending on the `IntegralType` -/// these are associated with different entities: -/// `IntegralType::cell`: cells -/// `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::interior_facet`: (cell, local_facet) -/// `IntegralType::vertex`: (cell, local_vertex) +/// @return List of integration entity data. For `integral_type.codim` is 0 this +/// is just a list of cells. For `codim>0` this is a list of tuples (cell, +/// local_entity_index) std::vector -compute_integration_domains(IntegralType integral_type, +compute_integration_domains(const IntegralType& integral_type, const mesh::Topology& topology, std::span entities); @@ -237,9 +229,11 @@ void build_sparsity_pattern(la::SparsityPattern& pattern, const Form& a) std::shared_ptr mesh1 = a.function_spaces().at(1)->mesh(); assert(mesh1); - const std::set types = a.integral_types(); - if (types.find(IntegralType::interior_facet) != types.end() - or types.find(IntegralType::facet) != types.end()) + const std::set> types + = a.integral_types(); + + if (types.find(IntegralType(1, 1)) != types.end() + or types.find(IntegralType(1, 2)) != types.end()) { // FIXME: cleanup these calls? Some of the happen internally again. int tdim = mesh->topology()->dim(); @@ -269,9 +263,8 @@ void build_sparsity_pattern(la::SparsityPattern& pattern, const Form& a) // Create and build sparsity pattern for (auto type : types) { - switch (type) + if (type == IntegralType(0, 1)) { - case IntegralType::cell: for (int i = 0; i < a.num_integrals(type, cell_type_idx); ++i) { sparsitybuild::cells(pattern, @@ -279,8 +272,9 @@ void build_sparsity_pattern(la::SparsityPattern& pattern, const Form& a) a.domain_arg(type, 1, i, cell_type_idx)}, {{dofmaps[0], dofmaps[1]}}); } - break; - case IntegralType::interior_facet: + } + else if (type == IntegralType(1, 2)) + { for (int i = 0; i < a.num_integrals(type, cell_type_idx); ++i) { sparsitybuild::interior_facets( @@ -289,8 +283,9 @@ void build_sparsity_pattern(la::SparsityPattern& pattern, const Form& a) extract_cells(a.domain_arg(type, 1, i, 0))}, {{dofmaps[0], dofmaps[1]}}); } - break; - case IntegralType::facet: + } + else if (type == IntegralType(1, 1)) + { for (int i = 0; i < a.num_integrals(type, cell_type_idx); ++i) { sparsitybuild::cells(pattern, @@ -298,10 +293,9 @@ void build_sparsity_pattern(la::SparsityPattern& pattern, const Form& a) extract_cells(a.domain_arg(type, 1, i, 0))}, {{dofmaps[0], dofmaps[1]}}); } - break; - default: - throw std::runtime_error("Unsupported integral type"); } + else + throw std::runtime_error("Unsupported integral type"); } } @@ -533,7 +527,7 @@ Form create_form_factory( std::span ids(ufcx_forms[0].get().form_integral_ids + integral_offsets[cell], num_integrals_type[cell]); - auto sd = subdomains.find(IntegralType::cell); + auto sd = subdomains.find(IntegralType(0, 1)); for (std::size_t form_idx = 0; form_idx < ufcx_forms.size(); ++form_idx) { const ufcx_form& ufcx_form = ufcx_forms[form_idx]; @@ -568,7 +562,7 @@ Form create_form_factory( default_cells.resize( topology->index_maps(tdim).at(form_idx)->size_local(), 0); std::iota(default_cells.begin(), default_cells.end(), 0); - integrals.insert({{IntegralType::cell, i, form_idx}, + integrals.insert({{IntegralType(0, 1), i, form_idx}, {k, default_cells, active_coeffs}}); } else if (sd != subdomains.end()) @@ -578,7 +572,7 @@ Form create_form_factory( [](auto& a) { return a.first; }); if (it != sd->second.end() and it->first == id) { - integrals.insert({{IntegralType::cell, i, form_idx}, + integrals.insert({{IntegralType(0, 1), i, form_idx}, {k, std::vector(it->second.begin(), it->second.end()), @@ -598,7 +592,7 @@ Form create_form_factory( std::span ids(ufcx_forms[0].get().form_integral_ids + integral_offsets[facet], num_integrals_type[facet]); - auto sd = subdomains.find(IntegralType::facet); + auto sd = subdomains.find(IntegralType(1, 1)); for (std::size_t form_idx = 0; form_idx < ufcx_forms.size(); ++form_idx) { const ufcx_form& ufcx_form = ufcx_forms[form_idx]; @@ -637,7 +631,7 @@ Form create_form_factory( default_facets_ext.insert(default_facets_ext.end(), pair.begin(), pair.end()); } - integrals.insert({{IntegralType::facet, i, form_idx}, + integrals.insert({{IntegralType(1, 1), i, form_idx}, {k, default_facets_ext, active_coeffs}}); } else if (sd != subdomains.end()) @@ -647,7 +641,7 @@ Form create_form_factory( [](auto& a) { return a.first; }); if (it != sd->second.end() and it->first == id) { - integrals.insert({{IntegralType::facet, i, form_idx}, + integrals.insert({{IntegralType(1, 1), i, form_idx}, {k, std::vector(it->second.begin(), it->second.end()), @@ -667,7 +661,7 @@ Form create_form_factory( std::span ids(ufcx_forms[0].get().form_integral_ids + integral_offsets[interior_facet], num_integrals_type[interior_facet]); - auto sd = subdomains.find(IntegralType::interior_facet); + auto sd = subdomains.find(IntegralType(1, 1)); for (std::size_t form_idx = 0; form_idx < ufcx_forms.size(); ++form_idx) { const ufcx_form& ufcx_form = ufcx_forms[form_idx]; @@ -733,7 +727,7 @@ Form create_form_factory( "mesh"); } } - integrals.insert({{IntegralType::interior_facet, i, form_idx}, + integrals.insert({{IntegralType(1, 2), i, form_idx}, {k, default_facets_int, active_coeffs}}); } else if (sd != subdomains.end()) @@ -742,7 +736,7 @@ Form create_form_factory( [](auto& a) { return a.first; }); if (it != sd->second.end() and it->first == id) { - integrals.insert({{IntegralType::interior_facet, i, form_idx}, + integrals.insert({{IntegralType(1, 2), i, form_idx}, {k, std::vector(it->second.begin(), it->second.end()), @@ -765,7 +759,7 @@ Form create_form_factory( std::span ids(form.form_integral_ids + integral_offsets[vertex], num_integrals_type[vertex]); - auto sd = subdomains.find(IntegralType::vertex); + auto sd = subdomains.find(IntegralType(tdim, 1)); for (int i = 0; i < num_integrals_type[vertex]; ++i) { const int id = ids[i]; @@ -816,7 +810,7 @@ Form create_form_factory( std::vector cells_and_vertices = get_cells_and_vertices( std::ranges::views::iota(0, num_vertices)); - integrals.insert({{IntegralType::vertex, i, form_idx}, + integrals.insert({{IntegralType(tdim, 1), i, form_idx}, {k, cells_and_vertices, active_coeffs}}); } else if (sd != subdomains.end()) @@ -826,7 +820,7 @@ Form create_form_factory( [](auto& a) { return a.first; }); if (it != sd->second.end() and it->first == id) { - integrals.insert({{IntegralType::vertex, i, form_idx}, + integrals.insert({{IntegralType(tdim, 1), i, form_idx}, {k, std::vector(it->second.begin(), it->second.end()), diff --git a/python/dolfinx/fem/forms.py b/python/dolfinx/fem/forms.py index 4e5da546039..9b33e4e59d5 100644 --- a/python/dolfinx/fem/forms.py +++ b/python/dolfinx/fem/forms.py @@ -157,13 +157,11 @@ def get_integration_domains( else: domains = [] if not isinstance(subdomain, list): - if integral_type in (IntegralType.facet, IntegralType.interior_facet): - tdim = subdomain.topology.dim + tdim = subdomain.topology.dim + if integral_type.codim == 1: subdomain._cpp_object.topology.create_connectivity(tdim - 1, tdim) subdomain._cpp_object.topology.create_connectivity(tdim, tdim - 1) - - if integral_type is IntegralType.vertex: - tdim = subdomain.topology.dim + if integral_type.codim == tdim: subdomain._cpp_object.topology.create_connectivity(0, tdim) subdomain._cpp_object.topology.create_connectivity(tdim, 0) @@ -216,10 +214,10 @@ def form_cpp_class( _ufl_to_dolfinx_domain = { - "cell": IntegralType.cell, - "exterior_facet": IntegralType.facet, - "interior_facet": IntegralType.interior_facet, - "vertex": IntegralType.vertex, + "cell": IntegralType(0, 1), + "exterior_facet": IntegralType(1, 1), + "interior_facet": IntegralType(1, 2), + "vertex": IntegralType(2, 1), } @@ -371,8 +369,9 @@ def _form(form): # Extract subdomain ids from ufcx_form subdomain_ids = {_ufl_to_dolfinx_domain[type]: [] for type in sd.get(domain).keys()} integral_offsets = [ufcx_form.form_integral_offsets[i] for i in range(5)] + integral_types = [_ufl_to_dolfinx_domain[type] for type in ["cell", "exterior_facet", "interior_facet", "vertex"]] for i in range(len(integral_offsets) - 1): - integral_type = IntegralType(i) + integral_type = integral_types[i] for j in range(integral_offsets[i], integral_offsets[i + 1]): subdomain_ids[integral_type].append(ufcx_form.form_integral_ids[j]) diff --git a/python/dolfinx/wrappers/fem.cpp b/python/dolfinx/wrappers/fem.cpp index 595c3b69341..968094e7985 100644 --- a/python/dolfinx/wrappers/fem.cpp +++ b/python/dolfinx/wrappers/fem.cpp @@ -34,6 +34,7 @@ #include #include #include +#include #include #include #include @@ -767,22 +768,24 @@ void declare_form(nb::module_& m, std::string type) dolfinx::fem::IntegralType type, int i) { std::span _d = self.domain(type, i, 0); - switch (type) + if (type.codim == 0) { - case dolfinx::fem::IntegralType::cell: return nb::ndarray(_d.data(), {_d.size()}); - case dolfinx::fem::IntegralType::facet: + } + else if ((type.codim == 1) and (type.num_cells == 2)) { return nb::ndarray( - _d.data(), {_d.size() / 2, 2}); + _d.data(), {_d.size() / 4, 2, 2}); } - case dolfinx::fem::IntegralType::interior_facet: + else if (type.num_cells == 1) { return nb::ndarray( - _d.data(), {_d.size() / 4, 2, 2}); + _d.data(), {_d.size() / 2, 2}); } - default: + + else + { throw ::std::runtime_error("Integral type unsupported."); } }, @@ -1267,12 +1270,12 @@ void fem(nb::module_& m) }, nb::rv_policy::reference_internal); - nb::enum_(m, "_IntegralType") - .value("cell", dolfinx::fem::IntegralType::cell, "cell integral") - .value("facet", dolfinx::fem::IntegralType::facet, "facet integral") - .value("interior_facet", dolfinx::fem::IntegralType::interior_facet, - "exterior facet integral") - .value("vertex", dolfinx::fem::IntegralType::vertex, "vertex integral"); + nb::class_(m, "_IntegralType") + .def(nb::init(), + nb::arg("codim"), nb::arg("num_cells")) + .def_ro("codim", &dolfinx::fem::IntegralType::codim) + .def_ro("num_cells", &dolfinx::fem::IntegralType::num_cells) + .def("__eq__", &dolfinx::fem::IntegralType::operator==); declare_function_space(m, "float32"); declare_function_space(m, "float64"); From 192469087f715991363590a77b7d3f9d34cc4c20 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= Date: Wed, 10 Sep 2025 15:42:49 +0000 Subject: [PATCH 08/30] Ruff formatting --- python/dolfinx/fem/forms.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/python/dolfinx/fem/forms.py b/python/dolfinx/fem/forms.py index 9b33e4e59d5..1cde6830b46 100644 --- a/python/dolfinx/fem/forms.py +++ b/python/dolfinx/fem/forms.py @@ -369,7 +369,10 @@ def _form(form): # Extract subdomain ids from ufcx_form subdomain_ids = {_ufl_to_dolfinx_domain[type]: [] for type in sd.get(domain).keys()} integral_offsets = [ufcx_form.form_integral_offsets[i] for i in range(5)] - integral_types = [_ufl_to_dolfinx_domain[type] for type in ["cell", "exterior_facet", "interior_facet", "vertex"]] + integral_types = [ + _ufl_to_dolfinx_domain[type] + for type in ["cell", "exterior_facet", "interior_facet", "vertex"] + ] for i in range(len(integral_offsets) - 1): integral_type = integral_types[i] for j in range(integral_offsets[i], integral_offsets[i + 1]): From b0ffb36fa28a451ac60e9b63eb235763f4706dd8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= Date: Wed, 10 Sep 2025 15:50:00 +0000 Subject: [PATCH 09/30] Replace remaining integrlal types --- cpp/dolfinx/fem/assemble_scalar_impl.h | 3 ++- python/demo/demo_static-condensation.py | 2 +- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/cpp/dolfinx/fem/assemble_scalar_impl.h b/cpp/dolfinx/fem/assemble_scalar_impl.h index 0a6c0c1e982..a7146eb9ac0 100644 --- a/cpp/dolfinx/fem/assemble_scalar_impl.h +++ b/cpp/dolfinx/fem/assemble_scalar_impl.h @@ -277,7 +277,8 @@ T assemble_scalar( perms, cdofs_b); } - IntegralType vertex_type = IntegralType(1, 2); + std::int32_t tdim = mesh->topology()->dim(); + IntegralType vertex_type = IntegralType(tdim, 1); for (int i = 0; i < M.num_integrals(vertex_type, 0); ++i) { auto fn = M.kernel(vertex_type, i, 0); diff --git a/python/demo/demo_static-condensation.py b/python/demo/demo_static-condensation.py index 5be91e0eec6..815d245268b 100644 --- a/python/demo/demo_static-condensation.py +++ b/python/demo/demo_static-condensation.py @@ -183,7 +183,7 @@ def tabulate_A(A_, w_, c_, coords_, entity_local_index, permutation=ffi.NULL, cu # Prepare a Form with a condensed tabulation kernel formtype = form_cpp_class(PETSc.ScalarType) # type: ignore cells = np.arange(msh.topology.index_map(msh.topology.dim).size_local) -integrals = {IntegralType.cell: [(0, tabulate_A.address, cells, np.array([], dtype=np.int8))]} +integrals = {IntegralType(0, 1): [(0, tabulate_A.address, cells, np.array([], dtype=np.int8))]} a_cond = Form( formtype([U._cpp_object, U._cpp_object], integrals, [], [], False, [], mesh=msh._cpp_object) ) From 7ac5ebfccee8bec9ddfcc3715dde420c4a41d8ef Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= Date: Wed, 10 Sep 2025 15:59:34 +0000 Subject: [PATCH 10/30] Fix last bug --- cpp/dolfinx/fem/utils.h | 2 +- .../unit/fem/test_assemble_mesh_independent_form.py | 4 ++-- python/test/unit/fem/test_assemble_submesh.py | 2 +- python/test/unit/fem/test_assembler.py | 8 ++++---- python/test/unit/fem/test_custom_jit_kernels.py | 10 +++++----- python/test/unit/fem/test_forms.py | 4 ++-- 6 files changed, 15 insertions(+), 15 deletions(-) diff --git a/cpp/dolfinx/fem/utils.h b/cpp/dolfinx/fem/utils.h index bedbe0208bb..facf3ae9aa5 100644 --- a/cpp/dolfinx/fem/utils.h +++ b/cpp/dolfinx/fem/utils.h @@ -661,7 +661,7 @@ Form create_form_factory( std::span ids(ufcx_forms[0].get().form_integral_ids + integral_offsets[interior_facet], num_integrals_type[interior_facet]); - auto sd = subdomains.find(IntegralType(1, 1)); + auto sd = subdomains.find(IntegralType(1, 2)); for (std::size_t form_idx = 0; form_idx < ufcx_forms.size(); ++form_idx) { const ufcx_form& ufcx_form = ufcx_forms[form_idx]; 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 5d20745bd31..e5baf3f2a2e 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(1, 1), mesh.topology, facets ) - subdomains = {dolfinx.fem.IntegralType.facet: [(subdomain_id, facet_entities)]} + subdomains = {dolfinx.fem.IntegralType(1, 1): [(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 0334eb70878..38d6ca94f5b 100644 --- a/python/test/unit/fem/test_assemble_submesh.py +++ b/python/test/unit/fem/test_assemble_submesh.py @@ -442,7 +442,7 @@ def compute_mapped_interior_facet_data(mesh, facet_tag, value, parent_to_sub_map mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim) assert facet_tag.dim == mesh.topology.dim - 1 integration_data = compute_integration_domains( - fem.IntegralType.interior_facet, mesh.topology, facet_tag.find(value) + fem.IntegralType(1, 2), mesh.topology, facet_tag.find(value) ) mapped_cell_0 = parent_to_sub_map[integration_data[0::4]] mapped_cell_1 = parent_to_sub_map[integration_data[2::4]] diff --git a/python/test/unit/fem/test_assembler.py b/python/test/unit/fem/test_assembler.py index 72de8437822..49d27310f81 100644 --- a/python/test/unit/fem/test_assembler.py +++ b/python/test/unit/fem/test_assembler.py @@ -1712,8 +1712,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(1, 1), msh.topology, vertices) + subdomains = {fem.IntegralType(1, 1): [(0, cell_vertex_pairs)]} compiled_form = fem.compile_form( comm, x[0] * ufl.dP, form_compiler_options={"scalar_type": dtype} @@ -1852,8 +1852,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(1, 1), msh.topology, vertices) + subdomains = {fem.IntegralType(1, 1): [(0, cell_vertex_pairs)]} compiled_form = fem.compile_form( comm, x[0] * v * ufl.dP, form_compiler_options={"scalar_type": dtype} diff --git a/python/test/unit/fem/test_custom_jit_kernels.py b/python/test/unit/fem/test_custom_jit_kernels.py index 0f0a0acdd2c..f70cdec3b00 100644 --- a/python/test/unit/fem/test_custom_jit_kernels.py +++ b/python/test/unit/fem/test_custom_jit_kernels.py @@ -93,7 +93,7 @@ def test_numba_assembly(dtype): cells = np.arange(mesh.topology.index_map(mesh.topology.dim).size_local, dtype=np.int32) active_coeffs = np.array([], dtype=np.int8) integrals = { - IntegralType.cell: [ + IntegralType(0, 1): [ (0, k2.address, cells, active_coeffs), (1, k2.address, np.arange(0), active_coeffs), (2, k2.address, np.arange(0), active_coeffs), @@ -105,7 +105,7 @@ def test_numba_assembly(dtype): [V._cpp_object, V._cpp_object], integrals, [], [], False, [], mesh=mesh._cpp_object ) ) - integrals = {IntegralType.cell: [(0, k1.address, cells, active_coeffs)]} + integrals = {IntegralType(0, 1): [(0, k1.address, cells, active_coeffs)]} L = Form(formtype([V._cpp_object], integrals, [], [], False, [], mesh=mesh._cpp_object)) A = dolfinx.fem.assemble_matrix(a) @@ -136,7 +136,7 @@ def test_coefficient(dtype): num_cells = mesh.topology.index_map(tdim).size_local + mesh.topology.index_map(tdim).num_ghosts active_coeffs = np.array([0], dtype=np.int8) integrals = { - IntegralType.cell: [(0, k1.address, np.arange(num_cells, dtype=np.int32), active_coeffs)] + IntegralType(0, 1): [(0, k1.address, np.arange(num_cells, dtype=np.int32), active_coeffs)] } formtype = form_cpp_class(dtype) L = Form( @@ -276,7 +276,7 @@ def test_cffi_assembly(): ptrA = ffi.cast("intptr_t", ffi.addressof(lib, "tabulate_tensor_poissonA")) active_coeffs = np.array([], dtype=np.int8) - integrals = {IntegralType.cell: [(0, ptrA, cells, active_coeffs)]} + integrals = {IntegralType(0, 1): [(0, ptrA, cells, active_coeffs)]} a = Form( _cpp.fem.Form_float64( [V._cpp_object, V._cpp_object], integrals, [], [], False, [], mesh=mesh._cpp_object @@ -284,7 +284,7 @@ def test_cffi_assembly(): ) ptrL = ffi.cast("intptr_t", ffi.addressof(lib, "tabulate_tensor_poissonL")) - integrals = {IntegralType.cell: [(0, ptrL, cells, active_coeffs)]} + integrals = {IntegralType(0, 1): [(0, ptrL, cells, active_coeffs)]} L = Form( _cpp.fem.Form_float64([V._cpp_object], integrals, [], [], False, [], mesh=mesh._cpp_object) ) diff --git a/python/test/unit/fem/test_forms.py b/python/test/unit/fem/test_forms.py index f6b2865984e..b0677741fab 100644 --- a/python/test/unit/fem/test_forms.py +++ b/python/test/unit/fem/test_forms.py @@ -93,7 +93,7 @@ def test_incorrect_element(): [space._cpp_object, space._cpp_object], [], [], - {IntegralType.cell: []}, + {IntegralType(0, 1): []}, [], mesh._cpp_object, ) @@ -105,7 +105,7 @@ def test_incorrect_element(): [incorrect_space._cpp_object, incorrect_space._cpp_object], [], [], - {IntegralType.cell: []}, + {IntegralType(0, 1): []}, [], mesh._cpp_object, ) From 1ec3cf933b85fd867ee5966e110d2a6a37ed4adb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= Date: Wed, 10 Sep 2025 16:15:02 +0000 Subject: [PATCH 11/30] Fix codim --- cpp/dolfinx/fem/Form.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cpp/dolfinx/fem/Form.h b/cpp/dolfinx/fem/Form.h index 8416a231e7b..4081fad27df 100644 --- a/cpp/dolfinx/fem/Form.h +++ b/cpp/dolfinx/fem/Form.h @@ -286,7 +286,7 @@ class Form std::vector e; if (type.codim == 0) e = emap.sub_topology_to_topology(itg.entities, inverse); - else if (type.codim = 1) + else if (type.codim == 1) { const mesh::Topology topology = *_mesh->topology(); int tdim = topology.dim(); From 187d0b5eea4291b2ceacfbb3f2d19662d6f90e62 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= Date: Wed, 10 Sep 2025 16:16:57 +0000 Subject: [PATCH 12/30] Remove nanobind operators --- python/dolfinx/wrappers/fem.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/python/dolfinx/wrappers/fem.cpp b/python/dolfinx/wrappers/fem.cpp index 968094e7985..8f7fe9669fc 100644 --- a/python/dolfinx/wrappers/fem.cpp +++ b/python/dolfinx/wrappers/fem.cpp @@ -34,7 +34,6 @@ #include #include #include -#include #include #include #include From af02ec6d185b51dcc54adb38d0b7d600587a17d2 Mon Sep 17 00:00:00 2001 From: jorgensd Date: Wed, 10 Sep 2025 17:48:05 +0000 Subject: [PATCH 13/30] Add missing fem:: declaration --- cpp/demo/mixed_poisson/main.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/cpp/demo/mixed_poisson/main.cpp b/cpp/demo/mixed_poisson/main.cpp index a6f330affe7..41b6e825549 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( - IntegralType(1, 1), *mesh->topology(), dfacets); + fem::IntegralType(1, 1), *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< std::int32_t, std::vector>>> - subdomain_data{{IntegralType(1, 1), {{1, domains}}}}; + subdomain_data{{fem::IntegralType(1, 1), {{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 From a805bce66f66a69deb06057ba117e412bd416264 Mon Sep 17 00:00:00 2001 From: jorgensd Date: Wed, 10 Sep 2025 18:00:07 +0000 Subject: [PATCH 14/30] Revert change in demo --- cpp/demo/mixed_poisson/main.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/cpp/demo/mixed_poisson/main.cpp b/cpp/demo/mixed_poisson/main.cpp index 41b6e825549..fd700dc5eb4 100644 --- a/cpp/demo/mixed_poisson/main.cpp +++ b/cpp/demo/mixed_poisson/main.cpp @@ -278,9 +278,9 @@ int main(int argc, char* argv[]) // in the map), and exterior facet domain marked as '1' in the UFL // file, and `domains` holds the necessary data to perform // integration of selected facets. - std::map< - std::int32_t, - std::vector>>> + std::map>>> subdomain_data{{fem::IntegralType(1, 1), {{1, domains}}}}; // Since we are doing a `ds(1)` integral on mesh and `u0` is defined From d6b9577a7a33108055dc467a557cc8bb60dc5a17 Mon Sep 17 00:00:00 2001 From: jorgensd Date: Wed, 10 Sep 2025 18:01:22 +0000 Subject: [PATCH 15/30] Another revert --- cpp/demo/mixed_poisson/main.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cpp/demo/mixed_poisson/main.cpp b/cpp/demo/mixed_poisson/main.cpp index fd700dc5eb4..63b8d2ba8fe 100644 --- a/cpp/demo/mixed_poisson/main.cpp +++ b/cpp/demo/mixed_poisson/main.cpp @@ -278,7 +278,7 @@ int main(int argc, char* argv[]) // in the map), and exterior facet domain marked as '1' in the UFL // file, and `domains` holds the necessary data to perform // integration of selected facets. - std::map>>> subdomain_data{{fem::IntegralType(1, 1), {{1, domains}}}}; From 471a3574eb1cfa12938d254fc40252b18879fe23 Mon Sep 17 00:00:00 2001 From: jorgensd Date: Wed, 10 Sep 2025 18:20:43 +0000 Subject: [PATCH 16/30] Fix demo properly --- cpp/demo/mixed_poisson/main.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/cpp/demo/mixed_poisson/main.cpp b/cpp/demo/mixed_poisson/main.cpp index 63b8d2ba8fe..22d156706a3 100644 --- a/cpp/demo/mixed_poisson/main.cpp +++ b/cpp/demo/mixed_poisson/main.cpp @@ -278,9 +278,9 @@ int main(int argc, char* argv[]) // in the map), and exterior facet domain marked as '1' in the UFL // file, and `domains` holds the necessary data to perform // integration of selected facets. - std::map>>> + std::map< + fem::IntegralType, + std::vector>>> subdomain_data{{fem::IntegralType(1, 1), {{1, domains}}}}; // Since we are doing a `ds(1)` integral on mesh and `u0` is defined From 765733edd9ee7e7b9d8a6cd3963cdded12811962 Mon Sep 17 00:00:00 2001 From: jorgensd Date: Wed, 10 Sep 2025 18:30:17 +0000 Subject: [PATCH 17/30] Move less than operator into struct --- cpp/dolfinx/fem/CMakeLists.txt | 1 - cpp/dolfinx/fem/Form.cpp | 24 ------------------------ cpp/dolfinx/fem/Form.h | 30 ++++++++++++++++++++++-------- 3 files changed, 22 insertions(+), 33 deletions(-) delete mode 100644 cpp/dolfinx/fem/Form.cpp diff --git a/cpp/dolfinx/fem/CMakeLists.txt b/cpp/dolfinx/fem/CMakeLists.txt index d72c37a5390..548df96a011 100644 --- a/cpp/dolfinx/fem/CMakeLists.txt +++ b/cpp/dolfinx/fem/CMakeLists.txt @@ -34,7 +34,6 @@ target_sources( ${CMAKE_CURRENT_SOURCE_DIR}/DofMap.cpp ${CMAKE_CURRENT_SOURCE_DIR}/ElementDofLayout.cpp ${CMAKE_CURRENT_SOURCE_DIR}/FiniteElement.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/Form.cpp ${CMAKE_CURRENT_SOURCE_DIR}/dofmapbuilder.cpp ${CMAKE_CURRENT_SOURCE_DIR}/petsc.cpp ${CMAKE_CURRENT_SOURCE_DIR}/sparsitybuild.cpp diff --git a/cpp/dolfinx/fem/Form.cpp b/cpp/dolfinx/fem/Form.cpp deleted file mode 100644 index c8aa0854b83..00000000000 --- a/cpp/dolfinx/fem/Form.cpp +++ /dev/null @@ -1,24 +0,0 @@ -// Copyright (C) 2025 Jørgen S. Dokken -// -// This file is part of DOLFINx (https://www.fenicsproject.org) -// -// SPDX-License-Identifier: LGPL-3.0-or-later - -#include "Form.h" - -bool dolfinx::fem::operator<(const dolfinx::fem::IntegralType& self, - const dolfinx::fem::IntegralType& other) -{ - // First, compare codim - if (self.codim < other.codim) - { - return true; - } - // If codims are equal, compare num_cells - if (self.codim == other.codim) - { - return self.num_cells < other.num_cells; - } - // Otherwise, this object is not less than 'other' - return false; -} diff --git a/cpp/dolfinx/fem/Form.h b/cpp/dolfinx/fem/Form.h index 4081fad27df..905ed48e79a 100644 --- a/cpp/dolfinx/fem/Form.h +++ b/cpp/dolfinx/fem/Form.h @@ -48,10 +48,24 @@ struct IntegralType { return codim == other.codim && num_cells == other.num_cells; } -}; -// Overload the less-than operator -bool operator<(const IntegralType& self, const IntegralType& other); + /// @brief Less than operator. + /// Required for set operations. + /// @param other The other IntegralType to compare with + bool operator<(const IntegralType& other) const + { + // First, compare codim + if (codim < other.codim) + return true; + + // If codims are equal, compare num_cells + if (codim == other.codim) + return num_cells < other.num_cells; + + // Otherwise, this object is not less than 'other' + return false; + } +}; /// @brief Represents integral data, containing the kernel, and a list /// of entities to integrate over and the indicies of the coefficient @@ -140,9 +154,9 @@ class Form /// trial function spaces. /// @param[in] integrals Integrals in the form, where /// `integrals[IntegralType, i, kernel index]` returns the `i`th integral - /// (`integral_data`) of type `IntegralType` with kernel index `kernel index`. - /// The `i`-index refers to the position of a kernel when flattened by - /// sorted subdomain ids, sorted by subdomain ids. The subdomain ids can + /// (`integral_data`) of type `IntegralType` with kernel index `kernel + /// index`. The `i`-index refers to the position of a kernel when flattened + /// by sorted subdomain ids, sorted by subdomain ids. The subdomain ids can /// contain duplicate entries referring to different kernels over the same /// subdomain. /// @param[in] coefficients Coefficients in the form. @@ -516,8 +530,8 @@ class Form /// is marked with -1. /// /// @param[in] type Integral type. - /// @param[in] rank Argument index, e.g. `0` for the test function space, `1` - /// for the trial function space. + /// @param[in] rank Argument index, e.g. `0` for the test function space, + /// `1` for the trial function space. /// @param[in] idx Integral identifier. /// @param[in] kernel_idx Kernel index (cell type). /// @return Entity indices in the argument function space mesh that is From f0b5bc216100cb47452d04a8e761c73713016913 Mon Sep 17 00:00:00 2001 From: jorgensd Date: Wed, 10 Sep 2025 18:50:52 +0000 Subject: [PATCH 18/30] Simplify interface, set default number of cells to be 1 --- cpp/demo/codim_0_assembly/main.cpp | 5 ++-- cpp/demo/custom_kernel/main.cpp | 4 ++-- cpp/demo/mixed_poisson/main.cpp | 4 ++-- cpp/dolfinx/fem/Form.h | 5 ++++ cpp/dolfinx/fem/assemble_matrix_impl.h | 4 ++-- cpp/dolfinx/fem/assemble_scalar_impl.h | 6 ++--- cpp/dolfinx/fem/assemble_vector_impl.h | 10 ++++---- cpp/dolfinx/fem/pack.h | 10 ++++---- cpp/dolfinx/fem/utils.h | 24 +++++++++---------- python/demo/demo_static-condensation.py | 2 +- python/dolfinx/fem/forms.py | 24 ++++++++++--------- python/dolfinx/wrappers/fem.cpp | 2 +- .../test_assemble_mesh_independent_form.py | 4 ++-- python/test/unit/fem/test_assembler.py | 8 +++---- .../test/unit/fem/test_custom_jit_kernels.py | 10 ++++---- python/test/unit/fem/test_forms.py | 4 ++-- 16 files changed, 66 insertions(+), 60 deletions(-) diff --git a/cpp/demo/codim_0_assembly/main.cpp b/cpp/demo/codim_0_assembly/main.cpp index d227c01f24f..f4f4f85194e 100644 --- a/cpp/demo/codim_0_assembly/main.cpp +++ b/cpp/demo/codim_0_assembly/main.cpp @@ -100,12 +100,11 @@ int main(int argc, char* argv[]) // domain `mesh` std::vector integration_entities = fem::compute_integration_domains( - fem::IntegralType(0, 1), *mesh->topology(), cell_marker.find(2)); + fem::IntegralType(0), *mesh->topology(), cell_marker.find(2)); std::map< fem::IntegralType, std::vector>>> - subdomain_data - = {{fem::IntegralType(0, 1), {{3, integration_entities}}}}; + subdomain_data = {{fem::IntegralType(0), {{3, integration_entities}}}}; // A mixed-domain form involves functions defined over multiple // meshes. The mesh passed to `create_form` is called the diff --git a/cpp/demo/custom_kernel/main.cpp b/cpp/demo/custom_kernel/main.cpp index 8bdb1836e21..2aab6d7bccc 100644 --- a/cpp/demo/custom_kernel/main.cpp +++ b/cpp/demo/custom_kernel/main.cpp @@ -75,7 +75,7 @@ double assemble_matrix0(std::shared_ptr> V, { // Kernel data (ID, kernel function, cell indices to execute over) std::map integrals{ - std::pair{std::tuple{fem::IntegralType(0, 1), 0, 0}, + std::pair{std::tuple{fem::IntegralType(0), 0, 0}, fem::integral_data(kernel, cells, std::vector{})}}; fem::Form a({V, V}, integrals, V->mesh(), {}, {}, false, {}); @@ -105,7 +105,7 @@ double assemble_vector0(std::shared_ptr> V, { auto mesh = V->mesh(); std::map integrals{ - std::pair{std::tuple{fem::IntegralType(0, 1), 0, 0}, + std::pair{std::tuple{fem::IntegralType(0), 0, 0}, fem::integral_data(kernel, cells, std::vector{})}}; fem::Form L({V}, integrals, mesh, {}, {}, false, {}); auto dofmap = V->dofmap(); diff --git a/cpp/demo/mixed_poisson/main.cpp b/cpp/demo/mixed_poisson/main.cpp index 22d156706a3..96acd340905 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(1, 1), *mesh->topology(), dfacets); + fem::IntegralType(1), *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(1, 1), {{1, domains}}}}; + subdomain_data{{fem::IntegralType(1), {{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 905ed48e79a..09591bce727 100644 --- a/cpp/dolfinx/fem/Form.h +++ b/cpp/dolfinx/fem/Form.h @@ -41,6 +41,11 @@ struct IntegralType std::int32_t codim; ///< Codimension of the integral std::int32_t num_cells; ///< Number of cells in the integral + IntegralType(std::int32_t codim, std::int32_t num_cells = 1) + : codim(codim), num_cells(num_cells) + { + } // Default constructor + /// @brief Equality operator for integral type /// @param other The other IntegralType to compare with /// @return True if the integral types are equal, false otherwise diff --git a/cpp/dolfinx/fem/assemble_matrix_impl.h b/cpp/dolfinx/fem/assemble_matrix_impl.h index 9ee0c1cbf79..31c1e37f490 100644 --- a/cpp/dolfinx/fem/assemble_matrix_impl.h +++ b/cpp/dolfinx/fem/assemble_matrix_impl.h @@ -588,7 +588,7 @@ void assemble_matrix( cell_info0 = std::span(mesh0->topology()->get_cell_permutation_info()); cell_info1 = std::span(mesh1->topology()->get_cell_permutation_info()); } - IntegralType cell_integral = IntegralType(0, 1); + IntegralType cell_integral = IntegralType(0); for (int i = 0; i < a.num_integrals(cell_integral, cell_type_idx); ++i) { auto fn = a.kernel(cell_integral, i, cell_type_idx); @@ -618,7 +618,7 @@ void assemble_matrix( num_facets_per_cell); } - IntegralType facet = IntegralType(1, 1); + IntegralType facet = IntegralType(1); for (int i = 0; i < a.num_integrals(facet, cell_type_idx); ++i) { if (num_cell_types > 1) diff --git a/cpp/dolfinx/fem/assemble_scalar_impl.h b/cpp/dolfinx/fem/assemble_scalar_impl.h index a7146eb9ac0..fbdaf0f31d8 100644 --- a/cpp/dolfinx/fem/assemble_scalar_impl.h +++ b/cpp/dolfinx/fem/assemble_scalar_impl.h @@ -201,7 +201,7 @@ T assemble_scalar( std::vector> cdofs_b(2 * 3 * x_dofmap.extent(1)); T value = 0; - IntegralType cell_integral = IntegralType(0, 1); + IntegralType cell_integral = IntegralType(0); for (int i = 0; i < M.num_integrals(cell_integral, 0); ++i) { auto fn = M.kernel(cell_integral, i, 0); @@ -227,7 +227,7 @@ T assemble_scalar( num_facets_per_cell); } - IntegralType facet_type = IntegralType(1, 1); + IntegralType facet_type = IntegralType(1); for (int i = 0; i < M.num_integrals(facet_type, 0); ++i) { auto fn = M.kernel(facet_type, i, 0); @@ -278,7 +278,7 @@ T assemble_scalar( } std::int32_t tdim = mesh->topology()->dim(); - IntegralType vertex_type = IntegralType(tdim, 1); + IntegralType vertex_type = IntegralType(tdim); for (int i = 0; i < M.num_integrals(vertex_type, 0); ++i) { auto fn = M.kernel(vertex_type, i, 0); diff --git a/cpp/dolfinx/fem/assemble_vector_impl.h b/cpp/dolfinx/fem/assemble_vector_impl.h index 77a86d9ca8d..624e1100183 100644 --- a/cpp/dolfinx/fem/assemble_vector_impl.h +++ b/cpp/dolfinx/fem/assemble_vector_impl.h @@ -1100,7 +1100,7 @@ void lift_bc(V&& b, const Form& a, mdspan2_t x_dofmap, = element1->template dof_transformation_right_fn( doftransform::transpose); - IntegralType cell_integral_type = IntegralType(0, 1); + IntegralType cell_integral_type = IntegralType(0); for (int i = 0; i < a.num_integrals(cell_integral_type, 0); ++i) { auto kernel = a.kernel(cell_integral_type, i, 0); @@ -1146,7 +1146,7 @@ void lift_bc(V&& b, const Form& a, mdspan2_t x_dofmap, num_facets_per_cell); } - IntegralType facet_integral_type = IntegralType(1, 1); + IntegralType facet_integral_type = IntegralType(1); for (int i = 0; i < a.num_integrals(facet_integral_type, 0); ++i) { auto kernel = a.kernel(facet_integral_type, i, 0); @@ -1342,7 +1342,7 @@ void assemble_vector( cell_info0 = std::span(mesh0->topology()->get_cell_permutation_info()); } - IntegralType cell_integral_type = IntegralType(0, 1); + IntegralType cell_integral_type = IntegralType(0); for (int i = 0; i < L.num_integrals(cell_integral_type, 0); ++i) { auto fn = L.kernel(cell_integral_type, i, cell_type_idx); @@ -1388,7 +1388,7 @@ void assemble_vector( = md::mdspan>; - IntegralType facet_integral_type = IntegralType(1, 1); + IntegralType facet_integral_type = IntegralType(1); for (int i = 0; i < L.num_integrals(facet_integral_type, 0); ++i) { auto fn = L.kernel(facet_integral_type, i, 0); @@ -1474,7 +1474,7 @@ void assemble_vector( } } - IntegralType vertex_integral_type = IntegralType(2, 1); + IntegralType vertex_integral_type = IntegralType(mesh->topology()->dim()); for (int i = 0; i < L.num_integrals(vertex_integral_type, 0); ++i) { auto fn = L.kernel(vertex_integral_type, i, 0); diff --git a/cpp/dolfinx/fem/pack.h b/cpp/dolfinx/fem/pack.h index 07b747cf373..2d0e0c24274 100644 --- a/cpp/dolfinx/fem/pack.h +++ b/cpp/dolfinx/fem/pack.h @@ -237,7 +237,7 @@ void pack_coefficients(const Form& form, if (integral_type.codim == 0) { // Iterate over coefficients that are active in cell integrals - for (int coeff : form.active_coeffs(IntegralType(0, 1), id)) + for (int coeff : form.active_coeffs(IntegralType(0), id)) { // Get coefficient mesh auto mesh = coefficients[coeff]->function_space()->mesh(); @@ -255,7 +255,7 @@ void pack_coefficients(const Form& form, } std::span cells_b - = form.domain_coeff(IntegralType(0, 1), id, coeff); + = form.domain_coeff(IntegralType(0), id, coeff); md::mdspan cells(cells_b.data(), cells_b.size()); std::span cell_info = impl::get_cell_orientation_info(*coefficients[coeff]); @@ -269,11 +269,11 @@ void pack_coefficients(const Form& form, // Iterate over coefficients coefficients that are active in // exterior facet integrals for (int coeff : - form.active_coeffs(IntegralType(integral_type.codim, 1), id)) + form.active_coeffs(IntegralType(integral_type.codim), id)) { auto mesh = coefficients[coeff]->function_space()->mesh(); - std::span entities_b = form.domain_coeff( - IntegralType(integral_type.codim, 1), id, coeff); + std::span entities_b + = form.domain_coeff(IntegralType(integral_type.codim), id, coeff); md::mdspan> entities(entities_b.data(), entities_b.size() / 2, 2); diff --git a/cpp/dolfinx/fem/utils.h b/cpp/dolfinx/fem/utils.h index facf3ae9aa5..3f47c23fef0 100644 --- a/cpp/dolfinx/fem/utils.h +++ b/cpp/dolfinx/fem/utils.h @@ -232,7 +232,7 @@ void build_sparsity_pattern(la::SparsityPattern& pattern, const Form& a) const std::set> types = a.integral_types(); - if (types.find(IntegralType(1, 1)) != types.end() + if (types.find(IntegralType(1)) != types.end() or types.find(IntegralType(1, 2)) != types.end()) { // FIXME: cleanup these calls? Some of the happen internally again. @@ -263,7 +263,7 @@ void build_sparsity_pattern(la::SparsityPattern& pattern, const Form& a) // Create and build sparsity pattern for (auto type : types) { - if (type == IntegralType(0, 1)) + if (type == IntegralType(0)) { for (int i = 0; i < a.num_integrals(type, cell_type_idx); ++i) { @@ -284,7 +284,7 @@ void build_sparsity_pattern(la::SparsityPattern& pattern, const Form& a) {{dofmaps[0], dofmaps[1]}}); } } - else if (type == IntegralType(1, 1)) + else if (type == IntegralType(1)) { for (int i = 0; i < a.num_integrals(type, cell_type_idx); ++i) { @@ -527,7 +527,7 @@ Form create_form_factory( std::span ids(ufcx_forms[0].get().form_integral_ids + integral_offsets[cell], num_integrals_type[cell]); - auto sd = subdomains.find(IntegralType(0, 1)); + auto sd = subdomains.find(IntegralType(0)); for (std::size_t form_idx = 0; form_idx < ufcx_forms.size(); ++form_idx) { const ufcx_form& ufcx_form = ufcx_forms[form_idx]; @@ -562,7 +562,7 @@ Form create_form_factory( default_cells.resize( topology->index_maps(tdim).at(form_idx)->size_local(), 0); std::iota(default_cells.begin(), default_cells.end(), 0); - integrals.insert({{IntegralType(0, 1), i, form_idx}, + integrals.insert({{IntegralType(0), i, form_idx}, {k, default_cells, active_coeffs}}); } else if (sd != subdomains.end()) @@ -572,7 +572,7 @@ Form create_form_factory( [](auto& a) { return a.first; }); if (it != sd->second.end() and it->first == id) { - integrals.insert({{IntegralType(0, 1), i, form_idx}, + integrals.insert({{IntegralType(0), i, form_idx}, {k, std::vector(it->second.begin(), it->second.end()), @@ -592,7 +592,7 @@ Form create_form_factory( std::span ids(ufcx_forms[0].get().form_integral_ids + integral_offsets[facet], num_integrals_type[facet]); - auto sd = subdomains.find(IntegralType(1, 1)); + auto sd = subdomains.find(IntegralType(1)); for (std::size_t form_idx = 0; form_idx < ufcx_forms.size(); ++form_idx) { const ufcx_form& ufcx_form = ufcx_forms[form_idx]; @@ -631,7 +631,7 @@ Form create_form_factory( default_facets_ext.insert(default_facets_ext.end(), pair.begin(), pair.end()); } - integrals.insert({{IntegralType(1, 1), i, form_idx}, + integrals.insert({{IntegralType(1), i, form_idx}, {k, default_facets_ext, active_coeffs}}); } else if (sd != subdomains.end()) @@ -641,7 +641,7 @@ Form create_form_factory( [](auto& a) { return a.first; }); if (it != sd->second.end() and it->first == id) { - integrals.insert({{IntegralType(1, 1), i, form_idx}, + integrals.insert({{IntegralType(1), i, form_idx}, {k, std::vector(it->second.begin(), it->second.end()), @@ -759,7 +759,7 @@ Form create_form_factory( std::span ids(form.form_integral_ids + integral_offsets[vertex], num_integrals_type[vertex]); - auto sd = subdomains.find(IntegralType(tdim, 1)); + auto sd = subdomains.find(IntegralType(tdim)); for (int i = 0; i < num_integrals_type[vertex]; ++i) { const int id = ids[i]; @@ -810,7 +810,7 @@ Form create_form_factory( std::vector cells_and_vertices = get_cells_and_vertices( std::ranges::views::iota(0, num_vertices)); - integrals.insert({{IntegralType(tdim, 1), i, form_idx}, + integrals.insert({{IntegralType(tdim), i, form_idx}, {k, cells_and_vertices, active_coeffs}}); } else if (sd != subdomains.end()) @@ -820,7 +820,7 @@ Form create_form_factory( [](auto& a) { return a.first; }); if (it != sd->second.end() and it->first == id) { - integrals.insert({{IntegralType(tdim, 1), i, form_idx}, + integrals.insert({{IntegralType(tdim), i, form_idx}, {k, std::vector(it->second.begin(), it->second.end()), diff --git a/python/demo/demo_static-condensation.py b/python/demo/demo_static-condensation.py index 815d245268b..1d296217015 100644 --- a/python/demo/demo_static-condensation.py +++ b/python/demo/demo_static-condensation.py @@ -183,7 +183,7 @@ def tabulate_A(A_, w_, c_, coords_, entity_local_index, permutation=ffi.NULL, cu # Prepare a Form with a condensed tabulation kernel formtype = form_cpp_class(PETSc.ScalarType) # type: ignore cells = np.arange(msh.topology.index_map(msh.topology.dim).size_local) -integrals = {IntegralType(0, 1): [(0, tabulate_A.address, cells, np.array([], dtype=np.int8))]} +integrals = {IntegralType(0): [(0, tabulate_A.address, cells, np.array([], dtype=np.int8))]} a_cond = Form( formtype([U._cpp_object, U._cpp_object], integrals, [], [], False, [], mesh=msh._cpp_object) ) diff --git a/python/dolfinx/fem/forms.py b/python/dolfinx/fem/forms.py index 1cde6830b46..51d698e65d5 100644 --- a/python/dolfinx/fem/forms.py +++ b/python/dolfinx/fem/forms.py @@ -213,11 +213,12 @@ def form_cpp_class( raise NotImplementedError(f"Type {dtype} not supported.") -_ufl_to_dolfinx_domain = { - "cell": IntegralType(0, 1), - "exterior_facet": IntegralType(1, 1), - "interior_facet": IntegralType(1, 2), - "vertex": IntegralType(2, 1), +def _ufl_to_dolfinx_domain(tdim:int) -> dict[str, IntegralType]: + return{ + "cell": IntegralType(0), + "exterior_facet": IntegralType(1), + "interior_facet": IntegralType(1, 2), + "vertex": IntegralType(tdim), } @@ -365,12 +366,13 @@ def _form(form): for i in range(ufcx_form.num_coefficients) ] constants = [c._cpp_object for c in form.constants()] - # Extract subdomain ids from ufcx_form - subdomain_ids = {_ufl_to_dolfinx_domain[type]: [] for type in sd.get(domain).keys()} + tdim = domain.topological_dimension() + ufl_to_itg = _ufl_to_dolfinx_domain(tdim) + subdomain_ids = {ufl_to_itg[type]: [] for type in sd.get(domain).keys()} integral_offsets = [ufcx_form.form_integral_offsets[i] for i in range(5)] integral_types = [ - _ufl_to_dolfinx_domain[type] + ufl_to_itg[type] for type in ["cell", "exterior_facet", "interior_facet", "vertex"] ] for i in range(len(integral_offsets) - 1): @@ -380,10 +382,10 @@ def _form(form): # Subdomain markers (possibly empty list for some integral types) subdomains = { - _ufl_to_dolfinx_domain[key]: get_integration_domains( - _ufl_to_dolfinx_domain[key], + ufl_to_itg[key]: get_integration_domains( + ufl_to_itg[key], subdomain_data[0], - subdomain_ids[_ufl_to_dolfinx_domain[key]], + subdomain_ids[ufl_to_itg[key]], ) for (key, subdomain_data) in sd.get(domain).items() } diff --git a/python/dolfinx/wrappers/fem.cpp b/python/dolfinx/wrappers/fem.cpp index 8f7fe9669fc..0e5bc43cf08 100644 --- a/python/dolfinx/wrappers/fem.cpp +++ b/python/dolfinx/wrappers/fem.cpp @@ -1271,7 +1271,7 @@ void fem(nb::module_& m) nb::class_(m, "_IntegralType") .def(nb::init(), - nb::arg("codim"), nb::arg("num_cells")) + nb::arg("codim"), nb::arg("num_cells") = 1) .def_ro("codim", &dolfinx::fem::IntegralType::codim) .def_ro("num_cells", &dolfinx::fem::IntegralType::num_cells) .def("__eq__", &dolfinx::fem::IntegralType::operator==); 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 e5baf3f2a2e..cde088b5a04 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(1, 1), mesh.topology, facets + dolfinx.fem.IntegralType(1), mesh.topology, facets ) - subdomains = {dolfinx.fem.IntegralType(1, 1): [(subdomain_id, facet_entities)]} + subdomains = {dolfinx.fem.IntegralType(1): [(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 49d27310f81..1c9cecfffbf 100644 --- a/python/test/unit/fem/test_assembler.py +++ b/python/test/unit/fem/test_assembler.py @@ -1712,8 +1712,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(1, 1), msh.topology, vertices) - subdomains = {fem.IntegralType(1, 1): [(0, cell_vertex_pairs)]} + fem.compute_integration_domains(fem.IntegralType(1), msh.topology, vertices) + subdomains = {fem.IntegralType(1): [(0, cell_vertex_pairs)]} compiled_form = fem.compile_form( comm, x[0] * ufl.dP, form_compiler_options={"scalar_type": dtype} @@ -1852,8 +1852,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(1, 1), msh.topology, vertices) - subdomains = {fem.IntegralType(1, 1): [(0, cell_vertex_pairs)]} + fem.compute_integration_domains(fem.IntegralType(1), msh.topology, vertices) + subdomains = {fem.IntegralType(1): [(0, cell_vertex_pairs)]} compiled_form = fem.compile_form( comm, x[0] * v * ufl.dP, form_compiler_options={"scalar_type": dtype} diff --git a/python/test/unit/fem/test_custom_jit_kernels.py b/python/test/unit/fem/test_custom_jit_kernels.py index f70cdec3b00..e0d83c34d50 100644 --- a/python/test/unit/fem/test_custom_jit_kernels.py +++ b/python/test/unit/fem/test_custom_jit_kernels.py @@ -93,7 +93,7 @@ def test_numba_assembly(dtype): cells = np.arange(mesh.topology.index_map(mesh.topology.dim).size_local, dtype=np.int32) active_coeffs = np.array([], dtype=np.int8) integrals = { - IntegralType(0, 1): [ + IntegralType(0): [ (0, k2.address, cells, active_coeffs), (1, k2.address, np.arange(0), active_coeffs), (2, k2.address, np.arange(0), active_coeffs), @@ -105,7 +105,7 @@ def test_numba_assembly(dtype): [V._cpp_object, V._cpp_object], integrals, [], [], False, [], mesh=mesh._cpp_object ) ) - integrals = {IntegralType(0, 1): [(0, k1.address, cells, active_coeffs)]} + integrals = {IntegralType(0): [(0, k1.address, cells, active_coeffs)]} L = Form(formtype([V._cpp_object], integrals, [], [], False, [], mesh=mesh._cpp_object)) A = dolfinx.fem.assemble_matrix(a) @@ -136,7 +136,7 @@ def test_coefficient(dtype): num_cells = mesh.topology.index_map(tdim).size_local + mesh.topology.index_map(tdim).num_ghosts active_coeffs = np.array([0], dtype=np.int8) integrals = { - IntegralType(0, 1): [(0, k1.address, np.arange(num_cells, dtype=np.int32), active_coeffs)] + IntegralType(0): [(0, k1.address, np.arange(num_cells, dtype=np.int32), active_coeffs)] } formtype = form_cpp_class(dtype) L = Form( @@ -276,7 +276,7 @@ def test_cffi_assembly(): ptrA = ffi.cast("intptr_t", ffi.addressof(lib, "tabulate_tensor_poissonA")) active_coeffs = np.array([], dtype=np.int8) - integrals = {IntegralType(0, 1): [(0, ptrA, cells, active_coeffs)]} + integrals = {IntegralType(0): [(0, ptrA, cells, active_coeffs)]} a = Form( _cpp.fem.Form_float64( [V._cpp_object, V._cpp_object], integrals, [], [], False, [], mesh=mesh._cpp_object @@ -284,7 +284,7 @@ def test_cffi_assembly(): ) ptrL = ffi.cast("intptr_t", ffi.addressof(lib, "tabulate_tensor_poissonL")) - integrals = {IntegralType(0, 1): [(0, ptrL, cells, active_coeffs)]} + integrals = {IntegralType(0): [(0, ptrL, cells, active_coeffs)]} L = Form( _cpp.fem.Form_float64([V._cpp_object], integrals, [], [], False, [], mesh=mesh._cpp_object) ) diff --git a/python/test/unit/fem/test_forms.py b/python/test/unit/fem/test_forms.py index b0677741fab..b218d47e6e6 100644 --- a/python/test/unit/fem/test_forms.py +++ b/python/test/unit/fem/test_forms.py @@ -93,7 +93,7 @@ def test_incorrect_element(): [space._cpp_object, space._cpp_object], [], [], - {IntegralType(0, 1): []}, + {IntegralType(0): []}, [], mesh._cpp_object, ) @@ -105,7 +105,7 @@ def test_incorrect_element(): [incorrect_space._cpp_object, incorrect_space._cpp_object], [], [], - {IntegralType(0, 1): []}, + {IntegralType(0): []}, [], mesh._cpp_object, ) From bbe4c137978b348bbf988f2704e0ed04827794e2 Mon Sep 17 00:00:00 2001 From: jorgensd Date: Wed, 10 Sep 2025 19:03:50 +0000 Subject: [PATCH 19/30] Add docstring in ufl_to_dolfinx_domain --- python/dolfinx/fem/forms.py | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/python/dolfinx/fem/forms.py b/python/dolfinx/fem/forms.py index 51d698e65d5..11a846fdd28 100644 --- a/python/dolfinx/fem/forms.py +++ b/python/dolfinx/fem/forms.py @@ -213,13 +213,16 @@ def form_cpp_class( raise NotImplementedError(f"Type {dtype} not supported.") -def _ufl_to_dolfinx_domain(tdim:int) -> dict[str, IntegralType]: - return{ +def _ufl_to_dolfinx_domain(tdim: int) -> dict[str, IntegralType]: + """Return a mapping from the UFL integral types to the corresponding + `dolfinx.fem.IntegralType` for a given topological dimension. + """ + return { "cell": IntegralType(0), "exterior_facet": IntegralType(1), "interior_facet": IntegralType(1, 2), "vertex": IntegralType(tdim), -} + } def mixed_topology_form( @@ -372,8 +375,7 @@ def _form(form): subdomain_ids = {ufl_to_itg[type]: [] for type in sd.get(domain).keys()} integral_offsets = [ufcx_form.form_integral_offsets[i] for i in range(5)] integral_types = [ - ufl_to_itg[type] - for type in ["cell", "exterior_facet", "interior_facet", "vertex"] + ufl_to_itg[type] for type in ["cell", "exterior_facet", "interior_facet", "vertex"] ] for i in range(len(integral_offsets) - 1): integral_type = integral_types[i] From 27709e673a3963d1fc1c3f940e2ba9cc9b5c4dc8 Mon Sep 17 00:00:00 2001 From: jorgensd Date: Wed, 10 Sep 2025 19:19:26 +0000 Subject: [PATCH 20/30] Add docs for default constructor --- cpp/dolfinx/fem/Form.h | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/cpp/dolfinx/fem/Form.h b/cpp/dolfinx/fem/Form.h index 09591bce727..3814d177020 100644 --- a/cpp/dolfinx/fem/Form.h +++ b/cpp/dolfinx/fem/Form.h @@ -41,10 +41,13 @@ struct IntegralType std::int32_t codim; ///< Codimension of the integral std::int32_t num_cells; ///< Number of cells in the integral + /// @brief Create an integral type. + /// @param codim Codimension of the integral. + /// @param num_cells Number of cells in the integral. Default is 1. IntegralType(std::int32_t codim, std::int32_t num_cells = 1) : codim(codim), num_cells(num_cells) { - } // Default constructor + } /// @brief Equality operator for integral type /// @param other The other IntegralType to compare with From bba5adc51a0431b9b795666bff2e3dceafb34d27 Mon Sep 17 00:00:00 2001 From: jorgensd Date: Wed, 10 Sep 2025 20:08:14 +0000 Subject: [PATCH 21/30] Deactivate vertex integrals tests in 1D as their default integration domain is strange is not consistent with ds or dS (it is the sum of both). --- cpp/dolfinx/fem/pack.h | 15 +++++++-------- python/test/unit/fem/test_assembler.py | 4 ++-- 2 files changed, 9 insertions(+), 10 deletions(-) diff --git a/cpp/dolfinx/fem/pack.h b/cpp/dolfinx/fem/pack.h index 2d0e0c24274..8cde37e387c 100644 --- a/cpp/dolfinx/fem/pack.h +++ b/cpp/dolfinx/fem/pack.h @@ -237,7 +237,7 @@ void pack_coefficients(const Form& form, if (integral_type.codim == 0) { // Iterate over coefficients that are active in cell integrals - for (int coeff : form.active_coeffs(IntegralType(0), id)) + for (int coeff : form.active_coeffs(integral_type, id)) { // Get coefficient mesh auto mesh = coefficients[coeff]->function_space()->mesh(); @@ -255,7 +255,7 @@ void pack_coefficients(const Form& form, } std::span cells_b - = form.domain_coeff(IntegralType(0), id, coeff); + = form.domain_coeff(integral_type, id, coeff); md::mdspan cells(cells_b.data(), cells_b.size()); std::span cell_info = impl::get_cell_orientation_info(*coefficients[coeff]); @@ -267,13 +267,12 @@ void pack_coefficients(const Form& form, else if (integral_type.num_cells == 1) { // Iterate over coefficients coefficients that are active in - // exterior facet integrals - for (int coeff : - form.active_coeffs(IntegralType(integral_type.codim), id)) + // facet or vertex integrals. + for (int coeff : form.active_coeffs(integral_type, id)) { auto mesh = coefficients[coeff]->function_space()->mesh(); std::span entities_b - = form.domain_coeff(IntegralType(integral_type.codim), id, coeff); + = form.domain_coeff(integral_type, id, coeff); md::mdspan> entities(entities_b.data(), entities_b.size() / 2, 2); @@ -290,11 +289,11 @@ void pack_coefficients(const Form& form, { // Iterate over coefficients that are active in interior // facet integrals - for (int coeff : form.active_coeffs(IntegralType(1, 2), id)) + for (int coeff : form.active_coeffs(integral_type, id)) { auto mesh = coefficients[coeff]->function_space()->mesh(); std::span facets_b - = form.domain_coeff(IntegralType(1, 2), id, coeff); + = form.domain_coeff(integral_type, id, coeff); md::mdspan> facets(facets_b.data(), facets_b.size() / 4, 4); diff --git a/python/test/unit/fem/test_assembler.py b/python/test/unit/fem/test_assembler.py index 1c9cecfffbf..eb8db84552b 100644 --- a/python/test/unit/fem/test_assembler.py +++ b/python/test/unit/fem/test_assembler.py @@ -1591,7 +1591,7 @@ def vertex_to_dof_map(V): @pytest.mark.parametrize( "cell_type", [ - mesh.CellType.interval, + #mesh.CellType.interval, mesh.CellType.triangle, mesh.CellType.quadrilateral, mesh.CellType.tetrahedron, @@ -1727,7 +1727,7 @@ def check_vertex_integral_against_sum(form, vertices, weighted=False): @pytest.mark.parametrize( "cell_type", [ - mesh.CellType.interval, + # mesh.CellType.interval, mesh.CellType.triangle, mesh.CellType.quadrilateral, mesh.CellType.tetrahedron, From 0930b945ae92bafce112b7aa609cc0e5c8356d8d Mon Sep 17 00:00:00 2001 From: jorgensd Date: Wed, 10 Sep 2025 20:09:52 +0000 Subject: [PATCH 22/30] Ruff formatting --- 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 eb8db84552b..8dc8251cc84 100644 --- a/python/test/unit/fem/test_assembler.py +++ b/python/test/unit/fem/test_assembler.py @@ -1591,7 +1591,7 @@ def vertex_to_dof_map(V): @pytest.mark.parametrize( "cell_type", [ - #mesh.CellType.interval, + # mesh.CellType.interval, mesh.CellType.triangle, mesh.CellType.quadrilateral, mesh.CellType.tetrahedron, From 396e1a655142fa435ae68a8a4ce11ac0d625586c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= Date: Thu, 11 Sep 2025 06:50:46 +0000 Subject: [PATCH 23/30] Fix vertex integrals --- cpp/dolfinx/fem/assemble_scalar_impl.h | 2 +- cpp/dolfinx/fem/utils.cpp | 2 +- cpp/dolfinx/fem/utils.h | 6 +++--- python/dolfinx/fem/forms.py | 6 +++--- python/test/unit/fem/test_assembler.py | 4 ++-- 5 files changed, 10 insertions(+), 10 deletions(-) diff --git a/cpp/dolfinx/fem/assemble_scalar_impl.h b/cpp/dolfinx/fem/assemble_scalar_impl.h index fbdaf0f31d8..04b9df59b11 100644 --- a/cpp/dolfinx/fem/assemble_scalar_impl.h +++ b/cpp/dolfinx/fem/assemble_scalar_impl.h @@ -278,7 +278,7 @@ T assemble_scalar( } std::int32_t tdim = mesh->topology()->dim(); - IntegralType vertex_type = IntegralType(tdim); + IntegralType vertex_type = IntegralType(-1); for (int i = 0; i < M.num_integrals(vertex_type, 0); ++i) { auto fn = M.kernel(vertex_type, i, 0); diff --git a/cpp/dolfinx/fem/utils.cpp b/cpp/dolfinx/fem/utils.cpp index 7ef7a8f9bb8..513d945754b 100644 --- a/cpp/dolfinx/fem/utils.cpp +++ b/cpp/dolfinx/fem/utils.cpp @@ -239,7 +239,7 @@ fem::compute_integration_domains(const fem::IntegralType& integral_type, } else { - throw std::runtime_error("Unsupported integral type with codim " + throw std::runtime_error("Cannot compute integration domains with codim " + std::to_string(integral_type.codim) + " and num_cells " + std::to_string(integral_type.num_cells)); diff --git a/cpp/dolfinx/fem/utils.h b/cpp/dolfinx/fem/utils.h index 3f47c23fef0..e4cc8f331c9 100644 --- a/cpp/dolfinx/fem/utils.h +++ b/cpp/dolfinx/fem/utils.h @@ -759,7 +759,7 @@ Form create_form_factory( std::span ids(form.form_integral_ids + integral_offsets[vertex], num_integrals_type[vertex]); - auto sd = subdomains.find(IntegralType(tdim)); + auto sd = subdomains.find(IntegralType(-1)); for (int i = 0; i < num_integrals_type[vertex]; ++i) { const int id = ids[i]; @@ -810,7 +810,7 @@ Form create_form_factory( std::vector cells_and_vertices = get_cells_and_vertices( std::ranges::views::iota(0, num_vertices)); - integrals.insert({{IntegralType(tdim), i, form_idx}, + integrals.insert({{IntegralType(-1), i, form_idx}, {k, cells_and_vertices, active_coeffs}}); } else if (sd != subdomains.end()) @@ -820,7 +820,7 @@ Form create_form_factory( [](auto& a) { return a.first; }); if (it != sd->second.end() and it->first == id) { - integrals.insert({{IntegralType(tdim), i, form_idx}, + integrals.insert({{IntegralType(-1), i, form_idx}, {k, std::vector(it->second.begin(), it->second.end()), diff --git a/python/dolfinx/fem/forms.py b/python/dolfinx/fem/forms.py index 11a846fdd28..51fc63f5b79 100644 --- a/python/dolfinx/fem/forms.py +++ b/python/dolfinx/fem/forms.py @@ -213,7 +213,7 @@ def form_cpp_class( raise NotImplementedError(f"Type {dtype} not supported.") -def _ufl_to_dolfinx_domain(tdim: int) -> dict[str, IntegralType]: +def _ufl_to_dolfinx_domain() -> dict[str, IntegralType]: """Return a mapping from the UFL integral types to the corresponding `dolfinx.fem.IntegralType` for a given topological dimension. """ @@ -221,7 +221,7 @@ def _ufl_to_dolfinx_domain(tdim: int) -> dict[str, IntegralType]: "cell": IntegralType(0), "exterior_facet": IntegralType(1), "interior_facet": IntegralType(1, 2), - "vertex": IntegralType(tdim), + "vertex": IntegralType(-1), } @@ -371,7 +371,7 @@ def _form(form): constants = [c._cpp_object for c in form.constants()] # Extract subdomain ids from ufcx_form tdim = domain.topological_dimension() - ufl_to_itg = _ufl_to_dolfinx_domain(tdim) + ufl_to_itg = _ufl_to_dolfinx_domain() subdomain_ids = {ufl_to_itg[type]: [] for type in sd.get(domain).keys()} integral_offsets = [ufcx_form.form_integral_offsets[i] for i in range(5)] integral_types = [ diff --git a/python/test/unit/fem/test_assembler.py b/python/test/unit/fem/test_assembler.py index 8dc8251cc84..1c9cecfffbf 100644 --- a/python/test/unit/fem/test_assembler.py +++ b/python/test/unit/fem/test_assembler.py @@ -1591,7 +1591,7 @@ def vertex_to_dof_map(V): @pytest.mark.parametrize( "cell_type", [ - # mesh.CellType.interval, + mesh.CellType.interval, mesh.CellType.triangle, mesh.CellType.quadrilateral, mesh.CellType.tetrahedron, @@ -1727,7 +1727,7 @@ def check_vertex_integral_against_sum(form, vertices, weighted=False): @pytest.mark.parametrize( "cell_type", [ - # mesh.CellType.interval, + mesh.CellType.interval, mesh.CellType.triangle, mesh.CellType.quadrilateral, mesh.CellType.tetrahedron, From e20c7f3db6e3db61913639660606faf235262d49 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= Date: Thu, 11 Sep 2025 06:54:02 +0000 Subject: [PATCH 24/30] Fix further vertex integrals --- cpp/dolfinx/fem/utils.cpp | 2 +- python/dolfinx/fem/forms.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/cpp/dolfinx/fem/utils.cpp b/cpp/dolfinx/fem/utils.cpp index 513d945754b..d784ac583b0 100644 --- a/cpp/dolfinx/fem/utils.cpp +++ b/cpp/dolfinx/fem/utils.cpp @@ -226,7 +226,7 @@ fem::compute_integration_domains(const fem::IntegralType& integral_type, } } } - else if (integral_type.codim == tdim and integral_type.num_cells == 1) + else if (integral_type.codim == -1 and integral_type.num_cells == 1) { auto [v_to_c, c_to_v] = get_connectivities(0); for (auto vertex : entities) diff --git a/python/dolfinx/fem/forms.py b/python/dolfinx/fem/forms.py index 51fc63f5b79..48da35fe7dd 100644 --- a/python/dolfinx/fem/forms.py +++ b/python/dolfinx/fem/forms.py @@ -161,7 +161,7 @@ def get_integration_domains( if integral_type.codim == 1: subdomain._cpp_object.topology.create_connectivity(tdim - 1, tdim) subdomain._cpp_object.topology.create_connectivity(tdim, tdim - 1) - if integral_type.codim == tdim: + if integral_type.codim == -1: subdomain._cpp_object.topology.create_connectivity(0, tdim) subdomain._cpp_object.topology.create_connectivity(tdim, 0) From d010a719112ab174c86cd02ce90fe1df788f433f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= Date: Thu, 11 Sep 2025 07:46:32 +0000 Subject: [PATCH 25/30] Fix last remaining vertex integral bugs --- cpp/dolfinx/fem/assemble_scalar_impl.h | 2 +- cpp/dolfinx/fem/utils.cpp | 3 ++- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/cpp/dolfinx/fem/assemble_scalar_impl.h b/cpp/dolfinx/fem/assemble_scalar_impl.h index 04b9df59b11..ce181096ce8 100644 --- a/cpp/dolfinx/fem/assemble_scalar_impl.h +++ b/cpp/dolfinx/fem/assemble_scalar_impl.h @@ -287,7 +287,7 @@ T assemble_scalar( auto& [coeffs, cstride] = coefficients.at({vertex_type, i}); std::span vertices = M.domain(vertex_type, i, 0); - assert(vertices.size() * cstride == coeffs.size()); + assert(vertices.size() / 2 * cstride == coeffs.size()); constexpr std::size_t num_adjacent_cells = 1; // Two values per adj. cell (cell index and local vertex index). diff --git a/cpp/dolfinx/fem/utils.cpp b/cpp/dolfinx/fem/utils.cpp index d784ac583b0..b252bc87d16 100644 --- a/cpp/dolfinx/fem/utils.cpp +++ b/cpp/dolfinx/fem/utils.cpp @@ -139,7 +139,7 @@ fem::compute_integration_domains(const fem::IntegralType& integral_type, { const int tdim = topology.dim(); - int dim = tdim - integral_type.codim; + int dim = integral_type.codim == -1 ? 0 : tdim - integral_type.codim; assert(dim >= 0); { @@ -228,6 +228,7 @@ fem::compute_integration_domains(const fem::IntegralType& integral_type, } else if (integral_type.codim == -1 and integral_type.num_cells == 1) { + std::cout << "HERE!" << std::endl; auto [v_to_c, c_to_v] = get_connectivities(0); for (auto vertex : entities) { From 5a70ea1cb04481c5dcf47cbeb21725124b6c8d28 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= Date: Thu, 11 Sep 2025 07:46:51 +0000 Subject: [PATCH 26/30] Remove debug print --- cpp/dolfinx/fem/utils.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/cpp/dolfinx/fem/utils.cpp b/cpp/dolfinx/fem/utils.cpp index b252bc87d16..34d8bb0d5d4 100644 --- a/cpp/dolfinx/fem/utils.cpp +++ b/cpp/dolfinx/fem/utils.cpp @@ -228,7 +228,6 @@ fem::compute_integration_domains(const fem::IntegralType& integral_type, } else if (integral_type.codim == -1 and integral_type.num_cells == 1) { - std::cout << "HERE!" << std::endl; auto [v_to_c, c_to_v] = get_connectivities(0); for (auto vertex : entities) { From 3b081e85ba2cda27033bc54cda1d92a128cbc835 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= Date: Thu, 11 Sep 2025 07:47:53 +0000 Subject: [PATCH 27/30] Remove unused variable --- python/dolfinx/fem/forms.py | 1 - 1 file changed, 1 deletion(-) diff --git a/python/dolfinx/fem/forms.py b/python/dolfinx/fem/forms.py index 48da35fe7dd..71ec5217fc9 100644 --- a/python/dolfinx/fem/forms.py +++ b/python/dolfinx/fem/forms.py @@ -370,7 +370,6 @@ def _form(form): ] constants = [c._cpp_object for c in form.constants()] # Extract subdomain ids from ufcx_form - tdim = domain.topological_dimension() ufl_to_itg = _ufl_to_dolfinx_domain() subdomain_ids = {ufl_to_itg[type]: [] for type in sd.get(domain).keys()} integral_offsets = [ufcx_form.form_integral_offsets[i] for i in range(5)] From 7c67f85c9a5868dc6d796d75a46be08e9ec91d95 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= Date: Thu, 11 Sep 2025 09:07:03 +0000 Subject: [PATCH 28/30] Fix last vertex integral --- 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 624e1100183..29ad253b577 100644 --- a/cpp/dolfinx/fem/assemble_vector_impl.h +++ b/cpp/dolfinx/fem/assemble_vector_impl.h @@ -1474,7 +1474,7 @@ void assemble_vector( } } - IntegralType vertex_integral_type = IntegralType(mesh->topology()->dim()); + IntegralType vertex_integral_type = IntegralType(-1); for (int i = 0; i < L.num_integrals(vertex_integral_type, 0); ++i) { auto fn = L.kernel(vertex_integral_type, i, 0); From 397f4b75a7b623493c6ca52f35d08a186045dda1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= Date: Thu, 11 Sep 2025 09:23:43 +0000 Subject: [PATCH 29/30] Fix another vertex bug --- 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 29ad253b577..59da99de3f6 100644 --- a/cpp/dolfinx/fem/assemble_vector_impl.h +++ b/cpp/dolfinx/fem/assemble_vector_impl.h @@ -1486,7 +1486,7 @@ void assemble_vector( auto& [coeffs, cstride] = coefficients.at({vertex_integral_type, i}); - assert(vertices.size() * cstride == coeffs.size()); + assert(vertices.size() / 2 * cstride == coeffs.size()); impl::assemble_vertices( P0, b, x_dofmap, x, From 5bc330513651fdf24bc0f51705c70ec367545b14 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= Date: Thu, 11 Sep 2025 09:51:27 +0000 Subject: [PATCH 30/30] Remove unused variable --- cpp/dolfinx/fem/assemble_scalar_impl.h | 1 - 1 file changed, 1 deletion(-) diff --git a/cpp/dolfinx/fem/assemble_scalar_impl.h b/cpp/dolfinx/fem/assemble_scalar_impl.h index ce181096ce8..2df70cf932d 100644 --- a/cpp/dolfinx/fem/assemble_scalar_impl.h +++ b/cpp/dolfinx/fem/assemble_scalar_impl.h @@ -277,7 +277,6 @@ T assemble_scalar( perms, cdofs_b); } - std::int32_t tdim = mesh->topology()->dim(); IntegralType vertex_type = IntegralType(-1); for (int i = 0; i < M.num_integrals(vertex_type, 0); ++i) {