conspire/geometry/mesh/partition/
mod.rs

1#[cfg(test)]
2mod test;
3
4mod agglomerate;
5mod part;
6mod rcb;
7
8pub use self::rcb::Bisection;
9
10use crate::{geometry::mesh::Mesh, math::Tensor};
11use std::array::from_fn;
12
13#[derive(Clone, Debug, PartialEq)]
14pub struct PartitionQuality {
15    pub sizes: Vec<usize>,
16    pub imbalance: f64,
17    pub interface_nodes: usize,
18    pub disconnected_parts: usize,
19}
20
21#[derive(Clone, Debug, Default, PartialEq)]
22pub struct Partition {
23    elements_parts: Vec<usize>,
24    parts_elements: Vec<Vec<usize>>,
25    parts_nodes: Vec<Vec<usize>>,
26    nodes_parts: Vec<Vec<usize>>,
27}
28
29impl Partition {
30    pub fn new<const D: usize>(mesh: &Mesh<D>, elements_parts: Vec<usize>) -> Self {
31        assert_eq!(
32            elements_parts.len(),
33            mesh.number_of_elements(),
34            "assignment must have one entry per element"
35        );
36        let number_of_parts = elements_parts
37            .iter()
38            .map(|&part| part + 1)
39            .max()
40            .unwrap_or(0);
41        let mut parts_elements = vec![Vec::new(); number_of_parts];
42        elements_parts
43            .iter()
44            .enumerate()
45            .for_each(|(element, &part)| parts_elements[part].push(element));
46        let mut parts_nodes = vec![Vec::new(); number_of_parts];
47        mesh.element_nodes()
48            .zip(&elements_parts)
49            .for_each(|(nodes, &part)| parts_nodes[part].extend(nodes));
50        parts_nodes.iter_mut().for_each(|nodes| {
51            nodes.sort_unstable();
52            nodes.dedup();
53        });
54        let mut nodes_parts = vec![Vec::new(); mesh.number_of_nodes()];
55        parts_nodes
56            .iter()
57            .enumerate()
58            .for_each(|(part, nodes)| nodes.iter().for_each(|&node| nodes_parts[node].push(part)));
59        Self {
60            elements_parts,
61            parts_elements,
62            parts_nodes,
63            nodes_parts,
64        }
65    }
66    pub fn number_of_parts(&self) -> usize {
67        self.parts_elements.len()
68    }
69    pub fn elements_parts(&self) -> &[usize] {
70        &self.elements_parts
71    }
72    pub fn part_of(&self, element: usize) -> usize {
73        self.elements_parts[element]
74    }
75    pub fn part_elements(&self, part: usize) -> &[usize] {
76        &self.parts_elements[part]
77    }
78    pub fn part_nodes(&self, part: usize) -> &[usize] {
79        &self.parts_nodes[part]
80    }
81    pub fn parts_nodes(&self) -> &[Vec<usize>] {
82        &self.parts_nodes
83    }
84    pub fn node_parts(&self, node: usize) -> &[usize] {
85        &self.nodes_parts[node]
86    }
87    pub fn interface_nodes(&self) -> Vec<usize> {
88        self.nodes_parts
89            .iter()
90            .enumerate()
91            .filter(|(_, parts)| parts.len() > 1)
92            .map(|(node, _)| node)
93            .collect()
94    }
95    pub fn quality<const D: usize>(&self, mesh: &Mesh<D>) -> PartitionQuality {
96        let sizes = self.parts_elements.iter().map(Vec::len).collect::<Vec<_>>();
97        let mean = self.elements_parts.len() as f64 / sizes.len() as f64;
98        let mut roots = (0..self.elements_parts.len()).collect::<Vec<_>>();
99        mesh.node_element_connectivity()
100            .iter()
101            .for_each(|elements| {
102                let mut firsts: Vec<(usize, usize)> = Vec::new();
103                elements.iter().for_each(|&element| {
104                    match firsts
105                        .iter()
106                        .find(|&&(part, _)| part == self.elements_parts[element])
107                    {
108                        Some(&(_, first)) => {
109                            let (a, b) = (find(&mut roots, first), find(&mut roots, element));
110                            roots[b] = a;
111                        }
112                        None => firsts.push((self.elements_parts[element], element)),
113                    }
114                })
115            });
116        let mut components = vec![0usize; sizes.len()];
117        (0..self.elements_parts.len())
118            .filter(|&element| find(&mut roots, element) == element)
119            .for_each(|element| components[self.elements_parts[element]] += 1);
120        PartitionQuality {
121            imbalance: sizes.iter().copied().max().unwrap_or(0) as f64 / mean,
122            interface_nodes: self.interface_nodes().len(),
123            disconnected_parts: components.iter().filter(|&&count| count > 1).count(),
124            sizes,
125        }
126    }
127}
128
129impl<const D: usize> Mesh<D> {
130    pub fn partition_box(&self, divisions: [usize; D]) -> Partition {
131        assert!(
132            divisions.iter().all(|&n| n > 0),
133            "divisions must be positive"
134        );
135        let points = self.element_points();
136        let mut lower = [f64::INFINITY; D];
137        let mut upper = [f64::NEG_INFINITY; D];
138        points.iter().for_each(|point| {
139            (0..D).for_each(|axis| {
140                lower[axis] = lower[axis].min(point[axis]);
141                upper[axis] = upper[axis].max(point[axis]);
142            })
143        });
144        let assignment = points
145            .iter()
146            .map(|point| {
147                (0..D).rev().fold(0, |part, axis| {
148                    let extent = upper[axis] - lower[axis];
149                    let cell = if extent > 0.0 {
150                        (((point[axis] - lower[axis]) / extent * divisions[axis] as f64) as usize)
151                            .min(divisions[axis] - 1)
152                    } else {
153                        0
154                    };
155                    part * divisions[axis] + cell
156                })
157            })
158            .collect();
159        Partition::new(self, assignment)
160    }
161    fn element_points(&self) -> Vec<[f64; D]> {
162        self.centroids()
163            .iter()
164            .map(|centroid| from_fn(|axis| centroid[axis].value()))
165            .collect()
166    }
167    fn element_nodes(&self) -> impl Iterator<Item = Vec<usize>> {
168        self.iter().flat_map(|block| {
169            block
170                .iter()
171                .map(move |element| block.element_nodes(element))
172        })
173    }
174}
175
176fn find(roots: &mut [usize], mut element: usize) -> usize {
177    while roots[element] != element {
178        roots[element] = roots[roots[element]];
179        element = roots[element];
180    }
181    element
182}
conspire/geometry/mesh/partition/mod.rs

conspire/geometry/mesh/partition/
mod.rs
