Dolfinx: Mesh refinement

I want to refine my mesh (multiple times) based on the location of the cell. The reason is, that I need to have a fine mesh towards the boundary of the domain. Since my problem is complex, I am using dolfinx.

However, I cannot get it to run. There are some examples for dolfin but I did not spot one for dolfinx. My code so far is. I would really appreciate any help on that.

import dolfinx
import numpy as np
from dolfinx import RectangleMesh
from dolfinx.fem import locate_dofs_geometrical, locate_dofs_topological
from dolfinx.io import XDMFFile
from dolfinx.mesh import locate_entities_boundary, refine

W = 380e-6  # width
H = 180e-6  # height

d_ele = 20e-6
delta = 1e-7

n_W = int(np.ceil(W/d_ele))  # number of elements along W
n_H = int(np.ceil(H/d_ele))  # number of elements along H

# Create mesh
mesh = RectangleMesh(MPI.COMM_WORLD,
                     [np.array([-W/2, -H/2, 0]), np.array([W/2, H/2, 0])],
                     [n_W, n_H],
                     CellType.triangle, dolfinx.cpp.mesh.GhostMode.none)


def inside_delta(xs):
    out = []
    for i in np.arange(xs.shape[1]):
        x, y = xs[0, i], xs[1, i]
        if W/2 - abs(x) <= 30*delta or H/2 - abs(y) <= 30*delta:
            out.append(True)
        else:
            out.append(False)
    return np.asarray(out)

def mark_all(xs):
    return np.ones(xs.shape[1])

# boundary layer
boundaries = [(1, lambda x: inside_delta(x)),
        (0, lambda x: mark_all(x))]

refine_indices, refine_markers = [], []
dim = mesh.topology.dim

for (marker, locator) in boundaries:
    facets = dolfinx.mesh.locate_entities(mesh, dim, locator)
    refine_indices.append(facets)
    refine_markers.append(np.full(len(facets), marker))

refine_tag = dolfinx.MeshTags(mesh, dim,
                              np.array(np.hstack(refine_indices), dtype=np.int8),
                              np.array(np.hstack(refine_markers), dtype=np.int8))

mesh.topology.create_entities(1)
mesh = refine(mesh, refine_tag)

# in order to export it you have to change dtype=np.int8 to dtype=np.int32
with XDMFFile(mesh.mpi_comm(), "mesh.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)

I could partly resolve my problem by myself. I thought I need to mark every cell; even the once which should not be refined. But this is not needed.

import dolfinx
import numpy as np
from dolfinx import RectangleMesh
from dolfinx.fem import locate_dofs_geometrical, locate_dofs_topological
from dolfinx.io import XDMFFile
from dolfinx.mesh import locate_entities_boundary, refine

W = 380e-6  # width
H = 180e-6  # height

d_ele = 20e-6
delta = 1e-6

n_W = int(np.ceil(W/d_ele))  # number of elements along W
n_H = int(np.ceil(H/d_ele))  # number of elements along H

# Create mesh
mesh = RectangleMesh(MPI.COMM_WORLD,
                     [np.array([-W/2, -H/2, 0]), np.array([W/2, H/2, 0])],
                     [n_W, n_H],
                     CellType.triangle, dolfinx.cpp.mesh.GhostMode.none)


def inside_delta(xs):
    out = []
    for i in np.arange(xs.shape[1]):
        x, y = xs[0, i], xs[1, i]
        if W/2 - abs(x) <= 30*delta or H/2 - abs(y) <= 30*delta:
            out.append(True)
        else:
            out.append(False)
    return np.asarray(out)

# boundary layer
boundaries = [(1, lambda x: inside_delta(x))]

for _ in np.arange(10):
    refine_indices, refine_markers = [], []
    dim = mesh.topology.dim

    for (marker, locator) in boundaries:
        facets = dolfinx.mesh.locate_entities(mesh, dim, locator)
        refine_indices.append(facets)
        refine_markers.append(np.full(len(facets), marker))

    refine_tag = dolfinx.MeshTags(mesh, dim,
                                 np.array(np.hstack(refine_indices), dtype=np.int8),
                                  np.array(np.hstack(refine_markers), dtype=np.int8))

   mesh.topology.create_entities(1)
   mesh = refine(mesh, refine_tag)

# in order to export it you have to change dtype=np.int8 to dtype=np.int32
with XDMFFile(mesh.mpi_comm(), "mesh.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)

However, since I want to do that process repeatedly I get now the following error code already in the 2nd loop

[0]PETSC ERROR: ------------------------------------------------------------------------
[0]PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation, probably memory access out of range
[0]PETSC ERROR: Try option -start_in_debugger or -on_error_attach_debugger
[0]PETSC ERROR: or see https://www.mcs.anl.gov/petsc/documentation/faq.html#valgrind
[0]PETSC ERROR: or try http://valgrind.org on GNU/linux and Apple Mac OS X to find memory corruption errors
[0]PETSC ERROR: configure using --with-debugging=yes, recompile, link, and run
[0]PETSC ERROR: to get more information on the crash.
[0]PETSC ERROR: Run with -malloc_debug to check if memory corruption is causing the crash.
application called MPI_Abort(MPI_COMM_WORLD, 50176059) - process 0
[unset]: write_line error; fd=-1 buf=:cmd=abort exitcode=50176059
:
system msg for write_line failure : Bad file descriptor

I appreciate any help regarding this issue.

This indeed looks like a bug. I’ll have a closer look tomorrow.

1 Like

So the issue is that

is wrong. here, only the markers should be int8, not the indices. Note that you should also sort the indices, to ensure that meshtags works as it should. I’ve added a fully working code below:


import numpy as np
from dolfinx import RectangleMesh, MeshTags
from dolfinx.io import XDMFFile
from dolfinx.mesh import refine, locate_entities
from dolfinx.cpp.mesh import GhostMode, CellType
from mpi4py import MPI


W = 380e-6  # width
H = 180e-6  # height

d_ele = 20e-6
delta = 1e-6


n_W = int(np.ceil(W/d_ele))  # number of elements along W
n_H = int(np.ceil(H/d_ele))  # number of elements along H

# Create mesh
mesh = RectangleMesh(MPI.COMM_WORLD,
                     [np.array([-W/2, -H/2, 0]), np.array([W/2, H/2, 0])],
                     [n_W, n_H],
                     CellType.triangle, GhostMode.none)


def inside_delta(xs):
    out = []
    for i in np.arange(xs.shape[1]):
        x, y = xs[0, i], xs[1, i]
        if W/2 - abs(x) <= 30*delta or H/2 - abs(y) <= 30*delta:
            out.append(True)
        else:
            out.append(False)
    return np.asarray(out)


def mark_all(xs):
    return np.ones(xs.shape[1])


# boundary layer
boundaries = [(1, lambda x: inside_delta(x))]

for _ in np.arange(5):
    refine_indices, refine_markers = [], []
    dim = mesh.topology.dim
    for (marker, locator) in boundaries:
        facets = locate_entities(mesh, dim, locator)
        refine_indices.append(facets)
        refine_markers.append(np.full(len(facets), marker))

    # Only markers should be int8
    refine_indices = np.array(np.hstack(refine_indices), dtype=np.int32)
    refine_markers = np.array(np.hstack(refine_markers), dtype=np.int8)
    # Indices sent into meshtags should be sorted
    sort = np.argsort(refine_indices)
    refine_tag = MeshTags(
        mesh, dim, refine_indices[sort], refine_markers[sort])

    out_tag = MeshTags(
        mesh, dim, refine_indices[sort], np.asarray(refine_markers[sort], dtype=np.int32))

    with XDMFFile(mesh.mpi_comm(), f"mesh_{_}.xdmf", "w") as xdmf:
        xdmf.write_mesh(mesh)
        xdmf.write_meshtags(out_tag)

    mesh.topology.create_entities(1)
    mesh = refine(mesh, refine_tag)

with XDMFFile(mesh.mpi_comm(), f"mesh_{_+1}.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
3 Likes

The posted fixed example by @dokken no longer works, because the dolfinx API has changed.

The following code works with newer dolfinx versions (was tested with dolfinx 0.5.1)

import numpy as np
import pyvista
from dolfinx.mesh import CellType, create_rectangle, locate_entities, refine
from dolfinx.plot import create_vtk_mesh
from mpi4py import MPI

W = 380e-6  # width
H = 180e-6  # height

d_ele = 20e-6
delta = 1e-6


n_W = int(np.ceil(W / d_ele))  # number of elements along W
n_H = int(np.ceil(H / d_ele))  # number of elements along H

# Create mesh
mesh = create_rectangle(
    MPI.COMM_WORLD,
    points=((-W / 2, -H / 2), (W / 2, H / 2)),
    n=(n_W, n_H),
    cell_type=CellType.triangle,
)


def inside_delta(xs):
    # refine around the boundary
    return np.logical_or(
        W / 2 - abs(xs[0, :]) <= 30 * delta,
        H / 2 - abs(xs[1, :]) <= 30 * delta,
    )


for _ in np.arange(5):
    dim = mesh.topology.dim
    edges = locate_entities(mesh, dim - 1, inside_delta)
    mesh.topology.create_entities(1)
    mesh = refine(mesh, edges, redistribute=False)

    vtkdata = create_vtk_mesh(mesh, mesh.topology.dim)
    plotter = pyvista.Plotter()
    grid = pyvista.UnstructuredGrid(*vtkdata)
    actor = plotter.add_mesh(grid, show_edges=True)
    plotter.view_xy()
    plotter.show()

and it generates plots of the refined mesh

4 Likes