conspire/geometry/mesh/tessellation/cut/tables/
mod.rs

1#[cfg(test)]
2mod test;
3
4use super::{
5    CROSSING_TOLERANCE, Class, DIRECTIONS, Sign, Tables, Vertex,
6    face::face_cut,
7    geometry::dedupe,
8    topology::{element_edges, element_faces, face_owners, oriented_element_faces},
9};
10use crate::{
11    geometry::{
12        Coordinate, DirectionsRef,
13        mesh::{Mesh, tessellation::D, tessellation::Tessellation},
14    },
15    math::{Quantity, Scalar, Tensor},
16    units::Length,
17};
18use std::{array::from_fn, collections::HashMap, collections::HashSet};
19
20pub(super) struct GenericTables {
21    signs: HashMap<usize, Sign>,
22    crossings: HashMap<[usize; 2], Vec<Coordinate<D>>>,
23    faces: HashMap<Vec<usize>, Vec<usize>>,
24    segments: HashMap<Vec<usize>, Vec<[Vertex; 2]>>,
25}
26
27impl GenericTables {
28    pub(super) fn signs(&self) -> &HashMap<usize, Sign> {
29        &self.signs
30    }
31    pub(super) fn crossings(&self) -> &HashMap<[usize; 2], Vec<Coordinate<D>>> {
32        &self.crossings
33    }
34    pub(super) fn faces(&self) -> &HashMap<Vec<usize>, Vec<usize>> {
35        &self.faces
36    }
37    pub(super) fn segments(&self) -> &HashMap<Vec<usize>, Vec<[Vertex; 2]>> {
38        &self.segments
39    }
40}
41
42impl Tessellation {
43    #[allow(clippy::type_complexity)]
44    fn edges_signs_crossings(
45        &self,
46        mesh: &Mesh<D>,
47        classes: &[Class],
48        snapped: &HashSet<usize>,
49    ) -> Result<
50        (
51            HashSet<[usize; 2]>,
52            HashMap<usize, Sign>,
53            HashMap<[usize; 2], Vec<Coordinate<D>>>,
54        ),
55        &'static str,
56    > {
57        let surface = self.mesh();
58        let surface_coordinates = surface.coordinates();
59        let elements: Vec<&[usize]> = surface.connectivities().iter().flatten().collect();
60        let normals: DirectionsRef<'_, D> = self.normals().iter().flatten().collect();
61        let directions = DIRECTIONS.map(|direction| direction.normalized());
62        let bvh = self.bvh();
63        let coordinates = mesh.coordinates();
64        let mut edges = HashSet::new();
65        let mut offset = 0;
66        mesh.iter().for_each(|block| {
67            block.iter().enumerate().for_each(|(local, element)| {
68                if classes[offset + local] == Class::Cut {
69                    element_edges(&element_faces(block, element))
70                        .into_iter()
71                        .for_each(|key| {
72                            edges.insert(key);
73                        });
74                }
75            });
76            offset += block.number_of_elements();
77        });
78        let mut signs = HashMap::new();
79        edges.iter().flatten().for_each(|&node| {
80            signs.entry(node).or_insert_with(|| {
81                if snapped.contains(&node) {
82                    Sign::On
83                } else if self.encloses(
84                    &coordinates[node],
85                    surface_coordinates,
86                    &elements,
87                    &normals,
88                    &directions,
89                ) {
90                    Sign::Inside
91                } else {
92                    Sign::Outside
93                }
94            });
95        });
96        let mut crossings = HashMap::<[usize; 2], Vec<Coordinate<D>>>::new();
97        edges.iter().try_for_each(|&[a, b]| {
98            let span = &coordinates[b] - &coordinates[a];
99            let length = span.norm();
100            let margin = CROSSING_TOLERANCE.max(length * super::GRAZING_TOLERANCE);
101            match (signs[&a], signs[&b]) {
102                (Sign::On, Sign::On) => Ok(()),
103                (Sign::On, _) | (_, Sign::On) => {
104                    let (from, along) = if signs[&a] == Sign::On {
105                        (b, &coordinates[a] - &coordinates[b])
106                    } else {
107                        (a, span)
108                    };
109                    let hits = bvh.intersect_all(
110                        &(coordinates[from].clone(), along.clone()).into(),
111                        surface_coordinates,
112                        &elements,
113                    );
114                    let distances: Vec<Quantity<Length>> = dedupe(hits, margin)
115                        .into_iter()
116                        .filter(|&distance| distance < length - margin)
117                        .collect();
118                    if !distances.is_empty() {
119                        let points: Vec<Coordinate<D>> = distances
120                            .iter()
121                            .map(|&distance| &coordinates[from] + &(&along * (distance / length)))
122                            .collect();
123                        let ordered = if from == a {
124                            points
125                        } else {
126                            points.into_iter().rev().collect()
127                        };
128                        crossings.insert([a, b], ordered);
129                    }
130                    Ok(())
131                }
132                (inside, outside) if inside != outside => {
133                    let hits = bvh.intersect_all(
134                        &(coordinates[a].clone(), span.clone()).into(),
135                        surface_coordinates,
136                        &elements,
137                    );
138                    let distances: Vec<Quantity<Length>> = dedupe(hits, margin)
139                        .into_iter()
140                        .filter(|&distance| distance <= length + margin)
141                        .collect();
142                    if distances.is_empty() {
143                        return Err("crossing missing on a sign-change edge");
144                    }
145                    if distances.len().is_multiple_of(2) {
146                        return Err(
147                            "edge crosses the tessellation an inconsistent number of times",
148                        );
149                    }
150                    crossings.insert(
151                        [a, b],
152                        distances
153                            .iter()
154                            .map(|&distance| &coordinates[a] + &(&span * (distance / length)))
155                            .collect(),
156                    );
157                    Ok(())
158                }
159                _ => {
160                    let hits = bvh.intersect_all(
161                        &(coordinates[a].clone(), span.clone()).into(),
162                        surface_coordinates,
163                        &elements,
164                    );
165                    let distances: Vec<Quantity<Length>> = dedupe(hits, margin)
166                        .into_iter()
167                        .filter(|&distance| distance <= length + margin)
168                        .collect();
169                    if distances.len() % 2 == 1 {
170                        return Err(
171                            "edge crosses the tessellation an inconsistent number of times",
172                        );
173                    }
174                    if !distances.is_empty() {
175                        crossings.insert(
176                            [a, b],
177                            distances
178                                .iter()
179                                .map(|&distance| &coordinates[a] + &(&span * (distance / length)))
180                                .collect(),
181                        );
182                    }
183                    Ok(())
184                }
185            }
186        })?;
187        Ok((edges, signs, crossings))
188    }
189    pub fn tables(
190        &self,
191        mesh: &Mesh<D>,
192        classes: &[Class],
193        snapped: &HashSet<usize>,
194    ) -> Result<Tables, &'static str> {
195        let (_, signs, crossings) = self.edges_signs_crossings(mesh, classes, snapped)?;
196        let coordinates = mesh.coordinates();
197        let surface_coordinates = self.mesh().coordinates();
198        let elements: Vec<&[usize]> = self.mesh().connectivities().iter().flatten().collect();
199        let normals: DirectionsRef<'_, D> = self.normals().iter().flatten().collect();
200        let directions = DIRECTIONS.map(|direction| direction.normalized());
201        let contains = |point: &Coordinate<D>| {
202            self.encloses(point, surface_coordinates, &elements, &normals, &directions)
203        };
204        let mut face_loops = HashMap::new();
205        let mut offset = 0;
206        mesh.iter().for_each(|block| {
207            let local_faces = block.local_faces();
208            block.iter().enumerate().for_each(|(local, element)| {
209                if classes[offset + local] == Class::Cut {
210                    local_faces.iter().for_each(|face| {
211                        let corners = from_fn::<_, 4, _>(|i| element[face[i]]);
212                        let mut key = corners;
213                        key.sort_unstable();
214                        face_loops.entry(key).or_insert(corners);
215                    })
216                }
217            });
218            offset += block.number_of_elements();
219        });
220        let mut segments = HashMap::new();
221        face_loops.iter().try_for_each(|(key, corners)| {
222            let cut = face_cut(corners, &signs, &crossings)?;
223            let count = cut.endpoints.len();
224            if count == 0 {
225                return Ok(());
226            }
227            if count % 2 == 1 {
228                return Err("refinement required at a face");
229            }
230            let target = if count == 2 {
231                cut.sides[0]
232            } else {
233                let center = corners
234                    .iter()
235                    .map(|&node| coordinates[node].clone())
236                    .sum::<Coordinate<D>>()
237                    / 4.0;
238                if contains(&center) {
239                    Sign::Outside
240                } else {
241                    Sign::Inside
242                }
243            };
244            let pairs: Vec<[Vertex; 2]> = (0..count)
245                .filter(|&arc| cut.sides[arc] == target)
246                .map(|arc| [cut.endpoints[arc], cut.endpoints[(arc + 1) % count]])
247                .collect();
248            if pairs.len() != count / 2 {
249                return Err("inconsistent crossings around a face");
250            }
251            segments.insert(*key, pairs);
252            Ok(())
253        })?;
254        Ok(Tables {
255            signs,
256            crossings,
257            faces: face_loops,
258            segments,
259        })
260    }
261    pub(super) fn tables_generic(
262        &self,
263        mesh: &Mesh<D>,
264        classes: &[Class],
265        snapped: &HashSet<usize>,
266    ) -> Result<GenericTables, &'static str> {
267        let (_, signs, crossings) = self.edges_signs_crossings(mesh, classes, snapped)?;
268        let coordinates = mesh.coordinates();
269        let surface_coordinates = self.mesh().coordinates();
270        let elements: Vec<&[usize]> = self.mesh().connectivities().iter().flatten().collect();
271        let normals: DirectionsRef<'_, D> = self.normals().iter().flatten().collect();
272        let directions = DIRECTIONS.map(|direction| direction.normalized());
273        let contains = |point: &Coordinate<D>| {
274            self.encloses(point, surface_coordinates, &elements, &normals, &directions)
275        };
276        let mut face_loops = HashMap::<Vec<usize>, Vec<usize>>::new();
277        let mut offset = 0;
278        mesh.iter().for_each(|block| {
279            let owners = face_owners(block);
280            block.iter().enumerate().for_each(|(local, element)| {
281                if classes[offset + local] == Class::Cut {
282                    oriented_element_faces(block, element, local, owners.as_deref())
283                        .into_iter()
284                        .for_each(|corners| {
285                            let mut key = corners.clone();
286                            key.sort_unstable();
287                            face_loops.entry(key).or_insert(corners);
288                        })
289                }
290            });
291            offset += block.number_of_elements();
292        });
293        let mut segments = HashMap::new();
294        face_loops.iter().try_for_each(|(key, corners)| {
295            let cut = face_cut(corners, &signs, &crossings)?;
296            let count = cut.endpoints.len();
297            if count == 0 {
298                return Ok(());
299            }
300            if count % 2 == 1 {
301                return Err("refinement required at a face");
302            }
303            let n = corners.len();
304            let target = if count == 2 {
305                cut.sides[0]
306            } else {
307                let center = corners
308                    .iter()
309                    .map(|&node| coordinates[node].clone())
310                    .sum::<Coordinate<D>>()
311                    / n as Scalar;
312                if contains(&center) {
313                    Sign::Outside
314                } else {
315                    Sign::Inside
316                }
317            };
318            let pairs: Vec<[Vertex; 2]> = (0..count)
319                .filter(|&arc| cut.sides[arc] == target)
320                .map(|arc| [cut.endpoints[arc], cut.endpoints[(arc + 1) % count]])
321                .collect();
322            if pairs.len() != count / 2 {
323                return Err("inconsistent crossings around a face");
324            }
325            segments.insert(key.clone(), pairs);
326            Ok(())
327        })?;
328        Ok(GenericTables {
329            signs,
330            crossings,
331            faces: face_loops,
332            segments,
333        })
334    }
335}
conspire/geometry/mesh/tessellation/cut/tables/mod.rs

conspire/geometry/mesh/tessellation/cut/tables/
mod.rs
