A question about Function numerical setting

dolfinx
,
3

Hi yonggangli,

As @francesco-ballarin points out, the marker functions you use here will be problematic if the cells don’t align with the split coordinate. A good first step when dealing with problems where things are marked is plotting the mesh with the meshtags of the markings. Here is what I get when using your code:

mesh_tags_missing_cells
mesh_tags_missing_cells1013×706 91.7 KB

I have tagged volume0 with 2 and volume1 with 3, leaving any untagged cells with the tag 1. As you can see, the middle cells of the mesh are not tagged.

If you instead compute the midpoints of each cell and use these in the locator function, you would get this:

mesh_tags_all_cells
mesh_tags_all_cells990×685 93.4 KB

Note that this still won’t produce a split that is a horizontal line, since the cell midpoints are not aligned with the split coordinate. To achieve an even split, you would have to address this issue while meshing, as pointed out by @francesco-ballarin.

In case it is of interest, here is your code with the changes I made to mark all the cells using their midpoints as reference:

import numpy as np
from dolfinx.io import XDMFFile, gmshio
from mpi4py import MPI
from dolfinx.mesh import locate_entities, meshtags, compute_midpoints
from dolfinx import fem, default_scalar_type


len_side = 50
lc = 8
split = len_side / 2


def get_device():
    # Initialize Gmsh
    gmsh.initialize()

    model = gmsh.model()

    # Create a new model
    model.add("Device")

    p1 = model.occ.addPoint(0, 0, 0, lc, 1)
    p2 = model.occ.addPoint(len_side, 0, 0, lc, 2)
    p3 = model.occ.addPoint(len_side, len_side, 0, lc, 3)
    p4 = model.occ.addPoint(0, len_side, 0, lc, 4)

    line1 = model.occ.addLine(p1, p2, 1)
    line2 = model.occ.addLine(p2, p3, 2)
    line3 = model.occ.addLine(p3, p4, 3)
    line4 = model.occ.addLine(p4, p1, 4)

    loop = model.occ.addCurveLoop([line1, line2, line3, line4], 5)
    surface = model.occ.addPlaneSurface([loop], 1)

    model.occ.synchronize()

    model.add_physical_group(dim=2, tags=[surface])

    # Generate the 3D mesh
    model.mesh.generate(2)
    return model


def volume0(x):
    return x[1] >= split


def volume1(x):
    return x[1] <= split


model = get_device()
gdim  = 2

# Interfacing with dolfinx
mesh, cell_markers, facet_markers = gmshio.model_to_mesh(model, MPI.COMM_WORLD, 0, gdim=gdim)

Q = fem.functionspace(mesh, ("DG", 0))
epsilon = fem.Function(Q)

cells0 = locate_entities(mesh, mesh.topology.dim, volume0)
print('cells1', len(cells0), MPI.COMM_WORLD.rank)
cells1 = locate_entities(mesh, mesh.topology.dim, volume1)
print('cells2', len(cells1), MPI.COMM_WORLD.rank)

num_cells = mesh.topology.index_map(gdim).size_local +  + mesh.topology.index_map(gdim).num_ghosts
mids = compute_midpoints(mesh, mesh.topology.dim, np.arange(num_cells, dtype=np.int32))
locator0 = volume0(mids.T)
locator1 = volume1(mids.T)
cell_tags = np.full(num_cells, 1, dtype=np.int32)
cell_tags[locator0] = 2
cell_tags[locator1] = 3
tags = meshtags(mesh, gdim, np.arange(num_cells), cell_tags)

with XDMFFile(MPI.COMM_WORLD, "tags.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    xdmf.write_meshtags(tags, mesh.geometry)

print('length of array', len(epsilon.x.array), MPI.COMM_WORLD.rank)
print('total', np.sum(locator0) + np.sum(locator1), MPI.COMM_WORLD.rank)

if np.sum(locator0) + np.sum(locator1) == len(epsilon.x.array):
    print("All cells are included.")
else:
    print("Some cells are missing or duplicated.")

dofs0 = fem.locate_dofs_topological(Q, gdim, np.where(cell_tags==2)[0])
dofs1 = fem.locate_dofs_topological(Q, gdim, np.where(cell_tags==3)[0])

epsilon.x.array[dofs0] = np.full(np.sum(locator0), 1, dtype=default_scalar_type)
epsilon.x.array[dofs1] = np.full(np.sum(locator1), 2, dtype=default_scalar_type)

if np.any(epsilon.x.array == 0):
    print("Warning: Some cells have not been assigned values.")```

Cheers
1 Like