Consider the following code, that avoids loops in Python by using vectorized numpy operations and DOLFINx functions (written in C++)
import dolfinx
from mpi4py import MPI
import numpy as np
def cell_criterion(x):
"""Given mesh coordinates, return if each point
satisfies x[0]<0.3 or x[1]>0.4
:param x: Input coordinates, shape (num_points, 3)
:returns: Boolean array of shape (num_points, )
"""
return (x[0] < 0.3) | (x[1] > 0.4)
mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)
# Compute midpoints for all cells on process
cells_local = np.arange(mesh.topology.index_map(
mesh.topology.dim).size_local, dtype=np.int32)
midpoints = dolfinx.mesh.compute_midpoints(
mesh, mesh.topology.dim, cells_local).T
# Check midpoint criterion and find edges connected to cells
should_refine = np.flatnonzero(cell_criterion(midpoints)).astype(np.int32)
mesh.topology.create_entities(1)
local_edges = dolfinx.mesh.compute_incident_entities(
mesh, should_refine, mesh.topology.dim, 1)
refined_mesh = dolfinx.mesh.refine(mesh, local_edges)
with dolfinx.io.XDMFFile(refined_mesh.comm, "refined_mesh.xdmf", "w") as xdmf:
xdmf.write_mesh(refined_mesh)
giving the mesh