From 0cc1f689340d622afec1de9d60c0fc4a5ac8436b Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Sat, 18 Apr 2026 13:58:41 +0200 Subject: [PATCH 1/6] Add maximum marking (c++) --- cpp/dolfinx/refinement/CMakeLists.txt | 1 + cpp/dolfinx/refinement/mark.h | 59 +++++++++++++++++++++++++++ cpp/test/CMakeLists.txt | 1 + cpp/test/mesh/refinement/mark.cpp | 53 ++++++++++++++++++++++++ 4 files changed, 114 insertions(+) create mode 100644 cpp/dolfinx/refinement/mark.h create mode 100644 cpp/test/mesh/refinement/mark.cpp diff --git a/cpp/dolfinx/refinement/CMakeLists.txt b/cpp/dolfinx/refinement/CMakeLists.txt index 0d9d997c490..6db168cf3a4 100644 --- a/cpp/dolfinx/refinement/CMakeLists.txt +++ b/cpp/dolfinx/refinement/CMakeLists.txt @@ -1,6 +1,7 @@ set(HEADERS_refinement ${CMAKE_CURRENT_SOURCE_DIR}/dolfinx_refinement.h ${CMAKE_CURRENT_SOURCE_DIR}/interval.h + ${CMAKE_CURRENT_SOURCE_DIR}/mark.h ${CMAKE_CURRENT_SOURCE_DIR}/plaza.h ${CMAKE_CURRENT_SOURCE_DIR}/refine.h ${CMAKE_CURRENT_SOURCE_DIR}/utils.h diff --git a/cpp/dolfinx/refinement/mark.h b/cpp/dolfinx/refinement/mark.h new file mode 100644 index 00000000000..a622457408d --- /dev/null +++ b/cpp/dolfinx/refinement/mark.h @@ -0,0 +1,59 @@ +// Copyright (C) 2026 Paul T. Kühner +// +// This file is part of DOLFINX (https://www.fenicsproject.org) +// +// SPDX-License-Identifier: LGPL-3.0-or-later + +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "dolfinx/common/MPI.h" + +namespace dolfinx::refinement +{ + +/// @brief Maximum marking of a marker. +/// +/// @param[in] marker Input marker (local) - usually an error indicator per +/// entity +/// @param[in] theta Cut off parameter, 0 ≤ θ ≤ 1 +/// @param[in] comm Communicator over which the maximum is computed. +/// @return Indices (local) of marker elements, which satisfy: marker_i ≥ θ +/// max(marker). +template +std::vector mark_maximum(std::span marker, T theta, + MPI_Comm comm) +{ + assert((0 <= theta) && (theta <= 1)); + + T max = marker.empty() ? std::numeric_limits::min() + : std::ranges::max(marker); + MPI_Allreduce(MPI_IN_PLACE, &max, 1, dolfinx::MPI::mpi_t, MPI_MAX, comm); + + auto mark = [=](auto 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); + } + + spdlog::info("Marking (max) {} / {} (local) entities.", indices.size(), + marker.size()); + + return indices; +} + +} // namespace dolfinx::refinement diff --git a/cpp/test/CMakeLists.txt b/cpp/test/CMakeLists.txt index f446631aa06..c64cb105b3d 100644 --- a/cpp/test/CMakeLists.txt +++ b/cpp/test/CMakeLists.txt @@ -97,6 +97,7 @@ add_executable( mesh/refinement/interval.cpp mesh/refinement/option.cpp mesh/refinement/rectangle.cpp + mesh/refinement/mark.cpp ${CMAKE_CURRENT_BINARY_DIR}/expr.c ${CMAKE_CURRENT_BINARY_DIR}/poisson.c ) diff --git a/cpp/test/mesh/refinement/mark.cpp b/cpp/test/mesh/refinement/mark.cpp new file mode 100644 index 00000000000..1b2c79ae185 --- /dev/null +++ b/cpp/test/mesh/refinement/mark.cpp @@ -0,0 +1,53 @@ +// Copyright (C) 2026 Paul T. Kühner +// +// This file is part of DOLFINX (https://www.fenicsproject.org) +// +// SPDX-License-Identifier: LGPL-3.0-or-later + +#include +#include +#include +#include +#include +#include + +using namespace dolfinx; +using namespace dolfinx::refinement; + +TEMPLATE_TEST_CASE("Mark maximum empty", "[refinement][mark][maximum]", double, + float) +{ + std::vector marker; + auto indices = mark_maximum(marker, .5, MPI_COMM_WORLD); + CHECK(indices.size() == 0); +} + +TEMPLATE_TEST_CASE("Mark maximum", "[refinement][mark][maximum]", double, float) +{ + MPI_Comm comm = MPI_COMM_WORLD; + + std::vector marker; + marker.reserve(10); + for (std::size_t i = 0; i < marker.size(); i++) + marker.push_back(10 * dolfinx::MPI::rank(comm) + i); + + TestType theta = 0.5; + auto indices = mark_maximum(marker, theta, MPI_COMM_WORLD); + + CHECK(std::ranges::all_of( + indices, [&](auto e) + { return (0 <= e) && (e <= static_cast(marker.size())); })); + + TestType max = dolfinx::MPI::size(comm) * 10 - 1; + auto mark = [=](auto e) { return e >= theta * max; }; + + 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 60086aaa2cd45a42651041d7058e7896305c3940 Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Sat, 18 Apr 2026 16:21:50 +0200 Subject: [PATCH 2/6] Add maximum marking (py) --- python/dolfinx/mesh.py | 2 ++ .../wrappers/dolfinx_wrappers/refinement.h | 13 ++++++++ python/test/unit/refinement/test_mark.py | 33 +++++++++++++++++++ 3 files changed, 48 insertions(+) create mode 100644 python/test/unit/refinement/test_mark.py diff --git a/python/dolfinx/mesh.py b/python/dolfinx/mesh.py index c326acc04b0..28946adc7f1 100644 --- a/python/dolfinx/mesh.py +++ b/python/dolfinx/mesh.py @@ -35,6 +35,7 @@ from dolfinx.cpp.refinement import ( IdentityPartitionerPlaceholder, RefinementOption, + mark_maximum, ) from dolfinx.cpp.refinement import ( uniform_refine as _uniform_refine, @@ -72,6 +73,7 @@ "exterior_facet_indices", "locate_entities", "locate_entities_boundary", + "mark_maximum", "meshtags", "meshtags_from_entities", "refine", diff --git a/python/dolfinx/wrappers/dolfinx_wrappers/refinement.h b/python/dolfinx/wrappers/dolfinx_wrappers/refinement.h index f14d574b6ee..7ba4fd8ea22 100644 --- a/python/dolfinx/wrappers/dolfinx_wrappers/refinement.h +++ b/python/dolfinx/wrappers/dolfinx_wrappers/refinement.h @@ -7,10 +7,13 @@ #pragma once +#include "MPICommWrapper.h" #include "array.h" +#include "caster_mpi.h" #include "mesh.h" #include #include +#include #include #include #include @@ -121,6 +124,16 @@ void declare_refinement(nanobind::module_& m) }, nb::arg("mesh"), nb::arg("edges").none(), nb::arg("partitioner").none(), nb::arg("option")); + + m.def( + "mark_maximum", + [](nb::ndarray, nb::c_contig> marker, T theta, + MPICommWrapper comm) + { + return dolfinx_wrappers::as_nbarray(dolfinx::refinement::mark_maximum( + 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 new file mode 100644 index 00000000000..7e4868cf47a --- /dev/null +++ b/python/test/unit/refinement/test_mark.py @@ -0,0 +1,33 @@ +# Copyright (C) 2026 Paul T. Kühner +# +# This file is part of DOLFINx (https://www.fenicsproject.org) +# +# SPDX-License-Identifier: LGPL-3.0-or-later + +from mpi4py import MPI + +import numpy as np +import pytest + +from dolfinx import mesh + + +@pytest.mark.parametrize("theta", np.linspace(0, 1, num=5, endpoint=True)) +@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) + + 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_maximum(marker, theta, comm) + + assert np.allclose( + marked_cells, + np.argwhere(marker >= theta * comm.allreduce(np.max(marker), MPI.MAX)).flatten(), + ) + + msh.topology.create_entities(1) + marked_edges = mesh.compute_incident_entities(msh.topology, marked_cells, tdim, 1) + mesh.refine(msh, marked_edges) From bfedde2dbdfb341d93af654ef185ee129f8e9ff2 Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Mon, 20 Apr 2026 12:20:08 +0200 Subject: [PATCH 3/6] Change assert to exception --- 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 a622457408d..8dcd095b81e 100644 --- a/cpp/dolfinx/refinement/mark.h +++ b/cpp/dolfinx/refinement/mark.h @@ -33,7 +33,8 @@ template std::vector mark_maximum(std::span marker, T theta, MPI_Comm comm) { - assert((0 <= theta) && (theta <= 1)); + if ((theta < 0) || (theta > 1)) + throw std::invalid_argument("Theta needs to fullfill 0 ≤ θ ≤ 1."); T max = marker.empty() ? std::numeric_limits::min() : std::ranges::max(marker); From be03e8141e5438307665ae9cc740e87dbc2be6d3 Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Mon, 20 Apr 2026 12:24:57 +0200 Subject: [PATCH 4/6] Use T over auto --- 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 8dcd095b81e..1b00d7089a1 100644 --- a/cpp/dolfinx/refinement/mark.h +++ b/cpp/dolfinx/refinement/mark.h @@ -40,7 +40,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 = [=](auto e) { return e >= theta * max; }; + auto mark = [=](T e) { return e >= theta * max; }; std::vector indices; indices.reserve(std::ranges::count_if(marker, mark)); From fef53a42370823a11f24418bdf20014987dfa7a0 Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Mon, 20 Apr 2026 12:34:01 +0200 Subject: [PATCH 5/6] Fix: vector empty --- cpp/test/mesh/refinement/mark.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/cpp/test/mesh/refinement/mark.cpp b/cpp/test/mesh/refinement/mark.cpp index 1b2c79ae185..aac1cfaef3e 100644 --- a/cpp/test/mesh/refinement/mark.cpp +++ b/cpp/test/mesh/refinement/mark.cpp @@ -8,6 +8,7 @@ #include #include #include +#include #include #include @@ -28,7 +29,7 @@ TEMPLATE_TEST_CASE("Mark maximum", "[refinement][mark][maximum]", double, float) std::vector marker; marker.reserve(10); - for (std::size_t i = 0; i < marker.size(); i++) + for (std::size_t i = 0; i < 10; i++) marker.push_back(10 * dolfinx::MPI::rank(comm) + i); TestType theta = 0.5; From 0424889a25be96a7a5a464a55f1ff87be62d7d0f Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Mon, 20 Apr 2026 15:02:21 +0200 Subject: [PATCH 6/6] Fix: min -> lowest --- 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 1b00d7089a1..8f61a80cd3a 100644 --- a/cpp/dolfinx/refinement/mark.h +++ b/cpp/dolfinx/refinement/mark.h @@ -36,7 +36,7 @@ std::vector mark_maximum(std::span marker, T theta, if ((theta < 0) || (theta > 1)) throw std::invalid_argument("Theta needs to fullfill 0 ≤ θ ≤ 1."); - T max = marker.empty() ? std::numeric_limits::min() + T max = marker.empty() ? std::numeric_limits::lowest() : std::ranges::max(marker); MPI_Allreduce(MPI_IN_PLACE, &max, 1, dolfinx::MPI::mpi_t, MPI_MAX, comm);