mesh_sieve/data/
coordinates.rs

1//! Geometry/coordinates storage for mesh points.
2//!
3//! Coordinates are stored in a `Section` with fixed topological and embedding
4//! dimensions per point. Optional
5//! higher-order geometry data can be stored per entity (typically per-cell).
6
7use crate::data::atlas::Atlas;
8use crate::data::section::Section;
9use crate::data::storage::Storage;
10use crate::mesh_error::MeshSieveError;
11use crate::topology::point::PointId;
12
13/// Higher-order coordinate storage (e.g., per-cell geometry DOFs).
14#[derive(Clone, Debug)]
15pub struct HighOrderCoordinates<V, S: Storage<V>> {
16    dimension: usize,
17    section: Section<V, S>,
18}
19
20impl<V, S> HighOrderCoordinates<V, S>
21where
22    S: Storage<V>,
23{
24    /// Returns the spatial dimension per coordinate tuple.
25    #[inline]
26    pub fn dimension(&self) -> usize {
27        self.dimension
28    }
29
30    /// Returns a read-only reference to the underlying section.
31    #[inline]
32    pub fn section(&self) -> &Section<V, S> {
33        &self.section
34    }
35
36    /// Returns a mutable reference to the underlying section.
37    #[inline]
38    pub fn section_mut(&mut self) -> &mut Section<V, S> {
39        &mut self.section
40    }
41}
42
43impl<V, S> HighOrderCoordinates<V, S>
44where
45    V: Clone + Default,
46    S: Storage<V> + Clone,
47{
48    /// Construct a new higher-order coordinate section with a fixed dimension.
49    ///
50    /// Each entry length must be a non-zero multiple of `dimension`.
51    pub fn try_new(dimension: usize, atlas: Atlas) -> Result<Self, MeshSieveError> {
52        validate_high_order_dimension(dimension, &atlas)?;
53        Ok(Self {
54            dimension,
55            section: Section::new(atlas),
56        })
57    }
58
59    /// Wrap an existing section as higher-order coordinates, validating slice lengths.
60    pub fn from_section(dimension: usize, section: Section<V, S>) -> Result<Self, MeshSieveError> {
61        validate_high_order_dimension(dimension, section.atlas())?;
62        Ok(Self { dimension, section })
63    }
64}
65
66/// Velocity storage aligned with coordinate dimensions.
67#[derive(Clone, Debug)]
68pub struct MeshVelocity<V, S: Storage<V>> {
69    dimension: usize,
70    section: Section<V, S>,
71}
72
73impl<V, S> MeshVelocity<V, S>
74where
75    S: Storage<V>,
76{
77    /// Returns the spatial dimension per velocity tuple.
78    #[inline]
79    pub fn dimension(&self) -> usize {
80        self.dimension
81    }
82
83    /// Returns a read-only reference to the underlying section.
84    #[inline]
85    pub fn section(&self) -> &Section<V, S> {
86        &self.section
87    }
88
89    /// Returns a mutable reference to the underlying section.
90    #[inline]
91    pub fn section_mut(&mut self) -> &mut Section<V, S> {
92        &mut self.section
93    }
94
95    /// Read-only view of the velocity slice for a point `p`.
96    #[inline]
97    pub fn try_restrict(&self, p: PointId) -> Result<&[V], MeshSieveError> {
98        self.section.try_restrict(p)
99    }
100
101    /// Mutable view of the velocity slice for a point `p`.
102    #[inline]
103    pub fn try_restrict_mut(&mut self, p: PointId) -> Result<&mut [V], MeshSieveError> {
104        self.section.try_restrict_mut(p)
105    }
106}
107
108impl<V, S> MeshVelocity<V, S>
109where
110    V: Clone + Default,
111    S: Storage<V> + Clone,
112{
113    /// Construct a new velocity section with a fixed dimension.
114    ///
115    /// The provided `atlas` must store slices of length `dimension` for all points.
116    pub fn try_new(dimension: usize, atlas: Atlas) -> Result<Self, MeshSieveError> {
117        validate_dimension(dimension, &atlas)?;
118        Ok(Self {
119            dimension,
120            section: Section::new(atlas),
121        })
122    }
123
124    /// Wrap an existing section as velocity data, validating slice lengths.
125    pub fn from_section(dimension: usize, section: Section<V, S>) -> Result<Self, MeshSieveError> {
126        validate_dimension(dimension, section.atlas())?;
127        Ok(Self { dimension, section })
128    }
129
130    /// Adds a new point with the configured velocity dimension.
131    pub fn try_add_point(&mut self, p: PointId) -> Result<(), MeshSieveError> {
132        self.section.try_add_point(p, self.dimension)
133    }
134}
135
136/// Coordinate storage with an attached topological and embedding dimension.
137#[derive(Clone, Debug)]
138pub struct Coordinates<V, S: Storage<V>> {
139    topological_dimension: usize,
140    embedding_dimension: usize,
141    section: Section<V, S>,
142    high_order: Option<HighOrderCoordinates<V, S>>,
143}
144
145impl<V, S> Coordinates<V, S>
146where
147    S: Storage<V>,
148{
149    /// Returns the embedding dimension per point.
150    #[inline]
151    pub fn dimension(&self) -> usize {
152        self.embedding_dimension
153    }
154
155    /// Returns the topological dimension of the mesh.
156    #[inline]
157    pub fn topological_dimension(&self) -> usize {
158        self.topological_dimension
159    }
160
161    /// Returns the embedding dimension per point.
162    #[inline]
163    pub fn embedding_dimension(&self) -> usize {
164        self.embedding_dimension
165    }
166
167    /// Returns a read-only reference to the underlying section.
168    #[inline]
169    pub fn section(&self) -> &Section<V, S> {
170        &self.section
171    }
172
173    /// Returns a mutable reference to the underlying section.
174    #[inline]
175    pub fn section_mut(&mut self) -> &mut Section<V, S> {
176        &mut self.section
177    }
178
179    /// Optional higher-order coordinate data (e.g., per-cell geometry DOFs).
180    #[inline]
181    pub fn high_order(&self) -> Option<&HighOrderCoordinates<V, S>> {
182        self.high_order.as_ref()
183    }
184
185    /// Mutable view of optional higher-order coordinate data.
186    #[inline]
187    pub fn high_order_mut(&mut self) -> Option<&mut HighOrderCoordinates<V, S>> {
188        self.high_order.as_mut()
189    }
190
191    /// Attach higher-order coordinate data, validating the dimension.
192    pub fn set_high_order(
193        &mut self,
194        high_order: HighOrderCoordinates<V, S>,
195    ) -> Result<(), MeshSieveError> {
196        if high_order.dimension != self.embedding_dimension {
197            return Err(MeshSieveError::InvalidGeometry(format!(
198                "higher-order coordinate dimension {} does not match base dimension {}",
199                high_order.dimension, self.embedding_dimension
200            )));
201        }
202        self.high_order = Some(high_order);
203        Ok(())
204    }
205
206    /// Consumes the wrapper and returns the underlying section.
207    #[inline]
208    pub fn into_section(self) -> Section<V, S> {
209        self.section
210    }
211
212    /// Read-only view of the coordinate slice for a point `p`.
213    #[inline]
214    pub fn try_restrict(&self, p: PointId) -> Result<&[V], MeshSieveError> {
215        self.section.try_restrict(p)
216    }
217
218    /// Mutable view of the coordinate slice for a point `p`.
219    #[inline]
220    pub fn try_restrict_mut(&mut self, p: PointId) -> Result<&mut [V], MeshSieveError> {
221        self.section.try_restrict_mut(p)
222    }
223}
224
225impl<S> Coordinates<f64, S>
226where
227    S: Storage<f64>,
228{
229    /// Advance coordinates using a velocity field and timestep.
230    pub fn advance_with_velocity<St>(
231        &mut self,
232        velocity: &MeshVelocity<f64, St>,
233        dt: f64,
234    ) -> Result<(), MeshSieveError>
235    where
236        St: Storage<f64>,
237    {
238        let dim = self.embedding_dimension;
239        let points: Vec<PointId> = self.section.atlas().points().collect();
240        for point in points {
241            let vel = velocity.try_restrict(point)?;
242            if vel.len() != dim {
243                return Err(MeshSieveError::SliceLengthMismatch {
244                    point,
245                    expected: dim,
246                    found: vel.len(),
247                });
248            }
249            let coord = self.try_restrict_mut(point)?;
250            if coord.len() != dim {
251                return Err(MeshSieveError::SliceLengthMismatch {
252                    point,
253                    expected: dim,
254                    found: coord.len(),
255                });
256            }
257            for (coord_value, vel_value) in coord.iter_mut().zip(vel.iter()) {
258                *coord_value += dt * vel_value;
259            }
260        }
261        Ok(())
262    }
263}
264
265impl<V, S> Coordinates<V, S>
266where
267    V: Clone + Default,
268    S: Storage<V> + Clone,
269{
270    /// Construct a new coordinates section with fixed topological and embedding dimensions.
271    ///
272    /// The provided `atlas` must store slices of length `embedding_dimension`
273    /// for all points.
274    pub fn try_new(
275        topological_dimension: usize,
276        embedding_dimension: usize,
277        atlas: Atlas,
278    ) -> Result<Self, MeshSieveError> {
279        validate_coordinate_dimensions(topological_dimension, embedding_dimension, &atlas)?;
280        Ok(Self {
281            topological_dimension,
282            embedding_dimension,
283            section: Section::new(atlas),
284            high_order: None,
285        })
286    }
287
288    /// Wrap an existing section as coordinates, validating slice lengths.
289    pub fn from_section(
290        topological_dimension: usize,
291        embedding_dimension: usize,
292        section: Section<V, S>,
293    ) -> Result<Self, MeshSieveError> {
294        validate_coordinate_dimensions(
295            topological_dimension,
296            embedding_dimension,
297            section.atlas(),
298        )?;
299        Ok(Self {
300            topological_dimension,
301            embedding_dimension,
302            section,
303            high_order: None,
304        })
305    }
306
307    /// Adds a new point with the configured coordinate dimension.
308    pub fn try_add_point(&mut self, p: PointId) -> Result<(), MeshSieveError> {
309        self.section.try_add_point(p, self.embedding_dimension)
310    }
311}
312
313fn validate_coordinate_dimensions(
314    topological_dimension: usize,
315    embedding_dimension: usize,
316    atlas: &Atlas,
317) -> Result<(), MeshSieveError> {
318    if embedding_dimension == 0 {
319        return Err(MeshSieveError::ZeroLengthSlice);
320    }
321    if topological_dimension > embedding_dimension {
322        return Err(MeshSieveError::InvalidGeometry(format!(
323            "topological dimension {topological_dimension} exceeds embedding dimension {embedding_dimension}"
324        )));
325    }
326    for (point, (_offset, len)) in atlas.iter_entries() {
327        if len != embedding_dimension {
328            return Err(MeshSieveError::SliceLengthMismatch {
329                point,
330                expected: embedding_dimension,
331                found: len,
332            });
333        }
334    }
335    Ok(())
336}
337
338fn validate_dimension(dimension: usize, atlas: &Atlas) -> Result<(), MeshSieveError> {
339    if dimension == 0 {
340        return Err(MeshSieveError::ZeroLengthSlice);
341    }
342    for (point, (_offset, len)) in atlas.iter_entries() {
343        if len != dimension {
344            return Err(MeshSieveError::SliceLengthMismatch {
345                point,
346                expected: dimension,
347                found: len,
348            });
349        }
350    }
351    Ok(())
352}
353
354fn validate_high_order_dimension(dimension: usize, atlas: &Atlas) -> Result<(), MeshSieveError> {
355    if dimension == 0 {
356        return Err(MeshSieveError::ZeroLengthSlice);
357    }
358    for (point, (_offset, len)) in atlas.iter_entries() {
359        if len == 0 || len % dimension != 0 {
360            return Err(MeshSieveError::SliceLengthMismatch {
361                point,
362                expected: dimension,
363                found: len,
364            });
365        }
366    }
367    Ok(())
368}
369
370#[cfg(test)]
371mod tests {
372    use super::{Coordinates, MeshVelocity};
373    use crate::data::atlas::Atlas;
374    use crate::data::storage::VecStorage;
375    use crate::topology::point::PointId;
376
377    #[test]
378    fn advance_coordinates_over_multiple_steps() {
379        let mut atlas = Atlas::default();
380        let p1 = PointId::new(1).unwrap();
381        let p2 = PointId::new(2).unwrap();
382        atlas.try_insert(p1, 3).unwrap();
383        atlas.try_insert(p2, 3).unwrap();
384
385        let mut coords = Coordinates::<f64, VecStorage<f64>>::try_new(3, 3, atlas.clone()).unwrap();
386        let mut velocity = MeshVelocity::<f64, VecStorage<f64>>::try_new(3, atlas).unwrap();
387
388        coords.section_mut().try_set(p1, &[0.0, 0.0, 0.0]).unwrap();
389        coords.section_mut().try_set(p2, &[1.0, 1.0, 1.0]).unwrap();
390        velocity
391            .section_mut()
392            .try_set(p1, &[1.0, 0.0, -1.0])
393            .unwrap();
394        velocity
395            .section_mut()
396            .try_set(p2, &[0.5, -0.5, 1.0])
397            .unwrap();
398
399        let dt = 0.25;
400        for _ in 0..4 {
401            coords.advance_with_velocity(&velocity, dt).unwrap();
402        }
403
404        assert_eq!(coords.try_restrict(p1).unwrap(), &[1.0, 0.0, -1.0]);
405        assert_eq!(coords.try_restrict(p2).unwrap(), &[1.5, 0.5, 2.0]);
406    }
407}
mesh_sieve/data/coordinates.rs

mesh_sieve/data/
coordinates.rs
