dexterior_core/mesh/
views.rs

1use super::SimplicialMesh;
2
3use fixedbitset as fb;
4use nalgebra as na;
5use nalgebra_sparse as nas;
6
7/// A view into a single simplex's data.
8///
9/// This type can also be used as an index into a primal [`Cochain`][crate::Cochain]
10/// of the corresponding dimension, providing a shorter syntax
11/// and a type check to ensure the dimensions match:
12/// ```
13/// # use dexterior_core::{mesh::tiny_mesh_3d, Cochain, Primal, Dual};
14/// # let mesh_3d = tiny_mesh_3d();
15/// let c_primal: Cochain<1, Primal> = mesh_3d.new_zero_cochain();
16/// let c_dual: Cochain<2, Dual> = mesh_3d.new_zero_cochain();
17/// for edge in mesh_3d.simplices::<1>() {
18///     let primal_val = c_primal[edge];
19///     let dual_val = c_dual[edge.dual()];
20///     // ..is a typechecked equivalent to
21///     let primal_val = c_primal.values[edge.index()];
22///     let dual_val = c_dual.values[edge.index()];
23/// }
24/// ```
25/// Something like this wouldn't compile, for instance:
26/// ```compile_fail
27/// # use dexterior_core::{mesh::tiny_mesh_3d, Cochain, Dual};
28/// # let mesh_3d = tiny_mesh_3d();
29/// let c: Cochain<1, Dual> = mesh_3d.new_zero_cochain();
30/// for edge in mesh_3d.simplices::<1>() {
31///     let val = c[edge];
32///     // ..but the index method would still work (almost certainly incorrectly):
33///     let val = c.values[edge.index()];
34/// }
35/// ```
36#[derive(Clone, Copy, Debug)]
37pub struct SimplexView<'a, SimplexDim, const MESH_DIM: usize> {
38    pub(super) mesh: &'a SimplicialMesh<MESH_DIM>,
39    pub(super) index: usize,
40    pub(super) indices: &'a [usize],
41    pub(super) _marker: std::marker::PhantomData<SimplexDim>,
42}
43
44impl<'a, SimplexDim, const MESH_DIM: usize> PartialEq for SimplexView<'a, SimplexDim, MESH_DIM> {
45    fn eq(&self, other: &Self) -> bool {
46        self.index == other.index
47    }
48}
49impl<'a, SimplexDim, const MESH_DIM: usize> Eq for SimplexView<'a, SimplexDim, MESH_DIM> {}
50
51impl<'a, SimplexDim, const MESH_DIM: usize> SimplexView<'a, SimplexDim, MESH_DIM>
52where
53    SimplexDim: na::DimName,
54    na::Const<MESH_DIM>: na::DimNameSub<SimplexDim>,
55{
56    /// Iterate over the vertices of this simplex.
57    #[inline]
58    pub fn vertices(&self) -> impl '_ + Iterator<Item = na::SVector<f64, MESH_DIM>> {
59        self.indices.iter().map(|i| self.mesh.vertices[*i])
60    }
61
62    /// Iterate over the vertex indices of this simplex.
63    #[inline]
64    pub fn vertex_indices(&self) -> impl '_ + Iterator<Item = usize> {
65        self.indices.iter().cloned()
66    }
67
68    /// Get the index of this simplex in the ordering of `DIM`-simplices.
69    #[inline]
70    pub fn index(&self) -> usize {
71        self.index
72    }
73
74    /// Get the unsigned volume of this simplex.
75    #[inline]
76    pub fn volume(&self) -> f64 {
77        self.mesh.simplices[SimplexDim::USIZE].volumes[self.index]
78    }
79
80    /// Get the unsigned volume of the dual cell corresponding to this simplex.
81    #[inline]
82    pub fn dual_volume(&self) -> f64 {
83        self.mesh.simplices[SimplexDim::USIZE].dual_volumes[self.index]
84    }
85
86    /// Get the circumcenter of this simplex.
87    #[inline]
88    pub fn circumcenter(&self) -> na::SVector<f64, MESH_DIM> {
89        self.mesh.simplices[SimplexDim::USIZE].circumcenters[self.index]
90    }
91
92    /// Get the barycenter of this simplex.
93    #[inline]
94    pub fn barycenter(&self) -> na::SVector<f64, MESH_DIM> {
95        self.mesh.simplices[SimplexDim::USIZE].barycenters[self.index]
96    }
97
98    /// Get the dual cell corresponding to this simplex.
99    #[inline]
100    pub fn dual(
101        self,
102    ) -> DualCellView<'a, na::DimNameDiff<na::Const<MESH_DIM>, SimplexDim>, MESH_DIM> {
103        DualCellView {
104            mesh: self.mesh,
105            index: self.index,
106            _marker: std::marker::PhantomData,
107        }
108    }
109}
110
111impl<'a, SimplexDim, const MESH_DIM: usize> SimplexView<'a, SimplexDim, MESH_DIM>
112where
113    SimplexDim: na::DimNameAdd<na::U1>,
114    na::Const<MESH_DIM>: na::DimNameSub<na::DimNameSum<SimplexDim, na::U1>>,
115{
116    /// Iterate over the `k+1`-dimensional simplices on whose boundary this simplex lies.
117    pub fn coboundary(self) -> BoundaryIter<'a, na::DimNameSum<SimplexDim, na::U1>, MESH_DIM> {
118        BoundaryIter {
119            mesh: self.mesh,
120            index: 0,
121            map_row: self.mesh.simplices[SimplexDim::USIZE]
122                .coboundary_map
123                .row(self.index),
124            _marker: std::marker::PhantomData,
125        }
126    }
127}
128
129impl<'a, SimplexDim, const MESH_DIM: usize> SimplexView<'a, SimplexDim, MESH_DIM>
130where
131    SimplexDim: na::DimNameSub<na::U1>,
132{
133    /// Iterate over the `k-1`-dimensional simplices on the boundary of this simplex.
134    pub fn boundary(self) -> BoundaryIter<'a, na::DimNameDiff<SimplexDim, na::U1>, MESH_DIM> {
135        BoundaryIter {
136            mesh: self.mesh,
137            index: 0,
138            map_row: self.mesh.simplices[SimplexDim::USIZE]
139                .boundary_map
140                .row(self.index),
141            _marker: std::marker::PhantomData,
142        }
143    }
144}
145
146/// A view into a dual cell's data.
147///
148/// This type can also be used to index into a dual [`Cochain`][crate::Cochain]
149/// of the corresponding dimension in a type-checked fashion.
150/// See [`SimplexView`] for an example.
151#[derive(Clone, Copy, Debug)]
152pub struct DualCellView<'a, CellDim, const MESH_DIM: usize> {
153    pub(super) mesh: &'a SimplicialMesh<MESH_DIM>,
154    pub(super) index: usize,
155    pub(super) _marker: std::marker::PhantomData<CellDim>,
156}
157
158impl<'a, CellDim, const MESH_DIM: usize> PartialEq for DualCellView<'a, CellDim, MESH_DIM> {
159    fn eq(&self, other: &Self) -> bool {
160        self.index == other.index
161    }
162}
163impl<'a, CellDim, const MESH_DIM: usize> Eq for DualCellView<'a, CellDim, MESH_DIM> {}
164
165impl<'a, CellDim, const MESH_DIM: usize> DualCellView<'a, CellDim, MESH_DIM>
166where
167    CellDim: na::DimName,
168    na::Const<MESH_DIM>: na::DimNameSub<CellDim>,
169{
170    /// Get the index of this cell in the ordering of dual `DIM`-cells.
171    #[inline]
172    pub fn index(&self) -> usize {
173        self.index
174    }
175
176    /// Get the unsigned volume of this cell.
177    #[inline]
178    pub fn volume(&self) -> f64 {
179        self.mesh.simplices[MESH_DIM - CellDim::USIZE].dual_volumes[self.index]
180    }
181
182    /// Get the primal simplex corresponding to this cell.
183    #[inline]
184    pub fn dual(&self) -> SimplexView<'_, na::DimNameDiff<na::Const<MESH_DIM>, CellDim>, MESH_DIM> {
185        self.mesh.get_simplex_by_index_impl(self.index)
186    }
187}
188
189/// An iterator over the boundary or coboundary of a simplex,
190/// obtained with [`SimplexView::coboundary`] or [`SimplexView::boundary`].
191///
192/// The iterator returns a pair `(orientation, simplex)`
193/// where orientation is the relative orientation
194/// between the boundary and coboundary simplices in question.
195pub struct BoundaryIter<'a, SimplexDim, const MESH_DIM: usize> {
196    mesh: &'a SimplicialMesh<MESH_DIM>,
197    index: usize,
198    map_row: nas::csr::CsrRow<'a, i8>,
199    _marker: std::marker::PhantomData<SimplexDim>,
200}
201
202impl<'a, SimplexDim, const MESH_DIM: usize> Iterator for BoundaryIter<'a, SimplexDim, MESH_DIM>
203where
204    SimplexDim: na::DimName,
205{
206    type Item = (i8, SimplexView<'a, SimplexDim, MESH_DIM>);
207
208    fn next(&mut self) -> Option<Self::Item> {
209        let cols = self.map_row.col_indices();
210        if self.index >= cols.len() {
211            return None;
212        }
213        let next_idx = cols[self.index];
214        let next_ori = self.map_row.values()[self.index];
215
216        self.index += 1;
217
218        Some((next_ori, self.mesh.get_simplex_by_index_impl(next_idx)))
219    }
220}
221
222/// Iterator over a set of `DIM`-simplices in a mesh.
223pub struct SimplexIter<'a, const DIM: usize, const MESH_DIM: usize> {
224    pub(super) mesh: &'a SimplicialMesh<MESH_DIM>,
225    pub(super) idx_iter: IndexIter<'a>,
226}
227
228impl<'a, const DIM: usize, const MESH_DIM: usize> Iterator for SimplexIter<'a, DIM, MESH_DIM>
229where
230    na::Const<MESH_DIM>: na::DimNameSub<na::Const<DIM>>,
231{
232    type Item = SimplexView<'a, na::Const<DIM>, MESH_DIM>;
233
234    fn next(&mut self) -> Option<Self::Item> {
235        let next_idx = self.idx_iter.next()?;
236        Some(self.mesh.get_simplex_by_index::<DIM>(next_idx))
237    }
238}
239
240/// Iterator over a set of `DIM`-dimensional dual cells in a mesh.
241pub struct DualCellIter<'a, const DIM: usize, const MESH_DIM: usize> {
242    pub(super) mesh: &'a SimplicialMesh<MESH_DIM>,
243    pub(super) idx_iter: IndexIter<'a>,
244}
245
246impl<'a, const DIM: usize, const MESH_DIM: usize> Iterator for DualCellIter<'a, DIM, MESH_DIM>
247where
248    na::Const<MESH_DIM>: na::DimNameSub<na::Const<DIM>>,
249{
250    type Item = DualCellView<'a, na::Const<DIM>, MESH_DIM>;
251
252    fn next(&mut self) -> Option<Self::Item> {
253        let next_idx = self.idx_iter.next()?;
254        Some(DualCellView {
255            mesh: self.mesh,
256            index: next_idx,
257            _marker: std::marker::PhantomData,
258        })
259    }
260}
261
262/// A set of indices to iterate over,
263/// defined either as a contiguous range
264/// or an arbitrary set represented by a bitset.
265pub(super) enum IndexIter<'a> {
266    All(std::ops::Range<usize>),
267    Subset(fb::Ones<'a>),
268}
269
270impl<'a> Iterator for IndexIter<'a> {
271    type Item = usize;
272
273    fn next(&mut self) -> Option<Self::Item> {
274        match self {
275            Self::All(range) => range.next(),
276            Self::Subset(indices) => indices.next(),
277        }
278    }
279}
280
281#[cfg(test)]
282mod tests {
283    use super::super::*;
284
285    #[test]
286    fn iterators() {
287        let mesh = tiny_mesh_2d();
288        // for a quick sanity test, check that every edge
289        // is on the coboundary of all its boundary simplices and vice versa
290        for edge in mesh.simplices::<1>() {
291            for (b_ori, b) in edge.boundary() {
292                let (c_ori, _c) = b
293                    .coboundary()
294                    .find(|(_, e)| e.index() == edge.index())
295                    .expect("edge must be on its boundary's coboundary");
296                assert_eq!(
297                    b_ori, c_ori,
298                    "boundary and coboundary orientations should agree"
299                );
300            }
301
302            for (b_ori, b) in edge.coboundary() {
303                let (c_ori, _c) = b
304                    .boundary()
305                    .find(|(_, e)| e.index() == edge.index())
306                    .expect("edge must be on its coboundary's boundary");
307                assert_eq!(
308                    b_ori, c_ori,
309                    "boundary and coboundary orientations should agree"
310                );
311            }
312        }
313
314        // ensure simplex views are accessing the correct data
315
316        // volume
317        itertools::assert_equal(
318            mesh.simplices::<0>().map(|s| s.volume()),
319            mesh.simplices[0].volumes.iter().cloned(),
320        );
321        itertools::assert_equal(
322            mesh.simplices::<1>().map(|s| s.volume()),
323            mesh.simplices[1].volumes.iter().cloned(),
324        );
325        itertools::assert_equal(
326            mesh.simplices::<2>().map(|s| s.volume()),
327            mesh.simplices[2].volumes.iter().cloned(),
328        );
329        // dual volume
330        itertools::assert_equal(
331            mesh.simplices::<0>().map(|s| s.dual_volume()),
332            mesh.simplices[0].dual_volumes.iter().cloned(),
333        );
334        itertools::assert_equal(
335            mesh.simplices::<1>().map(|s| s.dual_volume()),
336            mesh.simplices[1].dual_volumes.iter().cloned(),
337        );
338        itertools::assert_equal(
339            mesh.simplices::<2>().map(|s| s.dual_volume()),
340            mesh.simplices[2].dual_volumes.iter().cloned(),
341        );
342        // circumcenter
343        itertools::assert_equal(
344            mesh.simplices::<0>().map(|s| s.circumcenter()),
345            mesh.simplices[0].circumcenters.iter().cloned(),
346        );
347        itertools::assert_equal(
348            mesh.simplices::<1>().map(|s| s.circumcenter()),
349            mesh.simplices[1].circumcenters.iter().cloned(),
350        );
351        itertools::assert_equal(
352            mesh.simplices::<2>().map(|s| s.circumcenter()),
353            mesh.simplices[2].circumcenters.iter().cloned(),
354        );
355        // barycenter
356        itertools::assert_equal(
357            mesh.simplices::<0>().map(|s| s.barycenter()),
358            mesh.simplices[0].barycenters.iter().cloned(),
359        );
360        itertools::assert_equal(
361            mesh.simplices::<1>().map(|s| s.barycenter()),
362            mesh.simplices[1].barycenters.iter().cloned(),
363        );
364        itertools::assert_equal(
365            mesh.simplices::<2>().map(|s| s.barycenter()),
366            mesh.simplices[2].barycenters.iter().cloned(),
367        );
368    }
369
370    #[test]
371    fn dual_cell_views() {
372        let mesh = tiny_mesh_3d();
373
374        // check that the dual of a dual is the original simplex
375
376        for edge in mesh.simplices::<1>() {
377            let dual_face = edge.dual();
378            let edge_again = dual_face.dual();
379            assert_eq!(edge.index(), edge_again.index());
380            itertools::assert_equal(edge.vertices(), edge_again.vertices());
381            assert_eq!(edge.dual_volume(), dual_face.volume());
382        }
383
384        // check that dual cell and simplex iterators agree
385
386        for (vert, dual_vol) in izip!(mesh.simplices::<0>(), mesh.dual_cells::<3>()) {
387            itertools::assert_equal(vert.vertex_indices(), dual_vol.dual().vertex_indices());
388        }
389        for (edge, dual_face) in izip!(mesh.simplices::<1>(), mesh.dual_cells::<2>()) {
390            itertools::assert_equal(edge.vertex_indices(), dual_face.dual().vertex_indices());
391        }
392        for (face, dual_edge) in izip!(mesh.simplices::<2>(), mesh.dual_cells::<1>()) {
393            itertools::assert_equal(face.vertex_indices(), dual_edge.dual().vertex_indices());
394        }
395        for (vol, dual_vert) in izip!(mesh.simplices::<3>(), mesh.dual_cells::<0>()) {
396            itertools::assert_equal(vol.vertex_indices(), dual_vert.dual().vertex_indices());
397        }
398    }
399}
dexterior_core/mesh/views.rs

dexterior_core/mesh/
views.rs
