algebraeon_geometry/
simplicial_disjoint_union.rs

1use super::*;
2use crate::{
3    ambient_space::AffineSpace,
4    convex_hull::ConvexHull,
5    partial_simplicial_complex::LabelledPartialSimplicialComplex,
6    simplex::Simplex,
7    simplex_collection::{InteriorOrBoundarySimplexCollection, LabelledSimplexCollection},
8};
9use rayon::iter::{IntoParallelIterator, ParallelIterator};
10use std::collections::{HashMap, HashSet};
11
12/// A collection of disjoint simplices labelled by T with no additional structure
13#[derive(Clone)]
14pub struct LabelledSimplicialDisjointUnion<
15    'f,
16    FS: OrderedRingSignature + FieldSignature,
17    T: Eq + Clone + Send + Sync,
18> where
19    FS::Set: Hash,
20{
21    ambient_space: AffineSpace<'f, FS>,
22    simplexes: HashMap<Simplex<'f, FS>, T>,
23}
24
25pub type SimplicialDisjointUnion<'f, FS> = LabelledSimplicialDisjointUnion<'f, FS, ()>;
26
27impl<'f, FS: OrderedRingSignature + FieldSignature, T: Eq + Clone + Send + Sync>
28    LabelledSimplexCollection<'f, FS, T> for LabelledSimplicialDisjointUnion<'f, FS, T>
29where
30    FS::Set: Hash,
31{
32    type WithLabel<S: Eq + Clone + Send + Sync> = LabelledSimplicialDisjointUnion<'f, FS, S>;
33    type SubsetType = LabelledSimplicialDisjointUnion<'f, FS, T>;
34
35    fn try_new_labelled(
36        ambient_space: AffineSpace<'f, FS>,
37        simplexes: HashMap<Simplex<'f, FS>, T>,
38    ) -> Result<Self, &'static str> {
39        //todo: check simplexes are disjoint
40        Ok(Self {
41            ambient_space,
42            simplexes,
43        })
44    }
45
46    fn new_labelled_unchecked(
47        ambient_space: AffineSpace<'f, FS>,
48        simplexes: HashMap<Simplex<'f, FS>, T>,
49    ) -> Self {
50        let s = Self {
51            ambient_space,
52            simplexes,
53        };
54        #[cfg(debug_assertions)]
55        s.check();
56        s
57    }
58
59    fn ambient_space(&self) -> AffineSpace<'f, FS> {
60        self.ambient_space
61    }
62
63    fn labelled_simplexes(&self) -> HashMap<&Simplex<'f, FS>, &T> {
64        self.simplexes.iter().collect()
65    }
66
67    fn into_labelled_simplexes(self) -> HashMap<Simplex<'f, FS>, T> {
68        self.simplexes
69    }
70
71    fn into_partial_simplicial_complex(self) -> LabelledPartialSimplicialComplex<'f, FS, T> {
72        self.refine_into_partial_simplicial_complex()
73    }
74
75    fn to_partial_simplicial_complex(&self) -> LabelledPartialSimplicialComplex<'f, FS, T> {
76        self.clone().refine_into_partial_simplicial_complex()
77    }
78
79    fn into_simplicial_disjoint_union(self) -> LabelledSimplicialDisjointUnion<'f, FS, T> {
80        self
81    }
82
83    fn to_simplicial_disjoint_union(&self) -> LabelledSimplicialDisjointUnion<'f, FS, T> {
84        self.clone()
85    }
86}
87
88impl<'f, FS: OrderedRingSignature + FieldSignature, T: Eq + Clone + Send + Sync>
89    LabelledSimplicialDisjointUnion<'f, FS, T>
90where
91    FS::Set: Hash,
92{
93    #[allow(unused)]
94    pub(crate) fn check(&self) {
95        for spx_a in self.simplexes.keys() {
96            for spx_b in self.simplexes.keys() {
97                let bdry_a = spx_a
98                    .sub_simplices_not_null()
99                    .into_iter()
100                    .collect::<HashSet<_>>();
101                let bdry_b = spx_b
102                    .sub_simplices_not_null()
103                    .into_iter()
104                    .collect::<HashSet<_>>();
105                if !bdry_a.contains(spx_b) && !bdry_b.contains(spx_a) {
106                    let overlap = ConvexHull::intersect(
107                        &ConvexHull::from_simplex(spx_a.clone()),
108                        &ConvexHull::from_simplex(spx_b.clone()),
109                    );
110
111                    if !(overlap.affine_span_dimension() < spx_a.n()
112                        && overlap.affine_span_dimension() < spx_b.n())
113                    {
114                        println!("spx_a = {spx_a:?}");
115                        println!("spx_b = {spx_b:?}");
116                        panic!("simplicial complex simplex overlap");
117                    }
118                }
119            }
120        }
121    }
122
123    pub fn into_partial_simplicial_complex(self) -> LabelledPartialSimplicialComplex<'f, FS, T> {
124        self.refine_into_partial_simplicial_complex().simplify()
125    }
126
127    pub fn refine_into_partial_simplicial_complex(
128        mut self,
129    ) -> LabelledPartialSimplicialComplex<'f, FS, T> {
130        let ambient_space = self.ambient_space();
131
132        //maintain a list of pairs of simplexes which may intersect on their boundary
133        let mut pairs_todo: HashMap<Simplex<'f, FS>, HashSet<Simplex<'f, FS>>> = HashMap::new();
134        let simplexes = self.simplexes().into_iter().collect::<Vec<_>>();
135
136        for (spx_i, spx_j) in (0..simplexes.len())
137            .flat_map(|i| ((i + 1)..simplexes.len()).map(move |j| (i, j)))
138            .collect::<Vec<_>>()
139            .into_par_iter()
140            .filter_map(|(i, j)| {
141                let spx_i = simplexes[i];
142                let spx_j = simplexes[j];
143                if simplex_overlap::simplex_closure_overlap(spx_i, spx_j) {
144                    Some((spx_i, spx_j))
145                } else {
146                    None
147                }
148            })
149            .collect::<Vec<_>>()
150        {
151            pairs_todo
152                .entry(spx_i.clone())
153                .or_default()
154                .insert(spx_j.clone());
155            pairs_todo
156                .entry(spx_j.clone())
157                .or_default()
158                .insert(spx_i.clone());
159        }
160
161        while !pairs_todo.is_empty() {
162            let spx1 = pairs_todo.keys().next().unwrap().clone();
163            match pairs_todo.get(&spx1).unwrap().iter().next().cloned() {
164                None => {
165                    pairs_todo.remove(&spx1);
166                }
167                Some(spx2) => {
168                    //The pair (spx1, spx2) no longer needs to be checked
169                    pairs_todo.get_mut(&spx1).unwrap().remove(&spx2);
170                    pairs_todo.get_mut(&spx2).unwrap().remove(&spx1);
171
172                    debug_assert_ne!(spx1, spx2);
173                    debug_assert!(self.simplexes.contains_key(&spx1));
174                    debug_assert!(self.simplexes.contains_key(&spx2));
175
176                    let overlap = ConvexHull::intersect(
177                        &ConvexHull::from_simplex(spx1.clone()),
178                        &ConvexHull::from_simplex(spx2.clone()),
179                    );
180
181                    if !overlap.is_empty()
182                        && match ambient_space
183                            .simplex(overlap.defining_points().into_iter().collect())
184                        {
185                            Ok(overlap_spx) => {
186                                let spx1_points = spx1.points().iter().collect::<HashSet<_>>();
187                                let spx2_points = spx2.points().iter().collect::<HashSet<_>>();
188                                !overlap_spx
189                                    .points()
190                                    .iter()
191                                    .all(|pt| spx1_points.contains(pt) && spx2_points.contains(pt))
192                            }
193                            Err(_) => true,
194                        }
195                    {
196                        //there is a bad overlap between spx1 and spx2
197                        let mut spx1_replacement = overlap.clone();
198                        for pt in spx1.points() {
199                            spx1_replacement.extend_by_point(pt.clone());
200                        }
201                        let mut spx2_replacement = overlap.clone();
202                        for pt in spx2.points() {
203                            spx2_replacement.extend_by_point(pt.clone());
204                        }
205
206                        //remove any pairs containing spx1 or spx2
207                        let mut spx1_paired = vec![];
208                        for spx in pairs_todo.get(&spx1).unwrap_or(&HashSet::new()).clone() {
209                            debug_assert_ne!(spx, spx1);
210                            debug_assert_ne!(spx, spx2);
211                            debug_assert!(self.simplexes.contains_key(&spx));
212                            pairs_todo
213                                .get_mut(&spx)
214                                .unwrap_or(&mut HashSet::new())
215                                .remove(&spx1);
216                            spx1_paired.push(spx);
217                        }
218                        pairs_todo.remove(&spx1);
219
220                        let mut spx2_paired = vec![];
221                        for spx in pairs_todo.get(&spx2).unwrap_or(&HashSet::new()).clone() {
222                            debug_assert_ne!(spx, spx1);
223                            debug_assert_ne!(spx, spx2);
224                            debug_assert!(self.simplexes.contains_key(&spx));
225                            pairs_todo
226                                .get_mut(&spx)
227                                .unwrap_or(&mut HashSet::new())
228                                .remove(&spx2);
229                            spx2_paired.push(spx);
230                        }
231                        pairs_todo.remove(&spx2);
232
233                        let spx1_label = self.simplexes.get(&spx1).unwrap().clone();
234                        let spx2_label = self.simplexes.get(&spx2).unwrap().clone();
235
236                        //remove spx1 and spx2
237                        self.simplexes.remove(&spx1);
238                        self.simplexes.remove(&spx2);
239
240                        //add the refinements of spx1 and spx2 and update the pairs todo
241                        for spx1_repl in spx1_replacement
242                            .to_simplicial_complex()
243                            .interior()
244                            .into_simplexes()
245                        {
246                            for spx in &spx1_paired {
247                                pairs_todo
248                                    .entry(spx1_repl.clone())
249                                    .or_default()
250                                    .insert(spx.clone());
251                                pairs_todo
252                                    .entry(spx.clone())
253                                    .or_default()
254                                    .insert(spx1_repl.clone());
255                            }
256                            self.simplexes.insert(spx1_repl, spx1_label.clone());
257                        }
258
259                        for spx2_repl in spx2_replacement
260                            .to_simplicial_complex()
261                            .interior()
262                            .into_simplexes()
263                        {
264                            for spx in &spx2_paired {
265                                pairs_todo
266                                    .entry(spx2_repl.clone())
267                                    .or_default()
268                                    .insert(spx.clone());
269                                pairs_todo
270                                    .entry(spx.clone())
271                                    .or_default()
272                                    .insert(spx2_repl.clone());
273                            }
274                            self.simplexes.insert(spx2_repl, spx2_label.clone());
275                        }
276                    }
277                }
278            }
279        }
280
281        LabelledPartialSimplicialComplex::new_labelled_unchecked(ambient_space, self.simplexes)
282    }
283}
algebraeon_geometry/simplicial_disjoint_union.rs

algebraeon_geometry/
simplicial_disjoint_union.rs
