For every cell in a 2D mesh I am trying to obtain the transformation from the reference cell to that cell. For this I need the location of the vertices and the adjacency list between the cells and vertices. The following snippet of code extracts that information and saves it into a file. Running it with multiple MPI processes, it ignores ghost cells and stores the information for each processes separately.
save.py
from mpi4py import MPI
import dolfinx.mesh as mesh
import numpy as np
import basix
import sys
from pathlib import Path
if __name__ == "__main__":
assert (
len(sys.argv) == 2
), "Usage: python3 save.py <Number of cells in both x and y direction>"
n_cells: int = int(sys.argv[1])
width: float = 10.0
height: float = 8.0
cell_type: basix.CellType = basix.CellType.triangle
msh: mesh.Mesh = mesh.create_rectangle(
comm=MPI.COMM_WORLD,
points=[[0.0, 0.0], [width, height]],
n=[n_cells, n_cells],
dtype=np.float64,
)
tdim: int = msh.topology.dim
gdim: int = msh.geometry.dim
local_vertices: np.ndarray = np.asarray(msh.geometry.x.copy())[:, :gdim]
num_local_cells: int = msh.topology.index_map(tdim).size_local
msh.topology.create_connectivity(tdim, 0)
adj_list = msh.topology.connectivity(tdim, 0)
adj_list_arr = adj_list.array.reshape((-1, tdim + 1))[:num_local_cells]
vertex_max = np.max(adj_list_arr)
local_vertices = local_vertices[: vertex_max + 1]
np.save(Path(f"./vertices_{MPI.COMM_WORLD.rank}.npy"), local_vertices)
np.save(Path(f"./adjlist_{MPI.COMM_WORLD.rank}.npy"), adj_list_arr)
With the next snippet of code I reconstruct the mesh and plot it using matplotlib
view.py
from pathlib import Path
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.tri as tri
import re
if __name__ == "__main__":
par_dir = Path(__file__).parent
pattern = re.compile(r"_(\d+)\.npy")
sort_fnc = lambda f: int(pattern.search(f.name).group(1))
vertex_files: list[Path] = sorted(par_dir.glob("vertices_*.npy"), key=sort_fnc)
adjlist_files: list[Path] = sorted(par_dir.glob("adjlist_*.npy"), key=sort_fnc)
vertices: list[np.ndarray] = [np.load(v) for v in vertex_files]
adjlists: list[np.ndarray] = [np.load(a) for a in adjlist_files]
cmap = plt.get_cmap("winter")
colors = cmap(np.linspace(0, 1, len(vertex_files)))
fig, ax = plt.subplots()
for v, a, c in zip(vertices, adjlists, colors):
ax.scatter(v[:, 0], v[:, 1], color=c)
ax.triplot(tri.Triangulation(x=v[:, 0], y=v[:, 1], triangles=a), color=c)
plt.show()
Running mpiexec -np 2 python3 save.py 4 a rectangular mesh with 4 cells in each direction is created and subsequently running python3 view.py plots it, and everything looks as expected. Similarly, using a single MPI process everything works fine for arbitrary number of cells. However, running — for example — mpiexec -np 2 python3 save.py 5 and viewing again with python3 view.py the result is completely wrong near the “MPI boundary”. Here is an stdout of the latter, printing the vertices and adjacency list:
RANK: 0
VERTICES:
[[ 8. 0. ]
[10. 0. ]
[10. 1.6]
[ 8. 1.6]
[ 6. 0. ]
[10. 3.2]
[ 6. 1.6]
[ 8. 3.2]
[ 4. 0. ]
[10. 4.8]
[ 4. 1.6]
[ 6. 3.2]
[ 2. 0. ]
[ 2. 1.6]
[ 6. 4.8]
[ 0. 0. ]
[ 8. 4.8]
[10. 6.4]
[ 4. 3.2]
[ 0. 1.6]
[ 2. 3.2]]
ADJACENCY:
<AdjacencyList> with 30 nodes
0: [0 1 2 ]
1: [0 3 2 ]
2: [4 0 3 ]
3: [3 2 5 ]
4: [4 6 3 ]
5: [3 7 5 ]
6: [8 4 6 ]
7: [6 3 7 ]
8: [7 5 9 ]
9: [8 10 6 ]
10: [6 11 7 ]
11: [7 19 9 ]
12: [12 8 10 ]
13: [10 6 11 ]
14: [11 7 19 ]
15: [19 9 20 ]
16: [12 13 10 ]
17: [10 18 11 ]
18: [11 14 19 ]
19: [15 12 13 ]
20: [13 10 18 ]
21: [18 11 14 ]
22: [15 16 13 ]
23: [13 17 18 ]
24: [16 13 17 ]
25: [16 21 17 ]
26: [17 18 22 ]
27: [18 22 14 ]
28: [14 19 23 ]
29: [19 23 20 ]
RANK: 1
VERTICES:
[[ 0. 1.6]
[ 0. 3.2]
[ 2. 3.2]
[ 0. 6.4]
[ 0. 8. ]
[ 2. 8. ]
[ 0. 4.8]
[ 2. 4.8]
[ 2. 6.4]
[ 4. 4.8]
[ 4. 8. ]
[ 4. 6.4]
[ 4. 3.2]
[ 6. 6.4]
[ 6. 8. ]
[ 8. 6.4]
[ 8. 8. ]
[ 8. 4.8]
[10. 6.4]
[10. 8. ]
[ 6. 4.8]]
ADJACENCY:
<AdjacencyList> with 30 nodes
0: [0 1 2 ]
1: [3 4 5 ]
2: [1 6 7 ]
3: [6 3 8 ]
4: [1 2 7 ]
5: [3 8 5 ]
6: [6 7 8 ]
7: [2 7 9 ]
8: [8 5 10 ]
9: [7 8 11 ]
10: [2 12 9 ]
11: [7 9 11 ]
12: [8 11 10 ]
13: [12 9 20 ]
14: [9 11 13 ]
15: [11 10 14 ]
16: [9 20 13 ]
17: [11 13 14 ]
18: [20 13 15 ]
19: [13 14 16 ]
20: [20 17 15 ]
21: [13 15 16 ]
22: [17 15 18 ]
23: [15 16 19 ]
24: [15 18 19 ]
25: [0 21 2 ]
26: [21 2 12 ]
27: [12 22 20 ]
28: [22 20 17 ]
29: [17 23 18 ]
For example — from drawing it out by hand — the 11th cell on rank 0 has connectivity [7 19 9] with corresponding vertices (8, 3.2), (0, 1.6) and (10, 4.8) respectively. This triangle is not what I would
have expected. It should probably have been [7 16 9]. For reference, vertex 16 corresponds to (8, 4.8).
I installed the v0.9.0 version of dolfinx (tagged on GitHub) manually with Intel oneAPI compilers and MPI libraries.
Thanks in advance for helping to resolve this issue.