mesh_sieve/data/refine/
sieved_array.rs

1//! A generic array of values indexed by mesh points, supporting refine/assemble.
2//! (Extracted from section.rs)
3//!
4//! This module provides [`SievedArray`], a flexible structure for storing and
5//! manipulating per-point data in mesh refinement and assembly operations.
6
7use crate::data::atlas::Atlas;
8use crate::data::refine::delta::SliceDelta;
9use crate::topology::arrow::Polarity;
10use crate::topology::point::PointId;
11
12/// A generic array of values indexed by mesh points, supporting refinement and assembly.
13///
14/// # Type Parameters
15/// - `P`: Point identifier type (must convert to [`PointId`]).
16/// - `V`: Value type stored for each point.
17#[derive(Clone, Debug)]
18pub struct SievedArray<P, V> {
19    pub(crate) atlas: Atlas,
20    pub(crate) data: Vec<V>,
21    _phantom: std::marker::PhantomData<P>,
22}
23
24impl<P, V> SievedArray<P, V>
25where
26    P: Into<PointId> + Copy + Eq,
27{
28    /// Access the underlying atlas for point-to-slice mapping.
29    pub fn atlas(&self) -> &Atlas {
30        &self.atlas
31    }
32
33    /// Get a read-only slice for the given point, or an error if not present.
34    pub fn try_get(&self, p: PointId) -> Result<&[V], crate::mesh_error::MeshSieveError> {
35        let (off, len) = self
36            .atlas
37            .get(p)
38            .ok_or(crate::mesh_error::MeshSieveError::SievedArrayPointNotInAtlas(p))?;
39        Ok(&self.data[off..off + len])
40    }
41    /// Get a mutable slice for the given point, or an error if not present.
42    pub fn try_get_mut(
43        &mut self,
44        p: PointId,
45    ) -> Result<&mut [V], crate::mesh_error::MeshSieveError> {
46        let (off, len) = self
47            .atlas
48            .get(p)
49            .ok_or(crate::mesh_error::MeshSieveError::SievedArrayPointNotInAtlas(p))?;
50        Ok(&mut self.data[off..off + len])
51    }
52    /// Fallible iterator over `(PointId, &[V])` in atlas insertion order.
53    ///
54    /// # Complexity
55    /// **O(n)** to traverse, **O(1)** per element.
56    ///
57    /// # Determinism
58    /// Iteration strictly follows atlas insertion order.
59    #[inline]
60    pub fn try_iter_in_order(
61        &self,
62    ) -> impl Iterator<Item = Result<(PointId, &[V]), crate::mesh_error::MeshSieveError>> + '_ {
63        self.atlas
64            .points()
65            .map(move |pid| self.try_get(pid).map(|sl| (pid, sl)))
66    }
67
68    /// Non-fallible iterator over `(PointId, &[V])` in atlas insertion order.
69    ///
70    /// # Complexity
71    /// **O(n)** to traverse, **O(1)** per element.
72    ///
73    /// # Determinism
74    /// Iteration strictly follows atlas insertion order.
75    #[inline]
76    pub fn iter_in_order(&self) -> impl Iterator<Item = (PointId, &[V])> + '_ {
77        self.atlas.points().map(move |pid| {
78            let (off, len) = self.atlas.get(pid).expect("atlas missing point");
79            let sl = &self.data[off..off + len];
80            (pid, sl)
81        })
82    }
83}
84
85impl<P, V: Clone> SievedArray<P, V>
86where
87    P: Into<PointId> + Copy + Eq,
88{
89    /// Set the values for the given point from a slice, or return an error if lengths mismatch or point not found.
90    pub fn try_set(
91        &mut self,
92        p: PointId,
93        val: &[V],
94    ) -> Result<(), crate::mesh_error::MeshSieveError> {
95        let tgt = self.try_get_mut(p)?;
96        if tgt.len() != val.len() {
97            return Err(
98                crate::mesh_error::MeshSieveError::SievedArraySliceLengthMismatch {
99                    point: p,
100                    expected: tgt.len(),
101                    found: val.len(),
102                },
103            );
104        }
105        tgt.clone_from_slice(val);
106        Ok(())
107    }
108}
109
110impl<P, V: Clone + Default> SievedArray<P, V>
111where
112    P: Into<PointId> + Copy + Eq,
113{
114    /// Create a new `SievedArray` with the given atlas.
115    pub fn new(atlas: Atlas) -> Self {
116        let data = vec![V::default(); atlas.total_len()];
117        Self {
118            atlas,
119            data,
120            _phantom: std::marker::PhantomData,
121        }
122    }
123
124    /// Refine this array from a coarse array using a sifter (with orientations).
125    ///
126    /// The `refinement` mapping must be a partial function from fine points to
127    /// coarse points; any fine point appearing more than once results in
128    /// [`MeshSieveError::DuplicateRefinementTarget`].
129    ///
130    /// # Complexity
131    /// **O(m · k)**, where `m` is the number of fine targets and `k` is slice length.
132    ///
133    /// # Determinism
134    /// Serial: deterministic. Parallel: deterministic final state if refinement
135    /// has **no duplicate fine targets** (duplicates are rejected). Writes occur
136    /// after a read-only computation phase.
137    pub fn try_refine_with_sifter(
138        &mut self,
139        coarse: &SievedArray<P, V>,
140        refinement: &[(P, Vec<(P, Polarity)>)],
141    ) -> Result<(), crate::mesh_error::MeshSieveError> {
142        use crate::mesh_error::MeshSieveError;
143
144        let mut updates = Vec::<(P, Vec<V>)>::new();
145        for (coarse_pt, fine_pts) in refinement.iter() {
146            let cpid = (*coarse_pt).into();
147            let coarse_slice = coarse.try_get(cpid)?;
148            for (fine_pt, orient) in fine_pts.iter() {
149                let fpid = (*fine_pt).into();
150                let (_off, len) = self
151                    .atlas
152                    .get(fpid)
153                    .ok_or(MeshSieveError::SievedArrayPointNotInAtlas(fpid))?;
154                if coarse_slice.len() != len {
155                    return Err(MeshSieveError::SievedArraySliceLengthMismatch {
156                        point: fpid,
157                        expected: coarse_slice.len(),
158                        found: len,
159                    });
160                }
161                let mut data = vec![V::default(); len];
162                orient.apply(coarse_slice, &mut data)?;
163                updates.push((*fine_pt, data));
164            }
165        }
166
167        updates.sort_unstable_by_key(|(f, _)| (*f).into());
168        for w in updates.windows(2) {
169            let f0: PointId = (w[0].0).into();
170            let f1: PointId = (w[1].0).into();
171            if f0 == f1 {
172                return Err(MeshSieveError::DuplicateRefinementTarget { fine: f0 });
173            }
174        }
175
176        for (fine_pt, data) in updates {
177            let dst = self.try_get_mut(fine_pt.into())?;
178            debug_assert_eq!(dst.len(), data.len());
179            dst.clone_from_slice(&data);
180        }
181        Ok(())
182    }
183
184    /// Refine this array from a coarse array using a simple mapping (all forward), propagating errors.
185    ///
186    /// # Complexity
187    /// **O(m · k)**, where `m` is the number of fine targets and `k` is slice length.
188    ///
189    /// # Determinism
190    /// Serial: deterministic. Parallel: deterministic final state if refinement
191    /// has **no duplicate fine targets** (duplicates are rejected). Writes occur
192    /// after a read-only computation phase.
193    pub fn try_refine(
194        &mut self,
195        coarse: &SievedArray<P, V>,
196        refinement: &[(P, Vec<P>)],
197    ) -> Result<(), crate::mesh_error::MeshSieveError> {
198        let sifter: Vec<_> = refinement
199            .iter()
200            .map(|(c, fs)| (*c, fs.iter().map(|f| (*f, Polarity::Forward)).collect()))
201            .collect();
202        self.try_refine_with_sifter(coarse, &sifter)
203    }
204}
205
206impl<P, V> SievedArray<P, V>
207where
208    P: Into<PointId> + Copy + Eq,
209    V: num_traits::FromPrimitive
210        + std::ops::AddAssign
211        + std::ops::Div<Output = V>
212        + Clone
213        + Default,
214{
215    /// Assemble fine data into coarse by averaging over refinement, propagating errors.
216    ///
217    /// # Complexity
218    /// **O(m · k)**; performs element-wise reduction per coarse point.
219    ///
220    /// # Determinism
221    /// Determined by the reduction order, which is fixed by the refinement input order.
222    pub fn try_assemble(
223        &self,
224        coarse: &mut SievedArray<P, V>,
225        refinement: &[(P, Vec<P>)],
226    ) -> Result<(), crate::mesh_error::MeshSieveError> {
227        for (coarse_pt, fine_pts) in refinement.iter() {
228            let mut accum = {
229                let coarse_slice = coarse.try_get((*coarse_pt).into())?;
230                vec![V::default(); coarse_slice.len()]
231            };
232            let mut count = 0;
233            for fine_pt in fine_pts {
234                let slice = self.try_get((*fine_pt).into())?;
235                if slice.len() != accum.len() {
236                    return Err(
237                        crate::mesh_error::MeshSieveError::SievedArraySliceLengthMismatch {
238                            point: (*fine_pt).into(),
239                            expected: accum.len(),
240                            found: slice.len(),
241                        },
242                    );
243                }
244                for (a, v) in accum.iter_mut().zip(slice.iter()) {
245                    *a += v.clone();
246                }
247                count += 1;
248            }
249            if count > 0 {
250                let divisor: V = num_traits::FromPrimitive::from_usize(count).ok_or(
251                    crate::mesh_error::MeshSieveError::SievedArrayPrimitiveConversionFailure(count),
252                )?;
253                for a in accum.iter_mut() {
254                    *a = a.clone() / divisor.clone();
255                }
256                coarse.try_set((*coarse_pt).into(), &accum)?;
257            }
258        }
259        Ok(())
260    }
261}
262
263#[cfg(feature = "rayon")]
264use rayon::prelude::*;
265
266impl<P, V: Clone + Default + Send + Sync> SievedArray<P, V>
267where
268    P: Into<PointId> + Copy + Eq + Send + Sync,
269{
270    /// Parallel refinement using a sifter, enabled with the `rayon` feature.
271    ///
272    /// Computes slice updates in parallel, short-circuiting on the first error
273    /// and rejecting duplicate fine targets deterministically.
274    ///
275    /// # Complexity
276    /// **O(m · k)**, where `m` is the number of fine targets and `k` is slice length.
277    /// Parallel variant short-circuits on first error.
278    ///
279    /// # Determinism
280    /// Serial: deterministic. Parallel: deterministic final state if refinement
281    /// has **no duplicate fine targets** (duplicates are rejected). Writes occur
282    /// after a read-only computation phase.
283    #[cfg(feature = "rayon")]
284    pub fn try_refine_with_sifter_parallel(
285        &mut self,
286        coarse: &Self,
287        refinement: &[(P, Vec<(P, Polarity)>)],
288    ) -> Result<(), crate::mesh_error::MeshSieveError> {
289        use crate::mesh_error::MeshSieveError;
290        use std::collections::HashMap;
291
292        let fine_spans: HashMap<PointId, (usize, usize)> = self
293            .atlas
294            .iter_entries()
295            .map(|(pid, span)| (pid, span))
296            .collect();
297
298        let updates: Vec<(P, Vec<V>)> = refinement
299            .par_iter()
300            .try_fold(
301                || Vec::<(P, Vec<V>)>::new(),
302                |mut local, (coarse_pt, fine_pts)| -> Result<_, MeshSieveError> {
303                    let cpid = (*coarse_pt).into();
304                    let coarse_slice = coarse.try_get(cpid)?;
305                    for (fine_pt, orient) in fine_pts {
306                        let fpid = (*fine_pt).into();
307                        let (_off, len) = fine_spans
308                            .get(&fpid)
309                            .copied()
310                            .ok_or(MeshSieveError::SievedArrayPointNotInAtlas(fpid))?;
311                        if coarse_slice.len() != len {
312                            return Err(MeshSieveError::SievedArraySliceLengthMismatch {
313                                point: fpid,
314                                expected: coarse_slice.len(),
315                                found: len,
316                            });
317                        }
318                        let mut data = vec![V::default(); len];
319                        orient.apply(coarse_slice, &mut data)?;
320                        local.push((*fine_pt, data));
321                    }
322                    Ok(local)
323                },
324            )
325            .try_reduce(
326                || Vec::<(P, Vec<V>)>::new(),
327                |mut a, mut b| -> Result<_, MeshSieveError> {
328                    a.append(&mut b);
329                    Ok(a)
330                },
331            )?;
332
333        let mut updates = updates;
334        updates.sort_unstable_by_key(|(f, _)| (*f).into());
335        for w in updates.windows(2) {
336            let f0: PointId = (w[0].0).into();
337            let f1: PointId = (w[1].0).into();
338            if f0 == f1 {
339                return Err(MeshSieveError::DuplicateRefinementTarget { fine: f0 });
340            }
341        }
342
343        for (fine_pt, data) in updates {
344            let dst = self.try_get_mut(fine_pt.into())?;
345            debug_assert_eq!(dst.len(), data.len());
346            dst.clone_from_slice(&data);
347        }
348        Ok(())
349    }
350}
351
352#[cfg(test)]
353mod tests {
354    use crate::data::atlas::Atlas;
355    use crate::data::refine::sieved_array::SievedArray;
356    use crate::mesh_error::MeshSieveError;
357    use crate::topology::arrow::Polarity;
358    use crate::topology::point::PointId;
359
360    fn pt(i: u64) -> PointId {
361        PointId::new(i).unwrap()
362    }
363    fn make_sieved() -> SievedArray<PointId, i32> {
364        let mut atlas = Atlas::default();
365        atlas.try_insert(pt(1), 2).unwrap();
366        atlas.try_insert(pt(2), 2).unwrap();
367        atlas.try_insert(pt(3), 2).unwrap();
368        SievedArray::new(atlas)
369    }
370
371    #[test]
372    fn sieved_array_basic_get_set_iter() {
373        let mut atlas = Atlas::default();
374        atlas.try_insert(pt(1), 2).unwrap();
375        atlas.try_insert(pt(2), 1).unwrap();
376        let mut arr = SievedArray::<PointId, i32>::new(atlas);
377        arr.try_set(pt(1), &[1, 2]).unwrap();
378        arr.try_set(pt(2), &[3]).unwrap();
379        assert_eq!(arr.try_get(pt(1)).unwrap(), &[1, 2]);
380        assert_eq!(arr.try_get(pt(2)).unwrap(), &[3]);
381        let vals: Vec<_> = arr.try_iter_in_order().map(|r| r.unwrap().1[0]).collect();
382        assert_eq!(vals, vec![1, 3]);
383    }
384
385    #[test]
386    fn sieved_array_refine_with_sifter_forward_and_reverse() {
387        let mut cat = Atlas::default();
388        cat.try_insert(pt(1), 2).unwrap();
389        let mut fat = Atlas::default();
390        fat.try_insert(pt(2), 2).unwrap();
391        fat.try_insert(pt(3), 2).unwrap();
392        let mut coarse = SievedArray::new(cat);
393        let mut fine = SievedArray::new(fat);
394        coarse.try_set(pt(1), &[10, 20]).unwrap();
395        let refinement = vec![(
396            pt(1),
397            vec![(pt(2), Polarity::Forward), (pt(3), Polarity::Reverse)],
398        )];
399        fine.try_refine_with_sifter(&coarse, &refinement).unwrap();
400        assert_eq!(fine.try_get(pt(2)).unwrap(), &[10, 20]);
401        assert_eq!(fine.try_get(pt(3)).unwrap(), &[20, 10]);
402    }
403
404    #[test]
405    fn sieved_array_refine_forward_only() {
406        let mut coarse = make_sieved();
407        let mut fine = make_sieved();
408        coarse.try_set(pt(1), &[5, 6]).unwrap();
409        fine.try_refine(&coarse, &[(pt(1), vec![pt(2), pt(3)])])
410            .unwrap();
411        assert_eq!(fine.try_get(pt(2)).unwrap(), &[5, 6]);
412        assert_eq!(fine.try_get(pt(3)).unwrap(), &[5, 6]);
413    }
414
415    #[test]
416    fn sieved_array_assemble_average() {
417        let mut coarse = make_sieved();
418        let mut fine = make_sieved();
419        fine.try_set(pt(1), &[2, 4]).unwrap();
420        fine.try_set(pt(2), &[6, 8]).unwrap();
421        fine.try_assemble(&mut coarse, &[(pt(3), vec![pt(1), pt(2)])])
422            .unwrap();
423        assert_eq!(coarse.try_get(pt(3)).unwrap(), &[4, 6]);
424    }
425
426    #[test]
427    fn sieved_array_set_wrong_length_error() {
428        let mut arr = make_sieved();
429        let err = arr.try_set(pt(1), &[1]).unwrap_err();
430        match err {
431            MeshSieveError::SievedArraySliceLengthMismatch {
432                point,
433                expected,
434                found,
435            } => {
436                assert_eq!(point, pt(1));
437                assert_eq!(expected, 2);
438                assert_eq!(found, 1);
439            }
440            _ => panic!("wrong error variant: {err:?}"),
441        }
442    }
443
444    #[test]
445    fn sieved_array_assemble_mismatch_error() {
446        use crate::data::atlas::Atlas;
447        let mut coarse_atlas = Atlas::default();
448        let mut fine_atlas = Atlas::default();
449        // pt(1) has length 2 in coarse, 1 in fine
450        coarse_atlas.try_insert(pt(1), 2).unwrap();
451        fine_atlas.try_insert(pt(1), 1).unwrap();
452        let mut coarse = SievedArray::<PointId, i32>::new(coarse_atlas);
453        let fine = SievedArray::new(fine_atlas);
454        let err = fine
455            .try_assemble(&mut coarse, &[(pt(1), vec![pt(1)])])
456            .unwrap_err();
457        match err {
458            MeshSieveError::SievedArraySliceLengthMismatch {
459                point,
460                expected,
461                found,
462            } => {
463                assert_eq!(point, pt(1));
464                assert_eq!(expected, 2);
465                assert_eq!(found, 1);
466            }
467            _ => panic!("wrong error variant: {err:?}"),
468        }
469    }
470
471    #[test]
472    fn sieved_array_point_not_in_atlas_error() {
473        let arr = make_sieved();
474        let missing = pt(99);
475        let err = arr.try_get(missing).unwrap_err();
476        match err {
477            MeshSieveError::SievedArrayPointNotInAtlas(p) => assert_eq!(p, missing),
478            _ => panic!("wrong error variant: {err:?}"),
479        }
480    }
481
482    #[cfg(feature = "rayon")]
483    #[test]
484    fn sieved_array_refine_with_sifter_parallel_works() {
485        let mut coarse = make_sieved();
486        let mut fine = make_sieved();
487        coarse.try_set(pt(1), &[2, 3]).unwrap();
488        let refinement = vec![(pt(1), vec![(pt(2), Polarity::Forward)])];
489        fine.try_refine_with_sifter_parallel(&coarse, &refinement)
490            .expect("parallel refinement failed");
491        assert_eq!(fine.try_get(pt(2)).unwrap(), &[2, 3]);
492    }
493}
mesh_sieve/data/refine/sieved_array.rs

mesh_sieve/data/refine/
sieved_array.rs
