Hello everyone,
Recently, I am running my code parallelly, and find a question about the cell indices in each processor. Let’s use the subdomain case in FEniCSx tutorial as an example Defining subdomains for different materials . Here is the code. The only thing I did, is just print the cells_0 and cells_1 out with which processor they are in.
import gmsh
import numpy as np
import pyvistafrom dolfinx.fem import (Constant, dirichletbc, Function, FunctionSpace, assemble_scalar,
form, locate_dofs_geometrical, locate_dofs_topological)
from dolfinx.fem.petsc import LinearProblem
from dolfinx.io import XDMFFile
from dolfinx.mesh import create_unit_square, locate_entities
from dolfinx.plot import create_vtk_meshfrom ufl import (SpatialCoordinate, TestFunction, TrialFunction,
dx, grad, inner)from mpi4py import MPI
from petsc4py.PETSc import ScalarTypemesh = create_unit_square(MPI.COMM_WORLD, 10, 10)
Q = FunctionSpace(mesh, (“DG”, 0))def Omega_0(x):
return x[1] <= 0.5def Omega_1(x):
return x[1] >= 0.5kappa = Function(Q)
cells_0 = locate_entities(mesh, mesh.topology.dim, Omega_0)
cells_1 = locate_entities(mesh, mesh.topology.dim, Omega_1)kappa.x.array[cells_0] = np.full_like(cells_0, 1, dtype=ScalarType)
kappa.x.array[cells_1] = np.full_like(cells_1, 0.1, dtype=ScalarType)imap = kappa.function_space.dofmap.index_map
print(“global size”, imap.size_global)
print(“processor”, MPI.COMM_WORLD.rank, “local size”, imap.size_local)V = FunctionSpace(mesh, (“CG”, 1))
u, v = TrialFunction(V), TestFunction(V)
a = inner(kappa*grad(u), grad(v)) * dx
x = SpatialCoordinate(mesh)
L = Constant(mesh, ScalarType(1)) * v * dx
dofs = locate_dofs_geometrical(V, lambda x: np.isclose(x[0], 0))
bcs = [dirichletbc(ScalarType(1), dofs, V)]problem = LinearProblem(a, L, bcs=bcs, petsc_options={“ksp_type”: “preonly”, “pc_type”: “lu”})
uh = problem.solve()pyvista.set_jupyter_backend(“pythreejs”)
p = pyvista.Plotter(window_size=[800, 800], shape=(1,2))
Filter out ghosted cells
num_cells_local = mesh.topology.index_map(mesh.topology.dim).size_local
marker = np.zeros(num_cells_local, dtype=np.int32)
print(“processor”, MPI.COMM_WORLD.rank, “cells 0”, cells_0)
print(“processor”, MPI.COMM_WORLD.rank, “cells 1”, cells_1)
print(“processor”, MPI.COMM_WORLD.rank, “number of local cells”, num_cells_local)
cells_0 = cells_0[cells_0<num_cells_local]
cells_1 = cells_1[cells_1<num_cells_local]
marker[cells_0] = 1
marker[cells_1] = 2
topology, cell_types, x = create_vtk_mesh(mesh, mesh.topology.dim, np.arange(num_cells_local, dtype=np.int32))
grid = pyvista.UnstructuredGrid(topology, cell_types, x)
grid.cell_data[“Marker”] = marker
grid.set_active_scalars(“Marker”)
p.subplot(0,0)
actor0 = p.add_mesh(grid, show_edges=True)
p.subplot(0,1)
grid_uh = pyvista.UnstructuredGrid(*create_vtk_mesh(V))
grid_uh.point_data[“u”] = uh.x.array.real
grid_uh.set_active_scalars(“u”)
actor1 = p.add_mesh(grid_uh, show_edges=True)
if not pyvista.OFF_SCREEN:
p.show()
else:
pyvista.start_xvfb()
figure = p.screenshot(“subdomains_structured.png”)
Here is the result. We can see that, on each processor, the local size of cess are both 100, but in actual the numbers of the output cells are both 109.
My question is that, what are these cells whose indices are larger than 100? Are them ghost cells? But I think in this case, there should not be ghost cells. Or if we delete them, will it affect the accurancy of our result?
Thanks a lot for your time and help sincerely!