Find cell tags from two overlapped meshes with different resolutions

dolfinx
1

I try to mark all the overlapped region from two meshes with different resolutions. Here I use collision related functions in dolfinx.geometry and meshtag, but the number of overlapped elements is much less than the expected. The following is the code and screenshot with DOLFINx v0.6.0:

from mpi4py import MPI
import ufl
import numpy as np
from dolfinx import fem, io, mesh, geometry

def mark_cells(msh, cell_index):
    num_cells = msh.topology.index_map(
        msh.topology.dim).size_local + msh.topology.index_map(
        msh.topology.dim).num_ghosts
    cells  = np.arange(0, num_cells)
    values = np.full(cells.shape, 0, dtype=np.int32)
    values[cell_index] = np.full(len(cell_index), 1, dtype=np.int32)
    cell_tag = mesh.meshtags(msh, msh.topology.dim, cells, values)
    return cell_tag


mesh_big = mesh.create_unit_square(MPI.COMM_WORLD, 64, 64)
mesh_big.geometry.x[:, :2] -= 0.51
mesh_big.geometry.x[:, :2] *= 4

mesh_small = mesh.create_unit_square(MPI.COMM_WORLD, 1, 1)
mesh_small.geometry.x[:, :2] -= 0.5

bb_tree         = geometry.BoundingBoxTree(mesh_big, mesh_big.topology.dim)
cell_candidates = geometry.compute_collisions(bb_tree, mesh_small.geometry.x)
colliding_cells = geometry.compute_colliding_cells(mesh_big, cell_candidates, mesh_small.geometry.x)

cell_indices = []

for i, point in enumerate(mesh_small.geometry.x):
    if len(colliding_cells.links(i))>0:
        cell_indices.append(colliding_cells.links(i)[0])

cell_tags = mark_cells(mesh_big, cell_indices)

with io.XDMFFile(MPI.COMM_WORLD, "cell_tags.xdmf", "w") as file:
    file.write_mesh(mesh_big)
    file.write_meshtags(cell_tags)

with io.XDMFFile(MPI.COMM_WORLD, "mesh_small.xdmf", "w") as file:
    file.write_mesh(mesh_small)

# Validation with the area of the overlapped region
dxs = ufl.Measure("dx", domain=mesh_big, subdomain_data=cell_tags)
f = 1*dxs(1)
area = MPI.COMM_WORLD.allreduce(
    fem.assemble_scalar(fem.form(f)), op=MPI.SUM)
print("Area:", area)

collision
collision1976×995 381 KB

2

This is because you are computing the intersection between the nodes (mesh geometry, in your case the vertices), and the mesh, not the cells itself.
The following code runs in dolfinx v0.7.x (Only minor changes to adapt it to 0.6.x required):

from mpi4py import MPI
import ufl
import numpy as np
from dolfinx import fem, io, mesh, geometry, cpp


def mark_cells(msh, cell_index):
    num_cells = msh.topology.index_map(
        msh.topology.dim).size_local + msh.topology.index_map(
        msh.topology.dim).num_ghosts
    cells = np.arange(0, num_cells, dtype=np.int32)
    values = np.full(cells.shape, 0, dtype=np.int32)
    values[cell_index] = np.full(len(cell_index), 1, dtype=np.int32)
    cell_tag = mesh.meshtags(msh, msh.topology.dim, cells, values)
    return cell_tag


mesh_big = mesh.create_unit_square(MPI.COMM_WORLD, 64, 64)
mesh_big.geometry.x[:, :2] -= 0.51
mesh_big.geometry.x[:, :2] *= 4
num_big_cells = mesh_big.topology.index_map(mesh_big.topology.dim).size_local + \
    mesh_big.topology.index_map(mesh_big.topology.dim).num_ghosts


mesh_small = mesh.create_unit_square(MPI.COMM_WORLD, 1, 1)
mesh_small.geometry.x[:, :2] -= 0.5
num_small_cells = mesh_small.topology.index_map(mesh_small.topology.dim).size_local + \
    mesh_small.topology.index_map(mesh_small.topology.dim).num_ghosts

bb_tree = geometry.bb_tree(
    mesh_big, mesh_big.topology.dim, np.arange(num_big_cells, dtype=np.int32))
bb_small = geometry.bb_tree(
    mesh_small, mesh_small.topology.dim, np.arange(num_small_cells, dtype=np.int32))
collisions = geometry.compute_collisions_trees(bb_tree, bb_small)


def extract_cell_geometry(input_mesh, cell: int):
    mesh_nodes = cpp.mesh.entities_to_geometry(
        input_mesh._cpp_object, input_mesh.topology.dim, np.array([cell], dtype=np.int32), False)[0]

    return input_mesh.geometry.x[mesh_nodes]


tol = 1e-13
big_cells = []
small_cells = []
for i, (big_cell, small_cell) in enumerate(collisions):
    geom_small = extract_cell_geometry(mesh_small, small_cell)
    geom_big = extract_cell_geometry(mesh_big, big_cell)
    distance = geometry.compute_distance_gjk(geom_big, geom_small)
    if np.linalg.norm(distance) <= tol:
        big_cells.append(big_cell)
        small_cells.append(small_cell)
cell_tags_big = mark_cells(mesh_big, np.asarray(big_cells, dtype=np.int32))
cell_tags_small = mark_cells(
    mesh_small, np.asarray(small_cells, dtype=np.int32))

with io.XDMFFile(MPI.COMM_WORLD, "cell_tags.xdmf", "w") as file:
    file.write_mesh(mesh_big)
    file.write_meshtags(cell_tags_big, mesh_big.geometry)

with io.XDMFFile(MPI.COMM_WORLD, "mesh_small.xdmf", "w") as file:
    file.write_mesh(mesh_small)
    file.write_meshtags(cell_tags_small, mesh_small.geometry)

# Validation with the area of the overlapped region
dxs = ufl.Measure("dx", domain=mesh_big, subdomain_data=cell_tags_big)
f = 1*dxs(1)
area = MPI.COMM_WORLD.allreduce(
    fem.assemble_scalar(fem.form(f)), op=MPI.SUM)
print("Area:", area)

Note that to get this to work in parallel, you need to do slightly more work than this, as the bounding boxes only work on the given process, no parallel communication internally.
It is not the most complicated thing to add, but as I’m not sure you need it I won’t add it here.

1 Like
3

I get that you used two bounding box trees of cells to find colliding cells with compute_collisions_trees (or compute_collisions in v0.6.0). As you have already found the colliding cells in two meshes respectively, why you further compute the gjk distance between two colliding cells? The result looks equivalent if I replace the following code

with

big_cells = [x for x, y in collisions]
small_cells = [y for x, y in collisions]

The cell_tag shows only a part if I run in parallel as you said. I am wondering how to adapt it to the parallel version.

4

A off-topic question: The calling of input_mesh._cpp_object returns a dolfinx.cpp.mesh.Mesh_float64 object while input_mesh gives dolfinx.mesh.Mesh object. The API of cpp.mesh.entities_to_geometry in v0.7.x gives one more calling method than that in v0.6.x. What is the goal of this design?

5

This is because your grids are axis aligned. The first computation is between bounding boxes (not triangular cells), thus will in most cases give more results than computing the collisions between the cells. We first use the collisions between the trees as it is much faster than computing than the collision between all cells.

All core structures of DOLFINx are written in C++.
These are often templated (for instance over mesh geometry precision). In Python, you do not have the templating notion, a mesh is a mesh irregardless of the floating precision. Additionally, the mesh has to interact with other packages written in pure Python (ufl).

6

Thanks, and how can get it to work in parallel when the bounding box tree only works on the given process?

7

This is quite more elaborate, and I’ve published a code under MIT license that does exactly this.

Eventually, I might refine it and make an extension to DOLFINx that performs these operations.

2 Likes
8

In a parallel run with mpirun -n 2 python .. (v0.7.2), the program stucks with the following message. Looks like the MPI alltoall issue?

--------------------------------------------------------------------------
WARNING: Linux kernel CMA support was requested via the
btl_vader_single_copy_mechanism MCA variable, but CMA support is
not available due to restrictive ptrace settings.

The vader shared memory BTL will fall back on another single-copy
mechanism if one is available. This may result in lower performance.

  Local host: *****
--------------------------------------------------------------------------
Traceback (most recent call last):
  File "*****/chenyongxin/test/fenicsx/parallel_celltags.py", line 111, in <module>
    small_to_big_comm.Neighbor_alltoall(send_geom_size, recv_geom_size)
  File "mpi4py/MPI/Comm.pyx", line 2154, in mpi4py.MPI.Topocomm.Neighbor_alltoall
mpi4py.MPI.Exception: MPI_ERR_TRUNCATE: message truncated
[*****:1221478] 1 more process has sent help message help-btl-vader.txt / cma-permission-denied
[*****:1221478] Set MCA parameter "orte_base_help_aggregate" to 0 to see all help / error messages

Additionally, the Linux kernel warning pops up in every parallel run even if those programs complete successfully.

9

I cannot reproduce this with the docker images for v0.7.2 or v0.7.3
docker run -ti --network=host -e DISPLAY=$DISPLAY -v /tmp/.X11-unix:/tmp/.X11-unix -v $(pwd):/root/shared -w /root/shared --rm ghcr.io/fenics/dolfinx/dolfinx:v0.7.3
docker run -ti --network=host -e DISPLAY=$DISPLAY -v /tmp/.X11-unix:/tmp/.X11-unix -v $(pwd):/root/shared -w /root/shared --rm ghcr.io/fenics/dolfinx/dolfinx:v0.7.2

How did you install DOLFINx?

10

I installed dolfinx with conda:

conda install -c conda-forge fenics-dolfinx=0.7.2
11

This is actual an issue with openmpi. I’ve pushed some updates that should fix this issue.
Thanks to @minrk for discussing these issues.
Minimal failing example with openmpi (works with mpich)

from mpi4py import MPI

import numpy as np

comm = MPI.COMM_WORLD
assert (comm.size == 2)
dtype = np.int64
if comm.rank == 0:
    outgoing_edges = np.array([], dtype=dtype)
elif comm.rank == 1:
    outgoing_edges = np.array([0], dtype=dtype)
sub_comm = comm.Create_dist_graph(
    list([comm.rank]), [len(np.unique(outgoing_edges))], outgoing_edges, reorder=False)
source, dest, _ = sub_comm.Get_dist_neighbors()

send_data = np.full(len(dest), comm.rank, dtype=dtype)
recv_data = np.full(len(source), -1, dtype=dtype)
print(comm.rank, source, dest,
      send_data, recv_data, flush=True)
sub_comm.Neighbor_alltoall(send_data, recv_data)
comm.Barrier()
print(comm.rank, source, dest,
      send_data, recv_data, flush=True)

Corrected code to work with openmpi

from mpi4py import MPI

import numpy as np

comm = MPI.COMM_WORLD
assert (comm.size == 2)
dtype = np.int64
if comm.rank == 0:
    outgoing_edges = np.array([], dtype=dtype)
elif comm.rank == 1:
    outgoing_edges = np.array([0], dtype=dtype)
sub_comm = comm.Create_dist_graph(
    list([comm.rank]), [len(np.unique(outgoing_edges))], outgoing_edges, reorder=False)
source, dest, _ = sub_comm.Get_dist_neighbors()

send_data = np.full(max(len(dest), 1), comm.rank, dtype=dtype)
recv_data = np.full(max(len(source), 1), -1, dtype=dtype)
print(comm.rank, source, dest,
      send_data, recv_data, flush=True)
sub_comm.Neighbor_alltoall(send_data, recv_data)
send_data = send_data[:len(dest)]
recv_data = recv_data[:len(source)]
comm.Barrier()
print(comm.rank, source, dest,
      send_data, recv_data, flush=True)