Mesh refinement in dolfinx

dolfinx
1

Hi everyone,

I need help translating some legacy dolfin code to dolfinx (0.6.0). The code is about marking cells depending on the value of the cell midpoint of function cn1_ and then refine all cells that are marked. This was my approach in legacy dolfin:

cell_markers = MeshFunction("bool", mesh_initial, mesh_initial.topology().dim())
            cell_markers.set_all(False)
            for cell in cells(mesh_initial):
                if cn1_(cell.midpoint()) > 0.01 and cn1_(cell.midpoint()) < 0.99:
                    cell_markers[cell] = True
            mesh = refine(mesh_initial, cell_markers)

any help is greatly appreciated.

Best regards,
Sven

2

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
3

This approach gives the following error in version v.0.8.0:

Traceback (most recent call last):
  File "/root/shared/refinement_dolfinx_test.py", line 27, in <module>
    local_edges = dolfinx.mesh.compute_incident_entities(
  File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/mesh.py", line 216, in compute_incident_entities
    return _cpp.mesh.compute_incident_entities(topology, entities, d0, d1)
TypeError: compute_incident_entities(): incompatible function arguments. The following argument types are supported:
    1. compute_incident_entities(mesh: dolfinx.cpp.mesh.Topology, entities: ndarray[dtype=int32, writable=False, shape=(*), order='C'], d0: int, d1: int) -> numpy.ndarray[dtype=int32, ]

Invoked with types: dolfinx.mesh.Mesh, ndarray, int, int
4

Please read the error message carefully, as it was surely giving an hint on how to proceed: the first input argument is a mesh, while it was supposed to be a Topology object. mesh.topology is used all over the place in the code, so a simple guess would have been to replace mesh with mesh.topology.

The updated code is (changes highlighted with my initials FB)

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)
mesh.topology.create_connectivity(mesh.topology.dim, mesh.topology.dim)  # FB: ADDED HERE
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.topology, should_refine, mesh.topology.dim, 1)  # FB: CHANGED HERE
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)
3 Likes
5

I appreciate the help, thank you!