Skip to main content

conspire/math/tensor/rank_1/vec/
mod.rs

1#[cfg(test)]
2mod test;
3use crate::math::{Current, Reference};
4use crate::units::{Dimensionless, UnitDiv};
5
6use crate::math::{
7    Jacobian, Quantity, Solution, Tensor, TensorRank0, TensorRank1, TensorRank1List,
8    TensorRank2SparseVec2D, TensorRank2SparseVec2DSymmetric, TensorRank2Vec2D, TensorVec, Vector,
9    tensor::vec::TensorVector,
10};
11use std::{
12    array::from_fn,
13    mem::forget,
14    ops::{Div, Sub},
15};
16
17use crate::math::assert::FiniteDifference;
18
19/// A vector of rank-1 tensors.
20pub type TensorRank1Vec<const D: usize, I, U = Dimensionless> = TensorVector<TensorRank1<D, I, U>>;
21
22impl<const D: usize, I, U> TensorRank1Vec<D, I, U> {
23    pub fn bounding_box(&self) -> TensorRank1List<D, I, 2, U> {
24        self.iter()
25            .skip(1)
26            .fold(
27                [self[0].clone(), self[0].clone()],
28                |[mut min, mut max], entry| {
29                    entry
30                        .iter()
31                        .zip(min.iter_mut().zip(max.iter_mut()))
32                        .for_each(|(&entry_i, (min_i, max_i))| {
33                            *min_i = min_i.min(entry_i);
34                            *max_i = max_i.max(entry_i);
35                        });
36                    [min, max]
37                },
38            )
39            .into()
40    }
41    pub fn zero(len: usize) -> Self {
42        (0..len).map(|_| super::zero()).collect()
43    }
44}
45
46impl<const D: usize, I, const N: usize, U> From<[[TensorRank0; D]; N]> for TensorRank1Vec<D, I, U> {
47    fn from(array: [[TensorRank0; D]; N]) -> Self {
48        array.into_iter().map(TensorRank1::from).collect()
49    }
50}
51
52impl<const D: usize, I, U> From<Vec<[TensorRank0; D]>> for TensorRank1Vec<D, I, U> {
53    fn from(vec: Vec<[TensorRank0; D]>) -> Self {
54        let (length, capacity) = (vec.len(), vec.capacity());
55        let pointer = vec.as_ptr() as *mut TensorRank1<D, I, U>;
56        forget(vec);
57        unsafe { Self::from(Vec::from_raw_parts(pointer, length, capacity)) }
58    }
59}
60
61impl<const D: usize, I, U> From<TensorRank1Vec<D, I, U>> for Vec<[TensorRank0; D]> {
62    fn from(tensor_rank_1_vec: TensorRank1Vec<D, I, U>) -> Self {
63        let vec = Vec::<TensorRank1<D, I, U>>::from(tensor_rank_1_vec);
64        let (length, capacity) = (vec.len(), vec.capacity());
65        let pointer = vec.as_ptr() as *mut [TensorRank0; D];
66        forget(vec);
67        unsafe { Vec::from_raw_parts(pointer, length, capacity) }
68    }
69}
70
71impl<const D: usize, I, U> From<Vec<Vec<TensorRank0>>> for TensorRank1Vec<D, I, U> {
72    fn from(vec: Vec<Vec<TensorRank0>>) -> Self {
73        vec.into_iter()
74            .map(|tensor_rank_1| tensor_rank_1.into())
75            .collect()
76    }
77}
78
79impl<const D: usize, I, U> From<TensorRank1Vec<D, I, U>> for Vec<Vec<TensorRank0>> {
80    fn from(tensor_rank_1_vec: TensorRank1Vec<D, I, U>) -> Self {
81        tensor_rank_1_vec
82            .into_iter()
83            .map(|tensor_rank_1| tensor_rank_1.into())
84            .collect()
85    }
86}
87
88impl<const D: usize, I, U> TryFrom<[Vec<TensorRank0>; D]> for TensorRank1Vec<D, I, U> {
89    type Error = String;
90    fn try_from(vec_array: [Vec<TensorRank0>; D]) -> Result<Self, Self::Error> {
91        let length = vec_array[0].len();
92        if vec_array.iter().any(|vec| vec.len() != length) {
93            Err("Vector length mismatch in type conversion".to_string())
94        } else {
95            Ok((0..length)
96                .map(|j| TensorRank1::const_from(from_fn(|i| vec_array[i][j])))
97                .collect())
98        }
99    }
100}
101
102impl<const D: usize, I, U> From<TensorRank1Vec<D, I, U>> for [Vec<TensorRank0>; D] {
103    fn from(tensor_rank_1_vec: TensorRank1Vec<D, I, U>) -> Self {
104        let length = tensor_rank_1_vec.len();
105        let mut output = from_fn(|_| Vec::with_capacity(length));
106        tensor_rank_1_vec.into_iter().for_each(|tensor_rank_1| {
107            output
108                .iter_mut()
109                .zip(tensor_rank_1)
110                .for_each(|(entry, value)| entry.push(value.value()))
111        });
112        output
113    }
114}
115
116impl<const D: usize, I, U> From<&TensorRank1Vec<D, I, U>> for [Vec<TensorRank0>; D] {
117    fn from(tensor_rank_1_vec: &TensorRank1Vec<D, I, U>) -> Self {
118        let length = tensor_rank_1_vec.len();
119        let mut output = from_fn(|_| Vec::with_capacity(length));
120        tensor_rank_1_vec.iter().for_each(|tensor_rank_1| {
121            output
122                .iter_mut()
123                .zip(tensor_rank_1.iter())
124                .for_each(|(entry, &value)| entry.push(value.value()))
125        });
126        output
127    }
128}
129
130impl<const D: usize, U> From<TensorRank1Vec<D, Reference, U>> for TensorRank1Vec<D, Current, U> {
131    fn from(tensor_rank_1_vec: TensorRank1Vec<D, Reference, U>) -> Self {
132        let (length, capacity) = (tensor_rank_1_vec.len(), tensor_rank_1_vec.capacity());
133        let pointer = tensor_rank_1_vec.as_ptr() as *mut TensorRank1<D, Current, U>;
134        forget(tensor_rank_1_vec);
135        unsafe { Self::from(Vec::from_raw_parts(pointer, length, capacity)) }
136    }
137}
138
139impl<const D: usize, U> From<&TensorRank1Vec<D, Reference, U>> for TensorRank1Vec<D, Current, U> {
140    fn from(tensor_rank_1_vec: &TensorRank1Vec<D, Reference, U>) -> Self {
141        tensor_rank_1_vec
142            .iter()
143            .map(|tensor_rank_1| tensor_rank_1.into())
144            .collect()
145    }
146}
147
148impl<const D: usize, U> From<TensorRank1Vec<D, Current, U>> for TensorRank1Vec<D, Reference, U> {
149    fn from(tensor_rank_1_vec: TensorRank1Vec<D, Current, U>) -> Self {
150        let (length, capacity) = (tensor_rank_1_vec.len(), tensor_rank_1_vec.capacity());
151        let pointer = tensor_rank_1_vec.as_ptr() as *mut TensorRank1<D, Reference, U>;
152        forget(tensor_rank_1_vec);
153        unsafe { Self::from(Vec::from_raw_parts(pointer, length, capacity)) }
154    }
155}
156
157impl<const D: usize, U> From<&TensorRank1Vec<D, Current, U>> for TensorRank1Vec<D, Reference, U> {
158    fn from(tensor_rank_1_vec: &TensorRank1Vec<D, Current, U>) -> Self {
159        tensor_rank_1_vec
160            .iter()
161            .map(|tensor_rank_1| tensor_rank_1.into())
162            .collect()
163    }
164}
165
166impl<const D: usize, I, U> From<Vector> for TensorRank1Vec<D, I, U> {
167    fn from(vector: Vector) -> Self {
168        let n = vector.len();
169        if !n.is_multiple_of(D) {
170            panic!("Vector length mismatch.")
171        } else if vector.capacity().is_multiple_of(D) {
172            let (length, capacity) = (n / D, vector.capacity() / D);
173            let pointer = vector.as_ptr() as *mut TensorRank1<D, I, U>;
174            forget(vector);
175            unsafe { Self::from(Vec::from_raw_parts(pointer, length, capacity)) }
176        } else {
177            (0..n / D)
178                .map(|i| TensorRank1::const_from(from_fn(|j| vector[D * i + j])))
179                .collect()
180        }
181    }
182}
183
184impl<const D: usize, I, U> Jacobian for TensorRank1Vec<D, I, U> {
185    fn fill_into(&self, vector: &mut Vector) {
186        self.iter()
187            .flat_map(|entry| entry.iter())
188            .zip(vector.iter_mut())
189            .for_each(|(self_i, vector_i)| *vector_i = self_i.value())
190    }
191    fn fill_into_chained(self, other: Vector, vector: &mut Vector) {
192        self.into_iter()
193            .flatten()
194            .map(|entry| entry.value())
195            .chain(other)
196            .zip(vector.iter_mut())
197            .for_each(|(self_i, vector_i)| *vector_i = self_i)
198    }
199    fn retain_from(self, retained: &[bool]) -> Vector {
200        self.into_iter()
201            .flatten()
202            .zip(retained.iter())
203            .filter(|(_, retained)| **retained)
204            .map(|(entry, _)| entry.value())
205            .collect()
206    }
207    fn zero_out(&mut self, indices: &[usize]) {
208        indices
209            .iter()
210            .for_each(|index| self[index / D][index % D] = Quantity::new(0.0))
211    }
212}
213
214impl<const D: usize, I, U> Solution for TensorRank1Vec<D, I, U> {
215    fn decrement_from(&mut self, other: &Vector) {
216        self.iter_mut()
217            .flat_map(|x| x.iter_mut())
218            .zip(other.iter())
219            .for_each(|(self_i, vector_i)| *self_i -= Quantity::new(*vector_i))
220    }
221    fn decrement_from_chained(&mut self, other: &mut Vector, vector: &Vector) {
222        let mut values = vector.iter();
223        self.iter_mut()
224            .flat_map(|x| x.iter_mut())
225            .zip(values.by_ref())
226            .for_each(|(entry_i, vector_i)| *entry_i -= Quantity::new(*vector_i));
227        other
228            .iter_mut()
229            .zip(values)
230            .for_each(|(entry_i, vector_i)| *entry_i -= vector_i)
231    }
232    fn decrement_from_retained(&mut self, retained: &[bool], other: &Vector) {
233        self.iter_mut()
234            .flat_map(|x| x.iter_mut())
235            .zip(retained.iter())
236            .filter(|(_, retained_i)| **retained_i)
237            .zip(other.iter())
238            .for_each(|((self_i, _), vector_i)| *self_i -= Quantity::new(*vector_i))
239    }
240}
241
242impl<const D: usize, I, U> Sub<Vector> for TensorRank1Vec<D, I, U> {
243    type Output = Self;
244    fn sub(mut self, vector: Vector) -> Self::Output {
245        self.iter_mut().enumerate().for_each(|(a, self_a)| {
246            self_a
247                .iter_mut()
248                .enumerate()
249                .for_each(|(i, self_a_i)| *self_a_i -= Quantity::new(vector[D * a + i]))
250        });
251        self
252    }
253}
254
255impl<const D: usize, I, U> Sub<&Vector> for TensorRank1Vec<D, I, U> {
256    type Output = Self;
257    fn sub(mut self, vector: &Vector) -> Self::Output {
258        self.iter_mut().enumerate().for_each(|(a, self_a)| {
259            self_a
260                .iter_mut()
261                .enumerate()
262                .for_each(|(i, self_a_i)| *self_a_i -= Quantity::new(vector[D * a + i]))
263        });
264        self
265    }
266}
267
268impl<const D: usize, I, J, U> Div<TensorRank2Vec2D<D, I, J, U>> for &TensorRank1Vec<D, I, U> {
269    type Output = TensorRank1Vec<D, J, U>;
270    fn div(self, _tensor_rank_2_vec_2d: TensorRank2Vec2D<D, I, J, U>) -> Self::Output {
271        unimplemented!(
272            "A mesh-scale step wants the sparse solver the caller supplies, which a division has nowhere to hold."
273        )
274    }
275}
276
277impl<const D: usize, I, J, U, V> Div<TensorRank2SparseVec2D<D, I, J, V>>
278    for &TensorRank1Vec<D, I, U>
279where
280    U: UnitDiv<V>,
281{
282    type Output = TensorRank1Vec<D, J, <U as UnitDiv<V>>::Output>;
283    fn div(self, _tensor_rank_2_sparse_vec_2d: TensorRank2SparseVec2D<D, I, J, V>) -> Self::Output {
284        unimplemented!(
285            "A mesh-scale step wants the sparse solver the caller supplies, which a division has nowhere to hold."
286        )
287    }
288}
289
290impl<const D: usize, I, J, U, V> Div<TensorRank2SparseVec2DSymmetric<D, I, J, V>>
291    for &TensorRank1Vec<D, I, U>
292where
293    U: UnitDiv<V>,
294{
295    type Output = TensorRank1Vec<D, J, <U as UnitDiv<V>>::Output>;
296    fn div(
297        self,
298        _tensor_rank_2_sparse_symmetric_vec_2d: TensorRank2SparseVec2DSymmetric<D, I, J, V>,
299    ) -> Self::Output {
300        unimplemented!(
301            "A mesh-scale step wants the sparse solver the caller supplies, which a division has nowhere to hold."
302        )
303    }
304}
305
306impl<const D: usize, I, U> FiniteDifference for TensorRank1Vec<D, I, U> {
307    fn error_fd(&self, comparator: &Self, epsilon: TensorRank0) -> Option<(bool, usize)> {
308        let error_count = self
309            .iter()
310            .zip(comparator.iter())
311            .map(|(entry, comparator_entry)| {
312                entry
313                    .iter()
314                    .zip(comparator_entry.iter())
315                    .filter(|&(&entry_i, &comparator_entry_i)| {
316                        entry_i.differs(comparator_entry_i, epsilon)
317                    })
318                    .count()
319            })
320            .sum();
321        if error_count > 0 {
322            let auxiliary = self
323                .iter()
324                .zip(comparator.iter())
325                .map(|(entry, comparator_entry)| {
326                    entry
327                        .iter()
328                        .zip(comparator_entry.iter())
329                        .filter(|&(&entry_i, &comparator_entry_i)| {
330                            entry_i.differs_severely(comparator_entry_i, epsilon)
331                        })
332                        .count()
333                })
334                .sum::<usize>()
335                > 0;
336            Some((auxiliary, error_count))
337        } else {
338            None
339        }
340    }
341}