diff --git a/cpp/dolfinx/refinement/mark.h b/cpp/dolfinx/refinement/mark.h index 395f54fcd2..17fccbb484 100644 --- a/cpp/dolfinx/refinement/mark.h +++ b/cpp/dolfinx/refinement/mark.h @@ -1,4 +1,4 @@ -// Copyright (C) 2026 Paul T. Kühner +// Copyright (C) 2026 Paul T. Kühner, Jack S. Hale // // This file is part of DOLFINX (https://www.fenicsproject.org) // @@ -8,11 +8,15 @@ #include #include +#include #include #include #include #include +#include +#include #include +#include #include #include "dolfinx/common/MPI.h" @@ -20,38 +24,155 @@ namespace dolfinx::refinement { -/// @brief Maximum marking of a marker. +namespace impl +{ + +/// @brief Computes local threshold marking of indicators. +/// +/// Helper for other marking routines. +/// +/// Returns the indices \f$ i \f$ of the indicators \f$ \eta_i \f$ that satisfy +/// the threshold: \f$ \eta_i > \text{threshold} \f$. /// -/// @param[in] marker Input marker (local) - usually an error indicator per -/// entity -/// @param[in] theta Cut off parameter, 0 < θ < 1 -/// @param[in] comm Communicator over which the maximum is computed. -/// @return Indices (local) of marker elements, which satisfy: marker_i ≥ θ -/// max(marker). +/// @param[in] indicators Indicators (local) \f$ \eta_i \f$ - +/// usually an error indicator associated with mesh entity \f$ i \f$. +/// @param[in] threshold Threshold value; indicators greater than this are +/// marked. +/// @return Local indices of marked entities. template -std::vector mark_maximum(std::span marker, T theta, - MPI_Comm comm) +std::vector mark_threshold(std::span indicators, + T threshold) { - if ((theta <= 0) || (theta >= 1)) - throw std::invalid_argument("Theta needs to fullfill 0 < θ < 1."); - - T max = marker.empty() ? std::numeric_limits::lowest() - : std::ranges::max(marker); - MPI_Allreduce(MPI_IN_PLACE, &max, 1, dolfinx::MPI::mpi_t, MPI_MAX, comm); - - auto mark = [=](T e) { return e > theta * max; }; + auto mark = [threshold](T e) { return e > threshold; }; std::vector indices; - indices.reserve(std::ranges::count_if(marker, mark)); + indices.reserve(std::ranges::count_if(indicators, mark)); - for (std::int32_t i = 0; i < static_cast(marker.size()); ++i) + for (std::int32_t i = 0; i < static_cast(indicators.size()); + ++i) { - if (mark(marker[i])) + if (mark(indicators[i])) indices.push_back(i); } - spdlog::info("Marking (max) {} / {} (local) entities.", indices.size(), - marker.size()); + return indices; +} + +} // namespace impl + +/// @brief Computes maximum-based marking of indicators. +/// +/// Returns the indices \f$ i \f$ of the indicators \f$ \eta_i \f$ that satisfy +/// the maximum threshold: \f$ \eta_i > \theta \max_j \eta_j \f$. +/// +/// @param[in] comm Communicator to compute the maximum over. +/// @param[in] indicators Indicators (local) \f$ \eta_i \f$ - +/// usually an error indicator associated with mesh entity \f$ i \f$. +/// @param[in] theta Parameter, \f$ 0 < \theta < 1 \f$. +/// @return Local indices of marked entities. +template +std::vector mark_maximum(MPI_Comm comm, + std::span indicators, T theta) +{ + if ((theta <= 0) || (theta >= 1)) + throw std::invalid_argument("theta must fulfill 0 < theta < 1."); + + T max = indicators.empty() ? std::numeric_limits::lowest() + : std::ranges::max(indicators); + MPI_Allreduce(MPI_IN_PLACE, &max, 1, dolfinx::MPI::mpi_t, MPI_MAX, comm); + + std::vector indices + = impl::mark_threshold(indicators, theta * max); + + spdlog::info("Marking (maximum) {} / {} (local) entities.", indices.size(), + indicators.size()); + + return indices; +} + +/// @brief Computes equidistribution threshold marking of indicators. +/// +/// Returns the indices \f$ i \f$ of the indicators \f$ \eta_i \f$ that satisfy +/// the equidistribution threshold: \f$\eta_i > \theta +/// \frac{||\eta||}{\sqrt{N}} \f$ where \f$ N \f$ is the (global) number of +/// indicators. +/// +/// @warning `indicators` must contain owned entities only. Ghost values cause +/// double-counting in the global sum reductions. +/// +/// @param[in] comm Communicator over which the global equidistribution +/// threshold is computed. +/// @param[in] indicators Indicators for owned local entities \f$ \eta_i \f$ +/// - usually associated with mesh entity \f$ i \f$. +/// @param[in] theta Parameter, \f$ 0 < \theta < 1 \f$. +/// @return Local indices of indicators that satisfy the threshold. +template +std::vector +mark_equidistribution(MPI_Comm comm, std::span indicators, T theta) +{ + if ((theta <= 0) || (theta >= 1)) + throw std::invalid_argument("theta must fulfill 0 < theta < 1."); + + T norm = std::inner_product(indicators.begin(), indicators.end(), + indicators.begin(), T{0}); + + MPI_Allreduce(MPI_IN_PLACE, &norm, 1, dolfinx::MPI::mpi_t, MPI_SUM, comm); + + T sqrt_norm = std::sqrt(norm); + + // int64_t gives headroom for global sum across ranks. + std::int64_t count = indicators.size(); + MPI_Allreduce(MPI_IN_PLACE, &count, 1, dolfinx::MPI::mpi_t, + MPI_SUM, comm); + + std::vector indices = impl::mark_threshold( + indicators, theta * sqrt_norm / std::sqrt(static_cast(count))); + + spdlog::info("Marking (equidistribution) {} / {} (local) entities.", + indices.size(), indicators.size()); + + return indices; +} + +/// @brief Computes equidistribution threshold marking of a squared indicator. +/// +/// Returns the indices \f$i\f$ of the squared indicators \f$ \eta_i^2 \f$ that +/// satisfy the equidistribution threshold: \f$ \eta_i^2 > \theta^2 +/// \frac{||\eta||^2}{N} \f$ where \f$ N \f$ is the (global) number of +/// indicators. +/// +/// @warning `squared_indicators` must contain owned entities only. Ghost +/// values cause double-counting in the global sum reductions. +/// +/// @param[in] comm Communicator over which the global equidistribution +/// threshold is computed. +/// @param[in] squared_indicators Input squared indicators for owned local +/// entities \f$ \eta^2_i \f$ - usually associated with mesh entity \f$ i \f$. +/// @param[in] theta Parameter, \f$ 0 < \theta < 1 \f$. +/// @return Local indices of squared indicators that satisfy the threshold. +template +std::vector +mark_equidistribution_squared(MPI_Comm comm, + std::span squared_indicators, T theta) +{ + if ((theta <= 0) || (theta >= 1)) + throw std::invalid_argument("theta must fulfill 0 < theta < 1."); + + T norm = std::accumulate(squared_indicators.begin(), squared_indicators.end(), + T{0}); + + MPI_Allreduce(MPI_IN_PLACE, &norm, 1, dolfinx::MPI::mpi_t, MPI_SUM, comm); + + // int64_t gives headroom for global sum across ranks. + std::int64_t count = squared_indicators.size(); + MPI_Allreduce(MPI_IN_PLACE, &count, 1, dolfinx::MPI::mpi_t, + MPI_SUM, comm); + + std::vector indices = impl::mark_threshold( + squared_indicators, theta * theta * norm / static_cast(count)); + + spdlog::info("Marking (equidistribution) {} of {} local entities.", + indices.size(), squared_indicators.size()); return indices; } diff --git a/cpp/test/mesh/refinement/mark.cpp b/cpp/test/mesh/refinement/mark.cpp index b4bed6329d..9ebe7ca45e 100644 --- a/cpp/test/mesh/refinement/mark.cpp +++ b/cpp/test/mesh/refinement/mark.cpp @@ -17,8 +17,8 @@ using namespace dolfinx::refinement; TEMPLATE_TEST_CASE("Mark maximum empty", "[refinement][mark][maximum]", double, float) { - std::vector marker; - auto indices = mark_maximum(marker, .5, MPI_COMM_WORLD); + std::vector indicators; + auto indices = mark_maximum(MPI_COMM_WORLD, indicators, .5); CHECK(indices.size() == 0); } @@ -26,28 +26,86 @@ TEMPLATE_TEST_CASE("Mark maximum", "[refinement][mark][maximum]", double, float) { MPI_Comm comm = MPI_COMM_WORLD; - std::vector marker; - marker.reserve(10); + std::vector indicators; + indicators.reserve(10); for (std::size_t i = 0; i < 10; i++) - marker.push_back(10 * dolfinx::MPI::rank(comm) + i); + indicators.push_back(10 * dolfinx::MPI::rank(comm) + i); TestType theta = 0.5; - auto indices = mark_maximum(marker, theta, comm); + auto indices = mark_maximum(comm, indicators, theta); CHECK(std::ranges::all_of( - indices, [&](auto e) - { return (0 <= e) && (e <= static_cast(marker.size())); })); + indices, + [&](auto e) + { + return (0 <= e) && (e <= static_cast(indicators.size())); + })); TestType max = dolfinx::MPI::size(comm) * 10 - 1; auto mark = [=](auto e) { return e > theta * max; }; - CHECK(std::ranges::count_if(marker, mark) + CHECK(std::ranges::count_if(indicators, mark) == static_cast(indices.size())); - for (std::int32_t i = 0; i < static_cast(marker.size()); ++i) + for (std::int32_t i = 0; i < static_cast(indicators.size()); + ++i) { - bool expect_marked = mark(marker[i]); + bool expect_marked = mark(indicators[i]); bool marked = std::ranges::find(indices, i) != indices.end(); CHECK(expect_marked == marked); } } + +TEMPLATE_TEST_CASE("Mark equidistribution empty", + "[refinement][mark][equidistribution]", double, float) +{ + std::vector indicators; + auto indices + = mark_equidistribution(MPI_COMM_WORLD, indicators, .5); + CHECK(indices.size() == 0); +} + +TEMPLATE_TEST_CASE("Mark equidistribution", + "[refinement][mark][equidistribution]", double, float) +{ + MPI_Comm comm = MPI_COMM_WORLD; + + std::vector indicators; + indicators.reserve(10); + for (std::size_t i = 0; i < 10; i++) + indicators.push_back(10 * dolfinx::MPI::rank(comm) + i); + + TestType theta = 0.5; + auto indices = mark_equidistribution(comm, indicators, theta); + + CHECK(std::ranges::all_of( + indices, + [&](auto e) + { + return (0 <= e) && (e <= static_cast(indicators.size())); + })); + + std::int32_t n = dolfinx::MPI::size(comm) * 10 - 1; + std::int32_t count = dolfinx::MPI::size(comm) * 10; + TestType norm = std::sqrt(n * (n + 1) * (2 * n + 1) / 6); + + auto mark = [=](auto e) { return e > theta * norm / std::sqrt(count); }; + + CHECK(std::ranges::count_if(indicators, mark) + == static_cast(indices.size())); + + for (std::int32_t i = 0; i < static_cast(indicators.size()); + ++i) + { + bool expect_marked = mark(indicators[i]); + bool marked = std::ranges::find(indices, i) != indices.end(); + CHECK(expect_marked == marked); + } + + // squared input + auto indicators_sq = indicators; + std::ranges::for_each(indicators_sq, [](auto& e) { e = std::pow(e, 2); }); + + CHECK(std::ranges::equal(indices, mark_equidistribution_squared( + comm, indicators_sq, theta))); +} diff --git a/python/dolfinx/mesh.py b/python/dolfinx/mesh.py index 7ef875fda3..e3ec70ddf1 100644 --- a/python/dolfinx/mesh.py +++ b/python/dolfinx/mesh.py @@ -35,7 +35,15 @@ from dolfinx.cpp.refinement import ( IdentityPartitionerPlaceholder, RefinementOption, - mark_maximum, +) +from dolfinx.cpp.refinement import ( + mark_equidistribution as _mark_equidistribution, +) +from dolfinx.cpp.refinement import ( + mark_equidistribution_squared as _mark_equidistribution_squared, +) +from dolfinx.cpp.refinement import ( + mark_maximum as _mark_maximum, ) from dolfinx.cpp.refinement import ( uniform_refine as _uniform_refine, @@ -75,6 +83,8 @@ "exterior_facet_indices", "locate_entities", "locate_entities_boundary", + "mark_equidistribution", + "mark_equidistribution_squared", "mark_maximum", "meshtags", "meshtags_from_entities", @@ -760,6 +770,83 @@ def refine( return Mesh(mesh1, ufl_domain), parent_cell, parent_facet +def mark_maximum( + comm: _MPI.Comm, + indicators: npt.NDArray[Real], + theta: float, +) -> npt.NDArray[np.int32]: + r"""Compute maximum-based marking of indicators. + + Returns the indices :math:`i` of the indicators :math:`\eta_i` that + satisfy the maximum threshold: + :math:`\eta_i > \theta \max_j \eta_j`. + + Args: + comm: Communicator to compute the maximum over. + indicators: Indicators (local) :math:`\eta_i` - usually an error + indicator associated with mesh entity :math:`i`. + theta: Parameter, :math:`0 < \theta < 1`. + + Returns: + Local indices of marked entities. + """ + return _mark_maximum(comm, indicators, theta) + + +def mark_equidistribution( + comm: _MPI.Comm, + indicators: npt.NDArray[Real], + theta: float, +) -> npt.NDArray[np.int32]: + r"""Compute equidistribution threshold marking of indicators. + + Returns the indices :math:`i` of the indicators :math:`\eta_i` that + satisfy the equidistribution threshold: + :math:`\eta_i > \theta \frac{\|\eta\|}{\sqrt{N}}` where + :math:`N` is the (global) number of indicators. + + Args: + comm: Communicator over which the global equidistribution + threshold is computed. + indicators: Indicators for owned local entities :math:`\eta_i` - + usually associated with mesh entity :math:`i`. Must contain owned + entities only; ghost values cause double-counting in the global + sum reductions. + theta: Parameter, :math:`0 < \theta < 1`. + + Returns: + Local indices of indicators that satisfy the threshold. + """ + return _mark_equidistribution(comm, indicators, theta) + + +def mark_equidistribution_squared( + comm: _MPI.Comm, + squared_indicators: npt.NDArray[Real], + theta: float, +) -> npt.NDArray[np.int32]: + r"""Compute equidistribution threshold marking of a squared indicator. + + Returns the indices :math:`i` of the squared indicators + :math:`\eta_i^2` that satisfy the equidistribution threshold: + :math:`\eta_i^2 > \theta^2 \frac{\|\eta\|^2}{N}` where + :math:`N` is the (global) number of indicators. + + Args: + comm: Communicator over which the global equidistribution + threshold is computed. + squared_indicators: Input squared indicators for owned local + entities :math:`\eta_i^2` - usually associated with mesh + entity :math:`i`. Must contain owned entities only; ghost + values cause double-counting in the global sum reductions. + theta: Parameter, :math:`0 < \theta < 1`. + + Returns: + Local indices of squared indicators that satisfy the threshold. + """ + return _mark_equidistribution_squared(comm, squared_indicators, theta) + + def create_mesh( comm: _MPI.Comm, cells: npt.NDArray[np.int64], diff --git a/python/dolfinx/wrappers/dolfinx_wrappers/refinement.h b/python/dolfinx/wrappers/dolfinx_wrappers/refinement.h index 7ba4fd8ea2..7eac92d1ee 100644 --- a/python/dolfinx/wrappers/dolfinx_wrappers/refinement.h +++ b/python/dolfinx/wrappers/dolfinx_wrappers/refinement.h @@ -127,13 +127,40 @@ void declare_refinement(nanobind::module_& m) m.def( "mark_maximum", - [](nb::ndarray, nb::c_contig> marker, T theta, - MPICommWrapper comm) + [](MPICommWrapper comm, + nb::ndarray, nb::c_contig> indicators, T theta) { return dolfinx_wrappers::as_nbarray(dolfinx::refinement::mark_maximum( - std::span(marker.data(), marker.size()), theta, comm.get())); + comm.get(), std::span(indicators.data(), indicators.size()), + theta)); }, - nb::arg("marker"), nb::arg("theta"), nb::arg("comm")); + nb::arg("comm"), nb::arg("indicators"), nb::arg("theta")); + + m.def( + "mark_equidistribution", + [](MPICommWrapper comm, + nb::ndarray, nb::c_contig> indicators, T theta) + { + return dolfinx_wrappers::as_nbarray( + dolfinx::refinement::mark_equidistribution( + comm.get(), std::span(indicators.data(), indicators.size()), + theta)); + }, + nb::arg("comm"), nb::arg("indicators"), nb::arg("theta")); + + m.def( + "mark_equidistribution_squared", + [](MPICommWrapper comm, + nb::ndarray, nb::c_contig> squared_indicators, + T theta) + { + return dolfinx_wrappers::as_nbarray( + dolfinx::refinement::mark_equidistribution_squared( + comm.get(), + std::span(squared_indicators.data(), squared_indicators.size()), + theta)); + }, + nb::arg("comm"), nb::arg("squared_indicators"), nb::arg("theta")); } } // namespace dolfinx_wrappers diff --git a/python/test/unit/refinement/test_mark.py b/python/test/unit/refinement/test_mark.py index 07a4ae6189..22c9b6cc0b 100644 --- a/python/test/unit/refinement/test_mark.py +++ b/python/test/unit/refinement/test_mark.py @@ -15,19 +15,48 @@ @pytest.mark.parametrize("theta", [0.2, 0.4, 0.6, 0.8]) @pytest.mark.parametrize("dtype", [np.float32, np.float64]) def test_mark_maximum(theta: float, dtype: np.dtype) -> None: - msh = mesh.create_unit_square(comm := MPI.COMM_WORLD, n := 10, n, dtype=dtype) + comm = MPI.COMM_WORLD + n = 10 + msh = mesh.create_unit_square(comm, n, n, dtype=dtype) tdim = msh.topology.dim - cell_count = (cell_im := msh.topology.index_map(tdim)).size_local + cell_im.num_ghosts - marker = np.random.default_rng(0).random(cell_count) + cell_count = msh.topology.index_map(tdim).size_local + indicators = np.random.default_rng(0).random(cell_count, dtype=dtype) - marked_cells = mesh.mark_maximum(marker, theta, comm) + marked_cells = mesh.mark_maximum(comm, indicators, theta) assert np.allclose( marked_cells, - np.argwhere(marker > theta * comm.allreduce(np.max(marker), MPI.MAX)).flatten(), + np.argwhere(indicators > theta * comm.allreduce(np.max(indicators), MPI.MAX)).flatten(), ) msh.topology.create_entities(1) marked_edges = mesh.compute_incident_entities(msh.topology, marked_cells, tdim, 1) mesh.refine(msh, marked_edges) + + +@pytest.mark.parametrize("theta", [0.2, 0.4, 0.6, 0.8]) +@pytest.mark.parametrize("dtype", [np.float32, np.float64]) +def test_mark_equidistribution(theta: float, dtype: np.dtype) -> None: + comm = MPI.COMM_WORLD + n = 10 + msh = mesh.create_unit_square(comm, n, n, dtype=dtype) + + tdim = msh.topology.dim + cell_count = msh.topology.index_map(tdim).size_local + indicators = np.random.default_rng(0).random(cell_count, dtype=dtype) + + marked_cells = mesh.mark_equidistribution(comm, indicators, theta) + + norm = np.sqrt(comm.allreduce(np.sum(indicators**2), MPI.SUM)) + count = comm.allreduce(indicators.size) + assert np.allclose( + marked_cells, + np.argwhere(indicators > theta * norm / np.sqrt(count)).flatten(), + ) + + msh.topology.create_entities(1) + marked_edges = mesh.compute_incident_entities(msh.topology, marked_cells, tdim, 1) + mesh.refine(msh, marked_edges) + + assert np.all(marked_cells == mesh.mark_equidistribution_squared(comm, indicators**2, theta))