diff --git a/cpp/doc/source/api.rst b/cpp/doc/source/api.rst index f9b914c732b..8464594f7a2 100644 --- a/cpp/doc/source/api.rst +++ b/cpp/doc/source/api.rst @@ -14,6 +14,7 @@ generated documentation is `here `_. io la mesh + multigrid refinement * :ref:`genindex` diff --git a/cpp/doc/source/multigrid.rst b/cpp/doc/source/multigrid.rst new file mode 100644 index 00000000000..f0fe0e6ecb1 --- /dev/null +++ b/cpp/doc/source/multigrid.rst @@ -0,0 +1,5 @@ +Multigrid (``dolfinx::multigrid``) +==================================== + +.. doxygennamespace:: dolfinx::multigrid + :project: DOLFINx diff --git a/cpp/dolfinx/CMakeLists.txt b/cpp/dolfinx/CMakeLists.txt index 063be0e1eb3..405986f819e 100644 --- a/cpp/dolfinx/CMakeLists.txt +++ b/cpp/dolfinx/CMakeLists.txt @@ -17,6 +17,7 @@ set(DOLFINX_DIRS io la mesh + multigrid nls refinement ) diff --git a/cpp/dolfinx/multigrid/CMakeLists.txt b/cpp/dolfinx/multigrid/CMakeLists.txt new file mode 100644 index 00000000000..81450cdd7d3 --- /dev/null +++ b/cpp/dolfinx/multigrid/CMakeLists.txt @@ -0,0 +1,9 @@ +set(HEADERS_multigrid + ${CMAKE_CURRENT_SOURCE_DIR}/inclusion.h + PARENT_SCOPE +) + +target_sources( + dolfinx + PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/inclusion.cpp +) diff --git a/cpp/dolfinx/multigrid/dolfinx_multigrid.h b/cpp/dolfinx/multigrid/dolfinx_multigrid.h new file mode 100644 index 00000000000..562acb169ab --- /dev/null +++ b/cpp/dolfinx/multigrid/dolfinx_multigrid.h @@ -0,0 +1,11 @@ +#pragma once + +/// @brief Multigrid algorithms. +/// +namespace dolfinx::multigrid +{ +} + +// DOLFINx multigrid interface + +#include diff --git a/cpp/dolfinx/multigrid/inclusion.cpp b/cpp/dolfinx/multigrid/inclusion.cpp new file mode 100644 index 00000000000..9e6b66a9081 --- /dev/null +++ b/cpp/dolfinx/multigrid/inclusion.cpp @@ -0,0 +1,78 @@ +// 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 + +#include "inclusion.h" + +#include +#include +#include +#include +#include + +#include + +#include "dolfinx/common/IndexMap.h" +#include "dolfinx/la/SparsityPattern.h" +#include "dolfinx/mesh/Mesh.h" + +namespace dolfinx::multigrid +{ + +template +std::vector +inclusion_mapping(const dolfinx::mesh::Mesh& mesh_from, + const dolfinx::mesh::Mesh& mesh_to) +{ + { + // Check comms equal + int result; + MPI_Comm_compare(mesh_from.comm(), mesh_to.comm(), &result); + assert(result == MPI_CONGRUENT); + } + + const common::IndexMap& im_from = *mesh_from.topology()->index_map(0); + const common::IndexMap& im_to = *mesh_to.topology()->index_map(0); + + std::vector map(im_from.size_local() + im_from.num_ghosts(), + -1); + + std::span x_from = mesh_from.geometry().x(); + std::span x_to = mesh_to.geometry().x(); + + for (std::int32_t i = 0; i < im_from.size_local() + im_from.num_ghosts(); i++) + { + std::ranges::subrange vertex_from(std::next(x_from.begin(), 3 * i), + std::next(x_from.begin(), 3 * (i + 1))); + for (std::int32_t j = 0; j < im_to.size_local() + im_to.num_ghosts(); j++) + { + std::ranges::subrange vertex_to(std::next(x_to.begin(), 3 * j), + std::next(x_to.begin(), 3 * (j + 1))); + + if (std::ranges::equal( + vertex_from, vertex_to, [](auto a, auto b) + { return std::abs(a - b) <= std::numeric_limits::epsilon(); })) + { + assert(map[i] == -1); + map[i] = j; + break; + } + } + } + + return map; +} + +/// @cond +template std::vector +inclusion_mapping(const dolfinx::mesh::Mesh& mesh_from, + const dolfinx::mesh::Mesh& mesh_to); + +template std::vector +inclusion_mapping(const dolfinx::mesh::Mesh& mesh_from, + const dolfinx::mesh::Mesh& mesh_to); +/// @endcond + +} // namespace dolfinx::multigrid diff --git a/cpp/dolfinx/multigrid/inclusion.h b/cpp/dolfinx/multigrid/inclusion.h new file mode 100644 index 00000000000..75c9a4e7496 --- /dev/null +++ b/cpp/dolfinx/multigrid/inclusion.h @@ -0,0 +1,39 @@ +// 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 + +#include "dolfinx/mesh/Mesh.h" + +namespace dolfinx::multigrid +{ + +/// @brief Computes an inclusion map: a map between vertex indices from one mesh +/// to another. +/// +/// @param mesh_from Domain of the map +/// @param mesh_to Range of the map +/// +/// @return Inclusion map, the `i`-th component is the vertex index of the +/// vertex with the same coordinates in `mesh_to` and `-1` if it can not be +/// found (locally!) in `mesh_to`. If `map[i] != -1` it holds +/// `mesh_from.geometry.x()[i:i+3] == mesh_to.geometry.x()[map[i]:map[i]+3]`. +/// +/// @note Invoking `inclusion_map` on a `(mesh_coarse, mesh_fine)` tuple, where +/// `mesh_fine` is produced by refinement with +/// `IdentityPartitionerPlaceholder()` option, the returned `map` is guaranteed +/// to match all vertices for all locally owned vertices (not for the ghost +/// vertices). +template +std::vector +inclusion_mapping(const dolfinx::mesh::Mesh& mesh_from, + const dolfinx::mesh::Mesh& mesh_to); + +} // namespace dolfinx::multigrid diff --git a/cpp/test/CMakeLists.txt b/cpp/test/CMakeLists.txt index 3f812819a10..b91d7277757 100644 --- a/cpp/test/CMakeLists.txt +++ b/cpp/test/CMakeLists.txt @@ -67,6 +67,7 @@ add_executable( mesh/refinement/interval.cpp mesh/refinement/option.cpp mesh/refinement/rectangle.cpp + multigrid/inclusion.cpp ${CMAKE_CURRENT_BINARY_DIR}/expr.c ${CMAKE_CURRENT_BINARY_DIR}/poisson.c ) diff --git a/cpp/test/multigrid/inclusion.cpp b/cpp/test/multigrid/inclusion.cpp new file mode 100644 index 00000000000..e016849040d --- /dev/null +++ b/cpp/test/multigrid/inclusion.cpp @@ -0,0 +1,111 @@ +// 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 + +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +#include +#include +#include + +#include + +#include +#include +#include +#include +#include + +using namespace dolfinx; +using namespace Catch::Matchers; + +template +void CHECK_inclusion_map(const mesh::Mesh& mesh_coarse, + const mesh::Mesh& mesh_fine, + const std::vector& map) +{ + const common::IndexMap& im_from = *mesh_coarse.topology()->index_map(0); + const common::IndexMap& im_to = *mesh_fine.topology()->index_map(0); + CHECK( + map.size() + == static_cast(im_from.size_local() + im_from.num_ghosts())); + for (std::size_t i = 0; i < map.size(); i++) + { + if (map[i] == -1) + continue; + + CHECK(0 <= map[i]); + CHECK(map[i] <= im_to.size_local() + im_to.num_ghosts()); + + auto x = std::array{mesh_coarse.geometry().x()[3 * i], + mesh_coarse.geometry().x()[3 * i + 1], + mesh_coarse.geometry().x()[3 * i + 2]}; + auto global_x = std::array{mesh_fine.geometry().x()[3 * map[i]], + mesh_fine.geometry().x()[3 * map[i] + 1], + mesh_fine.geometry().x()[3 * map[i] + 2]}; + + CHECK_THAT(x, Catch::Matchers::RangeEquals(global_x)); + } +} + +template +void TEST_inclusion(dolfinx::mesh::Mesh&& mesh_coarse) +{ + mesh_coarse.topology()->create_entities(1); + + std::array, + 3> + ghost_modes + = {mesh::create_cell_partitioner(mesh::GhostMode::none), + mesh::create_cell_partitioner(mesh::GhostMode::shared_facet), + refinement::IdentityPartitionerPlaceholder()}; + for (const auto& ghost_mode : ghost_modes) + { + auto [mesh_fine, parent_cell, parent_facet] + = refinement::refine(mesh_coarse, std::nullopt, ghost_mode); + mesh_fine.topology()->create_connectivity(1, 0); + mesh_fine.topology()->create_connectivity(0, 1); + std::vector inclusion_map + = multigrid::inclusion_mapping(mesh_coarse, mesh_fine); + + if (std::holds_alternative( + ghost_mode)) + CHECK(std::ranges::all_of(inclusion_map, [](auto e) { return e >= 0; })); + + CHECK_inclusion_map(mesh_coarse, mesh_fine, inclusion_map); + } +} + +TEMPLATE_TEST_CASE("Inclusion (interval)", "[multigrid][inclusion]", double, + float) +{ + TEST_inclusion( + dolfinx::mesh::create_interval(MPI_COMM_WORLD, 10, {0.0, 1.0})); +} + +TEMPLATE_TEST_CASE("Inclusion (triangle)", "[multigrid][inclusion]", double, + float) +{ + TEST_inclusion(dolfinx::mesh::create_rectangle( + MPI_COMM_WORLD, {{{0, 0}, {1, 1}}}, {5, 5}, mesh::CellType::triangle)); +} + +TEMPLATE_TEST_CASE("Inclusion (tetrahedron)", "[multigrid][inclusion]", double, + float) +{ + TEST_inclusion(dolfinx::mesh::create_box( + MPI_COMM_WORLD, {{{0, 0, 0}, {1, 1, 1}}}, {5, 5, 5}, + mesh::CellType::tetrahedron)); +} diff --git a/python/CMakeLists.txt b/python/CMakeLists.txt index 5893048c3c5..790751e0e99 100644 --- a/python/CMakeLists.txt +++ b/python/CMakeLists.txt @@ -53,6 +53,7 @@ nanobind_add_module( dolfinx/wrappers/la.cpp dolfinx/wrappers/log.cpp dolfinx/wrappers/mesh.cpp + dolfinx/wrappers/multigrid.cpp dolfinx/wrappers/petsc.cpp dolfinx/wrappers/refinement.cpp ) diff --git a/python/dolfinx/multigrid.py b/python/dolfinx/multigrid.py new file mode 100644 index 00000000000..a66bd21cbb6 --- /dev/null +++ b/python/dolfinx/multigrid.py @@ -0,0 +1,43 @@ +# 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 + +import numpy as np +from numpy.typing import NDArray + +from dolfinx.cpp.multigrid import inclusion_mapping_float32 as _inclusion_mapping_float32 +from dolfinx.cpp.multigrid import inclusion_mapping_float64 as _inclusion_mapping_float64 +from dolfinx.mesh import Mesh + +__all__ = ["inclusion_mapping"] + + +def inclusion_mapping(mesh_from: Mesh, mesh_to: Mesh) -> NDArray[np.int64]: + """Computes an inclusion map: a map between vertex indices from one + mesh to another. + + Args: + mesh_from: Domain of the map. + mesh_to: Range of the map. + + Returns: + Inclusion map, the `i`-th component is the vertex index of the + vertex with the same coordinates in `mesh_to` and `-1` if it can + not be found (locally!) in `mesh_to`. If `map[i] != -1` it holds + `mesh_from.geometry.x[i] == mesh_to.geometry.x[map[i]]`. + + Note: + Invoking `inclusion_map` on a `(mesh_coarse, mesh_fine)` tuple, + where `mesh_fine` is produced by refinement with + `IdentityPartitionerPlaceholder()`option, the returned `map` is + guaranteed to match all vertices for all locally owned vertices + (not for the ghost vertices). + """ + if np.issubdtype(mesh_from.geometry.x.dtype, np.float32): + return _inclusion_mapping_float32(mesh_from._cpp_object, mesh_to._cpp_object) + elif np.issubdtype(mesh_from.geometry.x.dtype, np.float64): + return _inclusion_mapping_float64(mesh_from._cpp_object, mesh_to._cpp_object) + else: + raise RuntimeError("Unsupported mesh dtype.") diff --git a/python/dolfinx/wrappers/dolfinx.cpp b/python/dolfinx/wrappers/dolfinx.cpp index 9c0864f9ef1..21bcec93f4f 100644 --- a/python/dolfinx/wrappers/dolfinx.cpp +++ b/python/dolfinx/wrappers/dolfinx.cpp @@ -25,6 +25,8 @@ void la(nb::module_& m); void mesh(nb::module_& m); void nls(nb::module_& m); void refinement(nb::module_& m); +void multigrid(nb::module_& m); + } // namespace dolfinx_wrappers NB_MODULE(cpp, m) @@ -74,6 +76,10 @@ NB_MODULE(cpp, m) nb::module_ refinement = m.def_submodule("refinement", "Refinement module"); dolfinx_wrappers::refinement(refinement); + // Create multigrid submodule + nb::module_ multigrid = m.def_submodule("multigrid", "Multigrid module"); + dolfinx_wrappers::multigrid(multigrid); + #if defined(HAS_PETSC) && defined(HAS_PETSC4PY) // PETSc-specific wrappers nb::module_ nls = m.def_submodule("nls", "Nonlinear solver module"); diff --git a/python/dolfinx/wrappers/multigrid.cpp b/python/dolfinx/wrappers/multigrid.cpp new file mode 100644 index 00000000000..4913700886f --- /dev/null +++ b/python/dolfinx/wrappers/multigrid.cpp @@ -0,0 +1,42 @@ +// 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 + +#include + +#include +#include + +#include +#include + +#include "dolfinx_wrappers/array.h" + +namespace nb = nanobind; + +namespace dolfinx_wrappers +{ + +template +void declare_inlcusion_mapping(nb::module_& m, const std::string& type) +{ + m.def(("inclusion_mapping_" + type).c_str(), + [](const dolfinx::mesh::Mesh& mesh_from, + const dolfinx::mesh::Mesh& mesh_to) + { + return dolfinx_wrappers::as_nbarray( + dolfinx::multigrid::inclusion_mapping(mesh_from, mesh_to)); + }, + nb::arg("mesh_from"), nb::arg("mesh_to"), + "Computes inclusion mapping between two meshes"); +} + +void multigrid(nb::module_& m) +{ + declare_inlcusion_mapping(m, "float32"); + declare_inlcusion_mapping(m, "float64"); +} + +} // namespace dolfinx_wrappers diff --git a/python/test/unit/multigrid/test_inclusion.py b/python/test/unit/multigrid/test_inclusion.py new file mode 100644 index 00000000000..33728b4dce2 --- /dev/null +++ b/python/test/unit/multigrid/test_inclusion.py @@ -0,0 +1,66 @@ +# 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 + +from mpi4py import MPI + +import numpy as np +import pytest +from numpy import typing as npt + +from dolfinx.mesh import ( + GhostMode, + IdentityPartitionerPlaceholder, + Mesh, + create_cell_partitioner, + create_unit_cube, + create_unit_interval, + create_unit_square, + refine, +) +from dolfinx.multigrid import inclusion_mapping + + +def check_inclusion_map( + mesh_from: Mesh, mesh_to: Mesh, map: npt.NDArray[np.int32], expect_all: bool = False +): + if expect_all: + assert np.all(map[: mesh_from.topology.index_map(0).size_local] >= 0) + + for i in range(len(map)): + if map[i] == -1: + continue + assert np.allclose(mesh_from.geometry.x[i], mesh_to.geometry.x[map[i]]) + + +@pytest.mark.parametrize("gdim", [1, 2, 3]) +@pytest.mark.parametrize("dtype", [np.float32, np.float64]) +@pytest.mark.parametrize("ghost_mode_coarse", [GhostMode.none, GhostMode.shared_facet]) +@pytest.mark.parametrize( + "partitioner_fine", + [ + create_cell_partitioner(GhostMode.none), + create_cell_partitioner(GhostMode.shared_facet), + IdentityPartitionerPlaceholder(), + ], +) +def test_inclusion_map(gdim, dtype, ghost_mode_coarse, partitioner_fine): + comm = MPI.COMM_WORLD + if gdim == 1: + mesh = create_unit_interval(comm, 10, ghost_mode=ghost_mode_coarse, dtype=dtype) + elif gdim == 2: + mesh = create_unit_square(comm, 5, 5, ghost_mode=ghost_mode_coarse, dtype=dtype) + else: + mesh = create_unit_cube(comm, 5, 5, 5, ghost_mode=ghost_mode_coarse, dtype=dtype) + + mesh.topology.create_entities(1) + mesh_fine, _, _ = refine(mesh, partitioner=partitioner_fine) + map = inclusion_mapping(mesh, mesh_fine) + check_inclusion_map( + mesh, mesh_fine, map, type(partitioner_fine) is IdentityPartitionerPlaceholder + ) + + map = inclusion_mapping(mesh_fine, mesh) + check_inclusion_map(mesh_fine, mesh, map)