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    pub fn new(data: Vec<T>) -> Self {
40        Self { data }
41    }
42
43    pub fn uniform(n: usize, value: T) -> Self {
44        Self {
45            data: vec![value; n],
46        }
47    }
48
49    pub fn from_fn(n: usize, f: impl Fn(usize) -> T) -> Self {
50        Self {
51            data: (0..n).map(f).collect(),
52        }
53    }
54
55    pub fn len(&self) -> usize {
56        self.data.len()
57    }
58
59    pub fn is_empty(&self) -> bool {
60        self.data.is_empty()
61    }
62
63    pub fn as_slice(&self) -> &[T] {
64        &self.data
65    }
66
67    pub fn as_mut_slice(&mut self) -> &mut [T] {
68        &mut self.data
69    }
70
71    pub fn into_vec(self) -> Vec<T> {
72        self.data
73    }
74
75    pub fn iter(&self) -> std::slice::Iter<'_, T> {
76        self.data.iter()
77    }
78
79    pub fn iter_mut(&mut self) -> std::slice::IterMut<'_, T> {
80        self.data.iter_mut()
81    }
82
83    pub fn map<U: Clone>(&self, f: impl Fn(&T) -> U) -> Field<U> {
84        Field {
85            data: self.data.iter().map(f).collect(),
86        }
87    }
88}
89
90impl Field<f64> {
91    pub fn zeros(n: usize) -> Self {
92        Self::uniform(n, 0.0)
93    }
94
95    pub fn ones(n: usize) -> Self {
96        Self::uniform(n, 1.0)
97    }
98
99    pub fn sum(&self) -> f64 {
100        self.data.iter().copied().sum()
101    }
102
103    pub fn mean(&self) -> f64 {
104        if self.data.is_empty() {
105            0.0
106        } else {
107            self.sum() / self.data.len() as f64
108        }
109    }
110
111    pub fn min(&self) -> f64 {
112        self.data.iter().cloned().fold(f64::INFINITY, f64::min)
113    }
114
115    pub fn max(&self) -> f64 {
116        self.data.iter().cloned().fold(f64::NEG_INFINITY, f64::max)
117    }
118
119    pub fn l2_norm(&self) -> f64 {
120        self.data.iter().map(|x| x * x).sum::<f64>().sqrt()
121    }
122
123    /// Element-wise absolute value.
124    pub fn abs(&self) -> Self {
125        self.map(|x| x.abs())
126    }
127
128    /// Element-wise clamp.
129    pub fn clamp(&self, lo: f64, hi: f64) -> Self {
130        self.map(|x| x.clamp(lo, hi))
131    }
132
133    /// Element-wise product of two scalar fields.
134    pub fn pointwise_mul(&self, rhs: &Self) -> Self {
135        assert_eq!(
136            self.len(),
137            rhs.len(),
138            "Field length mismatch in pointwise_mul"
139        );
140        Field::from_fn(self.len(), |i| self.data[i] * rhs.data[i])
141    }
142
143    /// Element-wise division of two scalar fields.
144    pub fn pointwise_div(&self, rhs: &Self) -> Self {
145        assert_eq!(
146            self.len(),
147            rhs.len(),
148            "Field length mismatch in pointwise_div"
149        );
150        Field::from_fn(self.len(), |i| self.data[i] / rhs.data[i])
151    }
152
153    /// Weighted sum: sum(w[i] * x[i]).
154    pub fn weighted_sum(&self, weights: &Field<f64>) -> f64 {
155        assert_eq!(self.len(), weights.len());
156        self.data
157            .iter()
158            .zip(weights.data.iter())
159            .map(|(x, w)| x * w)
160            .sum()
161    }
162}
163
164impl Field<Vector3> {
165    pub fn zero_vec(n: usize) -> Self {
166        Self::uniform(n, Vector3::ZERO)
167    }
168
169    /// Element-wise dot product → scalar field.
170    pub fn dot_field(&self, rhs: &Field<Vector3>) -> Field<f64> {
171        assert_eq!(self.len(), rhs.len());
172        Field::from_fn(self.len(), |i| self.data[i].dot(rhs.data[i]))
173    }
174
175    /// Scale each element by the corresponding scalar field entry.
176    pub fn scale(&self, s: &Field<f64>) -> Self {
177        assert_eq!(self.len(), s.len());
178        Field::from_fn(self.len(), |i| self.data[i] * s[i])
179    }
180}
181
182// ── Index ─────────────────────────────────────────────────────────────────────
183
184impl<T> std::ops::Index<usize> for Field<T> {
185    type Output = T;
186    fn index(&self, i: usize) -> &T {
187        &self.data[i]
188    }
189}
190
191impl<T> std::ops::IndexMut<usize> for Field<T> {
192    fn index_mut(&mut self, i: usize) -> &mut T {
193        &mut self.data[i]
194    }
195}
196
197impl<T> AsRef<[T]> for Field<T> {
198    fn as_ref(&self) -> &[T] {
199        &self.data
200    }
201}
202
203impl<T> IntoIterator for Field<T> {
204    type Item = T;
205    type IntoIter = std::vec::IntoIter<T>;
206    fn into_iter(self) -> Self::IntoIter {
207        self.data.into_iter()
208    }
209}
210
211impl<'a, T> IntoIterator for &'a Field<T> {
212    type Item = &'a T;
213    type IntoIter = std::slice::Iter<'a, T>;
214    fn into_iter(self) -> Self::IntoIter {
215        self.data.iter()
216    }
217}
218
219impl<T: Clone> From<Vec<T>> for Field<T> {
220    fn from(v: Vec<T>) -> Self {
221        Self::new(v)
222    }
223}
224
225// ── Arithmetic — Field<T> op Field<T> ────────────────────────────────────────
226
227impl<T: Add<Output = T> + Clone> Add for Field<T> {
228    type Output = Self;
229    fn add(self, rhs: Self) -> Self {
230        assert_eq!(self.len(), rhs.len(), "Field length mismatch in add");
231        Field::from_fn(self.len(), |i| self.data[i].clone() + rhs.data[i].clone())
232    }
233}
234
235impl<T: Sub<Output = T> + Clone> Sub for Field<T> {
236    type Output = Self;
237    fn sub(self, rhs: Self) -> Self {
238        assert_eq!(self.len(), rhs.len(), "Field length mismatch in sub");
239        Field::from_fn(self.len(), |i| self.data[i].clone() - rhs.data[i].clone())
240    }
241}
242
243impl<T: Neg<Output = T> + Clone> Neg for Field<T> {
244    type Output = Self;
245    fn neg(self) -> Self {
246        self.map(|x| -x.clone())
247    }
248}
249
250impl<T: Add<Output = T> + Clone> AddAssign for Field<T> {
251    fn add_assign(&mut self, rhs: Self) {
252        assert_eq!(self.len(), rhs.len());
253        for (a, b) in self.data.iter_mut().zip(rhs.data) {
254            *a = a.clone() + b;
255        }
256    }
257}
258
259impl<T: Sub<Output = T> + Clone> SubAssign for Field<T> {
260    fn sub_assign(&mut self, rhs: Self) {
261        assert_eq!(self.len(), rhs.len());
262        for (a, b) in self.data.iter_mut().zip(rhs.data) {
263            *a = a.clone() - b;
264        }
265    }
266}
267
268// ── Arithmetic — Field<T> op f64 ─────────────────────────────────────────────
269
270impl<T: Mul<f64, Output = T> + Clone> Mul<f64> for Field<T> {
271    type Output = Self;
272    fn mul(self, rhs: f64) -> Self {
273        self.map(|x| x.clone() * rhs)
274    }
275}
276
277// f64 * Field<T>  (commutativity helper)
278impl Mul<Field<f64>> for f64 {
279    type Output = Field<f64>;
280    fn mul(self, rhs: Field<f64>) -> Field<f64> {
281        rhs.map(|x| self * x)
282    }
283}
284
285impl Mul<Field<Vector3>> for f64 {
286    type Output = Field<Vector3>;
287    fn mul(self, rhs: Field<Vector3>) -> Field<Vector3> {
288        rhs.map(|v| v.clone() * self)
289    }
290}
291
292impl<T: Mul<f64, Output = T> + Clone> Div<f64> for Field<T> {
293    type Output = Self;
294    fn div(self, rhs: f64) -> Self {
295        let inv = 1.0 / rhs;
296        self.map(|x| x.clone() * inv)
297    }
298}
299
300// ── Field<f64> specific binary ops ───────────────────────────────────────────
301
302impl Mul for Field<f64> {
303    type Output = Self;
304    fn mul(self, rhs: Self) -> Self {
305        self.pointwise_mul(&rhs)
306    }
307}
308
309impl Div for Field<f64> {
310    type Output = Self;
311    fn div(self, rhs: Self) -> Self {
312        self.pointwise_div(&rhs)
313    }
314}
315
316// Field<Vector3> scaled by Field<f64>
317impl Mul<Field<f64>> for Field<Vector3> {
318    type Output = Field<Vector3>;
319    fn mul(self, rhs: Field<f64>) -> Field<Vector3> {
320        self.scale(&rhs)
321    }
322}
323
324impl Mul<Field<Vector3>> for Field<f64> {
325    type Output = Field<Vector3>;
326    fn mul(self, rhs: Field<Vector3>) -> Field<Vector3> {
327        rhs.scale(&self)
328    }
329}
330
331// ── Default ───────────────────────────────────────────────────────────────────
332
333impl Default for Field<f64> {
334    fn default() -> Self {
335        Self::new(vec![])
336    }
337}
338
339impl Default for Field<Vector3> {
340    fn default() -> Self {
341        Self::new(vec![])
342    }
343}
344
345impl Default for Field<Tensor> {
346    fn default() -> Self {
347        Self::new(vec![])
348    }
349}
350
351impl Default for Field<SymmTensor> {
352    fn default() -> Self {
353        Self::new(vec![])
354    }
355}
356
357impl Default for Field<SphericalTensor> {
358    fn default() -> Self {
359        Self::new(vec![])
360    }
361}
362
363// ── Tests ─────────────────────────────────────────────────────────────────────
364
365#[cfg(test)]
366mod tests {
367    use super::*;
368
369    #[test]
370    fn add_f64_fields() {
371        let a = Field::new(vec![1.0, 2.0, 3.0]);
372        let b = Field::new(vec![4.0, 5.0, 6.0]);
373        let c = a + b;
374        assert_eq!(c.as_slice(), &[5.0, 7.0, 9.0]);
375    }
376
377    #[test]
378    fn mul_scalar() {
379        let a = Field::new(vec![1.0, 2.0, 3.0]);
380        let b = a * 2.0;
381        assert_eq!(b.as_slice(), &[2.0, 4.0, 6.0]);
382    }
383
384    #[test]
385    fn div_scalar() {
386        let a = Field::new(vec![2.0, 4.0, 6.0]);
387        let b = a / 2.0;
388        assert_eq!(b.as_slice(), &[1.0, 2.0, 3.0]);
389    }
390
391    #[test]
392    fn neg_field() {
393        let a = Field::new(vec![1.0, -2.0, 3.0]);
394        let b = -a;
395        assert_eq!(b.as_slice(), &[-1.0, 2.0, -3.0]);
396    }
397
398    #[test]
399    fn sum_mean() {
400        let a = Field::new(vec![1.0, 2.0, 3.0, 4.0]);
401        assert_eq!(a.sum(), 10.0);
402        assert_eq!(a.mean(), 2.5);
403    }
404
405    #[test]
406    fn pointwise_mul_fields() {
407        let a = Field::new(vec![2.0, 3.0]);
408        let b = Field::new(vec![4.0, 5.0]);
409        let c = a.pointwise_mul(&b);
410        assert_eq!(c.as_slice(), &[8.0, 15.0]);
411    }
412
413    #[test]
414    fn vector3_add_fields() {
415        let a = Field::new(vec![
416            Vector3::new(1.0, 0.0, 0.0),
417            Vector3::new(0.0, 1.0, 0.0),
418        ]);
419        let b = Field::new(vec![
420            Vector3::new(0.0, 0.0, 1.0),
421            Vector3::new(1.0, 0.0, 0.0),
422        ]);
423        let c = a + b;
424        assert_eq!(c[0], Vector3::new(1.0, 0.0, 1.0));
425        assert_eq!(c[1], Vector3::new(1.0, 1.0, 0.0));
426    }
427
428    #[test]
429    fn vector3_scale_by_scalar_field() {
430        let v = Field::new(vec![
431            Vector3::new(1.0, 2.0, 3.0),
432            Vector3::new(4.0, 5.0, 6.0),
433        ]);
434        let s = Field::new(vec![2.0, 0.5]);
435        let result = v.scale(&s);
436        assert_eq!(result[0], Vector3::new(2.0, 4.0, 6.0));
437        assert_eq!(result[1], Vector3::new(2.0, 2.5, 3.0));
438    }
439
440    #[test]
441    fn uniform_and_len() {
442        let f: Field<f64> = Field::uniform(5, 3.14);
443        assert_eq!(f.len(), 5);
444        assert!((f[0] - 3.14).abs() < 1e-15);
445    }
446}
outram_foam_basic_lib/fields/field.rs

outram_foam_basic_lib/fields/
field.rs
