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

outram_foam_basic_lib/fields/
field.rs
