parallel_mesh/
parallel_mesh.rs

1use itertools::izip;
2use mpi::{
3    collective::SystemOperation, environment::Universe, topology::Communicator,
4    traits::CommunicatorCollectives,
5};
6use ndelement::types::ReferenceCellType;
7use ndmesh::{
8    mesh::local_mesh::SingleElementMeshBuilder,
9    traits::{Builder, Entity, Mesh, ParallelBuilder, ParallelMesh},
10    types::{GraphPartitioner, Ownership},
11};
12
13/// Creating a distributed parallel mesh
14fn main() {
15    // The SingleElementMeshBuilder is used to create the mesh
16    let mut b = SingleElementMeshBuilder::<f64>::new(2, (ReferenceCellType::Quadrilateral, 1));
17
18    let universe: Universe = mpi::initialize().unwrap();
19    let comm = universe.world();
20    let rank = comm.rank();
21
22    // Add points and cells to the builder on process 0
23    let n = 10;
24    let mesh = if rank == 0 {
25        let mut i = 0;
26        for y in 0..n {
27            for x in 0..n {
28                b.add_point(i, &[x as f64 / (n - 1) as f64, y as f64 / (n - 1) as f64]);
29                i += 1;
30            }
31        }
32
33        let mut i = 0;
34        for y in 0..n - 1 {
35            for x in 0..n - 1 {
36                let sw = n * y + x;
37                b.add_cell(i, &[sw, sw + 1, sw + n, sw + n + 1]);
38                i += 1;
39            }
40        }
41
42        // Distribute the mesh
43        // In this example, we use Scotch to partition the mesh into pieces to be handles by each process
44        b.create_parallel_mesh_root(&comm, GraphPartitioner::Scotch)
45    } else {
46        // receice the mesh
47        b.create_parallel_mesh(&comm, 0)
48    };
49
50    // Check that owned cells are sorted ahead of ghost cells
51    let cell_count_owned = mesh
52        .local_mesh()
53        .entity_iter(ReferenceCellType::Quadrilateral)
54        .filter(|entity| entity.is_owned())
55        .count();
56
57    // Now check that the first `cell_count_owned` entities are actually owned.
58    for cell in mesh
59        .local_mesh()
60        .entity_iter(ReferenceCellType::Quadrilateral)
61        .take(cell_count_owned)
62    {
63        assert!(cell.is_owned())
64    }
65
66    // Now make sure that the indices of the global cells are in consecutive order
67    for (cell_global_count, cell) in izip!(
68        mesh.cell_layout().local_range().0..,
69        mesh.local_mesh()
70            .entity_iter(ReferenceCellType::Quadrilateral)
71            .take(cell_count_owned)
72    ) {
73        assert_eq!(cell.global_index(), cell_global_count);
74    }
75
76    // Get the global indices
77    let global_vertices = mesh
78        .local_mesh()
79        .entity_iter(ReferenceCellType::Point)
80        .filter(|e| matches!(e.ownership(), Ownership::Owned))
81        .map(|e| e.global_index())
82        .collect::<Vec<_>>();
83
84    let nvertices = global_vertices.len();
85
86    let global_cells = mesh
87        .local_mesh()
88        .entity_iter(ReferenceCellType::Quadrilateral)
89        .filter(|e| matches!(e.ownership(), Ownership::Owned))
90        .map(|e| e.global_index())
91        .collect::<Vec<_>>();
92
93    let ncells = global_cells.len();
94
95    let mut total_cells: usize = 0;
96    let mut total_vertices: usize = 0;
97
98    comm.all_reduce_into(&ncells, &mut total_cells, SystemOperation::sum());
99    comm.all_reduce_into(&nvertices, &mut total_vertices, SystemOperation::sum());
100
101    // Check that the total number of cells and vertices are correct
102    assert_eq!(total_cells, (n - 1) * (n - 1));
103    assert_eq!(total_vertices, n * n);
104}
parallel_mesh/parallel_mesh.rs

parallel_mesh/
parallel_mesh.rs
