Fictitious domain with two meshes in MPI

dolfinx
,
1

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.

img-serial
img-serial1501×782 13.5 KB

img-mpi
img-mpi1501×782 24.1 KB

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

There are several problems here.

  1. If the fictitious domain mesh is distributed across multiple processes, you need a global gathering to inform all processes about the extent of this mesh.
  2. You only work on local data in this case, thus excluding all ghosts.

One way to by-pass this, is to let the ficticious domain exist on all processes:

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_SELF, [[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()

# Evaluate indicator function only on local dofs
local_points = dof_coords
mask = indicator(local_points)

# Assign values
local_values = np.full(dof_coords.shape[0], 0.0, dtype=np.float64)
local_values[mask] = 1.0

func.x.array[:] = 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)

A more elegant solution would be to use non-matching interpolation, where both meshes can be distributed:

import dolfinx
import numpy as np
from mpi4py import MPI
from dolfinx import fem

# 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)
V = fem.functionspace(mesh, ('CG', 1))


# Create marker function on fictitious domain mesh
Q = fem.functionspace(fict_dom_mesh, ("CG", 1))
q = fem.Function(Q)
q.x.array[:] = 1

cell_map = mesh.topology.index_map(mesh.topology.dim)
num_cells_on_proc = cell_map.size_local + cell_map.num_ghosts
cells = np.arange(num_cells_on_proc, dtype=np.int32)
interpolation_data = dolfinx.fem.create_interpolation_data(V, Q, cells, padding=1e-8)

func = fem.Function(V)
func.interpolate_nonmatching(q, cells, interpolation_data=interpolation_data)

# Save the function to a file
with dolfinx.io.VTKFile(MPI.COMM_WORLD, 'output/indicator_field.pvd', 'w') as f:
    f.write_function(func)