outram_foam_basic_lib/fields/
field.rs

1// SPDX-License-Identifier: GPL-3.0-only
2// Copyright (C) 2026 OUTRAM PARK contributors
3// Derived from OpenFOAM (www.openfoam.com)
4// Copyright (C) 2004-2023 OpenFOAM Foundation
5// Copyright (C) 2016-2023 OpenCFD Ltd.
6//
7// This file is part of OUTRAM PARK.
8//
9// OUTRAM PARK is free software: you can redistribute it and/or modify it
10// under the terms of the GNU General Public License as published by the
11// Free Software Foundation, either version 3 of the License, or (at your
12// option) any later version.
13//
14// OUTRAM PARK is distributed in the hope that it will be useful, but
15// WITHOUT ANY WARRANTY; without even the implied warranty of
16// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
17// General Public License for more details.
18//
19// You should have received a copy of the GNU General Public License along
20// with OUTRAM PARK.  If not, see <https://www.gnu.org/licenses/>.
21
22use std::ops::{Add, AddAssign, Div, Mul, Neg, Sub, SubAssign};
23
24use crate::primitives::{SphericalTensor, SymmTensor, Tensor, Vector3};
25
26/// A flat array over all cells or faces, with element-wise arithmetic.
27///
28/// Mirrors `Foam::Field<Type>` from `src/OpenFOAM/fields/Fields/Field/Field.H`.
29/// The raw storage is `Vec<T>` with no dimension bookkeeping — that lives in
30/// the wrapping `VolField`/`SurfaceField`.
31#[derive(Debug, Clone, PartialEq)]
32pub struct Field<T> {
33    data: Vec<T>,
34}
35
36// ── Construction ─────────────────────────────────────────────────────────────
37
38impl<T: Clone> Field<T> {
39    /// Wrap an existing `Vec<T>` as a field (no copy).
40    pub fn new(data: Vec<T>) -> Self {
41        Self { data }
42    }
43
44    /// Field of `n` elements all equal to `value`.
45    pub fn uniform(n: usize, value: T) -> Self {
46        Self {
47            data: vec![value; n],
48        }
49    }
50
51    /// Field of `n` elements, element `i` set to `f(i)`.
52    pub fn from_fn(n: usize, f: impl Fn(usize) -> T) -> Self {
53        Self {
54            data: (0..n).map(f).collect(),
55        }
56    }
57
58    /// Number of elements in the field.
59    pub fn len(&self) -> usize {
60        self.data.len()
61    }
62
63    /// True if the field has no elements.
64    pub fn is_empty(&self) -> bool {
65        self.data.is_empty()
66    }
67
68    /// Borrow the underlying storage as a slice.
69    pub fn as_slice(&self) -> &[T] {
70        &self.data
71    }
72
73    /// Mutably borrow the underlying storage as a slice.
74    pub fn as_mut_slice(&mut self) -> &mut [T] {
75        &mut self.data
76    }
77
78    /// Consume the field, returning its underlying `Vec<T>`.
79    pub fn into_vec(self) -> Vec<T> {
80        self.data
81    }
82
83    /// Iterator over element references.
84    pub fn iter(&self) -> std::slice::Iter<'_, T> {
85        self.data.iter()
86    }
87
88    /// Iterator over mutable element references.
89    pub fn iter_mut(&mut self) -> std::slice::IterMut<'_, T> {
90        self.data.iter_mut()
91    }
92
93    /// Map `f` element-wise, returning a new `Field<U>` of the same length.
94    pub fn map<U: Clone>(&self, f: impl Fn(&T) -> U) -> Field<U> {
95        Field {
96            data: self.data.iter().map(f).collect(),
97        }
98    }
99}
100
101impl Field<f64> {
102    /// Scalar field of `n` zeros.
103    pub fn zeros(n: usize) -> Self {
104        Self::uniform(n, 0.0)
105    }
106
107    /// Scalar field of `n` ones.
108    pub fn ones(n: usize) -> Self {
109        Self::uniform(n, 1.0)
110    }
111
112    /// Sum of all elements.
113    pub fn sum(&self) -> f64 {
114        self.data.iter().copied().sum()
115    }
116
117    /// Arithmetic mean of all elements; returns `0.0` for an empty field.
118    pub fn mean(&self) -> f64 {
119        if self.data.is_empty() {
120            0.0
121        } else {
122            self.sum() / self.data.len() as f64
123        }
124    }
125
126    /// Smallest element (`+∞` for an empty field).
127    pub fn min(&self) -> f64 {
128        self.data.iter().cloned().fold(f64::INFINITY, f64::min)
129    }
130
131    /// Largest element (`−∞` for an empty field).
132    pub fn max(&self) -> f64 {
133        self.data.iter().cloned().fold(f64::NEG_INFINITY, f64::max)
134    }
135
136    /// Euclidean (L2) norm: `sqrt(sum(x_i²))`.
137    pub fn l2_norm(&self) -> f64 {
138        self.data.iter().map(|x| x * x).sum::<f64>().sqrt()
139    }
140
141    /// Element-wise absolute value.
142    pub fn abs(&self) -> Self {
143        self.map(|x| x.abs())
144    }
145
146    /// Element-wise clamp.
147    pub fn clamp(&self, lo: f64, hi: f64) -> Self {
148        self.map(|x| x.clamp(lo, hi))
149    }
150
151    /// Element-wise product of two scalar fields.
152    pub fn pointwise_mul(&self, rhs: &Self) -> Self {
153        assert_eq!(
154            self.len(),
155            rhs.len(),
156            "Field length mismatch in pointwise_mul"
157        );
158        Field::from_fn(self.len(), |i| self.data[i] * rhs.data[i])
159    }
160
161    /// Element-wise division of two scalar fields.
162    pub fn pointwise_div(&self, rhs: &Self) -> Self {
163        assert_eq!(
164            self.len(),
165            rhs.len(),
166            "Field length mismatch in pointwise_div"
167        );
168        Field::from_fn(self.len(), |i| self.data[i] / rhs.data[i])
169    }
170
171    /// Weighted sum: sum(w[i] * x[i]).
172    pub fn weighted_sum(&self, weights: &Field<f64>) -> f64 {
173        assert_eq!(self.len(), weights.len());
174        self.data
175            .iter()
176            .zip(weights.data.iter())
177            .map(|(x, w)| x * w)
178            .sum()
179    }
180}
181
182impl Field<Vector3> {
183    /// Vector field of `n` zero vectors.
184    pub fn zero_vec(n: usize) -> Self {
185        Self::uniform(n, Vector3::ZERO)
186    }
187
188    /// Element-wise dot product → scalar field.
189    pub fn dot_field(&self, rhs: &Field<Vector3>) -> Field<f64> {
190        assert_eq!(self.len(), rhs.len());
191        Field::from_fn(self.len(), |i| self.data[i].dot(rhs.data[i]))
192    }
193
194    /// Scale each element by the corresponding scalar field entry.
195    pub fn scale(&self, s: &Field<f64>) -> Self {
196        assert_eq!(self.len(), s.len());
197        Field::from_fn(self.len(), |i| self.data[i] * s[i])
198    }
199}
200
201// ── Index ─────────────────────────────────────────────────────────────────────
202
203impl<T> std::ops::Index<usize> for Field<T> {
204    type Output = T;
205    fn index(&self, i: usize) -> &T {
206        &self.data[i]
207    }
208}
209
210impl<T> std::ops::IndexMut<usize> for Field<T> {
211    fn index_mut(&mut self, i: usize) -> &mut T {
212        &mut self.data[i]
213    }
214}
215
216impl<T> AsRef<[T]> for Field<T> {
217    fn as_ref(&self) -> &[T] {
218        &self.data
219    }
220}
221
222impl<T> IntoIterator for Field<T> {
223    type Item = T;
224    type IntoIter = std::vec::IntoIter<T>;
225    fn into_iter(self) -> Self::IntoIter {
226        self.data.into_iter()
227    }
228}
229
230impl<'a, T> IntoIterator for &'a Field<T> {
231    type Item = &'a T;
232    type IntoIter = std::slice::Iter<'a, T>;
233    fn into_iter(self) -> Self::IntoIter {
234        self.data.iter()
235    }
236}
237
238impl<T: Clone> From<Vec<T>> for Field<T> {
239    fn from(v: Vec<T>) -> Self {
240        Self::new(v)
241    }
242}
243
244// ── Arithmetic — Field<T> op Field<T> ────────────────────────────────────────
245
246impl<T: Add<Output = T> + Clone> Add for Field<T> {
247    type Output = Self;
248    fn add(self, rhs: Self) -> Self {
249        assert_eq!(self.len(), rhs.len(), "Field length mismatch in add");
250        Field::from_fn(self.len(), |i| self.data[i].clone() + rhs.data[i].clone())
251    }
252}
253
254impl<T: Sub<Output = T> + Clone> Sub for Field<T> {
255    type Output = Self;
256    fn sub(self, rhs: Self) -> Self {
257        assert_eq!(self.len(), rhs.len(), "Field length mismatch in sub");
258        Field::from_fn(self.len(), |i| self.data[i].clone() - rhs.data[i].clone())
259    }
260}
261
262impl<T: Neg<Output = T> + Clone> Neg for Field<T> {
263    type Output = Self;
264    fn neg(self) -> Self {
265        self.map(|x| -x.clone())
266    }
267}
268
269impl<T: Add<Output = T> + Clone> AddAssign for Field<T> {
270    fn add_assign(&mut self, rhs: Self) {
271        assert_eq!(self.len(), rhs.len());
272        for (a, b) in self.data.iter_mut().zip(rhs.data) {
273            *a = a.clone() + b;
274        }
275    }
276}
277
278impl<T: Sub<Output = T> + Clone> SubAssign for Field<T> {
279    fn sub_assign(&mut self, rhs: Self) {
280        assert_eq!(self.len(), rhs.len());
281        for (a, b) in self.data.iter_mut().zip(rhs.data) {
282            *a = a.clone() - b;
283        }
284    }
285}
286
287// ── Arithmetic — Field<T> op f64 ─────────────────────────────────────────────
288
289impl<T: Mul<f64, Output = T> + Clone> Mul<f64> for Field<T> {
290    type Output = Self;
291    fn mul(self, rhs: f64) -> Self {
292        self.map(|x| x.clone() * rhs)
293    }
294}
295
296// f64 * Field<T>  (commutativity helper)
297impl Mul<Field<f64>> for f64 {
298    type Output = Field<f64>;
299    fn mul(self, rhs: Field<f64>) -> Field<f64> {
300        rhs.map(|x| self * x)
301    }
302}
303
304impl Mul<Field<Vector3>> for f64 {
305    type Output = Field<Vector3>;
306    fn mul(self, rhs: Field<Vector3>) -> Field<Vector3> {
307        rhs.map(|v| v.clone() * self)
308    }
309}
310
311impl<T: Mul<f64, Output = T> + Clone> Div<f64> for Field<T> {
312    type Output = Self;
313    fn div(self, rhs: f64) -> Self {
314        let inv = 1.0 / rhs;
315        self.map(|x| x.clone() * inv)
316    }
317}
318
319// ── Field<f64> specific binary ops ───────────────────────────────────────────
320
321impl Mul for Field<f64> {
322    type Output = Self;
323    fn mul(self, rhs: Self) -> Self {
324        self.pointwise_mul(&rhs)
325    }
326}
327
328impl Div for Field<f64> {
329    type Output = Self;
330    fn div(self, rhs: Self) -> Self {
331        self.pointwise_div(&rhs)
332    }
333}
334
335// Field<Vector3> scaled by Field<f64>
336impl Mul<Field<f64>> for Field<Vector3> {
337    type Output = Field<Vector3>;
338    fn mul(self, rhs: Field<f64>) -> Field<Vector3> {
339        self.scale(&rhs)
340    }
341}
342
343impl Mul<Field<Vector3>> for Field<f64> {
344    type Output = Field<Vector3>;
345    fn mul(self, rhs: Field<Vector3>) -> Field<Vector3> {
346        rhs.scale(&self)
347    }
348}
349
350// ── Default ───────────────────────────────────────────────────────────────────
351
352impl Default for Field<f64> {
353    fn default() -> Self {
354        Self::new(vec![])
355    }
356}
357
358impl Default for Field<Vector3> {
359    fn default() -> Self {
360        Self::new(vec![])
361    }
362}
363
364impl Default for Field<Tensor> {
365    fn default() -> Self {
366        Self::new(vec![])
367    }
368}
369
370impl Default for Field<SymmTensor> {
371    fn default() -> Self {
372        Self::new(vec![])
373    }
374}
375
376impl Default for Field<SphericalTensor> {
377    fn default() -> Self {
378        Self::new(vec![])
379    }
380}
381
382// ── Tests ─────────────────────────────────────────────────────────────────────
383
384#[cfg(test)]
385mod tests {
386    use super::*;
387
388    #[test]
389    fn add_f64_fields() {
390        let a = Field::new(vec![1.0, 2.0, 3.0]);
391        let b = Field::new(vec![4.0, 5.0, 6.0]);
392        let c = a + b;
393        assert_eq!(c.as_slice(), &[5.0, 7.0, 9.0]);
394    }
395
396    #[test]
397    fn mul_scalar() {
398        let a = Field::new(vec![1.0, 2.0, 3.0]);
399        let b = a * 2.0;
400        assert_eq!(b.as_slice(), &[2.0, 4.0, 6.0]);
401    }
402
403    #[test]
404    fn div_scalar() {
405        let a = Field::new(vec![2.0, 4.0, 6.0]);
406        let b = a / 2.0;
407        assert_eq!(b.as_slice(), &[1.0, 2.0, 3.0]);
408    }
409
410    #[test]
411    fn neg_field() {
412        let a = Field::new(vec![1.0, -2.0, 3.0]);
413        let b = -a;
414        assert_eq!(b.as_slice(), &[-1.0, 2.0, -3.0]);
415    }
416
417    #[test]
418    fn sum_mean() {
419        let a = Field::new(vec![1.0, 2.0, 3.0, 4.0]);
420        assert_eq!(a.sum(), 10.0);
421        assert_eq!(a.mean(), 2.5);
422    }
423
424    #[test]
425    fn pointwise_mul_fields() {
426        let a = Field::new(vec![2.0, 3.0]);
427        let b = Field::new(vec![4.0, 5.0]);
428        let c = a.pointwise_mul(&b);
429        assert_eq!(c.as_slice(), &[8.0, 15.0]);
430    }
431
432    #[test]
433    fn vector3_add_fields() {
434        let a = Field::new(vec![
435            Vector3::new(1.0, 0.0, 0.0),
436            Vector3::new(0.0, 1.0, 0.0),
437        ]);
438        let b = Field::new(vec![
439            Vector3::new(0.0, 0.0, 1.0),
440            Vector3::new(1.0, 0.0, 0.0),
441        ]);
442        let c = a + b;
443        assert_eq!(c[0], Vector3::new(1.0, 0.0, 1.0));
444        assert_eq!(c[1], Vector3::new(1.0, 1.0, 0.0));
445    }
446
447    #[test]
448    fn vector3_scale_by_scalar_field() {
449        let v = Field::new(vec![
450            Vector3::new(1.0, 2.0, 3.0),
451            Vector3::new(4.0, 5.0, 6.0),
452        ]);
453        let s = Field::new(vec![2.0, 0.5]);
454        let result = v.scale(&s);
455        assert_eq!(result[0], Vector3::new(2.0, 4.0, 6.0));
456        assert_eq!(result[1], Vector3::new(2.0, 2.5, 3.0));
457    }
458
459    #[test]
460    fn uniform_and_len() {
461        let f: Field<f64> = Field::uniform(5, 3.14);
462        assert_eq!(f.len(), 5);
463        assert!((f[0] - 3.14).abs() < 1e-15);
464    }
465}
outram_foam_basic_lib/fields/field.rs

outram_foam_basic_lib/fields/
field.rs
