Mesh refinement in dolfinx

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
image

2 Likes