How to make sure the subdomains are connected in built-in mesh?

General
1

Hi all,

I am modifying the this example to a geometry with an extra subdomain on the upper half. The subdomain definitions are

Omega = {
    0 : lambda x: x[1] <= 0.5, # lower half of the mesh
    # a block at the middle on the upper half of the mesh
    1 : lambda x : np.logical_and(
        np.logical_and(x[0] >= 1/3, x[0] <= 2/3),
        x[1] >= 0.75
    ),
}
# the remaining part in the upper mesh
Omega[2] = lambda x : np.logical_and(x[1] >= 0.5, np.logical_not(Omega[1](x)))

I assigned the same \kappa value to subdomain 1 and 2 so that the expect solution should be the same as if there are only 2 subdomains. However, the solution is Nan. I have tried chopping the upper mesh into half and solution is valid. So, what is the best practice to implement subdomains to make sure they are connected? I appreciate for the help. Below if the script to reproduce the issue.

import numpy as np
import ufl
from mpi4py import MPI
from dolfinx import mesh, fem
from dolfinx.fem import FunctionSpace
from petsc4py.PETSc import ScalarType

domain = mesh.create_unit_square(MPI.COMM_WORLD, 30, 30)
Q = FunctionSpace(domain, ('DG', 0))

Omega = {
    0 : lambda x: x[1] <= 0.5,
    # this division works
    # 1 : lambda x : np.logical_and(x[1] >= 0.5, x[0] <= 0.5), 
    # 2 : lambda x : np.logical_and(x[1] >= 0.5, x[0] >= 0.5),
    1 : lambda x : np.logical_and(
        np.logical_and(x[0] >= 1/3, x[0] <= 2/3),
        x[1] >= 0.75
    ),
}
Omega[2] = lambda x : np.logical_and(x[1] >= 0.5, np.logical_not(Omega[1](x)))

cells = {
    tag : mesh.locate_entities(domain, domain.topology.dim, omega) for tag, omega in Omega.items() 
}

domain_indices = [cell for _, cell in cells.items()]
domain_markers = [np.full_like(cell, tag) for tag, cell in cells.items()]
domain_indices = np.hstack(domain_indices).astype(np.int32)
domain_markers = np.hstack(domain_markers).astype(np.int32)
sorted_domain = np.argsort(domain_indices)
subdomain_tag = mesh.meshtags(domain, domain.topology.dim, domain_indices[sorted_domain], domain_markers[sorted_domain])
dx = ufl.Measure('dx', domain=domain, subdomain_data=subdomain_tag)


Q = FunctionSpace(domain, ('DG',0))
kappa = fem.Function(Q)

kappa_values = {
    0 : 1,
    1 : 0.1,
    2 : 0.1,
}

for tag, value in kappa_values.items():
    cell = cells[tag]
    kappa.x.array[cell] = np.full_like(cell, value, dtype=ScalarType)


V = FunctionSpace(domain, ('CG', 1))
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
x = ufl.SpatialCoordinate(domain)

s = fem.Constant(domain, ScalarType(1.))
a = kappa * ufl.inner(ufl.grad(u), ufl.grad(v)) * dx
L = ufl.inner(s, v) * dx

dofs = fem.locate_dofs_geometrical(V, lambda x : np.isclose(x[0], 0))
bcs = [fem.dirichletbc(ScalarType(1.), dofs, V)]

problem = fem.petsc.LinearProblem(a, L, bcs=bcs, petsc_options={'ksp_type': 'preonly', 'pc_type': 'lu'})
uh = problem.solve()

import pyvista
from dolfinx.plot import create_vtk_mesh
pyvista.start_xvfb()
print(pyvista.global_theme.jupyter_backend)


num_cells_local = domain.topology.index_map(domain.topology.dim).size_local
marker = np.zeros(num_cells_local, dtype=np.int32)

for tag, cell in cells.items():
    cell = cell[cell < num_cells_local]
    marker[cell] = tag

topology, cell_types, geometry = create_vtk_mesh(domain, domain.topology.dim, np.arange(num_cells_local, dtype=np.int32))


grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)
grid.cell_data['Marker'] = marker
grid.set_active_scalars('Marker')

plotter = pyvista.Plotter()
plotter.add_mesh(grid, show_edges=True)
plotter.view_xy()
if not pyvista.OFF_SCREEN:
    plotter.show()

u_topology, u_cell_types, u_geometry = create_vtk_mesh(V)
u_grid = pyvista.UnstructuredGrid(u_topology, u_cell_types, u_geometry)
u_grid.point_data['u'] = uh.x.array
u_grid.set_active_scalars('u')
u_plotter = pyvista.Plotter()
u_plotter.add_mesh(u_grid, show_edges=True)
u_plotter.view_xy()
if not pyvista.OFF_SCREEN:
    u_plotter.show()
2

So just to let me clarify, you want to split the domain into three separate regions, that doesn’t necessarily align with your grid. If you want to do so, I would start by marking all cells with a given value:

import numpy as np
from mpi4py import MPI
from dolfinx import mesh, fem
from dolfinx.fem import FunctionSpace
from petsc4py.PETSc import ScalarType

from dolfinx.io import XDMFFile


domain = mesh.create_unit_square(MPI.COMM_WORLD, 30, 30)
Q = FunctionSpace(domain, ('DG', 0))

Omega = {
    0 : lambda x: x[1] <= 0.5,
    1 : lambda x : np.logical_and(
        np.logical_and(x[0] >= 1/3, x[0] <= 2/3),
        x[1] >= 0.75
    ),
}
Omega[2] = lambda x : np.logical_and(x[1] >= 0.5, np.logical_not(Omega[1](x)))

cell_map = domain.topology.index_map(domain.topology.dim)
all_cells = np.arange(cell_map.size_local+cell_map.num_ghosts, dtype=np.int32)
# Start by marking all cells with 1
values = np.full_like(all_cells,1, dtype=np.int32)
# Mark 0-value cells
dofs_0 = mesh.locate_entities(domain, domain.topology.dim, Omega[0])
values[dofs_0] = np.full_like(dofs_0, 0, dtype=np.int32) 
# Mark 2-value cells
dofs_2 = mesh.locate_entities(domain, domain.topology.dim, Omega[2])
values[dofs_2] = np.full_like(dofs_2, 2, dtype=np.int32) 


subdomain_tag = mesh.meshtags(domain, domain.topology.dim, all_cells, values)

with XDMFFile(MPI.COMM_WORLD, "cell_tag.xdmf", "w") as xdmf:
    xdmf.write_mesh(domain)
    xdmf.write_meshtags(subdomain_tag, domain.geometry)


Q = FunctionSpace(domain, ('DG',0))
kappa = fem.Function(Q)
kappa_values = {
    0 : 1,
    1 : 0.1,
    2 : 0.3,
}

for tag, value in kappa_values.items():
    cell = subdomain_tag.find(tag)
    kappa.x.array[cell] = np.full_like(cell, value, dtype=ScalarType)
with XDMFFile(MPI.COMM_WORLD, "kappa.xdmf", "w") as xdmf:
    xdmf.write_mesh(domain)
    xdmf.write_function(kappa)

If you make the domains align with the grid, you won’t see a band of cells that doesn’t satisfy either criteria.

1 Like