A question about Function numerical setting

dolfinx
,
1

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

# 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?

2

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
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
4

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.

5

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