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/.github/workflows/windows.yml b/.github/workflows/windows.yml index f439017f100..0214ded3d17 100644 --- a/.github/workflows/windows.yml +++ b/.github/workflows/windows.yml @@ -50,6 +50,7 @@ jobs: - name: Load environment variables run: | cat dolfinx-src/.github/workflows/fenicsx-refs.env >> $env:GITHUB_ENV + - name: Set up Python uses: actions/setup-python@v6 with: @@ -95,6 +96,7 @@ 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: Checkout FFCx uses: actions/checkout@v5 with: diff --git a/cpp/demo/codim_0_assembly/main.cpp b/cpp/demo/codim_0_assembly/main.cpp index 5be7086e5e2..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::cell, *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::cell, {{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 6ddcfacabfb..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::cell, 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::cell, 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 65eeb8fa395..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::exterior_facet, *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::exterior_facet, {{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 8a361dcf0c8..3814d177020 100644 --- a/cpp/dolfinx/fem/Form.h +++ b/cpp/dolfinx/fem/Form.h @@ -35,12 +35,44 @@ template class Function; /// @brief Type of integral -enum class IntegralType : std::int8_t + +struct IntegralType { - cell = 0, ///< Cell - exterior_facet = 1, ///< Exterior 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 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) + { + } + + /// @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; + } + + /// @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 @@ -130,9 +162,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. @@ -274,10 +306,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::exterior_facet - or type == IntegralType::interior_facet) + else if (type.codim == 1) { const mesh::Topology topology = *_mesh->topology(); int tdim = topology.dim(); @@ -317,10 +348,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::exterior_facet - or type == IntegralType::interior_facet) + else if (type.codim == 1) { const mesh::Topology topology = *_mesh->topology(); int tdim = topology.dim(); @@ -454,16 +484,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::exterior_facet, returns a list with shape - /// `(num_facets, 2)`, where `[cell_index, 0]` is the cell index and - /// `[cell_index, 1]` is the local facet index relative to the cell. - /// - For IntegralType::interior_facet the shape is `(num_facets, 4)`, - /// 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. @@ -512,8 +538,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 diff --git a/cpp/dolfinx/fem/assemble_matrix_impl.h b/cpp/dolfinx/fem/assemble_matrix_impl.h index 90351a05b7a..31c1e37f490 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); + 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::exterior_facet, cell_type_idx); ++i) + IntegralType facet = IntegralType(1); + for (int i = 0; i < a.num_integrals(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(facet, i, 0); assert(fn); - auto& [coeffs, cstride] - = coefficients.at({IntegralType::exterior_facet, i}); + auto& [coeffs, cstride] = coefficients.at({facet, i}); - std::span f = a.domain(IntegralType::exterior_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::exterior_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::exterior_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( @@ -650,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) { @@ -666,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 1ade0a9792e..2df70cf932d 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); + 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,14 +227,14 @@ T assemble_scalar( num_facets_per_cell); } - for (int i = 0; i < M.num_integrals(IntegralType::exterior_facet, 0); ++i) + IntegralType facet_type = IntegralType(1); + for (int i = 0; i < M.num_integrals(facet_type, 0); ++i) { - auto fn = M.kernel(IntegralType::exterior_facet, i, 0); + auto fn = M.kernel(facet_type, i, 0); assert(fn); - auto& [coeffs, cstride] - = coefficients.at({IntegralType::exterior_facet, i}); + auto& [coeffs, cstride] = coefficients.at({facet_type, i}); - std::span facets = M.domain(IntegralType::exterior_facet, i, 0); + std::span facets = M.domain(facet_type, i, 0); // Two values per each adjacent cell (cell index and local facet // index) @@ -251,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). @@ -276,16 +277,16 @@ T assemble_scalar( perms, cdofs_b); } - for (int i = 0; i < M.num_integrals(IntegralType::vertex, 0); ++i) + IntegralType vertex_type = IntegralType(-1); + 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); - assert(vertices.size() * cstride == coeffs.size()); + std::span vertices = M.domain(vertex_type, i, 0); + 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/assemble_vector_impl.h b/cpp/dolfinx/fem/assemble_vector_impl.h index 3e8c1c2a13d..59da99de3f6 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); + 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,21 +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::exterior_facet, 0); ++i) + IntegralType facet_integral_type = IntegralType(1); + for (int i = 0; i < a.num_integrals(facet_integral_type, 0); ++i) { - auto kernel = a.kernel(IntegralType::exterior_facet, i, 0); + auto kernel = a.kernel(facet_integral_type, i, 0); assert(kernel); - auto& [coeffs, cstride] - = coefficients.at({IntegralType::exterior_facet, i}); + auto& [coeffs, cstride] = coefficients.at({facet_integral_type, i}); using mdspanx2_t = md::mdspan>; - std::span f = a.domain(IntegralType::exterior_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::exterior_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::exterior_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( @@ -1169,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, @@ -1340,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); + 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) { @@ -1385,15 +1388,15 @@ void assemble_vector( = md::mdspan>; - for (int i = 0; i < L.num_integrals(IntegralType::exterior_facet, 0); ++i) + IntegralType facet_integral_type = IntegralType(1); + for (int i = 0; i < L.num_integrals(facet_integral_type, 0); ++i) { - auto fn = L.kernel(IntegralType::exterior_facet, i, 0); + auto fn = L.kernel(facet_integral_type, 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({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::exterior_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) @@ -1419,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) { @@ -1470,18 +1474,19 @@ void assemble_vector( } } - for (int i = 0; i < L.num_integrals(IntegralType::vertex, 0); ++i) + IntegralType vertex_integral_type = IntegralType(-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()); + assert(vertices.size() / 2 * cstride == coeffs.size()); impl::assemble_vertices( P0, b, x_dofmap, x, diff --git a/cpp/dolfinx/fem/pack.h b/cpp/dolfinx/fem/pack.h index 9b08e845d81..8cde37e387c 100644 --- a/cpp/dolfinx/fem/pack.h +++ b/cpp/dolfinx/fem/pack.h @@ -179,11 +179,8 @@ allocate_coefficient_storage(const Form& form, IntegralType integral_type, const std::vector offsets = form.coefficient_offsets(); cstride = offsets.back(); num_entities = form.domain(integral_type, id, 0).size(); - if (integral_type == IntegralType::exterior_facet - or integral_type == IntegralType::interior_facet) - { + if (integral_type.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(integral_type, 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(integral_type, 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,20 @@ void pack_coefficients(const Form& form, *coefficients[coeff], cell_info, cells, offsets[coeff]); } - break; } - case IntegralType::exterior_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::exterior_facet, id)) + // facet or vertex integrals. + for (int coeff : form.active_coeffs(integral_type, id)) { auto mesh = coefficients[coeff]->function_space()->mesh(); - std::span facets_b - = form.domain_coeff(IntegralType::exterior_facet, id, coeff); + std::span entities_b + = form.domain_coeff(integral_type, 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 +284,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(integral_type, id)) { auto mesh = coefficients[coeff]->function_space()->mesh(); std::span facets_b - = form.domain_coeff(IntegralType::interior_facet, id, coeff); + = form.domain_coeff(integral_type, id, coeff); md::mdspan> facets(facets_b.data(), facets_b.size() / 4, 4); @@ -320,31 +313,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 ec7a8501ae3..34d8bb0d5d4 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::exterior_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 = integral_type.codim == -1 ? 0 : 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::exterior_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 == -1 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("Cannot compute integration domains 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 22436c8e50d..e4cc8f331c9 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::exterior_facet`: facets -/// `IntegralType::interior_facet`: facets -/// `IntegralType::vertex`: vertices -/// @return List of integration entity data, depending on the `IntegralType` the -/// data per entity has different layouts -/// `IntegralType::cell`: cell -/// `IntegralType::exterior_facet`: (cell, local_facet) -/// `IntegralType::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::exterior_facet) != types.end()) + const std::set> types + = a.integral_types(); + + 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. 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)) { - 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::exterior_facet: + } + else if (type == IntegralType(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"); } } @@ -502,8 +496,7 @@ Form create_form_factory( } // Create facets, if required - if (num_integrals_type[exterior_facet] > 0 - or num_integrals_type[interior_facet] > 0) + if (num_integrals_type[facet] > 0 or num_integrals_type[interior_facet] > 0) { mesh->topology_mutable()->create_entities(tdim - 1); mesh->topology_mutable()->create_connectivity(tdim - 1, tdim); @@ -534,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)); for (std::size_t form_idx = 0; form_idx < ufcx_forms.size(); ++form_idx) { const ufcx_form& ufcx_form = ufcx_forms[form_idx]; @@ -569,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), i, form_idx}, {k, default_cells, active_coeffs}}); } else if (sd != subdomains.end()) @@ -579,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), i, form_idx}, {k, std::vector(it->second.begin(), it->second.end()), @@ -597,17 +590,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]); - auto sd = subdomains.find(IntegralType::exterior_facet); + + integral_offsets[facet], + num_integrals_type[facet]); + 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]; - 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); @@ -638,7 +631,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(1), i, form_idx}, {k, default_facets_ext, active_coeffs}}); } else if (sd != subdomains.end()) @@ -648,7 +641,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(1), i, form_idx}, {k, std::vector(it->second.begin(), it->second.end()), @@ -668,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, 2)); for (std::size_t form_idx = 0; form_idx < ufcx_forms.size(); ++form_idx) { const ufcx_form& ufcx_form = ufcx_forms[form_idx]; @@ -734,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()) @@ -743,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()), @@ -766,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(-1)); for (int i = 0; i < num_integrals_type[vertex]; ++i) { const int id = ids[i]; @@ -817,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(-1), i, form_idx}, {k, cells_and_vertices, active_coeffs}}); } else if (sd != subdomains.end()) @@ -827,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(-1), 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 5be91e0eec6..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.cell: [(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 2674ca96590..f91870c7378 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.exterior_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 == -1: subdomain._cpp_object.topology.create_connectivity(0, tdim) subdomain._cpp_object.topology.create_connectivity(tdim, 0) @@ -215,12 +213,16 @@ def form_cpp_class( raise NotImplementedError(f"Type {dtype} not supported.") -_ufl_to_dolfinx_domain = { - "cell": IntegralType.cell, - "exterior_facet": IntegralType.exterior_facet, - "interior_facet": IntegralType.interior_facet, - "vertex": IntegralType.vertex, -} +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. + """ + return { + "cell": IntegralType(0), + "exterior_facet": IntegralType(1), + "interior_facet": IntegralType(1, 2), + "vertex": IntegralType(-1), + } def mixed_topology_form( @@ -367,19 +369,24 @@ 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 = {type: [] for type in sd.get(domain).keys()} + 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 = [ + ufl_to_itg[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.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_itg[key]: get_integration_domains( + ufl_to_itg[key], + subdomain_data[0], + 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 1acc71e53a2..0e5bc43cf08 100644 --- a/python/dolfinx/wrappers/fem.cpp +++ b/python/dolfinx/wrappers/fem.cpp @@ -767,22 +767,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::exterior_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,13 +1269,12 @@ void fem(nb::module_& m) }, nb::rv_policy::reference_internal); - nb::enum_(m, "_IntegralType") - .value("cell", dolfinx::fem::IntegralType::cell, "cell integral") - .value("exterior_facet", dolfinx::fem::IntegralType::exterior_facet, - "exterior 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") = 1) + .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"); 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..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.exterior_facet, mesh.topology, facets + dolfinx.fem.IntegralType(1), mesh.topology, facets ) - subdomains = {dolfinx.fem.IntegralType.exterior_facet: [(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_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 a9165008d82..fe511bc2b97 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(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.exterior_facet, msh.topology, vertices) - subdomains = {fem.IntegralType.exterior_facet: [(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 0f0a0acdd2c..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.cell: [ + 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.cell: [(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.cell: [(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.cell: [(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.cell: [(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 f6b2865984e..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.cell: []}, + {IntegralType(0): []}, [], mesh._cpp_object, ) @@ -105,7 +105,7 @@ def test_incorrect_element(): [incorrect_space._cpp_object, incorrect_space._cpp_object], [], [], - {IntegralType.cell: []}, + {IntegralType(0): []}, [], mesh._cpp_object, )