From e84da5d456cd98329ccadbabcce4ba58c10834a1 Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Fri, 24 Apr 2026 22:33:24 +0200 Subject: [PATCH 01/45] Add equidistribution marking (c++) --- cpp/dolfinx/refinement/mark.h | 44 +++++++++++++++++++++++++++++++ cpp/test/mesh/refinement/mark.cpp | 41 ++++++++++++++++++++++++++++ 2 files changed, 85 insertions(+) diff --git a/cpp/dolfinx/refinement/mark.h b/cpp/dolfinx/refinement/mark.h index 395f54fcd2..d156bc2406 100644 --- a/cpp/dolfinx/refinement/mark.h +++ b/cpp/dolfinx/refinement/mark.h @@ -56,4 +56,48 @@ std::vector mark_maximum(std::span marker, T theta, return indices; } +/// @brief Equidistribution marking of a marker. +/// +/// @param[in] marker Input marker (local) - usually an error indicator per +/// entity +/// @param[in] theta Parameter, 0 < θ < 1 +/// @param[in] comm Communicator over which the total marker is computed. +/// @return Indices (local) of marker elements, which satisfy: marker_i ≥ θ +/// (sum_i marker_i^2)^1/2 / N^1/2. +template +std::vector mark_equidistribution(std::span marker, + T theta, MPI_Comm comm) +{ + if ((theta <= 0) || (theta >= 1)) + throw std::invalid_argument("Theta needs to fullfill 0 < θ < 1."); + + T norm{0}; + for (T e : marker) + norm += std::pow(e, 2); + + MPI_Allreduce(MPI_IN_PLACE, &norm, 1, dolfinx::MPI::mpi_t, MPI_SUM, comm); + + norm = std::sqrt(norm); + + std::int32_t count = marker.size(); + MPI_Allreduce(MPI_IN_PLACE, &count, 1, dolfinx::MPI::mpi_t, + MPI_SUM, comm); + + auto mark = [=](T e) { return e > theta * norm / std::sqrt(count); }; + + std::vector indices; + indices.reserve(std::ranges::count_if(marker, mark)); + + for (std::int32_t i = 0; i < static_cast(marker.size()); ++i) + { + if (mark(marker[i])) + indices.push_back(i); + } + + spdlog::info("Marking (equi) {} / {} (local) entities.", indices.size(), + marker.size()); + + return indices; +} + } // namespace dolfinx::refinement diff --git a/cpp/test/mesh/refinement/mark.cpp b/cpp/test/mesh/refinement/mark.cpp index b4bed6329d..c227f21032 100644 --- a/cpp/test/mesh/refinement/mark.cpp +++ b/cpp/test/mesh/refinement/mark.cpp @@ -51,3 +51,44 @@ TEMPLATE_TEST_CASE("Mark maximum", "[refinement][mark][maximum]", double, float) CHECK(expect_marked == marked); } } + +TEMPLATE_TEST_CASE("Mark equidistribution empty", + "[refinement][mark][equidistribution]", double, float) +{ + std::vector marker; + auto indices = mark_equidistribution(marker, .5, MPI_COMM_WORLD); + CHECK(indices.size() == 0); +} + +TEMPLATE_TEST_CASE("Mark equidistribution", "[refinement][mark][equidistribution]", + double, float) +{ + MPI_Comm comm = MPI_COMM_WORLD; + + std::vector marker; + marker.reserve(10); + for (std::size_t i = 0; i < 10; i++) + marker.push_back(10 * dolfinx::MPI::rank(comm) + i); + + TestType theta = 0.5; + auto indices = mark_equidistribution(marker, theta, comm); + + CHECK(std::ranges::all_of( + indices, [&](auto e) + { return (0 <= e) && (e <= static_cast(marker.size())); })); + + std::int32_t n = dolfinx::MPI::size(comm) * 10 - 1; + TestType norm = std::sqrt(n * (n + 1) * (2 * n + 1) / 6); + + auto mark = [=](auto e) { return e > theta * norm / std::sqrt(n); }; + + CHECK(std::ranges::count_if(marker, mark) + == static_cast(indices.size())); + + for (std::int32_t i = 0; i < static_cast(marker.size()); ++i) + { + bool expect_marked = mark(marker[i]); + bool marked = std::ranges::find(indices, i) != indices.end(); + CHECK(expect_marked == marked); + } +} From 6945e6f5f63622406428ca6e81ff1448bc0b8395 Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Mon, 27 Apr 2026 22:28:43 +0200 Subject: [PATCH 02/45] Pull out common routine mark_threshold --- cpp/dolfinx/refinement/mark.h | 51 +++++++++++++++++++------------ cpp/test/mesh/refinement/mark.cpp | 4 +-- 2 files changed, 33 insertions(+), 22 deletions(-) diff --git a/cpp/dolfinx/refinement/mark.h b/cpp/dolfinx/refinement/mark.h index d156bc2406..9f430c1c35 100644 --- a/cpp/dolfinx/refinement/mark.h +++ b/cpp/dolfinx/refinement/mark.h @@ -20,6 +20,34 @@ namespace dolfinx::refinement { +namespace impl +{ + +/// @brief Threshold marking helper +/// +/// @param[in] marker Input marker +/// @param[in] threshold Lower bound for values to mark +/// +/// @returns indices i which satisfy e_i > threshold. +template +std::vector mark_threshold(std::span marker, T threshold) +{ + auto mark = [=](T e) { return e > threshold; }; + + std::vector indices; + indices.reserve(std::ranges::count_if(marker, mark)); + + for (std::int32_t i = 0; i < static_cast(marker.size()); ++i) + { + if (mark(marker[i])) + indices.push_back(i); + } + + return indices; +} + +} // namespace impl + /// @brief Maximum marking of a marker. /// /// @param[in] marker Input marker (local) - usually an error indicator per @@ -39,16 +67,7 @@ std::vector mark_maximum(std::span marker, T theta, : 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; }; - - std::vector indices; - indices.reserve(std::ranges::count_if(marker, mark)); - - for (std::int32_t i = 0; i < static_cast(marker.size()); ++i) - { - if (mark(marker[i])) - indices.push_back(i); - } + auto indices = impl::mark_threshold(marker, theta * max); spdlog::info("Marking (max) {} / {} (local) entities.", indices.size(), marker.size()); @@ -83,16 +102,8 @@ std::vector mark_equidistribution(std::span marker, MPI_Allreduce(MPI_IN_PLACE, &count, 1, dolfinx::MPI::mpi_t, MPI_SUM, comm); - auto mark = [=](T e) { return e > theta * norm / std::sqrt(count); }; - - std::vector indices; - indices.reserve(std::ranges::count_if(marker, mark)); - - for (std::int32_t i = 0; i < static_cast(marker.size()); ++i) - { - if (mark(marker[i])) - indices.push_back(i); - } + auto indices + = impl::mark_threshold(marker, theta * norm / std::sqrt(count)); spdlog::info("Marking (equi) {} / {} (local) entities.", indices.size(), marker.size()); diff --git a/cpp/test/mesh/refinement/mark.cpp b/cpp/test/mesh/refinement/mark.cpp index c227f21032..4da5ff81d4 100644 --- a/cpp/test/mesh/refinement/mark.cpp +++ b/cpp/test/mesh/refinement/mark.cpp @@ -60,8 +60,8 @@ TEMPLATE_TEST_CASE("Mark equidistribution empty", CHECK(indices.size() == 0); } -TEMPLATE_TEST_CASE("Mark equidistribution", "[refinement][mark][equidistribution]", - double, float) +TEMPLATE_TEST_CASE("Mark equidistribution", + "[refinement][mark][equidistribution]", double, float) { MPI_Comm comm = MPI_COMM_WORLD; From c4a2bdb6e8292f94782f08cadfe912fbb1b1e63f Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Mon, 27 Apr 2026 22:38:28 +0200 Subject: [PATCH 03/45] Add equidistribution marking (py) --- python/dolfinx/mesh.py | 2 ++ .../wrappers/dolfinx_wrappers/refinement.h | 11 +++++++++ python/test/unit/refinement/test_mark.py | 23 +++++++++++++++++++ 3 files changed, 36 insertions(+) diff --git a/python/dolfinx/mesh.py b/python/dolfinx/mesh.py index 4cfaeb39c5..0e25c88f35 100644 --- a/python/dolfinx/mesh.py +++ b/python/dolfinx/mesh.py @@ -35,6 +35,7 @@ from dolfinx.cpp.refinement import ( IdentityPartitionerPlaceholder, RefinementOption, + mark_equidistribution, mark_maximum, ) from dolfinx.cpp.refinement import ( @@ -73,6 +74,7 @@ "exterior_facet_indices", "locate_entities", "locate_entities_boundary", + "mark_equidistribution", "mark_maximum", "meshtags", "meshtags_from_entities", diff --git a/python/dolfinx/wrappers/dolfinx_wrappers/refinement.h b/python/dolfinx/wrappers/dolfinx_wrappers/refinement.h index 7ba4fd8ea2..e19dc646ee 100644 --- a/python/dolfinx/wrappers/dolfinx_wrappers/refinement.h +++ b/python/dolfinx/wrappers/dolfinx_wrappers/refinement.h @@ -134,6 +134,17 @@ void declare_refinement(nanobind::module_& m) std::span(marker.data(), marker.size()), theta, comm.get())); }, nb::arg("marker"), nb::arg("theta"), nb::arg("comm")); + + m.def( + "mark_equidistribution", + [](nb::ndarray, nb::c_contig> marker, T theta, + MPICommWrapper comm) + { + return dolfinx_wrappers::as_nbarray( + dolfinx::refinement::mark_equidistribution( + std::span(marker.data(), marker.size()), theta, comm.get())); + }, + nb::arg("marker"), nb::arg("theta"), nb::arg("comm")); } } // namespace dolfinx_wrappers diff --git a/python/test/unit/refinement/test_mark.py b/python/test/unit/refinement/test_mark.py index 07a4ae6189..3a9dbb15bd 100644 --- a/python/test/unit/refinement/test_mark.py +++ b/python/test/unit/refinement/test_mark.py @@ -31,3 +31,26 @@ def test_mark_maximum(theta: float, dtype: np.dtype) -> None: 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: + msh = mesh.create_unit_square(comm := MPI.COMM_WORLD, n := 10, 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) + + marked_cells = mesh.mark_equidistribution(marker, theta, comm) + + norm = np.sqrt(comm.allreduce(np.sum(marker**2), MPI.SUM)) + count = comm.allreduce(marker.size) + assert np.allclose( + marked_cells, + np.argwhere(marker > 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) From e997738584847e4929685b21070f07634aacf85c Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Thu, 30 Apr 2026 15:22:14 +0200 Subject: [PATCH 04/45] Use inner prodcut --- cpp/dolfinx/refinement/mark.h | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/cpp/dolfinx/refinement/mark.h b/cpp/dolfinx/refinement/mark.h index 9f430c1c35..dfd55b0531 100644 --- a/cpp/dolfinx/refinement/mark.h +++ b/cpp/dolfinx/refinement/mark.h @@ -90,9 +90,8 @@ std::vector mark_equidistribution(std::span marker, if ((theta <= 0) || (theta >= 1)) throw std::invalid_argument("Theta needs to fullfill 0 < θ < 1."); - T norm{0}; - for (T e : marker) - norm += std::pow(e, 2); + auto norm + = std::inner_product(marker.begin(), marker.end(), marker.begin(), T{0}); MPI_Allreduce(MPI_IN_PLACE, &norm, 1, dolfinx::MPI::mpi_t, MPI_SUM, comm); From b9831ba556bfb2b88749155983a5aa27406efe51 Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Thu, 30 Apr 2026 15:26:23 +0200 Subject: [PATCH 05/45] Improve docstring --- cpp/dolfinx/refinement/mark.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cpp/dolfinx/refinement/mark.h b/cpp/dolfinx/refinement/mark.h index dfd55b0531..a802bfb929 100644 --- a/cpp/dolfinx/refinement/mark.h +++ b/cpp/dolfinx/refinement/mark.h @@ -81,7 +81,7 @@ std::vector mark_maximum(std::span marker, T theta, /// entity /// @param[in] theta Parameter, 0 < θ < 1 /// @param[in] comm Communicator over which the total marker is computed. -/// @return Indices (local) of marker elements, which satisfy: marker_i ≥ θ +/// @return Local indices of marked entities, which satisfy: marker_i ≥ θ /// (sum_i marker_i^2)^1/2 / N^1/2. template std::vector mark_equidistribution(std::span marker, From 6e9b9c473caf0ff784f34e5884707e334c44b140 Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Thu, 30 Apr 2026 15:31:02 +0200 Subject: [PATCH 06/45] use latex --- cpp/dolfinx/refinement/mark.h | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/cpp/dolfinx/refinement/mark.h b/cpp/dolfinx/refinement/mark.h index a802bfb929..e2befcb4dd 100644 --- a/cpp/dolfinx/refinement/mark.h +++ b/cpp/dolfinx/refinement/mark.h @@ -28,7 +28,7 @@ namespace impl /// @param[in] marker Input marker /// @param[in] threshold Lower bound for values to mark /// -/// @returns indices i which satisfy e_i > threshold. +/// @returns indices i which satisfy \f$ e_i > \text{threshold} \f$. template std::vector mark_threshold(std::span marker, T threshold) { @@ -52,10 +52,10 @@ std::vector mark_threshold(std::span marker, T threshold) /// /// @param[in] marker Input marker (local) - usually an error indicator per /// entity -/// @param[in] theta Cut off parameter, 0 < θ < 1 +/// @param[in] theta Cut off parameter, \f$ 0 < \theta < 1 \f$ /// @param[in] comm Communicator over which the maximum is computed. -/// @return Indices (local) of marker elements, which satisfy: marker_i ≥ θ -/// max(marker). +/// @return Indices (local) of marker elements, which satisfy: \f$ +/// \text{marker}_i \geq \theta \max \text{marker} \f$. template std::vector mark_maximum(std::span marker, T theta, MPI_Comm comm) @@ -79,10 +79,10 @@ std::vector mark_maximum(std::span marker, T theta, /// /// @param[in] marker Input marker (local) - usually an error indicator per /// entity -/// @param[in] theta Parameter, 0 < θ < 1 +/// @param[in] theta Parameter, \f$ 0 < \theta < 1 \f$ /// @param[in] comm Communicator over which the total marker is computed. -/// @return Local indices of marked entities, which satisfy: marker_i ≥ θ -/// (sum_i marker_i^2)^1/2 / N^1/2. +/// @return Local indices of marked entities, which satisfy: \f$ +/// \text{marker}_i \geq \theta \frac{|\text{marker}_i|}{\sqrt{N}} \f$. template std::vector mark_equidistribution(std::span marker, T theta, MPI_Comm comm) From de0b6a06e8c07e8e6f43bf1c565872288d7226ac Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Fri, 1 May 2026 10:47:32 +0200 Subject: [PATCH 07/45] Docs: differentiate between eta and marker consistently --- cpp/dolfinx/refinement/mark.h | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/cpp/dolfinx/refinement/mark.h b/cpp/dolfinx/refinement/mark.h index e2befcb4dd..c7f4f0872b 100644 --- a/cpp/dolfinx/refinement/mark.h +++ b/cpp/dolfinx/refinement/mark.h @@ -25,10 +25,10 @@ namespace impl /// @brief Threshold marking helper /// -/// @param[in] marker Input marker +/// @param[in] marker Input marker \f$ \eta \f$ /// @param[in] threshold Lower bound for values to mark /// -/// @returns indices i which satisfy \f$ e_i > \text{threshold} \f$. +/// @returns indices \f$ i \f$ which satisfy \f$ \eta_i > \text{threshold} \f$. template std::vector mark_threshold(std::span marker, T threshold) { @@ -50,12 +50,12 @@ std::vector mark_threshold(std::span marker, T threshold) /// @brief Maximum marking of a marker. /// -/// @param[in] marker Input marker (local) - usually an error indicator per -/// entity +/// @param[in] marker Input marker (local) \f$ \eta \f$ - usually an error +/// indicator per entity /// @param[in] theta Cut off parameter, \f$ 0 < \theta < 1 \f$ /// @param[in] comm Communicator over which the maximum is computed. /// @return Indices (local) of marker elements, which satisfy: \f$ -/// \text{marker}_i \geq \theta \max \text{marker} \f$. +/// \eta_i \geq \theta \max_j \eta_j \f$. template std::vector mark_maximum(std::span marker, T theta, MPI_Comm comm) @@ -77,12 +77,12 @@ std::vector mark_maximum(std::span marker, T theta, /// @brief Equidistribution marking of a marker. /// -/// @param[in] marker Input marker (local) - usually an error indicator per -/// entity +/// @param[in] marker Input marker (local) \f$ \eta \f$ - usually an error +/// indicator per entity /// @param[in] theta Parameter, \f$ 0 < \theta < 1 \f$ /// @param[in] comm Communicator over which the total marker is computed. /// @return Local indices of marked entities, which satisfy: \f$ -/// \text{marker}_i \geq \theta \frac{|\text{marker}_i|}{\sqrt{N}} \f$. +/// \eta_i \geq \theta \frac{|\eta|_2}{\sqrt{N}} \f$. template std::vector mark_equidistribution(std::span marker, T theta, MPI_Comm comm) From 77d35fe64715f2c11e010006a743174a8149b03f Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Thu, 7 May 2026 13:09:06 +0200 Subject: [PATCH 08/45] Add mark_equidistribution_squared --- cpp/dolfinx/refinement/mark.h | 34 +++++++++++++++++++++++++++++ cpp/test/mesh/refinement/mark.cpp | 36 +++++++++++++++++++++++++++++++ 2 files changed, 70 insertions(+) diff --git a/cpp/dolfinx/refinement/mark.h b/cpp/dolfinx/refinement/mark.h index c7f4f0872b..d83874fd34 100644 --- a/cpp/dolfinx/refinement/mark.h +++ b/cpp/dolfinx/refinement/mark.h @@ -110,4 +110,38 @@ std::vector mark_equidistribution(std::span marker, return indices; } +/// @brief Equidistribution marking of a marker with squared input. +/// +/// @param[in] marker Input marker (local) \f$ \eta^2 \f$ - usually an error +/// indicator per entity +/// @param[in] theta Parameter, \f$ 0 < \theta < 1 \f$ +/// @param[in] comm Communicator over which the total marker is computed. +/// @return Local indices of marked entities, which satisfy: \f$ +/// \eta_i \geq \theta \frac{|\eta|_2}{\sqrt{N}} \f$. +template +std::vector +mark_equidistribution_squared(std::span marker, T theta, MPI_Comm comm) +{ + if ((theta <= 0) || (theta >= 1)) + throw std::invalid_argument("Theta needs to fullfill 0 < θ < 1."); + + auto norm = std::accumulate(marker.begin(), marker.end(), T{0}); + + MPI_Allreduce(MPI_IN_PLACE, &norm, 1, dolfinx::MPI::mpi_t, MPI_SUM, comm); + + norm = std::sqrt(norm); + + std::int32_t count = marker.size(); + MPI_Allreduce(MPI_IN_PLACE, &count, 1, dolfinx::MPI::mpi_t, + MPI_SUM, comm); + + auto indices + = impl::mark_threshold(marker, theta * norm / std::sqrt(count)); + + spdlog::info("Marking (equi) {} / {} (local) entities.", indices.size(), + marker.size()); + + return indices; +} + } // namespace dolfinx::refinement diff --git a/cpp/test/mesh/refinement/mark.cpp b/cpp/test/mesh/refinement/mark.cpp index 4da5ff81d4..65cd159649 100644 --- a/cpp/test/mesh/refinement/mark.cpp +++ b/cpp/test/mesh/refinement/mark.cpp @@ -92,3 +92,39 @@ TEMPLATE_TEST_CASE("Mark equidistribution", CHECK(expect_marked == marked); } } + +TEMPLATE_TEST_CASE("Mark equidistribution squared", + "[refinement][mark][equidistributionsq]", double, float) +{ + MPI_Comm comm = MPI_COMM_WORLD; + + std::vector marker; + marker.reserve(10); + for (std::size_t i = 0; i < 10; i++) + marker.push_back(10 * dolfinx::MPI::rank(comm) + i); + + // square input + std::ranges::for_each(marker, [](auto& e) { e = std::pow(e, 2); }); + + TestType theta = 0.5; + auto indices = mark_equidistribution_squared(marker, theta, comm); + + CHECK(std::ranges::all_of( + indices, [&](auto e) + { return (0 <= e) && (e <= static_cast(marker.size())); })); + + std::int32_t n = dolfinx::MPI::size(comm) * 10 - 1; + TestType norm = std::sqrt(n * (n + 1) * (2 * n + 1) / 6); + + auto mark = [=](auto e) { return e > theta * norm / std::sqrt(n); }; + + CHECK(std::ranges::count_if(marker, mark) + == static_cast(indices.size())); + + for (std::int32_t i = 0; i < static_cast(marker.size()); ++i) + { + bool expect_marked = mark(marker[i]); + bool marked = std::ranges::find(indices, i) != indices.end(); + CHECK(expect_marked == marked); + } +} From 425a6420a8a4d5b410ccd4aeeb7bcd4fbbc2ad86 Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Thu, 7 May 2026 13:38:03 +0200 Subject: [PATCH 09/45] Fix and python exports --- cpp/dolfinx/refinement/mark.h | 10 ++++---- cpp/test/mesh/refinement/mark.cpp | 6 +++-- python/dolfinx/mesh.py | 2 ++ .../wrappers/dolfinx_wrappers/refinement.h | 11 +++++++++ python/test/unit/refinement/test_mark.py | 24 +++++++++++++++++++ 5 files changed, 45 insertions(+), 8 deletions(-) diff --git a/cpp/dolfinx/refinement/mark.h b/cpp/dolfinx/refinement/mark.h index d83874fd34..1a8435d799 100644 --- a/cpp/dolfinx/refinement/mark.h +++ b/cpp/dolfinx/refinement/mark.h @@ -117,7 +117,7 @@ std::vector mark_equidistribution(std::span marker, /// @param[in] theta Parameter, \f$ 0 < \theta < 1 \f$ /// @param[in] comm Communicator over which the total marker is computed. /// @return Local indices of marked entities, which satisfy: \f$ -/// \eta_i \geq \theta \frac{|\eta|_2}{\sqrt{N}} \f$. +/// \eta_i^2 \geq \theta^2 \frac{|\eta|_2^2}{N} \f$. template std::vector mark_equidistribution_squared(std::span marker, T theta, MPI_Comm comm) @@ -125,18 +125,16 @@ mark_equidistribution_squared(std::span marker, T theta, MPI_Comm comm) if ((theta <= 0) || (theta >= 1)) throw std::invalid_argument("Theta needs to fullfill 0 < θ < 1."); - auto norm = std::accumulate(marker.begin(), marker.end(), T{0}); + auto norm_sq = std::accumulate(marker.begin(), marker.end(), T{0}); - MPI_Allreduce(MPI_IN_PLACE, &norm, 1, dolfinx::MPI::mpi_t, MPI_SUM, comm); - - norm = std::sqrt(norm); + MPI_Allreduce(MPI_IN_PLACE, &norm_sq, 1, dolfinx::MPI::mpi_t, MPI_SUM, comm); std::int32_t count = marker.size(); MPI_Allreduce(MPI_IN_PLACE, &count, 1, dolfinx::MPI::mpi_t, MPI_SUM, comm); auto indices - = impl::mark_threshold(marker, theta * norm / std::sqrt(count)); + = impl::mark_threshold(marker, std::pow(theta, 2) * norm_sq / count); spdlog::info("Marking (equi) {} / {} (local) entities.", indices.size(), marker.size()); diff --git a/cpp/test/mesh/refinement/mark.cpp b/cpp/test/mesh/refinement/mark.cpp index 65cd159649..51ec2e78d9 100644 --- a/cpp/test/mesh/refinement/mark.cpp +++ b/cpp/test/mesh/refinement/mark.cpp @@ -104,10 +104,12 @@ TEMPLATE_TEST_CASE("Mark equidistribution squared", marker.push_back(10 * dolfinx::MPI::rank(comm) + i); // square input - std::ranges::for_each(marker, [](auto& e) { e = std::pow(e, 2); }); + auto marker_sq = marker; + std::ranges::for_each(marker_sq, [](auto& e) { e = std::pow(e, 2); }); TestType theta = 0.5; - auto indices = mark_equidistribution_squared(marker, theta, comm); + auto indices + = mark_equidistribution_squared(marker_sq, theta, comm); CHECK(std::ranges::all_of( indices, [&](auto e) diff --git a/python/dolfinx/mesh.py b/python/dolfinx/mesh.py index 0cb5250f35..23f7c65562 100644 --- a/python/dolfinx/mesh.py +++ b/python/dolfinx/mesh.py @@ -36,6 +36,7 @@ IdentityPartitionerPlaceholder, RefinementOption, mark_equidistribution, + mark_equidistribution_squared, mark_maximum, ) from dolfinx.cpp.refinement import ( @@ -77,6 +78,7 @@ "locate_entities", "locate_entities_boundary", "mark_equidistribution", + "mark_equidistribution_squared", "mark_maximum", "meshtags", "meshtags_from_entities", diff --git a/python/dolfinx/wrappers/dolfinx_wrappers/refinement.h b/python/dolfinx/wrappers/dolfinx_wrappers/refinement.h index e19dc646ee..11db1688d3 100644 --- a/python/dolfinx/wrappers/dolfinx_wrappers/refinement.h +++ b/python/dolfinx/wrappers/dolfinx_wrappers/refinement.h @@ -145,6 +145,17 @@ void declare_refinement(nanobind::module_& m) std::span(marker.data(), marker.size()), theta, comm.get())); }, nb::arg("marker"), nb::arg("theta"), nb::arg("comm")); + + m.def( + "mark_equidistribution_squared", + [](nb::ndarray, nb::c_contig> marker, T theta, + MPICommWrapper comm) + { + return dolfinx_wrappers::as_nbarray( + dolfinx::refinement::mark_equidistribution_squared( + std::span(marker.data(), marker.size()), theta, comm.get())); + }, + nb::arg("marker"), nb::arg("theta"), nb::arg("comm")); } } // namespace dolfinx_wrappers diff --git a/python/test/unit/refinement/test_mark.py b/python/test/unit/refinement/test_mark.py index 3a9dbb15bd..9b7484b6e7 100644 --- a/python/test/unit/refinement/test_mark.py +++ b/python/test/unit/refinement/test_mark.py @@ -54,3 +54,27 @@ def test_mark_equidistribution(theta: float, dtype: np.dtype) -> None: 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_squared(theta: float, dtype: np.dtype) -> None: + msh = mesh.create_unit_square(comm := MPI.COMM_WORLD, n := 10, 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) + + marker_sq = np.square(marker) + marked_cells = mesh.mark_equidistribution_squared(marker_sq, theta, comm) + + norm = np.sqrt(comm.allreduce(marker_sq.sum(), MPI.SUM)) + count = comm.allreduce(marker.size) + assert np.allclose( + marked_cells, + np.argwhere(marker > 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) From 63ace3fe3c55276d4be2067501d774f38b7e6a99 Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Thu, 7 May 2026 13:41:00 +0200 Subject: [PATCH 10/45] Drop code duplication --- cpp/test/mesh/refinement/mark.cpp | 37 ++---------------------- python/test/unit/refinement/test_mark.py | 24 +-------------- 2 files changed, 4 insertions(+), 57 deletions(-) diff --git a/cpp/test/mesh/refinement/mark.cpp b/cpp/test/mesh/refinement/mark.cpp index 51ec2e78d9..f8f0947676 100644 --- a/cpp/test/mesh/refinement/mark.cpp +++ b/cpp/test/mesh/refinement/mark.cpp @@ -91,42 +91,11 @@ TEMPLATE_TEST_CASE("Mark equidistribution", bool marked = std::ranges::find(indices, i) != indices.end(); CHECK(expect_marked == marked); } -} - -TEMPLATE_TEST_CASE("Mark equidistribution squared", - "[refinement][mark][equidistributionsq]", double, float) -{ - MPI_Comm comm = MPI_COMM_WORLD; - std::vector marker; - marker.reserve(10); - for (std::size_t i = 0; i < 10; i++) - marker.push_back(10 * dolfinx::MPI::rank(comm) + i); - - // square input + // squared input auto marker_sq = marker; std::ranges::for_each(marker_sq, [](auto& e) { e = std::pow(e, 2); }); - TestType theta = 0.5; - auto indices - = mark_equidistribution_squared(marker_sq, theta, comm); - - CHECK(std::ranges::all_of( - indices, [&](auto e) - { return (0 <= e) && (e <= static_cast(marker.size())); })); - - std::int32_t n = dolfinx::MPI::size(comm) * 10 - 1; - TestType norm = std::sqrt(n * (n + 1) * (2 * n + 1) / 6); - - auto mark = [=](auto e) { return e > theta * norm / std::sqrt(n); }; - - CHECK(std::ranges::count_if(marker, mark) - == static_cast(indices.size())); - - for (std::int32_t i = 0; i < static_cast(marker.size()); ++i) - { - bool expect_marked = mark(marker[i]); - bool marked = std::ranges::find(indices, i) != indices.end(); - CHECK(expect_marked == marked); - } + CHECK(std::ranges::equal(indices, mark_equidistribution_squared( + marker_sq, theta, comm))); } diff --git a/python/test/unit/refinement/test_mark.py b/python/test/unit/refinement/test_mark.py index 9b7484b6e7..760155224e 100644 --- a/python/test/unit/refinement/test_mark.py +++ b/python/test/unit/refinement/test_mark.py @@ -55,26 +55,4 @@ def test_mark_equidistribution(theta: float, dtype: np.dtype) -> None: 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_squared(theta: float, dtype: np.dtype) -> None: - msh = mesh.create_unit_square(comm := MPI.COMM_WORLD, n := 10, 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) - - marker_sq = np.square(marker) - marked_cells = mesh.mark_equidistribution_squared(marker_sq, theta, comm) - - norm = np.sqrt(comm.allreduce(marker_sq.sum(), MPI.SUM)) - count = comm.allreduce(marker.size) - assert np.allclose( - marked_cells, - np.argwhere(marker > 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(marker**2, theta, comm)) From 8d2ae2993f5d14f46ce5da3060500c377a0a6a54 Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Thu, 7 May 2026 13:42:14 +0200 Subject: [PATCH 11/45] format --- cpp/dolfinx/refinement/mark.h | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/cpp/dolfinx/refinement/mark.h b/cpp/dolfinx/refinement/mark.h index 1a8435d799..b10d1b43ec 100644 --- a/cpp/dolfinx/refinement/mark.h +++ b/cpp/dolfinx/refinement/mark.h @@ -127,7 +127,8 @@ mark_equidistribution_squared(std::span marker, T theta, MPI_Comm comm) auto norm_sq = std::accumulate(marker.begin(), marker.end(), T{0}); - MPI_Allreduce(MPI_IN_PLACE, &norm_sq, 1, dolfinx::MPI::mpi_t, MPI_SUM, comm); + MPI_Allreduce(MPI_IN_PLACE, &norm_sq, 1, dolfinx::MPI::mpi_t, MPI_SUM, + comm); std::int32_t count = marker.size(); MPI_Allreduce(MPI_IN_PLACE, &count, 1, dolfinx::MPI::mpi_t, From c588fb712b8b172ebc409c2735b88990e012cdef Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Thu, 7 May 2026 13:46:20 +0200 Subject: [PATCH 12/45] Rename eta to indicator --- cpp/dolfinx/refinement/mark.h | 40 +++++++++---------- .../wrappers/dolfinx_wrappers/refinement.h | 13 +++--- 2 files changed, 27 insertions(+), 26 deletions(-) diff --git a/cpp/dolfinx/refinement/mark.h b/cpp/dolfinx/refinement/mark.h index b10d1b43ec..82efd6c96f 100644 --- a/cpp/dolfinx/refinement/mark.h +++ b/cpp/dolfinx/refinement/mark.h @@ -25,21 +25,21 @@ namespace impl /// @brief Threshold marking helper /// -/// @param[in] marker Input marker \f$ \eta \f$ +/// @param[in] marker Indicator \f$ \eta \f$ /// @param[in] threshold Lower bound for values to mark /// /// @returns indices \f$ i \f$ which satisfy \f$ \eta_i > \text{threshold} \f$. template -std::vector mark_threshold(std::span marker, T threshold) +std::vector mark_threshold(std::span indicator, T threshold) { auto mark = [=](T e) { return e > threshold; }; std::vector indices; - indices.reserve(std::ranges::count_if(marker, mark)); + indices.reserve(std::ranges::count_if(indicator, mark)); - for (std::int32_t i = 0; i < static_cast(marker.size()); ++i) + for (std::int32_t i = 0; i < static_cast(indicator.size()); ++i) { - if (mark(marker[i])) + if (mark(indicator[i])) indices.push_back(i); } @@ -50,62 +50,62 @@ std::vector mark_threshold(std::span marker, T threshold) /// @brief Maximum marking of a marker. /// -/// @param[in] marker Input marker (local) \f$ \eta \f$ - usually an error +/// @param[in] marker Indicator (local) \f$ \eta \f$ - usually an error /// indicator per entity /// @param[in] theta Cut off parameter, \f$ 0 < \theta < 1 \f$ /// @param[in] comm Communicator over which the maximum is computed. /// @return Indices (local) of marker elements, which satisfy: \f$ /// \eta_i \geq \theta \max_j \eta_j \f$. template -std::vector mark_maximum(std::span marker, T theta, +std::vector mark_maximum(std::span indicator, T theta, MPI_Comm comm) { 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); + T max = indicator.empty() ? std::numeric_limits::lowest() + : std::ranges::max(indicator); MPI_Allreduce(MPI_IN_PLACE, &max, 1, dolfinx::MPI::mpi_t, MPI_MAX, comm); - auto indices = impl::mark_threshold(marker, theta * max); + auto indices = impl::mark_threshold(indicator, theta * max); spdlog::info("Marking (max) {} / {} (local) entities.", indices.size(), - marker.size()); + indicator.size()); return indices; } /// @brief Equidistribution marking of a marker. /// -/// @param[in] marker Input marker (local) \f$ \eta \f$ - usually an error +/// @param[in] marker Indicator (local) \f$ \eta \f$ - usually an error /// indicator per entity /// @param[in] theta Parameter, \f$ 0 < \theta < 1 \f$ /// @param[in] comm Communicator over which the total marker is computed. /// @return Local indices of marked entities, which satisfy: \f$ /// \eta_i \geq \theta \frac{|\eta|_2}{\sqrt{N}} \f$. template -std::vector mark_equidistribution(std::span marker, +std::vector mark_equidistribution(std::span indicator, T theta, MPI_Comm comm) { if ((theta <= 0) || (theta >= 1)) throw std::invalid_argument("Theta needs to fullfill 0 < θ < 1."); auto norm - = std::inner_product(marker.begin(), marker.end(), marker.begin(), T{0}); + = std::inner_product(indicator.begin(), indicator.end(), indicator.begin(), T{0}); MPI_Allreduce(MPI_IN_PLACE, &norm, 1, dolfinx::MPI::mpi_t, MPI_SUM, comm); norm = std::sqrt(norm); - std::int32_t count = marker.size(); + std::int32_t count = indicator.size(); MPI_Allreduce(MPI_IN_PLACE, &count, 1, dolfinx::MPI::mpi_t, MPI_SUM, comm); auto indices - = impl::mark_threshold(marker, theta * norm / std::sqrt(count)); + = impl::mark_threshold(indicator, theta * norm / std::sqrt(count)); spdlog::info("Marking (equi) {} / {} (local) entities.", indices.size(), - marker.size()); + indicator.size()); return indices; } @@ -125,9 +125,9 @@ mark_equidistribution_squared(std::span marker, T theta, MPI_Comm comm) if ((theta <= 0) || (theta >= 1)) throw std::invalid_argument("Theta needs to fullfill 0 < θ < 1."); - auto norm_sq = std::accumulate(marker.begin(), marker.end(), T{0}); + auto norm = std::accumulate(marker.begin(), marker.end(), T{0}); - MPI_Allreduce(MPI_IN_PLACE, &norm_sq, 1, dolfinx::MPI::mpi_t, MPI_SUM, + MPI_Allreduce(MPI_IN_PLACE, &norm, 1, dolfinx::MPI::mpi_t, MPI_SUM, comm); std::int32_t count = marker.size(); @@ -135,7 +135,7 @@ mark_equidistribution_squared(std::span marker, T theta, MPI_Comm comm) MPI_SUM, comm); auto indices - = impl::mark_threshold(marker, std::pow(theta, 2) * norm_sq / count); + = impl::mark_threshold(marker, std::pow(theta, 2) * norm / count); spdlog::info("Marking (equi) {} / {} (local) entities.", indices.size(), marker.size()); diff --git a/python/dolfinx/wrappers/dolfinx_wrappers/refinement.h b/python/dolfinx/wrappers/dolfinx_wrappers/refinement.h index 11db1688d3..0bab27ecb2 100644 --- a/python/dolfinx/wrappers/dolfinx_wrappers/refinement.h +++ b/python/dolfinx/wrappers/dolfinx_wrappers/refinement.h @@ -127,24 +127,25 @@ void declare_refinement(nanobind::module_& m) m.def( "mark_maximum", - [](nb::ndarray, nb::c_contig> marker, T theta, + [](nb::ndarray, nb::c_contig> indicator, T theta, MPICommWrapper comm) { return dolfinx_wrappers::as_nbarray(dolfinx::refinement::mark_maximum( - std::span(marker.data(), marker.size()), theta, comm.get())); + std::span(indicator.data(), indicator.size()), theta, comm.get())); }, - nb::arg("marker"), nb::arg("theta"), nb::arg("comm")); + nb::arg("indicator"), nb::arg("theta"), nb::arg("comm")); m.def( "mark_equidistribution", - [](nb::ndarray, nb::c_contig> marker, T theta, + [](nb::ndarray, nb::c_contig> indicator, T theta, MPICommWrapper comm) { return dolfinx_wrappers::as_nbarray( dolfinx::refinement::mark_equidistribution( - std::span(marker.data(), marker.size()), theta, comm.get())); + std::span(indicator.data(), indicator.size()), theta, + comm.get())); }, - nb::arg("marker"), nb::arg("theta"), nb::arg("comm")); + nb::arg("indicator"), nb::arg("theta"), nb::arg("comm")); m.def( "mark_equidistribution_squared", From f56fae32108fff879e6a2ccea3fa4379f6563140 Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Thu, 7 May 2026 13:47:47 +0200 Subject: [PATCH 13/45] format --- cpp/dolfinx/refinement/mark.h | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/cpp/dolfinx/refinement/mark.h b/cpp/dolfinx/refinement/mark.h index 82efd6c96f..121b62289f 100644 --- a/cpp/dolfinx/refinement/mark.h +++ b/cpp/dolfinx/refinement/mark.h @@ -30,7 +30,8 @@ namespace impl /// /// @returns indices \f$ i \f$ which satisfy \f$ \eta_i > \text{threshold} \f$. template -std::vector mark_threshold(std::span indicator, T threshold) +std::vector mark_threshold(std::span indicator, + T threshold) { auto mark = [=](T e) { return e > threshold; }; @@ -64,7 +65,7 @@ std::vector mark_maximum(std::span indicator, T theta, throw std::invalid_argument("Theta needs to fullfill 0 < θ < 1."); T max = indicator.empty() ? std::numeric_limits::lowest() - : std::ranges::max(indicator); + : std::ranges::max(indicator); MPI_Allreduce(MPI_IN_PLACE, &max, 1, dolfinx::MPI::mpi_t, MPI_MAX, comm); auto indices = impl::mark_threshold(indicator, theta * max); @@ -90,8 +91,8 @@ std::vector mark_equidistribution(std::span indicator, if ((theta <= 0) || (theta >= 1)) throw std::invalid_argument("Theta needs to fullfill 0 < θ < 1."); - auto norm - = std::inner_product(indicator.begin(), indicator.end(), indicator.begin(), T{0}); + auto norm = std::inner_product(indicator.begin(), indicator.end(), + indicator.begin(), T{0}); MPI_Allreduce(MPI_IN_PLACE, &norm, 1, dolfinx::MPI::mpi_t, MPI_SUM, comm); @@ -127,8 +128,7 @@ mark_equidistribution_squared(std::span marker, T theta, MPI_Comm comm) auto norm = std::accumulate(marker.begin(), marker.end(), T{0}); - MPI_Allreduce(MPI_IN_PLACE, &norm, 1, dolfinx::MPI::mpi_t, MPI_SUM, - comm); + MPI_Allreduce(MPI_IN_PLACE, &norm, 1, dolfinx::MPI::mpi_t, MPI_SUM, comm); std::int32_t count = marker.size(); MPI_Allreduce(MPI_IN_PLACE, &count, 1, dolfinx::MPI::mpi_t, From b768685ae743db65b9b2a7cfcf672a5aa28e09cb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Paul=20T=2E=20K=C3=BChner?= <56360279+schnellerhase@users.noreply.github.com> Date: Thu, 7 May 2026 13:55:29 +0200 Subject: [PATCH 14/45] Apply suggestions from code review MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Paul T. Kühner <56360279+schnellerhase@users.noreply.github.com> --- cpp/dolfinx/refinement/mark.h | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/cpp/dolfinx/refinement/mark.h b/cpp/dolfinx/refinement/mark.h index 121b62289f..1191869f6d 100644 --- a/cpp/dolfinx/refinement/mark.h +++ b/cpp/dolfinx/refinement/mark.h @@ -25,7 +25,7 @@ namespace impl /// @brief Threshold marking helper /// -/// @param[in] marker Indicator \f$ \eta \f$ +/// @param[in] indicator Indicator \f$ \eta \f$ /// @param[in] threshold Lower bound for values to mark /// /// @returns indices \f$ i \f$ which satisfy \f$ \eta_i > \text{threshold} \f$. @@ -51,7 +51,7 @@ std::vector mark_threshold(std::span indicator, /// @brief Maximum marking of a marker. /// -/// @param[in] marker Indicator (local) \f$ \eta \f$ - usually an error +/// @param[in] indicator Indicator (local) \f$ \eta \f$ - usually an error /// indicator per entity /// @param[in] theta Cut off parameter, \f$ 0 < \theta < 1 \f$ /// @param[in] comm Communicator over which the maximum is computed. @@ -78,7 +78,7 @@ std::vector mark_maximum(std::span indicator, T theta, /// @brief Equidistribution marking of a marker. /// -/// @param[in] marker Indicator (local) \f$ \eta \f$ - usually an error +/// @param[in] indicator Indicator (local) \f$ \eta \f$ - usually an error /// indicator per entity /// @param[in] theta Parameter, \f$ 0 < \theta < 1 \f$ /// @param[in] comm Communicator over which the total marker is computed. From d3b4a9b8dd2a8f523f1fef40d9ab616caed536fe Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Mon, 11 May 2026 08:41:57 +0200 Subject: [PATCH 15/45] Explicit capture --- cpp/dolfinx/refinement/mark.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cpp/dolfinx/refinement/mark.h b/cpp/dolfinx/refinement/mark.h index 1191869f6d..b41955abab 100644 --- a/cpp/dolfinx/refinement/mark.h +++ b/cpp/dolfinx/refinement/mark.h @@ -33,7 +33,7 @@ template std::vector mark_threshold(std::span indicator, T threshold) { - auto mark = [=](T e) { return e > threshold; }; + auto mark = [threshold](T e) { return e > threshold; }; std::vector indices; indices.reserve(std::ranges::count_if(indicator, mark)); From 73ea789acd29b876b41157713c14a6f35d283ec6 Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Mon, 11 May 2026 09:02:35 +0200 Subject: [PATCH 16/45] Improve docs --- cpp/dolfinx/refinement/mark.h | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/cpp/dolfinx/refinement/mark.h b/cpp/dolfinx/refinement/mark.h index b41955abab..64473d0767 100644 --- a/cpp/dolfinx/refinement/mark.h +++ b/cpp/dolfinx/refinement/mark.h @@ -23,12 +23,15 @@ namespace dolfinx::refinement namespace impl { -/// @brief Threshold marking helper +/// @brief Computes a local marking based on thresholding - helper for other +/// marking routines. /// -/// @param[in] indicator Indicator \f$ \eta \f$ +/// @param[in] indicator Entity wise indicator \f$ \eta \f$ to compute threshold +/// marking of. /// @param[in] threshold Lower bound for values to mark /// -/// @returns indices \f$ i \f$ which satisfy \f$ \eta_i > \text{threshold} \f$. +/// @returns list of indices \f$ i \f$ of indicator entries which satisfy +/// \f$ \eta_i > \text{threshold} \f$. template std::vector mark_threshold(std::span indicator, T threshold) @@ -49,11 +52,12 @@ std::vector mark_threshold(std::span indicator, } // namespace impl -/// @brief Maximum marking of a marker. +/// @brief Compute maximum marking of a marker, i.e. all elements which are +/// greater equal the markers maximum times a cut off parameter. /// -/// @param[in] indicator Indicator (local) \f$ \eta \f$ - usually an error -/// indicator per entity -/// @param[in] theta Cut off parameter, \f$ 0 < \theta < 1 \f$ +/// @param[in] indicator Indicator (local) \f$ \eta \f$ - usually an entity wise +/// error indicator. +/// @param[in] theta Cut off parameter, \f$ 0 < \theta < 1 \f$. /// @param[in] comm Communicator over which the maximum is computed. /// @return Indices (local) of marker elements, which satisfy: \f$ /// \eta_i \geq \theta \max_j \eta_j \f$. From 6b2cca205001df4f04372e7a43e86e2ba65b3159 Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Tue, 12 May 2026 13:01:17 +0200 Subject: [PATCH 17/45] Apply docstring improvements --- cpp/dolfinx/refinement/mark.h | 61 +++++++++++++++++++---------------- 1 file changed, 33 insertions(+), 28 deletions(-) diff --git a/cpp/dolfinx/refinement/mark.h b/cpp/dolfinx/refinement/mark.h index 64473d0767..d96ff9aac2 100644 --- a/cpp/dolfinx/refinement/mark.h +++ b/cpp/dolfinx/refinement/mark.h @@ -23,14 +23,15 @@ namespace dolfinx::refinement namespace impl { -/// @brief Computes a local marking based on thresholding - helper for other -/// marking routines. +/// @brief Computes local threshold-based marking. /// -/// @param[in] indicator Entity wise indicator \f$ \eta \f$ to compute threshold -/// marking of. -/// @param[in] threshold Lower bound for values to mark +/// Helper for other marking routines. /// -/// @returns list of indices \f$ i \f$ of indicator entries which satisfy +/// @param[in] indicator Entity-wise indicator \f$ \eta \f$ used for marking. +/// @param[in] threshold Threshold value; indicators greater than this are +/// marked. +/// +/// @returns list of indices \f$ i \f$ of indicator entries that satisfy /// \f$ \eta_i > \text{threshold} \f$. template std::vector mark_threshold(std::span indicator, @@ -52,14 +53,16 @@ std::vector mark_threshold(std::span indicator, } // namespace impl -/// @brief Compute maximum marking of a marker, i.e. all elements which are -/// greater equal the markers maximum times a cut off parameter. +/// @brief Computes maximum-based marking of an indicator. +/// +/// Marks all local entries whose indicators are greater than or equal to +/// \f$ \theta \f$ times the maximum indicator value across the communicator. /// /// @param[in] indicator Indicator (local) \f$ \eta \f$ - usually an entity wise /// error indicator. /// @param[in] theta Cut off parameter, \f$ 0 < \theta < 1 \f$. -/// @param[in] comm Communicator over which the maximum is computed. -/// @return Indices (local) of marker elements, which satisfy: \f$ +/// @param[in] comm Communicator to compute the maximum over. +/// @return Indices (local) of marker elements, that satisfy: \f$ /// \eta_i \geq \theta \max_j \eta_j \f$. template std::vector mark_maximum(std::span indicator, T theta, @@ -74,20 +77,21 @@ std::vector mark_maximum(std::span indicator, T theta, auto indices = impl::mark_threshold(indicator, theta * max); - spdlog::info("Marking (max) {} / {} (local) entities.", indices.size(), + spdlog::info("Marking (maximum) {} / {} (local) entities.", indices.size(), indicator.size()); return indices; } -/// @brief Equidistribution marking of a marker. +/// @brief Computes equidistribution threshold marking of an indicator. /// /// @param[in] indicator Indicator (local) \f$ \eta \f$ - usually an error -/// indicator per entity -/// @param[in] theta Parameter, \f$ 0 < \theta < 1 \f$ -/// @param[in] comm Communicator over which the total marker is computed. -/// @return Local indices of marked entities, which satisfy: \f$ -/// \eta_i \geq \theta \frac{|\eta|_2}{\sqrt{N}} \f$. +/// indicator per entity. +/// @param[in] theta Parameter, \f$ 0 < \theta < 1 \f$. +/// @param[in] comm Communicator over which the equidistribution threshold is +/// computed. +/// @return Local indices of marked entities, that satisfy: \f$ +/// \eta_i \geq \theta \frac{||\eta||}{\sqrt{N}} \f$. template std::vector mark_equidistribution(std::span indicator, T theta, MPI_Comm comm) @@ -109,20 +113,21 @@ std::vector mark_equidistribution(std::span indicator, auto indices = impl::mark_threshold(indicator, theta * norm / std::sqrt(count)); - spdlog::info("Marking (equi) {} / {} (local) entities.", indices.size(), - indicator.size()); + spdlog::info("Marking (equidistribution) {} / {} (local) entities.", + indices.size(), indicator.size()); return indices; } -/// @brief Equidistribution marking of a marker with squared input. +/// @brief Computes equidistribution threshold marking of a squared indicator. /// -/// @param[in] marker Input marker (local) \f$ \eta^2 \f$ - usually an error -/// indicator per entity -/// @param[in] theta Parameter, \f$ 0 < \theta < 1 \f$ -/// @param[in] comm Communicator over which the total marker is computed. -/// @return Local indices of marked entities, which satisfy: \f$ -/// \eta_i^2 \geq \theta^2 \frac{|\eta|_2^2}{N} \f$. +/// @param[in] marker Input marker (local) \f$ \eta^2 \f$ - usually a squared +/// error indicator per entity. +/// @param[in] theta Parameter, \f$ 0 < \theta < 1 \f$. +/// @param[in] comm Communicator over which the equidistribution threshold is +/// computed. +/// @return Local indices of marked entities, that satisfy: \f$ +/// \eta_i^2 \geq \theta^2 \frac{||\eta||^2}{N} \f$. template std::vector mark_equidistribution_squared(std::span marker, T theta, MPI_Comm comm) @@ -141,8 +146,8 @@ mark_equidistribution_squared(std::span marker, T theta, MPI_Comm comm) auto indices = impl::mark_threshold(marker, std::pow(theta, 2) * norm / count); - spdlog::info("Marking (equi) {} / {} (local) entities.", indices.size(), - marker.size()); + spdlog::info("Marking (equidistribution) {} / {} (local) entities.", + indices.size(), marker.size()); return indices; } From 625957b45a0f42478c6ca67abefa5a21542b65cc Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Tue, 12 May 2026 13:09:48 +0200 Subject: [PATCH 18/45] Remove walrus --- python/test/unit/refinement/test_mark.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/python/test/unit/refinement/test_mark.py b/python/test/unit/refinement/test_mark.py index 760155224e..4676845064 100644 --- a/python/test/unit/refinement/test_mark.py +++ b/python/test/unit/refinement/test_mark.py @@ -15,7 +15,9 @@ @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 @@ -36,7 +38,9 @@ def test_mark_maximum(theta: float, dtype: np.dtype) -> None: @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: - 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 From d9656db0a9f80cb10b578eda5b5ce2049e7cad0f Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Tue, 12 May 2026 13:13:59 +0200 Subject: [PATCH 19/45] Pass dtype to random --- python/test/unit/refinement/test_mark.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/python/test/unit/refinement/test_mark.py b/python/test/unit/refinement/test_mark.py index 4676845064..f6b1be725c 100644 --- a/python/test/unit/refinement/test_mark.py +++ b/python/test/unit/refinement/test_mark.py @@ -21,7 +21,7 @@ def test_mark_maximum(theta: float, dtype: np.dtype) -> None: 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) + marker = np.random.default_rng(0).random(cell_count, dtype=dtype) marked_cells = mesh.mark_maximum(marker, theta, comm) @@ -44,7 +44,7 @@ def test_mark_equidistribution(theta: float, dtype: np.dtype) -> None: 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) + marker = np.random.default_rng(0).random(cell_count, dtype=dtype) marked_cells = mesh.mark_equidistribution(marker, theta, comm) From cf0121f94cdd54e51e4286e192b044a99b2066e0 Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Tue, 12 May 2026 13:21:00 +0200 Subject: [PATCH 20/45] Comm first arg --- cpp/dolfinx/refinement/mark.h | 10 +++++----- cpp/test/mesh/refinement/mark.cpp | 10 +++++----- .../wrappers/dolfinx_wrappers/refinement.h | 20 +++++++++---------- python/test/unit/refinement/test_mark.py | 6 +++--- 4 files changed, 23 insertions(+), 23 deletions(-) diff --git a/cpp/dolfinx/refinement/mark.h b/cpp/dolfinx/refinement/mark.h index d96ff9aac2..98bdf9f4e5 100644 --- a/cpp/dolfinx/refinement/mark.h +++ b/cpp/dolfinx/refinement/mark.h @@ -65,8 +65,8 @@ std::vector mark_threshold(std::span indicator, /// @return Indices (local) of marker elements, that satisfy: \f$ /// \eta_i \geq \theta \max_j \eta_j \f$. template -std::vector mark_maximum(std::span indicator, T theta, - MPI_Comm comm) +std::vector mark_maximum(MPI_Comm comm, + std::span indicator, T theta) { if ((theta <= 0) || (theta >= 1)) throw std::invalid_argument("Theta needs to fullfill 0 < θ < 1."); @@ -93,8 +93,8 @@ std::vector mark_maximum(std::span indicator, T theta, /// @return Local indices of marked entities, that satisfy: \f$ /// \eta_i \geq \theta \frac{||\eta||}{\sqrt{N}} \f$. template -std::vector mark_equidistribution(std::span indicator, - T theta, MPI_Comm comm) +std::vector +mark_equidistribution(MPI_Comm comm, std::span indicator, T theta) { if ((theta <= 0) || (theta >= 1)) throw std::invalid_argument("Theta needs to fullfill 0 < θ < 1."); @@ -130,7 +130,7 @@ std::vector mark_equidistribution(std::span indicator, /// \eta_i^2 \geq \theta^2 \frac{||\eta||^2}{N} \f$. template std::vector -mark_equidistribution_squared(std::span marker, T theta, MPI_Comm comm) +mark_equidistribution_squared(MPI_Comm comm, std::span marker, T theta) { if ((theta <= 0) || (theta >= 1)) throw std::invalid_argument("Theta needs to fullfill 0 < θ < 1."); diff --git a/cpp/test/mesh/refinement/mark.cpp b/cpp/test/mesh/refinement/mark.cpp index f8f0947676..fdebe0184b 100644 --- a/cpp/test/mesh/refinement/mark.cpp +++ b/cpp/test/mesh/refinement/mark.cpp @@ -18,7 +18,7 @@ TEMPLATE_TEST_CASE("Mark maximum empty", "[refinement][mark][maximum]", double, float) { std::vector marker; - auto indices = mark_maximum(marker, .5, MPI_COMM_WORLD); + auto indices = mark_maximum(MPI_COMM_WORLD, marker, .5); CHECK(indices.size() == 0); } @@ -32,7 +32,7 @@ TEMPLATE_TEST_CASE("Mark maximum", "[refinement][mark][maximum]", double, float) marker.push_back(10 * dolfinx::MPI::rank(comm) + i); TestType theta = 0.5; - auto indices = mark_maximum(marker, theta, comm); + auto indices = mark_maximum(comm, marker, theta); CHECK(std::ranges::all_of( indices, [&](auto e) @@ -56,7 +56,7 @@ TEMPLATE_TEST_CASE("Mark equidistribution empty", "[refinement][mark][equidistribution]", double, float) { std::vector marker; - auto indices = mark_equidistribution(marker, .5, MPI_COMM_WORLD); + auto indices = mark_equidistribution(MPI_COMM_WORLD, marker, .5); CHECK(indices.size() == 0); } @@ -71,7 +71,7 @@ TEMPLATE_TEST_CASE("Mark equidistribution", marker.push_back(10 * dolfinx::MPI::rank(comm) + i); TestType theta = 0.5; - auto indices = mark_equidistribution(marker, theta, comm); + auto indices = mark_equidistribution(comm, marker, theta); CHECK(std::ranges::all_of( indices, [&](auto e) @@ -97,5 +97,5 @@ TEMPLATE_TEST_CASE("Mark equidistribution", std::ranges::for_each(marker_sq, [](auto& e) { e = std::pow(e, 2); }); CHECK(std::ranges::equal(indices, mark_equidistribution_squared( - marker_sq, theta, comm))); + comm, marker_sq, theta))); } diff --git a/python/dolfinx/wrappers/dolfinx_wrappers/refinement.h b/python/dolfinx/wrappers/dolfinx_wrappers/refinement.h index 0bab27ecb2..b05436aff4 100644 --- a/python/dolfinx/wrappers/dolfinx_wrappers/refinement.h +++ b/python/dolfinx/wrappers/dolfinx_wrappers/refinement.h @@ -127,34 +127,34 @@ void declare_refinement(nanobind::module_& m) m.def( "mark_maximum", - [](nb::ndarray, nb::c_contig> indicator, T theta, - MPICommWrapper comm) + [](MPICommWrapper comm, + nb::ndarray, nb::c_contig> indicator, T theta) { return dolfinx_wrappers::as_nbarray(dolfinx::refinement::mark_maximum( - std::span(indicator.data(), indicator.size()), theta, comm.get())); + comm.get(), std::span(indicator.data(), indicator.size()), theta)); }, nb::arg("indicator"), nb::arg("theta"), nb::arg("comm")); m.def( "mark_equidistribution", - [](nb::ndarray, nb::c_contig> indicator, T theta, - MPICommWrapper comm) + [](MPICommWrapper comm, + nb::ndarray, nb::c_contig> indicator, T theta) { return dolfinx_wrappers::as_nbarray( dolfinx::refinement::mark_equidistribution( - std::span(indicator.data(), indicator.size()), theta, - comm.get())); + comm.get(), std::span(indicator.data(), indicator.size()), + theta)); }, nb::arg("indicator"), nb::arg("theta"), nb::arg("comm")); m.def( "mark_equidistribution_squared", - [](nb::ndarray, nb::c_contig> marker, T theta, - MPICommWrapper comm) + [](MPICommWrapper comm, + nb::ndarray, nb::c_contig> marker, T theta) { return dolfinx_wrappers::as_nbarray( dolfinx::refinement::mark_equidistribution_squared( - std::span(marker.data(), marker.size()), theta, comm.get())); + comm.get(), std::span(marker.data(), marker.size()), theta)); }, nb::arg("marker"), nb::arg("theta"), nb::arg("comm")); } diff --git a/python/test/unit/refinement/test_mark.py b/python/test/unit/refinement/test_mark.py index f6b1be725c..e3f5f5740d 100644 --- a/python/test/unit/refinement/test_mark.py +++ b/python/test/unit/refinement/test_mark.py @@ -23,7 +23,7 @@ def test_mark_maximum(theta: float, dtype: np.dtype) -> None: cell_count = (cell_im := msh.topology.index_map(tdim)).size_local + cell_im.num_ghosts marker = np.random.default_rng(0).random(cell_count, dtype=dtype) - marked_cells = mesh.mark_maximum(marker, theta, comm) + marked_cells = mesh.mark_maximum(comm, marker, theta) assert np.allclose( marked_cells, @@ -46,7 +46,7 @@ def test_mark_equidistribution(theta: float, dtype: np.dtype) -> None: cell_count = (cell_im := msh.topology.index_map(tdim)).size_local + cell_im.num_ghosts marker = np.random.default_rng(0).random(cell_count, dtype=dtype) - marked_cells = mesh.mark_equidistribution(marker, theta, comm) + marked_cells = mesh.mark_equidistribution(comm, marker, theta) norm = np.sqrt(comm.allreduce(np.sum(marker**2), MPI.SUM)) count = comm.allreduce(marker.size) @@ -59,4 +59,4 @@ def test_mark_equidistribution(theta: float, dtype: np.dtype) -> None: 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(marker**2, theta, comm)) + assert np.all(marked_cells == mesh.mark_equidistribution_squared(comm, marker**2, theta)) From a65c3c1a211ac6d0c99273231a782ef5fcd5b892 Mon Sep 17 00:00:00 2001 From: "Jack S. Hale" Date: Tue, 19 May 2026 19:08:24 +0200 Subject: [PATCH 21/45] WIP - Docstring adjustments --- cpp/dolfinx/refinement/mark.h | 78 +++++++++++++++++++---------------- 1 file changed, 42 insertions(+), 36 deletions(-) diff --git a/cpp/dolfinx/refinement/mark.h b/cpp/dolfinx/refinement/mark.h index 98bdf9f4e5..af34bc89d0 100644 --- a/cpp/dolfinx/refinement/mark.h +++ b/cpp/dolfinx/refinement/mark.h @@ -27,24 +27,24 @@ namespace impl /// /// Helper for other marking routines. /// -/// @param[in] indicator Entity-wise indicator \f$ \eta \f$ used for marking. +/// @param[in] indicators Entity-wise indicator \f$ \eta \f$ used for marking. /// @param[in] threshold Threshold value; indicators greater than this are /// marked. /// /// @returns list of indices \f$ i \f$ of indicator entries that satisfy /// \f$ \eta_i > \text{threshold} \f$. template -std::vector mark_threshold(std::span indicator, +std::vector mark_threshold(std::span indicators, T threshold) { auto mark = [threshold](T e) { return e > threshold; }; std::vector indices; - indices.reserve(std::ranges::count_if(indicator, mark)); + indices.reserve(std::ranges::count_if(indicators, mark)); - for (std::int32_t i = 0; i < static_cast(indicator.size()); ++i) + for (std::int32_t i = 0; i < static_cast(indicators.size()); ++i) { - if (mark(indicator[i])) + if (mark(indicators[i])) indices.push_back(i); } @@ -55,18 +55,18 @@ std::vector mark_threshold(std::span indicator, /// @brief Computes maximum-based marking of an indicator. /// -/// Marks all local entries whose indicators are greater than or equal to -/// \f$ \theta \f$ times the maximum indicator value across the communicator. +/// Returns the indices \f$ i \f$ of the indicators $eta_i$ that satisfy the +/// maximum threshold: \f$ \eta_i \geq \theta \max_j \eta_j \f$. /// -/// @param[in] indicator Indicator (local) \f$ \eta \f$ - usually an entity wise -/// error indicator. -/// @param[in] theta Cut off parameter, \f$ 0 < \theta < 1 \f$. /// @param[in] comm Communicator to compute the maximum over. -/// @return Indices (local) of marker elements, that satisfy: \f$ -/// \eta_i \geq \theta \max_j \eta_j \f$. +/// @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 Indices (local) of marker elements, that satisfy the maximum +/// threshold: template std::vector mark_maximum(MPI_Comm comm, - std::span indicator, T theta) + std::span indicators, T theta) { if ((theta <= 0) || (theta >= 1)) throw std::invalid_argument("Theta needs to fullfill 0 < θ < 1."); @@ -85,57 +85,63 @@ std::vector mark_maximum(MPI_Comm comm, /// @brief Computes equidistribution threshold marking of an indicator. /// -/// @param[in] indicator Indicator (local) \f$ \eta \f$ - usually an error -/// indicator per entity. +/// Returns the indices \f$i\f$ of the indicators $eta_i$ that satisfy the +/// equidistribution threshold: \f$\eta_i > \theta \frac{||\eta||}{\sqrt{N}} \f$ +/// where \f$ N \f$ is the (global) number of indicators. +/// +/// @param[in] comm Communicator over which the global equidistribution +/// threshold is computed. +/// @param[in] indicators Indicators (local) \f$ \eta_i \f$ - usually +/// associated with mesh entity \f$ i \f$. /// @param[in] theta Parameter, \f$ 0 < \theta < 1 \f$. -/// @param[in] comm Communicator over which the equidistribution threshold is -/// computed. -/// @return Local indices of marked entities, that satisfy: \f$ -/// \eta_i \geq \theta \frac{||\eta||}{\sqrt{N}} \f$. +/// @return Local indices of marked entities. template std::vector -mark_equidistribution(MPI_Comm comm, std::span indicator, T theta) +mark_equidistribution(MPI_Comm comm, std::span indicators, T theta) { if ((theta <= 0) || (theta >= 1)) throw std::invalid_argument("Theta needs to fullfill 0 < θ < 1."); - auto norm = std::inner_product(indicator.begin(), indicator.end(), - indicator.begin(), T{0}); + auto 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); norm = std::sqrt(norm); - std::int32_t count = indicator.size(); + std::int32_t count = indicators.size(); MPI_Allreduce(MPI_IN_PLACE, &count, 1, dolfinx::MPI::mpi_t, MPI_SUM, comm); auto indices - = impl::mark_threshold(indicator, theta * norm / std::sqrt(count)); + = impl::mark_threshold(indicators, theta * norm / std::sqrt(count)); spdlog::info("Marking (equidistribution) {} / {} (local) entities.", - indices.size(), indicator.size()); + indices.size(), indicators.size()); return indices; } /// @brief Computes equidistribution threshold marking of a squared indicator. /// -/// @param[in] marker Input marker (local) \f$ \eta^2 \f$ - usually a squared -/// error indicator per entity. -/// @param[in] theta Parameter, \f$ 0 < \theta < 1 \f$. -/// @param[in] comm Communicator over which the equidistribution threshold is +/// Returns the indices \f$i\f$ of the squared indicators $eta_i^2$ 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. +/// +/// @param[in] comm Communicator over which the global equidistribution threshold is /// computed. -/// @return Local indices of marked entities, that satisfy: \f$ -/// \eta_i^2 \geq \theta^2 \frac{||\eta||^2}{N} \f$. +/// @param[in] indicators Input indicators (local) \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 marked entities. template std::vector -mark_equidistribution_squared(MPI_Comm comm, std::span marker, T theta) +mark_equidistribution_squared(MPI_Comm comm, std::span indicators, T theta) { if ((theta <= 0) || (theta >= 1)) throw std::invalid_argument("Theta needs to fullfill 0 < θ < 1."); - auto norm = std::accumulate(marker.begin(), marker.end(), T{0}); + auto norm = std::accumulate(indicators.begin(), indicators.end(), T{0}); MPI_Allreduce(MPI_IN_PLACE, &norm, 1, dolfinx::MPI::mpi_t, MPI_SUM, comm); @@ -144,10 +150,10 @@ mark_equidistribution_squared(MPI_Comm comm, std::span marker, T theta) MPI_SUM, comm); auto indices - = impl::mark_threshold(marker, std::pow(theta, 2) * norm / count); + = impl::mark_threshold(indicators, std::pow(theta, 2) * norm / count); - spdlog::info("Marking (equidistribution) {} / {} (local) entities.", - indices.size(), marker.size()); + spdlog::info("Marking (equidistribution) {} of {} local entities.", + indices.size(), indicators.size()); return indices; } From 08bb59658ccf115b3c99dd61421cd8f38b777394 Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Tue, 19 May 2026 21:34:37 +0200 Subject: [PATCH 22/45] fixup docs --- cpp/dolfinx/refinement/mark.h | 32 +++++++++++++++++--------------- 1 file changed, 17 insertions(+), 15 deletions(-) diff --git a/cpp/dolfinx/refinement/mark.h b/cpp/dolfinx/refinement/mark.h index af34bc89d0..a34dd86f15 100644 --- a/cpp/dolfinx/refinement/mark.h +++ b/cpp/dolfinx/refinement/mark.h @@ -42,7 +42,8 @@ std::vector mark_threshold(std::span indicators, std::vector indices; indices.reserve(std::ranges::count_if(indicators, mark)); - for (std::int32_t i = 0; i < static_cast(indicators.size()); ++i) + for (std::int32_t i = 0; i < static_cast(indicators.size()); + ++i) { if (mark(indicators[i])) indices.push_back(i); @@ -63,7 +64,7 @@ std::vector mark_threshold(std::span indicators, /// usually an error indicator associated with mesh entity \f$ i \f$. /// @param[in] theta Parameter, \f$ 0 < \theta < 1 \f$. /// @return Indices (local) of marker elements, that satisfy the maximum -/// threshold: +/// threshold: template std::vector mark_maximum(MPI_Comm comm, std::span indicators, T theta) @@ -71,14 +72,14 @@ std::vector mark_maximum(MPI_Comm comm, if ((theta <= 0) || (theta >= 1)) throw std::invalid_argument("Theta needs to fullfill 0 < θ < 1."); - T max = indicator.empty() ? std::numeric_limits::lowest() - : std::ranges::max(indicator); + 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); - auto indices = impl::mark_threshold(indicator, theta * max); + auto indices = impl::mark_threshold(indicators, theta * max); spdlog::info("Marking (maximum) {} / {} (local) entities.", indices.size(), - indicator.size()); + indicators.size()); return indices; } @@ -89,9 +90,9 @@ std::vector mark_maximum(MPI_Comm comm, /// equidistribution threshold: \f$\eta_i > \theta \frac{||\eta||}{\sqrt{N}} \f$ /// where \f$ N \f$ is the (global) number of indicators. /// -/// @param[in] comm Communicator over which the global equidistribution +/// @param[in] comm Communicator over which the global equidistribution /// threshold is computed. -/// @param[in] indicators Indicators (local) \f$ \eta_i \f$ - usually +/// @param[in] indicators Indicators (local) \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 marked entities. @@ -124,19 +125,20 @@ mark_equidistribution(MPI_Comm comm, std::span indicators, T theta) /// @brief Computes equidistribution threshold marking of a squared indicator. /// -/// Returns the indices \f$i\f$ of the squared indicators $eta_i^2$ 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. +/// Returns the indices \f$i\f$ of the squared indicators $eta_i^2$ 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. /// -/// @param[in] comm Communicator over which the global equidistribution threshold is -/// computed. +/// @param[in] comm Communicator over which the global equidistribution +/// threshold is computed. /// @param[in] indicators Input indicators (local) \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 marked entities. template std::vector -mark_equidistribution_squared(MPI_Comm comm, std::span indicators, T theta) +mark_equidistribution_squared(MPI_Comm comm, std::span indicators, + T theta) { if ((theta <= 0) || (theta >= 1)) throw std::invalid_argument("Theta needs to fullfill 0 < θ < 1."); @@ -145,7 +147,7 @@ mark_equidistribution_squared(MPI_Comm comm, std::span indicators, T th MPI_Allreduce(MPI_IN_PLACE, &norm, 1, dolfinx::MPI::mpi_t, MPI_SUM, comm); - std::int32_t count = marker.size(); + std::int32_t count = indicators.size(); MPI_Allreduce(MPI_IN_PLACE, &count, 1, dolfinx::MPI::mpi_t, MPI_SUM, comm); From 6f3733c825e102a8d1a033dd78fb0fe5d7a3219f Mon Sep 17 00:00:00 2001 From: "Jack S. Hale" Date: Wed, 20 May 2026 10:30:33 +0200 Subject: [PATCH 23/45] Fixup docs and remove potential overflow on summing std::int32_t --- cpp/dolfinx/refinement/mark.h | 58 ++++++++++++++++++----------------- 1 file changed, 30 insertions(+), 28 deletions(-) diff --git a/cpp/dolfinx/refinement/mark.h b/cpp/dolfinx/refinement/mark.h index a34dd86f15..f0b2354c88 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) // @@ -23,16 +23,17 @@ namespace dolfinx::refinement namespace impl { -/// @brief Computes local threshold-based marking. +/// @brief Computes local threshold marking of indicators. /// /// Helper for other marking routines. /// -/// @param[in] indicators Entity-wise indicator \f$ \eta \f$ used for marking. +/// Returns the indices \f$ i \f$ of the indicators $eta_i$ that satisfy the +/// threshold: \f$ \eta_i > \text{threshold} \f$. +/// +/// @param[in] indicators Indicators \f$ \eta_i \f$ used for marking. /// @param[in] threshold Threshold value; indicators greater than this are /// marked. -/// -/// @returns list of indices \f$ i \f$ of indicator entries that satisfy -/// \f$ \eta_i > \text{threshold} \f$. +/// @return Local indices of marked entities. template std::vector mark_threshold(std::span indicators, T threshold) @@ -54,7 +55,7 @@ std::vector mark_threshold(std::span indicators, } // namespace impl -/// @brief Computes maximum-based marking of an indicator. +/// @brief Computes maximum-based marking of indicators. /// /// Returns the indices \f$ i \f$ of the indicators $eta_i$ that satisfy the /// maximum threshold: \f$ \eta_i \geq \theta \max_j \eta_j \f$. @@ -63,8 +64,7 @@ std::vector mark_threshold(std::span indicators, /// @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 Indices (local) of marker elements, that satisfy the maximum -/// threshold: +/// @return Local indices of marked entities. template std::vector mark_maximum(MPI_Comm comm, std::span indicators, T theta) @@ -76,7 +76,8 @@ std::vector mark_maximum(MPI_Comm comm, : std::ranges::max(indicators); MPI_Allreduce(MPI_IN_PLACE, &max, 1, dolfinx::MPI::mpi_t, MPI_MAX, comm); - auto indices = impl::mark_threshold(indicators, theta * max); + std::vector indices + = impl::mark_threshold(indicators, theta * max); spdlog::info("Marking (maximum) {} / {} (local) entities.", indices.size(), indicators.size()); @@ -84,38 +85,39 @@ std::vector mark_maximum(MPI_Comm comm, return indices; } -/// @brief Computes equidistribution threshold marking of an indicator. +/// @brief Computes equidistribution threshold marking of indicators. /// -/// Returns the indices \f$i\f$ of the indicators $eta_i$ that satisfy the -/// equidistribution threshold: \f$\eta_i > \theta \frac{||\eta||}{\sqrt{N}} \f$ -/// where \f$ N \f$ is the (global) number 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. /// /// @param[in] comm Communicator over which the global equidistribution /// threshold is computed. /// @param[in] indicators Indicators (local) \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 marked entities. +/// @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 needs to fullfill 0 < θ < 1."); + throw std::invalid_argument("Theta must fullfill 0 < θ < 1."); - auto norm = std::inner_product(indicators.begin(), indicators.end(), - indicators.begin(), T{0}); + 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); - norm = std::sqrt(norm); + T sqrt_norm = std::sqrt(norm); std::int32_t count = indicators.size(); - MPI_Allreduce(MPI_IN_PLACE, &count, 1, dolfinx::MPI::mpi_t, + MPI_Allreduce(MPI_IN_PLACE, &count, 1, dolfinx::MPI::mpi_t, MPI_SUM, comm); - auto indices - = impl::mark_threshold(indicators, theta * norm / std::sqrt(count)); + 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()); @@ -134,25 +136,25 @@ mark_equidistribution(MPI_Comm comm, std::span indicators, T theta) /// @param[in] indicators Input indicators (local) \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 marked entities. +/// @return Local indices of indicators that satisfy the threshold. template std::vector mark_equidistribution_squared(MPI_Comm comm, std::span indicators, T theta) { if ((theta <= 0) || (theta >= 1)) - throw std::invalid_argument("Theta needs to fullfill 0 < θ < 1."); + throw std::invalid_argument("Theta must fullfill 0 < θ < 1."); - auto norm = std::accumulate(indicators.begin(), indicators.end(), T{0}); + T norm = std::accumulate(indicators.begin(), indicators.end(), T{0}); MPI_Allreduce(MPI_IN_PLACE, &norm, 1, dolfinx::MPI::mpi_t, MPI_SUM, comm); std::int32_t count = indicators.size(); - MPI_Allreduce(MPI_IN_PLACE, &count, 1, dolfinx::MPI::mpi_t, + MPI_Allreduce(MPI_IN_PLACE, &count, 1, dolfinx::MPI::mpi_t, MPI_SUM, comm); - auto indices - = impl::mark_threshold(indicators, std::pow(theta, 2) * norm / count); + std::vector indices = impl::mark_threshold( + indicators, std::pow(theta, 2) * norm / static_cast(count)); spdlog::info("Marking (equidistribution) {} of {} local entities.", indices.size(), indicators.size()); From 477c05c3a03f53231b691912d9c31ba9a9fb1d2a Mon Sep 17 00:00:00 2001 From: "Jack S. Hale" Date: Wed, 20 May 2026 10:32:13 +0200 Subject: [PATCH 24/45] size_t --- cpp/dolfinx/refinement/mark.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/cpp/dolfinx/refinement/mark.h b/cpp/dolfinx/refinement/mark.h index f0b2354c88..93766619f5 100644 --- a/cpp/dolfinx/refinement/mark.h +++ b/cpp/dolfinx/refinement/mark.h @@ -112,7 +112,7 @@ mark_equidistribution(MPI_Comm comm, std::span indicators, T theta) T sqrt_norm = std::sqrt(norm); - std::int32_t count = indicators.size(); + std::size_t count = indicators.size(); MPI_Allreduce(MPI_IN_PLACE, &count, 1, dolfinx::MPI::mpi_t, MPI_SUM, comm); @@ -149,7 +149,7 @@ mark_equidistribution_squared(MPI_Comm comm, std::span indicators, MPI_Allreduce(MPI_IN_PLACE, &norm, 1, dolfinx::MPI::mpi_t, MPI_SUM, comm); - std::int32_t count = indicators.size(); + std::size_t count = indicators.size(); MPI_Allreduce(MPI_IN_PLACE, &count, 1, dolfinx::MPI::mpi_t, MPI_SUM, comm); From a498c3b311d73a97ad1b0e3dde0b0f4a4cfc2feb Mon Sep 17 00:00:00 2001 From: "Jack S. Hale" Date: Wed, 20 May 2026 10:42:21 +0200 Subject: [PATCH 25/45] This is optimal unless the expected mark fraction is very small (< 5%). --- cpp/dolfinx/refinement/mark.h | 13 +++++-------- 1 file changed, 5 insertions(+), 8 deletions(-) diff --git a/cpp/dolfinx/refinement/mark.h b/cpp/dolfinx/refinement/mark.h index 93766619f5..e0a16a1cbe 100644 --- a/cpp/dolfinx/refinement/mark.h +++ b/cpp/dolfinx/refinement/mark.h @@ -38,15 +38,12 @@ template std::vector mark_threshold(std::span indicators, T threshold) { - auto mark = [threshold](T e) { return e > threshold; }; - std::vector indices; - indices.reserve(std::ranges::count_if(indicators, mark)); - + indices.reserve(indicators.size()); for (std::int32_t i = 0; i < static_cast(indicators.size()); ++i) { - if (mark(indicators[i])) + if (indicators[i] > threshold) indices.push_back(i); } @@ -70,7 +67,7 @@ std::vector mark_maximum(MPI_Comm comm, std::span indicators, T theta) { if ((theta <= 0) || (theta >= 1)) - throw std::invalid_argument("Theta needs to fullfill 0 < θ < 1."); + throw std::invalid_argument("theta must fullfill 0 < theta < 1."); T max = indicators.empty() ? std::numeric_limits::lowest() : std::ranges::max(indicators); @@ -103,7 +100,7 @@ std::vector mark_equidistribution(MPI_Comm comm, std::span indicators, T theta) { if ((theta <= 0) || (theta >= 1)) - throw std::invalid_argument("Theta must fullfill 0 < θ < 1."); + throw std::invalid_argument("theta must fullfill 0 < theta < 1."); T norm = std::inner_product(indicators.begin(), indicators.end(), indicators.begin(), T{0}); @@ -143,7 +140,7 @@ mark_equidistribution_squared(MPI_Comm comm, std::span indicators, T theta) { if ((theta <= 0) || (theta >= 1)) - throw std::invalid_argument("Theta must fullfill 0 < θ < 1."); + throw std::invalid_argument("theta must fullfill 0 < theta < 1."); T norm = std::accumulate(indicators.begin(), indicators.end(), T{0}); From e533aed386b6a90258a3b80f50b14784c789af6f Mon Sep 17 00:00:00 2001 From: "Jack S. Hale" Date: Wed, 20 May 2026 10:52:18 +0200 Subject: [PATCH 26/45] Fixes --- cpp/dolfinx/refinement/mark.h | 30 ++++++++++++++++++------------ 1 file changed, 18 insertions(+), 12 deletions(-) diff --git a/cpp/dolfinx/refinement/mark.h b/cpp/dolfinx/refinement/mark.h index e0a16a1cbe..41e0e9768c 100644 --- a/cpp/dolfinx/refinement/mark.h +++ b/cpp/dolfinx/refinement/mark.h @@ -8,11 +8,14 @@ #include #include +#include #include #include #include #include +#include #include +#include #include #include "dolfinx/common/MPI.h" @@ -27,8 +30,8 @@ namespace impl /// /// Helper for other marking routines. /// -/// Returns the indices \f$ i \f$ of the indicators $eta_i$ that satisfy the -/// threshold: \f$ \eta_i > \text{threshold} \f$. +/// 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] indicators Indicators \f$ \eta_i \f$ used for marking. /// @param[in] threshold Threshold value; indicators greater than this are @@ -54,8 +57,8 @@ std::vector mark_threshold(std::span indicators, /// @brief Computes maximum-based marking of indicators. /// -/// Returns the indices \f$ i \f$ of the indicators $eta_i$ that satisfy the -/// maximum threshold: \f$ \eta_i \geq \theta \max_j \eta_j \f$. +/// 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$ - @@ -67,7 +70,7 @@ std::vector mark_maximum(MPI_Comm comm, std::span indicators, T theta) { if ((theta <= 0) || (theta >= 1)) - throw std::invalid_argument("theta must fullfill 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); @@ -100,7 +103,7 @@ std::vector mark_equidistribution(MPI_Comm comm, std::span indicators, T theta) { if ((theta <= 0) || (theta >= 1)) - throw std::invalid_argument("theta must fullfill 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}); @@ -109,7 +112,8 @@ mark_equidistribution(MPI_Comm comm, std::span indicators, T theta) T sqrt_norm = std::sqrt(norm); - std::size_t count = indicators.size(); + // Caution with headroom of global sum across ranks. + std::int64_t count = indicators.size(); MPI_Allreduce(MPI_IN_PLACE, &count, 1, dolfinx::MPI::mpi_t, MPI_SUM, comm); @@ -124,9 +128,10 @@ mark_equidistribution(MPI_Comm comm, std::span indicators, T theta) /// @brief Computes equidistribution threshold marking of a squared indicator. /// -/// Returns the indices \f$i\f$ of the squared indicators $eta_i^2$ 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. +/// 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. /// /// @param[in] comm Communicator over which the global equidistribution /// threshold is computed. @@ -140,13 +145,14 @@ mark_equidistribution_squared(MPI_Comm comm, std::span indicators, T theta) { if ((theta <= 0) || (theta >= 1)) - throw std::invalid_argument("theta must fullfill 0 < theta < 1."); + throw std::invalid_argument("theta must fulfill 0 < theta < 1."); T norm = std::accumulate(indicators.begin(), indicators.end(), T{0}); MPI_Allreduce(MPI_IN_PLACE, &norm, 1, dolfinx::MPI::mpi_t, MPI_SUM, comm); - std::size_t count = indicators.size(); + // Caution with headroom of global sum across ranks. + std::int64_t count = indicators.size(); MPI_Allreduce(MPI_IN_PLACE, &count, 1, dolfinx::MPI::mpi_t, MPI_SUM, comm); From 5fd1a9612d42771bcdc2798410ed44eaef5c9505 Mon Sep 17 00:00:00 2001 From: "Jack S. Hale" Date: Wed, 20 May 2026 10:55:59 +0200 Subject: [PATCH 27/45] More fixes --- cpp/dolfinx/refinement/mark.h | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/cpp/dolfinx/refinement/mark.h b/cpp/dolfinx/refinement/mark.h index 41e0e9768c..e980652e67 100644 --- a/cpp/dolfinx/refinement/mark.h +++ b/cpp/dolfinx/refinement/mark.h @@ -13,6 +13,7 @@ #include #include #include +#include #include #include #include @@ -30,7 +31,7 @@ namespace impl /// /// Helper for other marking routines. /// -/// Returns the indices \f$ i \f$ of the indicators \f$ eta_i \f$ that satisfy +/// 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] indicators Indicators \f$ \eta_i \f$ used for marking. @@ -57,7 +58,7 @@ std::vector mark_threshold(std::span indicators, /// @brief Computes maximum-based marking of indicators. /// -/// Returns the indices \f$ i \f$ of the indicators \f$ eta_i \f$ that satisfy +/// 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. @@ -87,7 +88,7 @@ std::vector mark_maximum(MPI_Comm comm, /// @brief Computes equidistribution threshold marking of indicators. /// -/// Returns the indices \f$ i \f$ of the indicators \f$eta_i \f$ that satisfy +/// 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. @@ -112,7 +113,7 @@ mark_equidistribution(MPI_Comm comm, std::span indicators, T theta) T sqrt_norm = std::sqrt(norm); - // Caution with headroom of global sum across ranks. + // 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); @@ -128,7 +129,7 @@ mark_equidistribution(MPI_Comm comm, std::span indicators, T theta) /// @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 +/// 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. @@ -151,7 +152,7 @@ mark_equidistribution_squared(MPI_Comm comm, std::span indicators, MPI_Allreduce(MPI_IN_PLACE, &norm, 1, dolfinx::MPI::mpi_t, MPI_SUM, comm); - // Caution with headroom of global sum across ranks. + // 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); From 5455f0061fd99595dc2243775fc48f6d4f4006b0 Mon Sep 17 00:00:00 2001 From: "Jack S. Hale" Date: Wed, 20 May 2026 11:02:55 +0200 Subject: [PATCH 28/45] Fix wrappers nbarg and consistent variable names --- .../wrappers/dolfinx_wrappers/refinement.h | 19 ++++++++++--------- 1 file changed, 10 insertions(+), 9 deletions(-) diff --git a/python/dolfinx/wrappers/dolfinx_wrappers/refinement.h b/python/dolfinx/wrappers/dolfinx_wrappers/refinement.h index b05436aff4..dc36db6422 100644 --- a/python/dolfinx/wrappers/dolfinx_wrappers/refinement.h +++ b/python/dolfinx/wrappers/dolfinx_wrappers/refinement.h @@ -128,35 +128,36 @@ void declare_refinement(nanobind::module_& m) m.def( "mark_maximum", [](MPICommWrapper comm, - nb::ndarray, nb::c_contig> indicator, T theta) + nb::ndarray, nb::c_contig> indicators, T theta) { return dolfinx_wrappers::as_nbarray(dolfinx::refinement::mark_maximum( - comm.get(), std::span(indicator.data(), indicator.size()), theta)); + comm.get(), std::span(indicators.data(), indicators.size()), theta)); }, - nb::arg("indicator"), 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> indicator, T theta) + nb::ndarray, nb::c_contig> indicators, T theta) { return dolfinx_wrappers::as_nbarray( dolfinx::refinement::mark_equidistribution( - comm.get(), std::span(indicator.data(), indicator.size()), + comm.get(), std::span(indicators.data(), indicators.size()), theta)); }, - nb::arg("indicator"), nb::arg("theta"), nb::arg("comm")); + nb::arg("comm"), nb::arg("indicators"), nb::arg("theta")); m.def( "mark_equidistribution_squared", [](MPICommWrapper comm, - nb::ndarray, nb::c_contig> marker, T theta) + nb::ndarray, nb::c_contig> indicators, T theta) { return dolfinx_wrappers::as_nbarray( dolfinx::refinement::mark_equidistribution_squared( - comm.get(), std::span(marker.data(), marker.size()), theta)); + 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")); } } // namespace dolfinx_wrappers From f17e463987371353b9d53d34c2afc587d1214db5 Mon Sep 17 00:00:00 2001 From: "Jack S. Hale" Date: Wed, 20 May 2026 11:03:48 +0200 Subject: [PATCH 29/45] clang-format --- python/dolfinx/wrappers/dolfinx_wrappers/refinement.h | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/python/dolfinx/wrappers/dolfinx_wrappers/refinement.h b/python/dolfinx/wrappers/dolfinx_wrappers/refinement.h index dc36db6422..d62d221af8 100644 --- a/python/dolfinx/wrappers/dolfinx_wrappers/refinement.h +++ b/python/dolfinx/wrappers/dolfinx_wrappers/refinement.h @@ -131,7 +131,8 @@ void declare_refinement(nanobind::module_& m) nb::ndarray, nb::c_contig> indicators, T theta) { return dolfinx_wrappers::as_nbarray(dolfinx::refinement::mark_maximum( - comm.get(), std::span(indicators.data(), indicators.size()), theta)); + comm.get(), std::span(indicators.data(), indicators.size()), + theta)); }, nb::arg("comm"), nb::arg("indicators"), nb::arg("theta")); From ddc3eb454d881cbffcfaca379c87101c61ffc55e Mon Sep 17 00:00:00 2001 From: "Jack S. Hale" Date: Wed, 20 May 2026 11:08:16 +0200 Subject: [PATCH 30/45] markers -> indicators in tests --- python/test/unit/refinement/test_mark.py | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/python/test/unit/refinement/test_mark.py b/python/test/unit/refinement/test_mark.py index e3f5f5740d..27e4a37b15 100644 --- a/python/test/unit/refinement/test_mark.py +++ b/python/test/unit/refinement/test_mark.py @@ -21,13 +21,13 @@ def test_mark_maximum(theta: float, dtype: np.dtype) -> None: 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, dtype=dtype) + indicators = np.random.default_rng(0).random(cell_count, dtype=dtype) - marked_cells = mesh.mark_maximum(comm, marker, theta) + 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) @@ -44,19 +44,19 @@ def test_mark_equidistribution(theta: float, dtype: np.dtype) -> None: 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, dtype=dtype) + indicators = np.random.default_rng(0).random(cell_count, dtype=dtype) - marked_cells = mesh.mark_equidistribution(comm, marker, theta) + marked_cells = mesh.mark_equidistribution(comm, indicators, theta) - norm = np.sqrt(comm.allreduce(np.sum(marker**2), MPI.SUM)) - count = comm.allreduce(marker.size) + norm = np.sqrt(comm.allreduce(np.sum(indicators**2), MPI.SUM)) + count = comm.allreduce(indicators.size) assert np.allclose( marked_cells, - np.argwhere(marker > theta * norm / np.sqrt(count)).flatten(), + 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, marker**2, theta)) + assert np.all(marked_cells == mesh.mark_equidistribution_squared(comm, indicators**2, theta)) From 15e04385d742ad2b7f181bff9a44e4a4d5ca7405 Mon Sep 17 00:00:00 2001 From: "Jack S. Hale" Date: Wed, 20 May 2026 11:14:53 +0200 Subject: [PATCH 31/45] Improve docstring --- cpp/dolfinx/refinement/mark.h | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/cpp/dolfinx/refinement/mark.h b/cpp/dolfinx/refinement/mark.h index e980652e67..984f525888 100644 --- a/cpp/dolfinx/refinement/mark.h +++ b/cpp/dolfinx/refinement/mark.h @@ -34,7 +34,8 @@ namespace impl /// 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] indicators Indicators \f$ \eta_i \f$ used for marking. +/// @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. From fa679a2da1fd2d07170d029419ca2a0bcac2fab4 Mon Sep 17 00:00:00 2001 From: "Jack S. Hale" Date: Wed, 20 May 2026 11:20:16 +0200 Subject: [PATCH 32/45] Add Python wrappers for mark_* refinement functions Wrap the C++ bindings so the public Python API carries proper type hints and docstrings instead of re-exporting the raw nanobind bindings. Co-Authored-By: Claude Opus 4.7 --- python/dolfinx/mesh.py | 85 ++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 82 insertions(+), 3 deletions(-) diff --git a/python/dolfinx/mesh.py b/python/dolfinx/mesh.py index c1fff7049f..c5d032987d 100644 --- a/python/dolfinx/mesh.py +++ b/python/dolfinx/mesh.py @@ -35,9 +35,15 @@ from dolfinx.cpp.refinement import ( IdentityPartitionerPlaceholder, RefinementOption, - mark_equidistribution, - mark_equidistribution_squared, - 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, @@ -764,6 +770,79 @@ def refine( return Mesh(mesh1, ufl_domain), parent_cell, parent_facet +def mark_maximum( + comm: _MPI.Comm, + indicators: npt.NDArray[np.floating], + theta: float, +) -> npt.NDArray[np.int32]: + """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[np.floating], + theta: float, +) -> npt.NDArray[np.int32]: + """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 (local) :math:`\\eta_i` - usually + associated with mesh entity :math:`i`. + 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, + indicators: npt.NDArray[np.floating], + theta: float, +) -> npt.NDArray[np.int32]: + """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. + indicators: Input indicators (local) :math:`\\eta_i^2` - usually + associated with mesh entity :math:`i`. + theta: Parameter, :math:`0 < \\theta < 1`. + + Returns: + Local indices of indicators that satisfy the threshold. + """ + return _mark_equidistribution_squared(comm, indicators, theta) + + def create_mesh( comm: _MPI.Comm, cells: npt.NDArray[np.int64], From c475ca0a0b5289847ee7f83d92522247905ebad2 Mon Sep 17 00:00:00 2001 From: "Jack S. Hale" Date: Wed, 20 May 2026 11:24:50 +0200 Subject: [PATCH 33/45] Use Real --- python/dolfinx/mesh.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/python/dolfinx/mesh.py b/python/dolfinx/mesh.py index c5d032987d..d01768a0e8 100644 --- a/python/dolfinx/mesh.py +++ b/python/dolfinx/mesh.py @@ -772,7 +772,7 @@ def refine( def mark_maximum( comm: _MPI.Comm, - indicators: npt.NDArray[np.floating], + indicators: npt.NDArray[Real], theta: float, ) -> npt.NDArray[np.int32]: """Compute maximum-based marking of indicators. @@ -795,7 +795,7 @@ def mark_maximum( def mark_equidistribution( comm: _MPI.Comm, - indicators: npt.NDArray[np.floating], + indicators: npt.NDArray[Real], theta: float, ) -> npt.NDArray[np.int32]: """Compute equidistribution threshold marking of indicators. @@ -820,7 +820,7 @@ def mark_equidistribution( def mark_equidistribution_squared( comm: _MPI.Comm, - indicators: npt.NDArray[np.floating], + indicators: npt.NDArray[Real], theta: float, ) -> npt.NDArray[np.int32]: """Compute equidistribution threshold marking of a squared indicator. From 43b52039a9bfbac8f0c6577024aa944f0a9d304b Mon Sep 17 00:00:00 2001 From: "Jack S. Hale" Date: Wed, 20 May 2026 11:29:50 +0200 Subject: [PATCH 34/45] markers -> indicators in mark.cpp test Co-Authored-By: Claude Sonnet 4.6 --- cpp/test/mesh/refinement/mark.cpp | 59 ++++++++++++++++++------------- 1 file changed, 34 insertions(+), 25 deletions(-) diff --git a/cpp/test/mesh/refinement/mark.cpp b/cpp/test/mesh/refinement/mark.cpp index fdebe0184b..af2b070d7c 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(MPI_COMM_WORLD, marker, .5); + std::vector indicators; + auto indices = mark_maximum(MPI_COMM_WORLD, indicators, .5); CHECK(indices.size() == 0); } @@ -26,27 +26,31 @@ 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(comm, marker, theta); + 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); } @@ -55,8 +59,9 @@ TEMPLATE_TEST_CASE("Mark maximum", "[refinement][mark][maximum]", double, float) TEMPLATE_TEST_CASE("Mark equidistribution empty", "[refinement][mark][equidistribution]", double, float) { - std::vector marker; - auto indices = mark_equidistribution(MPI_COMM_WORLD, marker, .5); + std::vector indicators; + auto indices + = mark_equidistribution(MPI_COMM_WORLD, indicators, .5); CHECK(indices.size() == 0); } @@ -65,37 +70,41 @@ TEMPLATE_TEST_CASE("Mark equidistribution", { 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_equidistribution(comm, marker, theta); + auto indices = mark_equidistribution(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())); + })); std::int32_t n = dolfinx::MPI::size(comm) * 10 - 1; TestType norm = std::sqrt(n * (n + 1) * (2 * n + 1) / 6); auto mark = [=](auto e) { return e > theta * norm / std::sqrt(n); }; - 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); } // squared input - auto marker_sq = marker; - std::ranges::for_each(marker_sq, [](auto& e) { e = std::pow(e, 2); }); + 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, marker_sq, theta))); + comm, indicators_sq, theta))); } From ffc96abe60ec668d8f8353e9bf84efb0b0e11538 Mon Sep 17 00:00:00 2001 From: "Jack S. Hale" Date: Wed, 20 May 2026 11:36:56 +0200 Subject: [PATCH 35/45] Make squared_indicators very explicit --- cpp/dolfinx/refinement/mark.h | 17 +++++++++-------- python/dolfinx/mesh.py | 8 ++++---- .../wrappers/dolfinx_wrappers/refinement.h | 9 ++++++--- 3 files changed, 19 insertions(+), 15 deletions(-) diff --git a/cpp/dolfinx/refinement/mark.h b/cpp/dolfinx/refinement/mark.h index 984f525888..f1961d0cb0 100644 --- a/cpp/dolfinx/refinement/mark.h +++ b/cpp/dolfinx/refinement/mark.h @@ -137,32 +137,33 @@ mark_equidistribution(MPI_Comm comm, std::span indicators, T theta) /// /// @param[in] comm Communicator over which the global equidistribution /// threshold is computed. -/// @param[in] indicators Input indicators (local) \f$ \eta^2_i \f$ - +/// @param[in] squared_indicators Input squared indicators (local) \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 indicators that satisfy the threshold. +/// @return Local indices of squared indicators that satisfy the threshold. template std::vector -mark_equidistribution_squared(MPI_Comm comm, std::span indicators, - T theta) +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(indicators.begin(), indicators.end(), T{0}); + 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 = indicators.size(); + 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( - indicators, std::pow(theta, 2) * norm / static_cast(count)); + squared_indicators, std::pow(theta, 2) * norm / static_cast(count)); spdlog::info("Marking (equidistribution) {} of {} local entities.", - indices.size(), indicators.size()); + indices.size(), squared_indicators.size()); return indices; } diff --git a/python/dolfinx/mesh.py b/python/dolfinx/mesh.py index d01768a0e8..c9b763063a 100644 --- a/python/dolfinx/mesh.py +++ b/python/dolfinx/mesh.py @@ -820,7 +820,7 @@ def mark_equidistribution( def mark_equidistribution_squared( comm: _MPI.Comm, - indicators: npt.NDArray[Real], + squared_indicators: npt.NDArray[Real], theta: float, ) -> npt.NDArray[np.int32]: """Compute equidistribution threshold marking of a squared indicator. @@ -833,14 +833,14 @@ def mark_equidistribution_squared( Args: comm: Communicator over which the global equidistribution threshold is computed. - indicators: Input indicators (local) :math:`\\eta_i^2` - usually + squared_indicators: Input squared indicators (local) :math:`\\eta_i^2` - usually associated with mesh entity :math:`i`. theta: Parameter, :math:`0 < \\theta < 1`. Returns: - Local indices of indicators that satisfy the threshold. + Local indices of squared indicators that satisfy the threshold. """ - return _mark_equidistribution_squared(comm, indicators, theta) + return _mark_equidistribution_squared(comm, squared_indicators, theta) def create_mesh( diff --git a/python/dolfinx/wrappers/dolfinx_wrappers/refinement.h b/python/dolfinx/wrappers/dolfinx_wrappers/refinement.h index d62d221af8..d1b3026bd6 100644 --- a/python/dolfinx/wrappers/dolfinx_wrappers/refinement.h +++ b/python/dolfinx/wrappers/dolfinx_wrappers/refinement.h @@ -151,14 +151,17 @@ void declare_refinement(nanobind::module_& m) m.def( "mark_equidistribution_squared", [](MPICommWrapper comm, - nb::ndarray, nb::c_contig> indicators, T theta) + nb::ndarray, nb::c_contig> squared_indicators, + T theta) { return dolfinx_wrappers::as_nbarray( dolfinx::refinement::mark_equidistribution_squared( - comm.get(), std::span(indicators.data(), indicators.size()), + comm.get(), + std::span(squared_indicators.data(), + squared_indicators.size()), theta)); }, - nb::arg("comm"), nb::arg("indicators"), nb::arg("theta")); + nb::arg("comm"), nb::arg("squared_indicators"), nb::arg("theta")); } } // namespace dolfinx_wrappers From 5dbcea44c5d53eed0c7f08cc7f7260d9765e5fb9 Mon Sep 17 00:00:00 2001 From: "Jack S. Hale" Date: Wed, 20 May 2026 11:39:10 +0200 Subject: [PATCH 36/45] Fix up --- cpp/dolfinx/refinement/mark.h | 8 ++++---- python/dolfinx/mesh.py | 4 ++-- python/dolfinx/wrappers/dolfinx_wrappers/refinement.h | 3 +-- 3 files changed, 7 insertions(+), 8 deletions(-) diff --git a/cpp/dolfinx/refinement/mark.h b/cpp/dolfinx/refinement/mark.h index f1961d0cb0..f16ee32e5e 100644 --- a/cpp/dolfinx/refinement/mark.h +++ b/cpp/dolfinx/refinement/mark.h @@ -137,8 +137,8 @@ mark_equidistribution(MPI_Comm comm, std::span indicators, T theta) /// /// @param[in] comm Communicator over which the global equidistribution /// threshold is computed. -/// @param[in] squared_indicators Input squared indicators (local) \f$ \eta^2_i \f$ - -/// usually associated with mesh entity \f$ i \f$. +/// @param[in] squared_indicators Input squared indicators (local) \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 @@ -149,8 +149,8 @@ mark_equidistribution_squared(MPI_Comm comm, 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}); + 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); diff --git a/python/dolfinx/mesh.py b/python/dolfinx/mesh.py index c9b763063a..e31eccb1ec 100644 --- a/python/dolfinx/mesh.py +++ b/python/dolfinx/mesh.py @@ -833,8 +833,8 @@ def mark_equidistribution_squared( Args: comm: Communicator over which the global equidistribution threshold is computed. - squared_indicators: Input squared indicators (local) :math:`\\eta_i^2` - usually - associated with mesh entity :math:`i`. + squared_indicators: Input squared indicators (local) :math:`\\eta_i^2` + - usually associated with mesh entity :math:`i`. theta: Parameter, :math:`0 < \\theta < 1`. Returns: diff --git a/python/dolfinx/wrappers/dolfinx_wrappers/refinement.h b/python/dolfinx/wrappers/dolfinx_wrappers/refinement.h index d1b3026bd6..7eac92d1ee 100644 --- a/python/dolfinx/wrappers/dolfinx_wrappers/refinement.h +++ b/python/dolfinx/wrappers/dolfinx_wrappers/refinement.h @@ -157,8 +157,7 @@ void declare_refinement(nanobind::module_& m) return dolfinx_wrappers::as_nbarray( dolfinx::refinement::mark_equidistribution_squared( comm.get(), - std::span(squared_indicators.data(), - squared_indicators.size()), + std::span(squared_indicators.data(), squared_indicators.size()), theta)); }, nb::arg("comm"), nb::arg("squared_indicators"), nb::arg("theta")); From 86f97650e542a432141f0d0544cd792c9e5398c4 Mon Sep 17 00:00:00 2001 From: "Jack S. Hale" Date: Wed, 20 May 2026 11:39:47 +0200 Subject: [PATCH 37/45] Ruff manual fix --- python/dolfinx/mesh.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/python/dolfinx/mesh.py b/python/dolfinx/mesh.py index e31eccb1ec..a8b3b8b77c 100644 --- a/python/dolfinx/mesh.py +++ b/python/dolfinx/mesh.py @@ -833,8 +833,9 @@ def mark_equidistribution_squared( Args: comm: Communicator over which the global equidistribution threshold is computed. - squared_indicators: Input squared indicators (local) :math:`\\eta_i^2` - - usually associated with mesh entity :math:`i`. + squared_indicators: Input squared indicators (local) + :math:`\\eta_i^2` - usually associated with mesh + entity :math:`i`. theta: Parameter, :math:`0 < \\theta < 1`. Returns: From 4936a9441af856f6327e0de913d36d359c277e49 Mon Sep 17 00:00:00 2001 From: "Jack S. Hale" Date: Wed, 20 May 2026 11:45:40 +0200 Subject: [PATCH 38/45] Raw strings --- python/dolfinx/mesh.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/python/dolfinx/mesh.py b/python/dolfinx/mesh.py index a8b3b8b77c..2beae6f557 100644 --- a/python/dolfinx/mesh.py +++ b/python/dolfinx/mesh.py @@ -775,7 +775,7 @@ def mark_maximum( indicators: npt.NDArray[Real], theta: float, ) -> npt.NDArray[np.int32]: - """Compute maximum-based marking of indicators. + r"""Compute maximum-based marking of indicators. Returns the indices :math:`i` of the indicators :math:`\\eta_i` that satisfy the maximum threshold: @@ -798,7 +798,7 @@ def mark_equidistribution( indicators: npt.NDArray[Real], theta: float, ) -> npt.NDArray[np.int32]: - """Compute equidistribution threshold marking of indicators. + r"""Compute equidistribution threshold marking of indicators. Returns the indices :math:`i` of the indicators :math:`\\eta_i` that satisfy the equidistribution threshold: @@ -823,7 +823,7 @@ def mark_equidistribution_squared( squared_indicators: npt.NDArray[Real], theta: float, ) -> npt.NDArray[np.int32]: - """Compute equidistribution threshold marking of a squared indicator. + 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: From 254c711e0d9f721bed51a93391b765731157b38d Mon Sep 17 00:00:00 2001 From: "Jack S. Hale" Date: Wed, 20 May 2026 11:52:49 +0200 Subject: [PATCH 39/45] Fix double backslashes in raw string docstrings Co-Authored-By: Claude Sonnet 4.6 --- python/dolfinx/mesh.py | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/python/dolfinx/mesh.py b/python/dolfinx/mesh.py index 2beae6f557..35f8ae6654 100644 --- a/python/dolfinx/mesh.py +++ b/python/dolfinx/mesh.py @@ -777,15 +777,15 @@ def mark_maximum( ) -> npt.NDArray[np.int32]: r"""Compute maximum-based marking of indicators. - Returns the indices :math:`i` of the indicators :math:`\\eta_i` that + Returns the indices :math:`i` of the indicators :math:`\eta_i` that satisfy the maximum threshold: - :math:`\\eta_i > \\theta \\max_j \\eta_j`. + :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 + indicators: Indicators (local) :math:`\eta_i` - usually an error indicator associated with mesh entity :math:`i`. - theta: Parameter, :math:`0 < \\theta < 1`. + theta: Parameter, :math:`0 < \theta < 1`. Returns: Local indices of marked entities. @@ -800,17 +800,17 @@ def mark_equidistribution( ) -> npt.NDArray[np.int32]: r"""Compute equidistribution threshold marking of indicators. - Returns the indices :math:`i` of the indicators :math:`\\eta_i` that + 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:`\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 (local) :math:`\\eta_i` - usually + indicators: Indicators (local) :math:`\eta_i` - usually associated with mesh entity :math:`i`. - theta: Parameter, :math:`0 < \\theta < 1`. + theta: Parameter, :math:`0 < \theta < 1`. Returns: Local indices of indicators that satisfy the threshold. @@ -826,17 +826,17 @@ def mark_equidistribution_squared( 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:`\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 (local) - :math:`\\eta_i^2` - usually associated with mesh + :math:`\eta_i^2` - usually associated with mesh entity :math:`i`. - theta: Parameter, :math:`0 < \\theta < 1`. + theta: Parameter, :math:`0 < \theta < 1`. Returns: Local indices of squared indicators that satisfy the threshold. From 7c5f5a6b06b9213a3f28075900100e7d648dad28 Mon Sep 17 00:00:00 2001 From: "Jack S. Hale" Date: Wed, 20 May 2026 12:38:42 +0200 Subject: [PATCH 40/45] Remove std::pow --- cpp/dolfinx/refinement/mark.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cpp/dolfinx/refinement/mark.h b/cpp/dolfinx/refinement/mark.h index f16ee32e5e..36993d211f 100644 --- a/cpp/dolfinx/refinement/mark.h +++ b/cpp/dolfinx/refinement/mark.h @@ -160,7 +160,7 @@ mark_equidistribution_squared(MPI_Comm comm, MPI_SUM, comm); std::vector indices = impl::mark_threshold( - squared_indicators, std::pow(theta, 2) * norm / static_cast(count)); + squared_indicators, theta * theta * norm / static_cast(count)); spdlog::info("Marking (equidistribution) {} of {} local entities.", indices.size(), squared_indicators.size()); From 6caed8d0e5259ced96bff3cec0ac2736814510f1 Mon Sep 17 00:00:00 2001 From: "Jack S. Hale" Date: Wed, 20 May 2026 12:56:35 +0200 Subject: [PATCH 41/45] Add warning to docstring --- cpp/dolfinx/refinement/mark.h | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/cpp/dolfinx/refinement/mark.h b/cpp/dolfinx/refinement/mark.h index 36993d211f..909bdadb81 100644 --- a/cpp/dolfinx/refinement/mark.h +++ b/cpp/dolfinx/refinement/mark.h @@ -94,10 +94,14 @@ std::vector mark_maximum(MPI_Comm comm, /// \frac{||\eta||}{\sqrt{N}} \f$ where \f$ N \f$ is the (global) number of /// indicators. /// +/// @warning `indicators` must contain values for owned entities only (sliced +/// to `index_map.size_local()`). Ghost values cause double-counting in the +/// global reductions and an incorrect threshold. +/// /// @param[in] comm Communicator over which the global equidistribution /// threshold is computed. -/// @param[in] indicators Indicators (local) \f$ \eta_i \f$ - usually -/// associated with mesh entity \f$ i \f$. +/// @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 @@ -135,10 +139,14 @@ mark_equidistribution(MPI_Comm comm, std::span indicators, T theta) /// \frac{||\eta||^2}{N} \f$ where \f$ N \f$ is the (global) number of /// indicators. /// +/// @warning `squared_indicators` must contain values for owned entities only +/// (sliced to `index_map.size_local()`). Ghost values cause double-counting in +/// the global reductions and an incorrect threshold. +/// /// @param[in] comm Communicator over which the global equidistribution /// threshold is computed. -/// @param[in] squared_indicators Input squared indicators (local) \f$ \eta^2_i -/// \f$ - usually associated with mesh entity \f$ i \f$. +/// @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 From bcc5bcd3367bd9844ec3e3c03e07a6d60fa77367 Mon Sep 17 00:00:00 2001 From: "Jack S. Hale" Date: Wed, 20 May 2026 12:59:48 +0200 Subject: [PATCH 42/45] Add to Python docstring --- python/dolfinx/mesh.py | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/python/dolfinx/mesh.py b/python/dolfinx/mesh.py index 35f8ae6654..736fec6bfb 100644 --- a/python/dolfinx/mesh.py +++ b/python/dolfinx/mesh.py @@ -808,8 +808,10 @@ def mark_equidistribution( Args: comm: Communicator over which the global equidistribution threshold is computed. - indicators: Indicators (local) :math:`\eta_i` - usually - associated with mesh entity :math:`i`. + indicators: Indicators for owned local entities :math:`\eta_i` - + usually associated with mesh entity :math:`i`. Ghost values + cause double-counting in the global reductions and an incorrect + threshold. theta: Parameter, :math:`0 < \theta < 1`. Returns: @@ -833,9 +835,10 @@ def mark_equidistribution_squared( Args: comm: Communicator over which the global equidistribution threshold is computed. - squared_indicators: Input squared indicators (local) - :math:`\eta_i^2` - usually associated with mesh - entity :math:`i`. + squared_indicators: Input squared indicators for owned local + entities :math:`\eta_i^2` - usually associated with mesh + entity :math:`i`. Ghost values cause double-counting in the + global reductions and an incorrect threshold. theta: Parameter, :math:`0 < \theta < 1`. Returns: From 3b464e1c538727fd3b990788711926d7c1d7e76f Mon Sep 17 00:00:00 2001 From: "Jack S. Hale" Date: Wed, 20 May 2026 13:13:06 +0200 Subject: [PATCH 43/45] Also fix up tests to show correct usage. --- cpp/test/mesh/refinement/mark.cpp | 3 ++- python/test/unit/refinement/test_mark.py | 4 ++-- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/cpp/test/mesh/refinement/mark.cpp b/cpp/test/mesh/refinement/mark.cpp index af2b070d7c..9ebe7ca45e 100644 --- a/cpp/test/mesh/refinement/mark.cpp +++ b/cpp/test/mesh/refinement/mark.cpp @@ -86,9 +86,10 @@ TEMPLATE_TEST_CASE("Mark equidistribution", })); 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(n); }; + auto mark = [=](auto e) { return e > theta * norm / std::sqrt(count); }; CHECK(std::ranges::count_if(indicators, mark) == static_cast(indices.size())); diff --git a/python/test/unit/refinement/test_mark.py b/python/test/unit/refinement/test_mark.py index 27e4a37b15..22c9b6cc0b 100644 --- a/python/test/unit/refinement/test_mark.py +++ b/python/test/unit/refinement/test_mark.py @@ -20,7 +20,7 @@ def test_mark_maximum(theta: float, dtype: np.dtype) -> None: 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 + 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(comm, indicators, theta) @@ -43,7 +43,7 @@ def test_mark_equidistribution(theta: float, dtype: np.dtype) -> None: 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 + 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) From 573b8daf8ae8658fe394c3b57b903dae9cef9457 Mon Sep 17 00:00:00 2001 From: "Jack S. Hale" Date: Wed, 20 May 2026 13:20:11 +0200 Subject: [PATCH 44/45] Simplify docstrings --- cpp/dolfinx/refinement/mark.h | 10 ++++------ python/dolfinx/mesh.py | 10 +++++----- 2 files changed, 9 insertions(+), 11 deletions(-) diff --git a/cpp/dolfinx/refinement/mark.h b/cpp/dolfinx/refinement/mark.h index 909bdadb81..bdee2b169b 100644 --- a/cpp/dolfinx/refinement/mark.h +++ b/cpp/dolfinx/refinement/mark.h @@ -94,9 +94,8 @@ std::vector mark_maximum(MPI_Comm comm, /// \frac{||\eta||}{\sqrt{N}} \f$ where \f$ N \f$ is the (global) number of /// indicators. /// -/// @warning `indicators` must contain values for owned entities only (sliced -/// to `index_map.size_local()`). Ghost values cause double-counting in the -/// global reductions and an incorrect threshold. +/// @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. @@ -139,9 +138,8 @@ mark_equidistribution(MPI_Comm comm, std::span indicators, T theta) /// \frac{||\eta||^2}{N} \f$ where \f$ N \f$ is the (global) number of /// indicators. /// -/// @warning `squared_indicators` must contain values for owned entities only -/// (sliced to `index_map.size_local()`). Ghost values cause double-counting in -/// the global reductions and an incorrect threshold. +/// @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. diff --git a/python/dolfinx/mesh.py b/python/dolfinx/mesh.py index 736fec6bfb..e3ec70ddf1 100644 --- a/python/dolfinx/mesh.py +++ b/python/dolfinx/mesh.py @@ -809,9 +809,9 @@ def mark_equidistribution( 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`. Ghost values - cause double-counting in the global reductions and an incorrect - threshold. + 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: @@ -837,8 +837,8 @@ def mark_equidistribution_squared( threshold is computed. squared_indicators: Input squared indicators for owned local entities :math:`\eta_i^2` - usually associated with mesh - entity :math:`i`. Ghost values cause double-counting in the - global reductions and an incorrect threshold. + 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: From e5e056450be4d013f924a463cd53f1da0a15f242 Mon Sep 17 00:00:00 2001 From: "Jack S. Hale" Date: Wed, 20 May 2026 13:59:49 +0200 Subject: [PATCH 45/45] Revert to count_if --- cpp/dolfinx/refinement/mark.h | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/cpp/dolfinx/refinement/mark.h b/cpp/dolfinx/refinement/mark.h index bdee2b169b..17fccbb484 100644 --- a/cpp/dolfinx/refinement/mark.h +++ b/cpp/dolfinx/refinement/mark.h @@ -43,12 +43,15 @@ template std::vector mark_threshold(std::span indicators, T threshold) { + auto mark = [threshold](T e) { return e > threshold; }; + std::vector indices; - indices.reserve(indicators.size()); + indices.reserve(std::ranges::count_if(indicators, mark)); + for (std::int32_t i = 0; i < static_cast(indicators.size()); ++i) { - if (indicators[i] > threshold) + if (mark(indicators[i])) indices.push_back(i); }