diff --git a/.github/workflows/ccpp.yml b/.github/workflows/ccpp.yml index 07cfcd5e47..087519a8de 100644 --- a/.github/workflows/ccpp.yml +++ b/.github/workflows/ccpp.yml @@ -97,9 +97,9 @@ jobs: - name: Install FEniCS Python components run: | - pip install --no-build-isolation git+https://github.com/fenics/ufl.git@${{ env.ufl_ref }} - pip install --no-build-isolation git+https://github.com/fenics/basix.git@${{ env.basix_ref }} - pip install --no-build-isolation git+https://github.com/fenics/ffcx.git@${{ env.ffcx_ref }} + pip install --no-build-isolation git+https://github.com/${{ env.ufl_repository }}.git@${{ env.ufl_ref }} + pip install --no-build-isolation git+https://github.com/${{ env.basix_repository }}.git@${{ env.basix_ref }} + pip install --no-build-isolation git+https://github.com/${{ env.ffcx_repository }}.git@${{ env.ffcx_ref }} - name: Configure (C++) working-directory: cpp @@ -189,9 +189,9 @@ jobs: - name: Install FEniCS Python components run: | - pip install --no-build-isolation git+https://github.com/FEniCS/ufl.git@${{ env.ufl_ref }} - pip install --no-build-isolation git+https://github.com/FEniCS/basix.git@${{ env.basix_ref }} - pip install --no-build-isolation git+https://github.com/FEniCS/ffcx.git@${{ env.ffcx_ref }} + pip install --no-build-isolation git+https://github.com/${{ env.ufl_repository }}.git@${{ env.ufl_ref }} + pip install --no-build-isolation git+https://github.com/${{ env.basix_repository }}.git@${{ env.basix_ref }} + pip install --no-build-isolation git+https://github.com/${{ env.ffcx_repository }}.git@${{ env.ffcx_ref }} - name: Configure (C++) working-directory: cpp @@ -293,9 +293,9 @@ jobs: - name: Install FEniCS Python components run: | pip install -r python/build-requirements.txt - pip install --no-build-isolation git+https://github.com/FEniCS/ufl.git@${{ env.ufl_ref }} - pip install --no-build-isolation git+https://github.com/FEniCS/basix.git@${{ env.basix_ref }} - pip install --no-build-isolation git+https://github.com/FEniCS/ffcx.git@${{ env.ffcx_ref }} + pip install --no-build-isolation git+https://github.com/${{ env.ufl_repository }}.git@${{ env.ufl_ref }} + pip install --no-build-isolation git+https://github.com/${{ env.basix_repository }}.git@${{ env.basix_ref }} + pip install --no-build-isolation git+https://github.com/${{ env.ffcx_repository }}.git@${{ env.ffcx_ref }} - name: Configure C++ run: | diff --git a/.github/workflows/ci-spack.yml b/.github/workflows/ci-spack.yml index fe549cb7b1..178f63ce11 100644 --- a/.github/workflows/ci-spack.yml +++ b/.github/workflows/ci-spack.yml @@ -86,9 +86,9 @@ jobs: run: | . ./spack-src/share/spack/setup-env.sh spack env activate . - pip install git+https://github.com/fenics/ufl.git@${{ env.ufl_ref }} - pip install git+https://github.com/fenics/basix.git@${{ env.basix_ref }} - pip install git+https://github.com/fenics/ffcx.git@${{ env.ffcx_ref }} + pip install git+https://github.com/${{ env.ufl_repository }}.git@${{ env.ufl_ref }} + pip install git+https://github.com/${{ env.basix_repository }}.git@${{ env.basix_ref }} + pip install git+https://github.com/${{ env.ffcx_repository }}.git@${{ env.ffcx_ref }} - name: Configure and build C++ run: | diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index bb2e85119a..d4dc16dee9 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -5,8 +5,8 @@ on: push: branches: - "**" - # pull_request: - # branches: [ "main" ] + pull_request: + branches: [ "main" ] # merge_group: # branches: # - main @@ -39,6 +39,28 @@ jobs: uses: actions/checkout@v5 with: path: dolfinx-src + + - name: Load dev branch environment variables + run: cat dolfinx-src/.github/workflows/fenicsx-refs.env >> $GITHUB_ENV + + - name: Get UFL code + uses: actions/checkout@v5 + with: + repository: ${{ env.ufl_repository }} + ref: ${{ env.ufl_ref }} + path: ufl-src + - name: Get FFCx code + uses: actions/checkout@v5 + with: + repository: ${{ env.ffcx_repository }} + ref: ${{ env.ffcx_ref }} + path: ffcx-src + - name: Get BASIx code + uses: actions/checkout@v5 + with: + repository: ${{ env.basix_repository }} + ref: ${{ env.basix_ref }} + path: basix-src - name: Get Spack uses: actions/checkout@v5 with: @@ -52,9 +74,6 @@ jobs: ref: fenics-testing path: spack-packages - - name: Load dev branch environment variables - run: cat dolfinx-src/.github/workflows/fenicsx-refs.env >> $GITHUB_ENV - - name: Add repo run: | . $GITHUB_WORKSPACE/spack-src/share/spack/setup-env.sh @@ -65,12 +84,12 @@ jobs: . $GITHUB_WORKSPACE/spack-src/share/spack/setup-env.sh spack info fenics-dolfinx spack env create cxx dolfinx-src/.github/workflows/spack-config/gh-actions-env-test.yml - spack -e cxx develop --path $GITHUB_WORKSPACE/dolfinx-src fenics-dolfinx@main - spack -e cxx develop fenics-basix@git.${{ env.basix_ref }}=main - spack -e cxx develop fenics-ufcx@git.${{ env.ffcx_ref }}=main - spack -e cxx develop py-fenics-basix@git.${{ env.basix_ref }}=main - spack -e cxx develop py-fenics-ffcx@git.${{ env.ffcx_ref }}=main - spack -e cxx develop py-fenics-ufl@git.${{ env.ufl_ref }}=main + spack -e cxx develop --no-clone --path $GITHUB_WORKSPACE/dolfinx-src fenics-dolfinx@main + spack -e cxx develop --no-clone --path $GITHUB_WORKSPACE/ffcx-src fenics-ufcx@main + spack -e cxx develop --no-clone --path $GITHUB_WORKSPACE/basix-src fenics-basix@main + spack -e cxx develop --no-clone --path $GITHUB_WORKSPACE/basix-src py-fenics-basix@main + spack -e cxx develop --no-clone --path $GITHUB_WORKSPACE/ffcx-src py-fenics-ffcx@main + spack -e cxx develop --no-clone --path $GITHUB_WORKSPACE/ufl-src py-fenics-ufl@main - name: Build library (C++) run: | @@ -115,13 +134,13 @@ jobs: . $GITHUB_WORKSPACE/spack-src/share/spack/setup-env.sh spack info py-fenics-dolfinx spack env create py dolfinx-src/.github/workflows/spack-config/gh-actions-env-test.yml - spack -e py develop --path $GITHUB_WORKSPACE/dolfinx-src fenics-dolfinx@main - spack -e py develop fenics-basix@git.${{ env.basix_ref }}=main - spack -e py develop fenics-ufcx@git.${{ env.ffcx_ref }}=main - spack -e py develop --path $GITHUB_WORKSPACE/dolfinx-src py-fenics-dolfinx@main - spack -e py develop py-fenics-basix@git.${{ env.basix_ref }}=main - spack -e py develop py-fenics-ffcx@git.${{ env.ffcx_ref }}=main - spack -e py develop py-fenics-ufl@git.${{ env.ufl_ref }}=main + spack -e py develop --no-clone --path $GITHUB_WORKSPACE/dolfinx-src fenics-dolfinx@main + spack -e py develop --no-clone --path $GITHUB_WORKSPACE/dolfinx-src py-fenics-dolfinx@main + spack -e py develop --no-clone --path $GITHUB_WORKSPACE/ffcx-src fenics-ufcx@main + spack -e py develop --no-clone --path $GITHUB_WORKSPACE/basix-src fenics-basix@main + spack -e py develop --no-clone --path $GITHUB_WORKSPACE/basix-src py-fenics-basix@main + spack -e py develop --no-clone --path $GITHUB_WORKSPACE/ffcx-src py-fenics-ffcx@main + spack -e py develop --no-clone --path $GITHUB_WORKSPACE/ufl-src py-fenics-ufl@main - name: Build library (Python) run: | diff --git a/.github/workflows/fenicsx-refs.env b/.github/workflows/fenicsx-refs.env index 8e958e229c..22c7e44a7e 100644 --- a/.github/workflows/fenicsx-refs.env +++ b/.github/workflows/fenicsx-refs.env @@ -1,4 +1,6 @@ +basix_repository=FEniCS/basix basix_ref=main +ufl_repository=FEniCS/ufl ufl_ref=main +ffcx_repository=FEniCS/ffcx ffcx_ref=main - diff --git a/.github/workflows/macos.yml b/.github/workflows/macos.yml index cd8b7f155e..f073088068 100644 --- a/.github/workflows/macos.yml +++ b/.github/workflows/macos.yml @@ -84,9 +84,9 @@ jobs: - name: Install FEniCSx dependencies run: | - pip install git+https://github.com/fenics/ufl.git@${{ env.ufl_ref }} - pip install git+https://github.com/fenics/basix.git@${{ env.basix_ref }} - pip install git+https://github.com/fenics/ffcx.git@${{ env.ffcx_ref }} + pip install git+https://github.com/${{ env.ufl_repository }}.git@${{ env.ufl_ref }} + pip install git+https://github.com/${{ env.basix_repository }}.git@${{ env.basix_ref }} + pip install git+https://github.com/${{ env.ffcx_repository }}.git@${{ env.ffcx_ref }} - name: Configure (C++) working-directory: cpp diff --git a/.github/workflows/sonarcloud.yml b/.github/workflows/sonarcloud.yml index 99cb683543..eb06e5a37f 100644 --- a/.github/workflows/sonarcloud.yml +++ b/.github/workflows/sonarcloud.yml @@ -62,9 +62,9 @@ jobs: echo "$HOME/.sonar/build-wrapper-linux-x86" >> $GITHUB_PATH - name: Install FEniCS Python components run: | - python -m pip install git+https://github.com/fenics/ufl.git@${{ env.ufl_ref }} - python -m pip install git+https://github.com/fenics/basix.git@${{ env.basix_ref }} - python -m pip install git+https://github.com/fenics/ffcx.git@${{ env.ffcx_ref }} + pip install git+https://github.com/${{ env.ufl_repository }}.git@${{ env.ufl_ref }} + pip install git+https://github.com/${{ env.basix_repository }}.git@${{ env.basix_ref }} + pip install git+https://github.com/${{ env.ffcx_repository }}.git@${{ env.ffcx_ref }} - name: Run build-wrapper run: | mkdir build diff --git a/cpp/dolfinx/fem/CMakeLists.txt b/cpp/dolfinx/fem/CMakeLists.txt index 8016d43187..548df96a01 100644 --- a/cpp/dolfinx/fem/CMakeLists.txt +++ b/cpp/dolfinx/fem/CMakeLists.txt @@ -18,6 +18,7 @@ set(HEADERS_fem ${CMAKE_CURRENT_SOURCE_DIR}/dofmapbuilder.h ${CMAKE_CURRENT_SOURCE_DIR}/dolfinx_fem.h ${CMAKE_CURRENT_SOURCE_DIR}/interpolate.h + ${CMAKE_CURRENT_SOURCE_DIR}/kernel.h ${CMAKE_CURRENT_SOURCE_DIR}/petsc.h ${CMAKE_CURRENT_SOURCE_DIR}/pack.h ${CMAKE_CURRENT_SOURCE_DIR}/sparsitybuild.h diff --git a/cpp/dolfinx/fem/assemble_scalar_impl.h b/cpp/dolfinx/fem/assemble_scalar_impl.h index 85c1e2cd6e..3f060e7853 100644 --- a/cpp/dolfinx/fem/assemble_scalar_impl.h +++ b/cpp/dolfinx/fem/assemble_scalar_impl.h @@ -1,4 +1,4 @@ -// Copyright (C) 2019-2025 Garth N. Wells +// Copyright (C) 2019-2025 Garth N. Wells and Paul T. Kühner // // This file is part of DOLFINx (https://www.fenicsproject.org) // @@ -150,6 +150,43 @@ T assemble_interior_facets( return value; } +/// Assemble functional over vertices +template +T assemble_vertices(mdspan2_t x_dofmap, + md::mdspan, + md::extents> + x, + md::mdspan> + vertices, + FEkernel auto fn, std::span constants, + md::mdspan> coeffs) +{ + T value(0); + if (vertices.empty()) + return value; + + // Create data structures used in assembly + std::vector> cdofs(3 * x_dofmap.extent(1)); + + // Iterate over all cells + for (std::size_t index = 0; index < vertices.extent(0); ++index) + { + std::int32_t cell = vertices(index, 0); + std::int32_t local_vertex_index = vertices(index, 1); + + // Get cell coordinates/geometry + auto x_dofs = md::submdspan(x_dofmap, cell, md::full_extent); + for (std::size_t i = 0; i < x_dofs.size(); ++i) + std::copy_n(&x(x_dofs[i], 0), 3, std::next(cdofs.begin(), 3 * i)); + + fn(&value, &coeffs(index, 0), constants.data(), cdofs.data(), + &local_vertex_index, nullptr, nullptr); + } + + return value; +} + /// Assemble functional into an scalar with provided mesh geometry. template T assemble_scalar( @@ -198,14 +235,19 @@ T assemble_scalar( = coefficients.at({IntegralType::exterior_facet, i}); std::span facets = M.domain(IntegralType::exterior_facet, i, 0); + + constexpr std::size_t num_adjacent_cells = 1; + // Two values per each adj. cell (cell index and local facet index). + constexpr std::size_t shape1 = 2 * num_adjacent_cells; + assert((facets.size() / 2) * cstride == coeffs.size()); value += impl::assemble_exterior_facets( x_dofmap, x, md::mdspan>( - facets.data(), facets.size() / 2, 2), - fn, constants, md::mdspan(coeffs.data(), facets.size() / 2, cstride), - perms); + facets.data(), facets.size() / shape1, 2), + fn, constants, + md::mdspan(coeffs.data(), facets.size() / shape1, cstride), perms); } for (int i = 0; i < M.num_integrals(IntegralType::interior_facet, 0); ++i) @@ -215,19 +257,48 @@ T assemble_scalar( auto& [coeffs, cstride] = coefficients.at({IntegralType::interior_facet, i}); std::span facets = M.domain(IntegralType::interior_facet, i, 0); - assert((facets.size() / 4) * 2 * cstride == coeffs.size()); + + constexpr std::size_t num_adjacent_cells = 2; + // Two values per each adj. cell (cell index and local facet index). + constexpr std::size_t shape1 = 2 * num_adjacent_cells; + + assert((facets.size() / shape1) * 2 * cstride == coeffs.size()); value += impl::assemble_interior_facets( x_dofmap, x, md::mdspan>( - facets.data(), facets.size() / 4, 2, 2), + facets.data(), facets.size() / shape1, 2, 2), fn, constants, md::mdspan>( - coeffs.data(), facets.size() / 4, 2, cstride), + coeffs.data(), facets.size() / shape1, 2, cstride), perms); } + for (int i = 0; i < M.num_integrals(IntegralType::vertex, 0); ++i) + { + auto fn = M.kernel(IntegralType::vertex, i, 0); + assert(fn); + + auto& [coeffs, cstride] = coefficients.at({IntegralType::vertex, i}); + + std::span vertices + = M.domain(IntegralType::vertex, i, 0); + assert(vertices.size() * cstride == coeffs.size()); + + constexpr std::size_t num_adjacent_cells = 1; + // Two values per adj. cell (cell index and local vertex index). + constexpr std::size_t shape1 = 2 * num_adjacent_cells; + + value += impl::assemble_vertices( + x_dofmap, x, + md::mdspan>( + vertices.data(), vertices.size() / shape1, shape1), + fn, constants, + md::mdspan(coeffs.data(), vertices.size() / shape1, cstride)); + } + return value; } diff --git a/cpp/dolfinx/fem/assemble_vector_impl.h b/cpp/dolfinx/fem/assemble_vector_impl.h index 855aedccd4..3e8c1c2a13 100644 --- a/cpp/dolfinx/fem/assemble_vector_impl.h +++ b/cpp/dolfinx/fem/assemble_vector_impl.h @@ -1,4 +1,4 @@ -// Copyright (C) 2018-2021 Garth N. Wells +// Copyright (C) 2018-2025 Garth N. Wells and Paul T. Kühner // // This file is part of DOLFINx (https://www.fenicsproject.org) // @@ -941,6 +941,94 @@ void assemble_interior_facets( } } +/// @brief Execute kernel over a set of vertices and accumulate result in +/// vector. +/// +/// @tparam T Scalar type +/// @tparam _bs Block size of the form test function dof map. If less +/// than zero the block size is determined at runtime. If `_bs` is +/// positive the block size is used as a compile-time constant, which +/// has performance benefits. +/// @param[in] P0 Function that applies transformation `P0.b` in-place +/// to `b` to transform test degrees-of-freedom. +/// @param[in,out] b Array to accumulate into. +/// @param[in] x_dofmap Dofmap for the mesh geometry. +/// @param[in] x Mesh geometry (coordinates). +/// @param[in] vertices Vertex indices `(vertices.size(), 2)` - first entry +/// holds the index of the cell adjacent to the vertex, and the second +/// stores the local index of the vertex within the cell. +/// @param[in] dofmap Test function (row) degree-of-freedom data holding +/// the (0) dofmap, (1) dofmap block size and (2) dofmap cell indices. +/// @param[in] kernel Kernel function to execute over each cell. +/// @param[in] constants Constant coefficient data in the kernel. +/// @param[in] coeffs Coefficient data in the kernel. It has shape +/// `(vertices.size(), num_cell_coeffs)`. `coeffs(i, j)` is the `j`th +/// coefficient for cell `i`. +/// @param[in] cell_info0 Cell permutation information for the test +/// function mesh. +template +void assemble_vertices( + fem::DofTransformKernel auto P0, std::span b, mdspan2_t x_dofmap, + md::mdspan, + md::extents> + x, + md::mdspan> + vertices, + std::tuple>> + dofmap, + FEkernel auto kernel, std::span constants, + md::mdspan> coeffs, + std::span cell_info0) +{ + if (vertices.empty()) + return; + + const auto [dmap, bs, vertices0] = dofmap; + assert(_bs < 0 or _bs == bs); + + // Create data structures used in assembly + std::vector> cdofs(3 * x_dofmap.extent(1)); + std::vector be(bs * dmap.extent(1)); + + // Iterate over active vertices + for (std::size_t index = 0; index < vertices.extent(0); ++index) + { + // Integration domain cell, local index, and test function cell + std::int32_t cell = vertices(index, 0); + std::int32_t local_index = vertices(index, 1); + std::int32_t c0 = vertices0(index, 0); + + // Get cell coordinates/geometry + auto x_dofs = md::submdspan(x_dofmap, cell, md::full_extent); + for (std::size_t i = 0; i < x_dofs.size(); ++i) + std::copy_n(&x(x_dofs[i], 0), 3, std::next(cdofs.begin(), 3 * i)); + + // Tabulate vector for vertex + std::ranges::fill(be, 0); + kernel(be.data(), &coeffs(index, 0), constants.data(), cdofs.data(), + &local_index, nullptr, nullptr); + P0(be, cell_info0, c0, 1); + + // Scatter vertex vector to 'global' vector array + auto dofs = md::submdspan(dmap, c0, md::full_extent); + if constexpr (_bs > 0) + { + for (std::size_t i = 0; i < dofs.size(); ++i) + for (int k = 0; k < _bs; ++k) + b[_bs * dofs[i] + k] += be[_bs * i + k]; + } + else + { + for (std::size_t i = 0; i < dofs.size(); ++i) + for (int k = 0; k < bs; ++k) + b[bs * dofs[i] + k] += be[bs * i + k]; + } + } +} + /// Modify RHS vector to account for boundary condition such that: /// /// b <- b - alpha * A.(x_bc - x0) @@ -1381,6 +1469,32 @@ void assemble_vector( cell_info0, perms); } } + + for (int i = 0; i < L.num_integrals(IntegralType::vertex, 0); ++i) + { + auto fn = L.kernel(IntegralType::vertex, i, 0); + assert(fn); + + std::span vertices = L.domain(IntegralType::vertex, i, cell_type_idx); + std::span vertices0 + = L.domain_arg(IntegralType::vertex, 0, i, cell_type_idx); + + auto& [coeffs, cstride] = coefficients.at({IntegralType::vertex, i}); + + assert(vertices.size() * cstride == coeffs.size()); + + impl::assemble_vertices( + P0, b, x_dofmap, x, + md::mdspan>( + vertices.data(), vertices.size() / 2, 2), + {dofs, bs, + md::mdspan>( + vertices0.data(), vertices0.size() / 2, 2)}, + fn, constants, + md::mdspan(coeffs.data(), vertices.size() / 2, cstride), cell_info0); + } } } diff --git a/cpp/dolfinx/fem/kernel.h b/cpp/dolfinx/fem/kernel.h new file mode 100644 index 0000000000..653611cb2f --- /dev/null +++ b/cpp/dolfinx/fem/kernel.h @@ -0,0 +1,57 @@ +// Copyright (C) 2025 Paul T. Kühner +// +// This file is part of DOLFINx (https://www.fenicsproject.org) +// +// SPDX-License-Identifier: LGPL-3.0-or-later + +#pragma once + +#include +#include +#include + +namespace dolfinx::fem::impl +{ +/// @brief Kernel C-pointer type. +/// @tparam T scalar type. +/// @tparam U geometry type. +template > +using kernelptr_t = void (*)(T*, const T*, const T*, const U*, const int*, + const std::uint8_t*, void*); + +/// @brief Kernel callback type. +/// @tparam T scalar type. +/// @tparam U geometry type. +template > +using kernel_t = std::function; + +/// @brief Extract correct kernel by type from UFCx integral. +/// @tparam T scalar type of kernel to extract. +/// @tparam U geometry type of kernel to extract. +/// @param integral UFCx integral to retrieve the kernel from. +/// @return Kernel callback. +template > +constexpr kernel_t extract_kernel(const ufcx_integral* integral) +{ + if constexpr (std::is_same_v) + return integral->tabulate_tensor_float32; + else if constexpr (std::is_same_v) + return integral->tabulate_tensor_float64; + else if constexpr (std::is_same_v> + && has_complex_ufcx_kernels()) + { + return reinterpret_cast>( + integral->tabulate_tensor_complex64); + } + else if constexpr (std::is_same_v> + && has_complex_ufcx_kernels()) + { + return reinterpret_cast>( + integral->tabulate_tensor_complex128); + } + else + throw std::runtime_error("Could not extract kernel from ufcx integral."); +} + +} // namespace dolfinx::fem::impl \ No newline at end of file diff --git a/cpp/dolfinx/fem/pack.h b/cpp/dolfinx/fem/pack.h index 0172101447..9b08e845d8 100644 --- a/cpp/dolfinx/fem/pack.h +++ b/cpp/dolfinx/fem/pack.h @@ -322,6 +322,28 @@ void pack_coefficients(const Form& form, } break; } + case IntegralType::vertex: + { + // 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 c3099306a6..ec7a8501ae 100644 --- a/cpp/dolfinx/fem/utils.cpp +++ b/cpp/dolfinx/fem/utils.cpp @@ -22,6 +22,7 @@ #include #include #include +#include #include #include @@ -137,7 +138,26 @@ fem::compute_integration_domains(fem::IntegralType integral_type, std::span entities) { const int tdim = topology.dim(); - const int dim = integral_type == IntegralType::cell ? tdim : tdim - 1; + + 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."); + } { // Create span of the owned entities (leaves off any ghosts) @@ -147,6 +167,30 @@ fem::compute_integration_domains(fem::IntegralType integral_type, entities = entities.first(std::distance(entities.begin(), it1)); } + auto get_connectivities = [tdim, &topology](int entity_dim) + -> std::pair>, + std::shared_ptr>> + { + auto e_to_c = topology.connectivity(entity_dim, tdim); + if (!e_to_c) + { + throw std::runtime_error( + std::format("Topology entity-to-cell connectivity has not been " + "computed for entity dim {}.", + entity_dim)); + } + + auto e_to_f = topology.connectivity(tdim, entity_dim); + if (!e_to_f) + { + throw std::runtime_error( + std::format("Topology cell-to-entity connectivity has not been " + "computed for entity dim {}.", + entity_dim)); + } + return {e_to_c, e_to_f}; + }; + std::vector entity_data; switch (integral_type) { @@ -155,74 +199,64 @@ fem::compute_integration_domains(fem::IntegralType integral_type, entity_data.insert(entity_data.begin(), entities.begin(), entities.end()); break; } - default: + case IntegralType::exterior_facet: { - auto f_to_c = topology.connectivity(tdim - 1, tdim); - if (!f_to_c) + auto [f_to_c, c_to_f] = get_connectivities(tdim - 1); + // Create list of tagged boundary facets + const std::vector bfacets = mesh::exterior_facet_indices(topology); + std::vector facets; + std::ranges::set_intersection(entities, bfacets, + std::back_inserter(facets)); + for (auto f : facets) { - throw std::runtime_error( - "Topology facet-to-cell connectivity has not been computed."); - } - - auto c_to_f = topology.connectivity(tdim, tdim - 1); - if (!c_to_f) - { - throw std::runtime_error( - "Topology cell-to-facet connectivity has not been computed."); + // Get the facet as a pair of (cell, local facet) + auto facet = impl::get_cell_facet_pairs<1>(f, f_to_c->links(f), *c_to_f); + entity_data.insert(entity_data.end(), facet.begin(), facet.end()); } + break; + } + case IntegralType::interior_facet: + { + auto [f_to_c, c_to_f] = get_connectivities(tdim - 1); - switch (integral_type) + // Create indicator for interprocess facets + assert(topology.index_map(tdim - 1)); + const std::vector& interprocess_facets + = topology.interprocess_facets(); + std::vector interprocess_marker( + topology.index_map(tdim - 1)->size_local() + + topology.index_map(tdim - 1)->num_ghosts(), + 0); + std::ranges::for_each(interprocess_facets, [&interprocess_marker](auto f) + { interprocess_marker[f] = 1; }); + for (auto f : entities) { - case IntegralType::exterior_facet: - { - // Create list of tagged boundary facets - const std::vector bfacets = mesh::exterior_facet_indices(topology); - std::vector facets; - std::ranges::set_intersection(entities, bfacets, - std::back_inserter(facets)); - for (auto f : facets) + if (f_to_c->num_links(f) == 2) { - // Get the facet as a pair of (cell, local facet) - auto facet - = impl::get_cell_facet_pairs<1>(f, f_to_c->links(f), *c_to_f); - entity_data.insert(entity_data.end(), facet.begin(), facet.end()); + // Get the facet as a pair of (cell, local facet) pairs, one + // for each cell + auto facets + = impl::get_cell_facet_pairs<2>(f, f_to_c->links(f), *c_to_f); + entity_data.insert(entity_data.end(), facets.begin(), facets.end()); } - } - break; - case IntegralType::interior_facet: - { - // Create indicator for interprocess facets - assert(topology.index_map(tdim - 1)); - const std::vector& interprocess_facets - = topology.interprocess_facets(); - std::vector interprocess_marker( - topology.index_map(tdim - 1)->size_local() - + topology.index_map(tdim - 1)->num_ghosts(), - 0); - std::ranges::for_each(interprocess_facets, [&interprocess_marker](auto f) - { interprocess_marker[f] = 1; }); - for (auto f : entities) + else if (interprocess_marker[f]) { - if (f_to_c->num_links(f) == 2) - { - // Get the facet as a pair of (cell, local facet) pairs, one - // for each cell - auto facets - = impl::get_cell_facet_pairs<2>(f, f_to_c->links(f), *c_to_f); - entity_data.insert(entity_data.end(), facets.begin(), facets.end()); - } - else if (interprocess_marker[f]) - { - throw std::runtime_error( - "Cannot compute interior facet integral over interprocess facet. " - "Use \"shared facet\" ghost mode when creating the mesh."); - } + throw std::runtime_error( + "Cannot compute interior facet integral over interprocess facet. " + "Use \"shared facet\" ghost mode when creating the mesh."); } } break; - default: - throw std::runtime_error( - "Cannot compute integration domains. Integral type not supported."); + } + case IntegralType::vertex: + { + auto [v_to_c, c_to_v] = get_connectivities(0); + for (auto vertex : entities) + { + std::array pair = impl::get_cell_vertex_pairs<1>( + vertex, v_to_c->links(vertex), *c_to_v); + + entity_data.insert(entity_data.end(), pair.begin(), pair.end()); } } } diff --git a/cpp/dolfinx/fem/utils.h b/cpp/dolfinx/fem/utils.h index a0e67af34a..22436c8e50 100644 --- a/cpp/dolfinx/fem/utils.h +++ b/cpp/dolfinx/fem/utils.h @@ -1,4 +1,5 @@ -// Copyright (C) 2013-2020 Johan Hake, Jan Blechta and Garth N. Wells +// Copyright (C) 2013-2025 Johan Hake, Jan Blechta, Garth N. Wells and Paul T. +// Kühner // // This file is part of DOLFINx (https://www.fenicsproject.org) // @@ -14,10 +15,13 @@ #include "Form.h" #include "Function.h" #include "FunctionSpace.h" +#include "kernel.h" #include "sparsitybuild.h" #include #include #include +#include +#include #include #include #include @@ -28,6 +32,7 @@ #include #include #include +#include #include #include #include @@ -86,6 +91,35 @@ get_cell_facet_pairs(std::int32_t f, std::span cells, return cell_local_facet_pairs; } + +/// Helper function to get an array of of (cell, local_vertex) pairs +/// corresponding to a given vertex index. +/// @note If the vertex is connected to multiple cells, the first one is picked. +/// @param[in] v vertex index +/// @param[in] cells List of cells incident to the vertex +/// @param[in] c_to_v Cell to vertex connectivity +/// @return Vector of (cell, local_vertex) pairs +template +std::array +get_cell_vertex_pairs(std::int32_t v, std::span cells, + const graph::AdjacencyList& c_to_v) +{ + static_assert(num_cells == 1); // Patch assembly not supported. + + assert(cells.size() > 0); + + // Use first cell for assembly over by default + std::int32_t cell = cells[0]; + + // Find local index of vertex within cell + auto cell_vertices = c_to_v.links(cell); + auto it = std::ranges::find(cell_vertices, v); + assert(it != cell_vertices.end()); + std::int32_t local_index = std::distance(cell_vertices.begin(), it); + + return {cell, local_index}; +} + } // namespace impl /// @brief Given an integral type and a set of entities, computes and @@ -109,11 +143,18 @@ get_cell_facet_pairs(std::int32_t f, std::span cells, /// /// @param[in] integral_type Integral type. /// @param[in] topology Mesh topology. -/// @param[in] entities List of mesh entities. For -/// `integral_type==IntegralType::cell`, `entities` should be cell -/// indices. For other `IntegralType`, `entities` should be facet -/// indices. -/// @return List of integration entity data. +/// @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) std::vector compute_integration_domains(IntegralType integral_type, const mesh::Topology& topology, @@ -449,10 +490,17 @@ Form create_form_factory( // integral offsets. Since the UFL forms for each type of cell should be // the same, I think this assumption is OK. const int* integral_offsets = ufcx_forms[0].get().form_integral_offsets; - std::vector num_integrals_type(3); - for (int i = 0; i < 3; ++i) + std::array num_integrals_type; + for (std::size_t i = 0; i < num_integrals_type.size(); ++i) num_integrals_type[i] = integral_offsets[i + 1] - integral_offsets[i]; + // Create vertices, if required + if (num_integrals_type[vertex] > 0) + { + mesh->topology_mutable()->create_connectivity(0, tdim); + mesh->topology_mutable()->create_connectivity(tdim, 0); + } + // Create facets, if required if (num_integrals_type[exterior_facet] > 0 or num_integrals_type[interior_facet] > 0) @@ -464,8 +512,6 @@ Form create_form_factory( // Get list of integral IDs, and load tabulate tensor into memory for // each - using kern_t = std::function; std::map, integral_data> integrals; auto check_geometry_hash @@ -508,30 +554,7 @@ Form create_form_factory( active_coeffs.push_back(j); } - kern_t k = nullptr; - if constexpr (std::is_same_v) - k = integral->tabulate_tensor_float32; -#ifndef DOLFINX_NO_STDC_COMPLEX_KERNELS - else if constexpr (std::is_same_v>) - { - k = reinterpret_cast*, const int*, - const unsigned char*, void*)>( - integral->tabulate_tensor_complex64); - } -#endif // DOLFINX_NO_STDC_COMPLEX_KERNELS - else if constexpr (std::is_same_v) - k = integral->tabulate_tensor_float64; -#ifndef DOLFINX_NO_STDC_COMPLEX_KERNELS - else if constexpr (std::is_same_v>) - { - k = reinterpret_cast*, const int*, - const unsigned char*, void*)>( - integral->tabulate_tensor_complex128); - } -#endif // DOLFINX_NO_STDC_COMPLEX_KERNELS - + impl::kernel_t k = impl::extract_kernel(integral); if (!k) { throw std::runtime_error( @@ -595,30 +618,7 @@ Form create_form_factory( active_coeffs.push_back(j); } - kern_t k = nullptr; - if constexpr (std::is_same_v) - k = integral->tabulate_tensor_float32; -#ifndef DOLFINX_NO_STDC_COMPLEX_KERNELS - else if constexpr (std::is_same_v>) - { - k = reinterpret_cast*, const int*, - const unsigned char*, void*)>( - integral->tabulate_tensor_complex64); - } -#endif // DOLFINX_NO_STDC_COMPLEX_KERNELS - else if constexpr (std::is_same_v) - k = integral->tabulate_tensor_float64; -#ifndef DOLFINX_NO_STDC_COMPLEX_KERNELS - else if constexpr (std::is_same_v>) - { - k = reinterpret_cast*, const int*, - const unsigned char*, void*)>( - integral->tabulate_tensor_complex128); - } -#endif // DOLFINX_NO_STDC_COMPLEX_KERNELS - assert(k); + impl::kernel_t k = impl::extract_kernel(integral); // Build list of entities to assembler over const std::vector bfacets = mesh::exterior_facet_indices(*topology); @@ -703,29 +703,7 @@ Form create_form_factory( active_coeffs.push_back(j); } - kern_t k = nullptr; - if constexpr (std::is_same_v) - k = integral->tabulate_tensor_float32; -#ifndef DOLFINX_NO_STDC_COMPLEX_KERNELS - else if constexpr (std::is_same_v>) - { - k = reinterpret_cast*, const int*, - const unsigned char*, void*)>( - integral->tabulate_tensor_complex64); - } -#endif // DOLFINX_NO_STDC_COMPLEX_KERNELS - else if constexpr (std::is_same_v) - k = integral->tabulate_tensor_float64; -#ifndef DOLFINX_NO_STDC_COMPLEX_KERNELS - else if constexpr (std::is_same_v>) - { - k = reinterpret_cast*, const int*, - const unsigned char*, void*)>( - integral->tabulate_tensor_complex128); - } -#endif // DOLFINX_NO_STDC_COMPLEX_KERNELS + impl::kernel_t k = impl::extract_kernel(integral); assert(k); // Build list of entities to assembler over @@ -779,6 +757,87 @@ Form create_form_factory( } } + // Attach vertex kernels + { + for (std::size_t form_idx = 0; form_idx < ufcx_forms.size(); ++form_idx) + { + const ufcx_form& form = ufcx_forms[form_idx]; + + std::span ids(form.form_integral_ids + + integral_offsets[vertex], + num_integrals_type[vertex]); + auto sd = subdomains.find(IntegralType::vertex); + for (int i = 0; i < num_integrals_type[vertex]; ++i) + { + const int id = ids[i]; + ufcx_integral* integral + = form.form_integrals[integral_offsets[vertex] + i]; + assert(integral); + check_geometry_hash(*integral, form_idx); + + std::vector active_coeffs; + for (int j = 0; j < form.num_coefficients; ++j) + { + if (integral->enabled_coefficients[j]) + active_coeffs.push_back(j); + } + + impl::kernel_t k = impl::extract_kernel(integral); + assert(k); + + // Build list of entities to assembler over + auto v_to_c = topology->connectivity(0, tdim); + assert(v_to_c); + auto c_to_v = topology->connectivity(tdim, 0); + assert(c_to_v); + + // pack for a range of vertices a flattened list of cell index c_i and + // local vertex index l_i: + // [c_0, l_0, ..., c_n, l_n] + auto get_cells_and_vertices = [v_to_c, c_to_v](auto vertices_range) + { + std::vector cell_and_vertex; + cell_and_vertex.reserve(2 * vertices_range.size()); + for (std::int32_t vertex : vertices_range) + { + std::array pair = impl::get_cell_vertex_pairs<1>( + vertex, v_to_c->links(vertex), *c_to_v); + + cell_and_vertex.insert(cell_and_vertex.end(), pair.begin(), + pair.end()); + } + assert(cell_and_vertex.size() == 2 * vertices_range.size()); + return cell_and_vertex; + }; + + if (id == -1) + { + // Default vertex kernel operates on all (owned) vertices + std::int32_t num_vertices = topology->index_map(0)->size_local(); + std::vector cells_and_vertices = get_cells_and_vertices( + std::ranges::views::iota(0, num_vertices)); + + integrals.insert({{IntegralType::vertex, i, form_idx}, + {k, cells_and_vertices, active_coeffs}}); + } + else if (sd != subdomains.end()) + { + // NOTE: This requires that pairs are sorted + auto it = std::ranges::lower_bound(sd->second, id, std::less<>{}, + [](auto& a) { return a.first; }); + if (it != sd->second.end() and it->first == id) + { + integrals.insert({{IntegralType::vertex, i, form_idx}, + {k, + std::vector(it->second.begin(), + it->second.end()), + active_coeffs}}); + } + } + } + } + } + return Form(spaces, std::move(integrals), mesh, coefficients, constants, needs_facet_permutations, entity_maps); } diff --git a/python/dolfinx/fem/forms.py b/python/dolfinx/fem/forms.py index f373eaf370..9967e3eadf 100644 --- a/python/dolfinx/fem/forms.py +++ b/python/dolfinx/fem/forms.py @@ -163,6 +163,12 @@ def get_integration_domains( tdim = subdomain.topology.dim 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 + subdomain._cpp_object.topology.create_connectivity(0, tdim) + subdomain._cpp_object.topology.create_connectivity(tdim, 0) + # Compute integration domains only for each subdomain id in # the integrals. If a process has no integral entities, # insert an empty array. @@ -363,8 +369,8 @@ def _form(form): # Extract subdomain ids from ufcx_form subdomain_ids = {type: [] for type in sd.get(domain).keys()} - integral_offsets = [ufcx_form.form_integral_offsets[i] for i in range(4)] - for i in range(3): + 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]) diff --git a/python/test/unit/fem/test_assembler.py b/python/test/unit/fem/test_assembler.py index 8c8caecc67..e5a87ac496 100644 --- a/python/test/unit/fem/test_assembler.py +++ b/python/test/unit/fem/test_assembler.py @@ -1,4 +1,4 @@ -# Copyright (C) 2018-2025 Garth N. Wells, Jørgen S. Dokken +# Copyright (C) 2018-2025 Garth N. Wells, Jørgen S. Dokken and Paul T. Kühner # # This file is part of DOLFINx (https://www.fenicsproject.org) # @@ -14,10 +14,11 @@ import scipy.sparse import basix +import dolfinx.cpp import ufl from basix.ufl import element, mixed_element from dolfinx import cpp as _cpp -from dolfinx import default_real_type, default_scalar_type, fem, graph, la +from dolfinx import default_real_type, default_scalar_type, fem, graph, la, mesh from dolfinx.fem import ( Constant, Function, @@ -1560,3 +1561,271 @@ def test_mixed_quadrature(dtype, method): local_sum = assemble_scalar(compiled_form) global_sum = mesh.comm.allreduce(local_sum, op=MPI.SUM) assert np.isclose(global_contribution, global_sum, rtol=tol, atol=tol) + + +def vertex_to_dof_map(V): + """Create a map from the vertices of the mesh to the corresponding degree of freedom.""" + mesh = V.mesh + num_vertices_per_cell = dolfinx.cpp.mesh.cell_num_entities(mesh.topology.cell_type, 0) + + dof_layout2 = np.empty((num_vertices_per_cell,), dtype=np.int32) + for i in range(num_vertices_per_cell): + var = V.dofmap.dof_layout.entity_dofs(0, i) + assert len(var) == 1 + dof_layout2[i] = var[0] + + num_vertices = mesh.topology.index_map(0).size_local + mesh.topology.index_map(0).num_ghosts + + c_to_v = mesh.topology.connectivity(mesh.topology.dim, 0) + assert ( + c_to_v.num_nodes == 0 + or (c_to_v.offsets[1:] - c_to_v.offsets[:-1] == c_to_v.offsets[1]).all() + ), "Single cell type supported" + + vertex_to_dof_map = np.empty(num_vertices, dtype=np.int32) + vertex_to_dof_map[c_to_v.array] = V.dofmap.list[:, dof_layout2].reshape(-1) + return vertex_to_dof_map + + +@pytest.mark.parametrize( + "cell_type", + [ + mesh.CellType.interval, + mesh.CellType.triangle, + mesh.CellType.quadrilateral, + mesh.CellType.tetrahedron, + # mesh.CellType.pyramid, + mesh.CellType.prism, + mesh.CellType.hexahedron, + ], +) +@pytest.mark.parametrize("ghost_mode", [mesh.GhostMode.none, mesh.GhostMode.shared_facet]) +@pytest.mark.parametrize("dtype", [np.float32, np.float64, np.complex64, np.complex128]) +def test_vertex_integral_rank_0(cell_type, ghost_mode, dtype): + comm = MPI.COMM_WORLD + rdtype = np.real(dtype(0)).dtype + + msh = None + cell_dim = mesh.cell_dim(cell_type) + if cell_dim == 1: + msh = mesh.create_unit_interval(comm, 4, dtype=rdtype, ghost_mode=ghost_mode) + elif cell_dim == 2: + msh = mesh.create_unit_square( + comm, 4, 4, cell_type=cell_type, dtype=rdtype, ghost_mode=ghost_mode + ) + elif cell_dim == 3: + msh = mesh.create_unit_cube( + comm, 4, 4, 4, cell_type=cell_type, dtype=rdtype, ghost_mode=ghost_mode + ) + else: + raise RuntimeError("Bad dimension") + + vertex_map = msh.topology.index_map(0) + + def check_vertex_integral_against_sum(form, vertices, weighted=False): + """Weighting assumes the vertex integral to be weighted by a P1 function, each vertex value + corresponding to its global index.""" + weights = vertex_map.local_to_global(vertices) if weighted else np.ones_like(vertices) + expected_value_l = np.sum(msh.geometry.x[vertices, 0] * weights) + value_l = fem.assemble_scalar(fem.form(form, dtype=dtype)) + assert expected_value_l == pytest.approx(value_l, abs=5e4 * np.finfo(rdtype).eps) + + expected_value = comm.allreduce(expected_value_l) + value = comm.allreduce(value_l) + assert expected_value == pytest.approx(value, abs=5e4 * np.finfo(rdtype).eps) + + num_vertices = vertex_map.size_local + x = ufl.SpatialCoordinate(msh) + + # Full domain + check_vertex_integral_against_sum(x[0] * ufl.dP, np.arange(num_vertices)) + + # Split domain into left half of vertices (1) and right half of vertices (2) + vertices_1 = mesh.locate_entities(msh, 0, lambda x: x[0] <= 0.5) + vertices_1 = vertices_1[vertices_1 < num_vertices] + vertices_2 = mesh.locate_entities(msh, 0, lambda x: x[0] > 0.5) + vertices_2 = vertices_2[vertices_2 < num_vertices] + + tags = np.full(num_vertices, 1) + tags[vertices_2] = 2 + vertices = np.arange(0, num_vertices, dtype=np.int32) + meshtags = mesh.meshtags(msh, 0, vertices, tags) + + dP = ufl.Measure("dP", domain=msh, subdomain_data=meshtags) + + # Combinations of sub domains + check_vertex_integral_against_sum(x[0] * dP(1), vertices_1) + check_vertex_integral_against_sum(x[0] * dP(2), vertices_2) + check_vertex_integral_against_sum(x[0] * (dP(1) + dP(2)), np.arange(num_vertices)) + + V = fem.functionspace(msh, ("P", 1)) + u = fem.Function(V, dtype=dtype) + vertex_to_dof = vertex_to_dof_map(V) + vertices = np.arange(num_vertices + vertex_map.num_ghosts) + u.x.array[vertex_to_dof[vertices]] = vertex_map.local_to_global(vertices) + + check_vertex_integral_against_sum(u * x[0] * ufl.dP, np.arange(num_vertices), True) + check_vertex_integral_against_sum(u * x[0] * dP(1), vertices_1, True) + check_vertex_integral_against_sum(u * x[0] * dP(2), vertices_2, True) + check_vertex_integral_against_sum(u * x[0] * (dP(1) + dP(2)), np.arange(num_vertices), True) + + # Check custom packing + if cell_type is mesh.CellType.prism: + return + + msh.topology.create_entities(1) + msh.topology.create_connectivity(cell_dim - 1, cell_dim) + + v_to_c = msh.topology.connectivity(0, cell_dim) + c_to_v = msh.topology.connectivity(cell_dim, 0) + + cell_vertex_pairs = np.array([], dtype=np.int32) + for v in range(num_vertices): + c = v_to_c.links(v)[0] + v_l = np.where(c_to_v.links(c) == v)[0] + cell_vertex_pairs = np.append(cell_vertex_pairs, [c, *v_l]) + + # a) With subdomain_data + check_vertex_integral_against_sum( + x[0] * ufl.dP(domain=msh, subdomain_data=[(1, cell_vertex_pairs)], subdomain_id=1), + np.arange(num_vertices), + ) + + # 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)]} + + compiled_form = fem.compile_form( + comm, x[0] * ufl.dP, form_compiler_options={"scalar_type": dtype} + ) + form = fem.create_form(compiled_form, [], msh, subdomains, {}, {}, []) + expected_value_l = np.sum(msh.geometry.x[vertices, 0]) + value_l = fem.assemble_scalar(form) + assert expected_value_l == pytest.approx(value_l, abs=5e4 * np.finfo(rdtype).eps) + + +@pytest.mark.parametrize( + "cell_type", + [ + mesh.CellType.interval, + mesh.CellType.triangle, + mesh.CellType.quadrilateral, + mesh.CellType.tetrahedron, + # mesh.CellType.pyramid, + mesh.CellType.prism, + mesh.CellType.hexahedron, + ], +) +@pytest.mark.parametrize("ghost_mode", [mesh.GhostMode.none, mesh.GhostMode.shared_facet]) +@pytest.mark.parametrize("dtype", [np.float32, np.float64, np.complex64, np.complex128]) +def test_vertex_integral_rank_1(cell_type, ghost_mode, dtype): + comm = MPI.COMM_WORLD + rdtype = np.real(dtype(0)).dtype + + msh = None + cell_dim = mesh.cell_dim(cell_type) + if cell_dim == 1: + msh = mesh.create_unit_interval(comm, 4, ghost_mode=ghost_mode, dtype=rdtype) + elif cell_dim == 2: + msh = mesh.create_unit_square( + comm, 4, 4, cell_type=cell_type, ghost_mode=ghost_mode, dtype=rdtype + ) + elif cell_dim == 3: + msh = mesh.create_unit_cube( + comm, 4, 4, 4, cell_type=cell_type, ghost_mode=ghost_mode, dtype=rdtype + ) + else: + raise RuntimeError("Bad dimension") + + vertex_map = msh.topology.index_map(0) + num_vertices = vertex_map.size_local + + def check_vertex_integral_against_sum(form, vertices, weighted=False): + """Weighting assumes the vertex integral to be weighted by a P1 function, each vertex value + corresponding to its global index.""" + weights = vertex_map.local_to_global(vertices) if weighted else np.ones_like(vertices) + expected_value_l = np.zeros(num_vertices, dtype=rdtype) + expected_value_l[vertices] = msh.geometry.x[vertices, 0] * weights + value_l = fem.assemble_vector(fem.form(form, dtype=dtype)) + equal_l = np.allclose( + expected_value_l, np.real(value_l.array[:num_vertices]), atol=1e3 * np.finfo(rdtype).eps + ) + assert equal_l + assert comm.allreduce(equal_l, MPI.BAND) + + x = ufl.SpatialCoordinate(msh) + V = fem.functionspace(msh, ("P", 1)) + v = ufl.conj(ufl.TestFunction(V)) + + # Full domain + check_vertex_integral_against_sum(x[0] * v * ufl.dP, np.arange(num_vertices)) + + # Split domain into left half of vertices (1) and right half of vertices (2) + vertices_1 = mesh.locate_entities(msh, 0, lambda x: x[0] <= 0.5) + vertices_1 = vertices_1[vertices_1 < num_vertices] + vertices_2 = mesh.locate_entities(msh, 0, lambda x: x[0] > 0.5) + vertices_2 = vertices_2[vertices_2 < num_vertices] + + tags = np.full(num_vertices, 1) + tags[vertices_2] = 2 + vertices = np.arange(0, num_vertices, dtype=np.int32) + meshtags = mesh.meshtags(msh, 0, vertices, tags) + + dP = ufl.Measure("dP", domain=msh, subdomain_data=meshtags) + + check_vertex_integral_against_sum(x[0] * v * dP(1), vertices_1) + check_vertex_integral_against_sum(x[0] * v * dP(2), vertices_2) + check_vertex_integral_against_sum(x[0] * v * (dP(1) + dP(2)), np.arange(num_vertices)) + + V = fem.functionspace(msh, ("P", 1)) + u = fem.Function(V, dtype=dtype) + u.x.array[:] = vertex_map.local_to_global(np.arange(num_vertices + vertex_map.num_ghosts)) + vertex_to_dof = vertex_to_dof_map(V) + vertices = np.arange(num_vertices + vertex_map.num_ghosts) + u.x.array[vertex_to_dof[vertices]] = vertex_map.local_to_global(vertices) + + check_vertex_integral_against_sum(u * x[0] * v * ufl.dP, np.arange(num_vertices), True) + check_vertex_integral_against_sum(u * x[0] * v * dP(1), vertices_1, True) + check_vertex_integral_against_sum(u * x[0] * v * dP(2), vertices_2, True) + check_vertex_integral_against_sum(u * x[0] * v * (dP(1) + dP(2)), np.arange(num_vertices), True) + + # Check custom packing + if cell_type is mesh.CellType.prism: + return + + msh.topology.create_entities(1) + msh.topology.create_connectivity(cell_dim - 1, cell_dim) + + v_to_c = msh.topology.connectivity(0, cell_dim) + c_to_v = msh.topology.connectivity(cell_dim, 0) + + cell_vertex_pairs = np.array([], dtype=np.int32) + for v in range(num_vertices): + c = v_to_c.links(v)[0] + v_l = np.where(c_to_v.links(c) == v)[0] + cell_vertex_pairs = np.append(cell_vertex_pairs, [c, *v_l]) + + # a) With subdomain_data + v = ufl.conj(ufl.TestFunction(V)) + check_vertex_integral_against_sum( + x[0] * v * ufl.dP(domain=msh, subdomain_data=[(1, cell_vertex_pairs)], subdomain_id=1), + np.arange(num_vertices), + ) + + # 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)]} + + compiled_form = fem.compile_form( + comm, x[0] * v * ufl.dP, form_compiler_options={"scalar_type": dtype} + ) + form = fem.create_form(compiled_form, [V], msh, subdomains, {}, {}, []) + expected_value_l = np.sum(msh.geometry.x[vertices, 0]) + expected_value_l = np.zeros(num_vertices, dtype=rdtype) + expected_value_l[vertices] = msh.geometry.x[vertices, 0] + value_l = fem.assemble_vector(form) + assert expected_value_l == pytest.approx( + value_l.array[: expected_value_l.size], abs=5e4 * np.finfo(rdtype).eps + )