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
19pub 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}