Hi Dolfinx community,
We are attempting to create a function that defines where a domain mesh and a indicator mesh intersect for use a fictitious domain formulation. We have been create this function by building a mask and assigning a value base of the bool value of the mask. This works well in serial but is failing in parallel.
I assume that since both the domain mesh and the fictitious domain mesh are both divided in amongst the processes, some points are getting marked as colliding correctly. I have scattered the function that we assign using the mesh forward, but I suspect the problem lies in creating the mask itself (within the indicator
function).
Thanks!
MWE:
import dolfinx
import numpy as np
from mpi4py import MPI
from dolfinx import fem, geometry
# Define rect mesh
mesh = dolfinx.mesh.create_rectangle(MPI.COMM_WORLD, [[0.0, 0.0], [10.0, 10.0]], [100, 100], dolfinx.mesh.CellType.triangle)
fict_dom_mesh = dolfinx.mesh.create_rectangle(MPI.COMM_WORLD, [[0.0, 0.0], [6.0, 6.0]], [20, 20], dolfinx.mesh.CellType.triangle)
# Create a function space on the mesh
V = fem.functionspace(mesh, ('CG', 1))
func = fem.Function(V)
# Define a mask for the indicator function
bb_tree = geometry.bb_tree(fict_dom_mesh, fict_dom_mesh.topology.dim)
def indicator(points):
pts = np.ascontiguousarray(points, dtype=np.float64)
pts.setflags(write=False)
# Compute candidates and then actual colliding cells
candidates = geometry.compute_collisions_points(bb_tree, pts)
colliding = geometry.compute_colliding_cells(fict_dom_mesh, candidates, pts)
# Determine which points lie inside any cell
mask = np.array([len(colliding.links(i)) > 0 for i in range(points.shape[0])], dtype=bool)
return mask
# Local dof coordinates on this MPI rank
dof_coords = V.tabulate_dof_coordinates()
local_size = V.dofmap.index_map.size_local
# Evaluate indicator function only on local dofs
local_points = dof_coords[:local_size]
mask = indicator(local_points)
# Assign values
local_values = np.full(local_size, 0.0, dtype=np.float64)
local_values[mask] = 1.0
func.x.array[:local_size] = local_values
func.x.scatter_forward()
# Save the function to a file
with dolfinx.io.VTKFile(MPI.COMM_WORLD, 'output/indicator_field.pvd', 'w') as f:
f.write_function(func)