Add maximum marking

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

cpp/dolfinx/refinement/mark.h

@@ -0,0 +1,61 @@
+// 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 <algorithm>
+#include <cassert>
+#include <concepts>
+#include <cstdint>
+#include <iterator>
+#include <limits>
+#include <mpi.h>
+#include <spdlog/spdlog.h>
+#include <vector>
+
+#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::floating_point T>
+std::vector<std::int32_t> mark_maximum(std::span<const T> marker, 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<T>::lowest()
+                         : std::ranges::max(marker);
+  MPI_Allreduce(MPI_IN_PLACE, &max, 1, dolfinx::MPI::mpi_t<T>, MPI_MAX, comm);
+
+  auto mark = [=](T e)
+  { return e + std::numeric_limits<T>::epsilon() * 1e2 > theta * max; };
+
+  std::vector<std::int32_t> indices;
+  indices.reserve(std::ranges::count_if(marker, mark));
+
+  for (std::int32_t i = 0; i < static_cast<std::int32_t>(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

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
 )

cpp/test/mesh/refinement/mark.cpp

@@ -0,0 +1,61 @@
+// 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 <algorithm>
+#include <catch2/catch_template_test_macros.hpp>
+#include <dolfinx/common/MPI.h>
+#include <dolfinx/refinement/mark.h>
+#include <mpi.h>
+#include <vector>
+
+using namespace dolfinx;
+using namespace dolfinx::refinement;
+
+TEMPLATE_TEST_CASE("Mark maximum empty", "[refinement][mark][maximum]", double,
+                   float)
+{
+  std::vector<TestType> marker;
+  auto indices = mark_maximum<TestType>(marker, .5, MPI_COMM_WORLD);
+  CHECK(indices.size() == 0);
+}
+
+TEMPLATE_TEST_CASE("Mark maximum ones", "[refinement][mark][maximum]", double,
+                   float)
+{
+  std::vector<TestType> marker(10, 1.0);
+  auto indices = mark_maximum<TestType>(marker, 1.0, MPI_COMM_WORLD);
+  CHECK(indices.size() == 10);
+}
+
+TEMPLATE_TEST_CASE("Mark maximum", "[refinement][mark][maximum]", double, float)
+{
+  MPI_Comm comm = MPI_COMM_WORLD;
+
+  std::vector<TestType> 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_maximum<TestType>(marker, theta, comm);
+
+  CHECK(std::ranges::all_of(
+      indices, [&](auto e)
+      { return (0 <= e) && (e <= static_cast<std::int32_t>(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<std::int32_t>(indices.size()));
+
+  for (std::int32_t i = 0; i < static_cast<std::int32_t>(marker.size()); ++i)
+  {
+    bool expect_marked = mark(marker[i]);
+    bool marked = std::ranges::find(indices, i) != indices.end();
+    CHECK(expect_marked == marked);
+  }
+}

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",

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 <concepts>
 #include <dolfinx/mesh/Mesh.h>
+#include <dolfinx/refinement/mark.h>
 #include <dolfinx/refinement/option.h>
 #include <dolfinx/refinement/refine.h>
 #include <dolfinx/refinement/uniform.h>
@@ -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<const T, nb::ndim<1>, 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

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)