Dolfinx: Mesh refinement

mesh
4

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