# A question about Function numerical setting

Hello everyone, while studying “Defining subdomains for different materials” in Defining subdomains for different materials — FEniCSx tutorial, I encountered an issue. The code is as follows:

``````from dolfinx import default_scalar_type
from dolfinx.fem import (Constant, dirichletbc, Function, functionspace, assemble_scalar,
form, locate_dofs_geometrical, locate_dofs_topological)
from dolfinx.mesh import create_unit_square, locate_entities, create_unit_cube
from mpi4py import MPI
import numpy as np

num_mesh = 10
mesh = create_unit_square(MPI.COMM_WORLD, num_mesh, num_mesh)
Q = functionspace(mesh, ("DG", 0))

def Omega_0(x):
return x[1] <= 0.5

def Omega_1(x):
return x[1] >= 0.5

kappa = Function(Q)
cells_0 = locate_entities(mesh, mesh.topology.dim, Omega_0)
print('cells_0', len(cells_0))
cells_1 = locate_entities(mesh, mesh.topology.dim, Omega_1)
print('cells_1', len(cells_1))
# print("kappa array before assignment:", kappa.x.array)
kappa.x.array[cells_0] = np.full_like(cells_0, 1, dtype=default_scalar_type)
kappa.x.array[cells_1] = np.full_like(cells_1, 0.1, dtype=default_scalar_type)
# print("kappa array after assignment:", kappa.x.array)

print('length of the array', len(kappa.x.array))
print('total number of cells', len(cells_0) + len(cells_1))
if np.any(kappa.x.array == 0):
print("Warning: Some cells have not been assigned values.")
``````

I found that when num_mesh is even, the values in kappa.x.array are all correctly set. However, when num_mesh is odd, some values in kappa.x.array remain 0 and are not set correctly. Then, I used gmsh to generate the mesh for testing, and the problem still persists. The elements in kappa.x.array are not all set to the correct values, which may cause errors in subsequent calculations. My test code is as follows:

``````import gmsh
import numpy as np
from dolfinx.io import XDMFFile, gmshio
from mpi4py import MPI
from dolfinx.mesh import locate_entities
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

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()

# 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()

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

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)

print('length of array', len(epsilon.x.array), MPI.COMM_WORLD.rank)
print('total', len(cells0) + len(cells1), MPI.COMM_WORLD.rank)

if len(cells0) + len(cells1) == len(epsilon.x.array):
print("All cells are included.")
else:
print("Some cells are missing or duplicated.")

epsilon.x.array[cells0] = np.full_like(cells0, 1, dtype=default_scalar_type)
epsilon.x.array[cells1] = np.full_like(cells1, 2, dtype=default_scalar_type)

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

Could you please help me identify the issue in my code?

`locate_entities` checks if all vertices satisfy the condition, e.g. the python function `volume0`.

If you are not careful in meshing, i.e. if you do not draw an horizontal line at `x[1] == split`, there will be a layer of cells which are neither in `volume0` nor in `volume1`, because some vertices satisfy `volume0` but not `volume1`, and the other vertices satisfy `volume1` but not `volume1`. This is easily seen with your odd vs even setting.

If you use gmsh, be sure to set the subdomains there to ensure that the mesh generates the desired horizontal line.

1 Like

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:

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:

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

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()

# 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

Thank you both for your responses. I now understand where my problem lies. In my previous code, some cells were determined to be neither in volume0 nor in volume1. Therefore, if I want to divide the space into two parts, I need to set a boundary at the split when generating the mesh.

You are correct, my issue has been resolved. Thank you very much for your reply.