Add maximum marking

[schnellerhase]

Implements a maximum marking criterion for AFEM schemes, i.e. for $\theta \in [0, 1]$ and a marker $\eta \in \mathbb{R}^N$, computes

$$ \mathcal{I} = \{ i \ | \ \eta_i \geq \theta \max_j \eta_j \}. $$

In particular this introduces/defines the interface for other marking routines (to follow in the future): based on a process local list of (entity-)indicators the (local) marked index set of it is computed.

Test program to visualise results (parallel ready).

from mpi4py import MPI

import numpy as np

from dolfinx import mesh, plot
import pyvista as pv

mesh_0 = mesh.create_unit_square(comm := MPI.COMM_WORLD, n := 10, n)

cell_im = mesh_0.topology.index_map(tdim := mesh_0.topology.dim)
cell_count = cell_im.size_local + cell_im.num_ghosts

marker = np.random.random(cell_count)
marked_cells = mesh.mark_maximum(marker, theta := 0.5, comm)

mesh_0.topology.create_entities(1)
marked_edges = mesh.compute_incident_entities(mesh_0.topology, marked_cells, tdim, 1)

mesh_1, _, _ = mesh.refine(mesh_0, marked_edges)

grid_0 = pv.UnstructuredGrid(*plot.vtk_mesh(mesh_0))
grid_1 = pv.UnstructuredGrid(*plot.vtk_mesh(mesh_1))

plotter = pv.Plotter(shape=(1, 2))
plotter.subplot(0, 0)
plotter.add_mesh(grid_0, show_edges=True)
plotter.view_xy()
plotter.add_text("mesh_0", position="upper_left")

plotter.subplot(0, 1)
plotter.add_mesh(grid_1, show_edges=True)
plotter.view_xy()
plotter.add_text("mesh_1", position="upper_left")

plotter.link_views()
plotter.show()

Work towards #4141.

[jhale]

Use const T instead of auto e?

Also remove inequality, pointless for floating point arithmetic.

[schnellerhase]

Changed to T - think we do not follow const correctness for trivial types.

[schnellerhase]

Also remove inequality, pointless for floating point arithmetic.

Don't think this is correct: for $\theta=1.0$ and $\eta = [ 1, \dots, 1]$ that would then mark nothing (tested).

[jhale]

And through what astronomically unlikely event would eta be binary 1.0 for every cell?

[schnellerhase]

Unlikely but valid input to the function. We can either leave as is or change to $e + \epsilon > \theta \max \eta$ for some close to machine epsilon $\epsilon$, which should ensure this edge case passes. Fine with both.

outdated

[jhale]

My comments on the floating point comparisons are not resolved.

[chrisrichardson]

@jhale - I think @schnellerhase did address it. Did you mean that it should be a simple > rather than >=?

[jhale]

Yes - the equality is completely redundant in a floating point context and the example given is a contrived one.