outram_foam_basic_lib/fields/
vol_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, Div, Mul, Neg, Sub};
23use std::sync::Arc;
24
25use super::boundary::bc::PatchField;
26use super::field::Field;
27use crate::mesh::fv_mesh::FvMesh;
28use crate::primitives::{SymmTensor, Tensor, Vector3};
29
30/// A volume field: one value per cell in the internal field, plus one
31/// `PatchField` per boundary patch.
32///
33/// Mirrors `Foam::volScalarField` / `Foam::VolumeField<Type>`.
34/// The internal field has length `mesh.n_cells`.
35#[derive(Debug, Clone)]
36pub struct VolField<T: Clone> {
37    /// Field name (diagnostic label, e.g. `"p"`, `"U"`, `"T"`).
38    pub name: String,
39    /// Mesh this field is defined on.
40    pub mesh: Arc<FvMesh>,
41    /// Cell-centred values; length == `mesh.n_cells`.
42    pub internal: Field<T>,
43    /// One entry per boundary patch; `boundary[i].values` has length
44    /// `mesh.patches[i].size`.
45    pub boundary: Vec<PatchField<T>>,
46}
47
48// ── Type aliases (matching OpenFOAM names) ────────────────────────────────────
49
50/// Scalar volume field: one `f64` per cell (e.g. pressure `[Pa]`, temperature `[K]`).
51pub type VolScalarField = VolField<f64>;
52/// Vector volume field: one `Vector3` per cell (e.g. velocity `[m/s]`).
53pub type VolVectorField = VolField<Vector3>;
54/// General (rank-2) tensor volume field: one `Tensor` per cell.
55pub type VolTensorField = VolField<Tensor>;
56/// Symmetric (rank-2) tensor volume field: one `SymmTensor` per cell (e.g. stress).
57pub type VolSymmTensorField = VolField<SymmTensor>;
58
59// ── Construction ─────────────────────────────────────────────────────────────
60
61impl<T: Clone + Default> VolField<T> {
62    /// Assemble a volume field from its cell-centred values and per-patch
63    /// boundary fields. In debug builds, asserts that `internal` has length
64    /// `mesh.n_cells` and that `boundary` has one entry per patch.
65    pub fn new(
66        name: impl Into<String>,
67        mesh: Arc<FvMesh>,
68        internal: Field<T>,
69        boundary: Vec<PatchField<T>>,
70    ) -> Self {
71        debug_assert_eq!(
72            internal.len(),
73            mesh.n_cells,
74            "VolField internal field length must equal n_cells"
75        );
76        debug_assert_eq!(
77            boundary.len(),
78            mesh.patches.len(),
79            "VolField boundary length must equal number of patches"
80        );
81        Self {
82            name: name.into(),
83            mesh,
84            internal,
85            boundary,
86        }
87    }
88}
89
90impl VolScalarField {
91    /// Uniform scalar field over the entire domain.
92    pub fn uniform(name: impl Into<String>, mesh: Arc<FvMesh>, value: f64) -> Self {
93        let n_cells = mesh.n_cells;
94        let boundary = mesh
95            .patches
96            .iter()
97            .map(|p| PatchField::zero_gradient(p.size))
98            .collect();
99        Self::new(name, mesh, Field::uniform(n_cells, value), boundary)
100    }
101
102    /// Zero-valued scalar field over the entire domain (zero-gradient boundaries).
103    pub fn zeros(name: impl Into<String>, mesh: Arc<FvMesh>) -> Self {
104        Self::uniform(name, mesh, 0.0)
105    }
106}
107
108impl VolVectorField {
109    /// Uniform vector field over the entire domain.
110    pub fn uniform(name: impl Into<String>, mesh: Arc<FvMesh>, value: Vector3) -> Self {
111        let n_cells = mesh.n_cells;
112        let boundary = mesh
113            .patches
114            .iter()
115            .map(|p| PatchField::zero_gradient_vec(p.size))
116            .collect();
117        Self::new(name, mesh, Field::uniform(n_cells, value), boundary)
118    }
119
120    /// Zero-valued vector field over the entire domain (zero-gradient boundaries).
121    pub fn zero(name: impl Into<String>, mesh: Arc<FvMesh>) -> Self {
122        Self::uniform(name, mesh, Vector3::ZERO)
123    }
124}
125
126// ── Arithmetic — borrows rhs so we don't move the mesh Arc ───────────────────
127
128impl<T> Add for VolField<T>
129where
130    T: Add<Output = T> + Clone + Default,
131{
132    type Output = VolField<T>;
133    fn add(mut self, rhs: Self) -> Self::Output {
134        // NB: deliberately keep `self.name` unchanged. Rebuilding a compound
135        // name here (e.g. "(rho + div(phi))") makes the `name` String grow —
136        // and in a solver loop where a field is repeatedly reassigned from an
137        // expression containing itself, it grows *exponentially* (doubling per
138        // step), exhausting memory within ~25 steps. The name is only a
139        // diagnostic label; the left operand's name is the useful one.
140        self.internal = self.internal + rhs.internal;
141        for (l, r) in self.boundary.iter_mut().zip(rhs.boundary) {
142            l.values = l.values.clone() + r.values;
143        }
144        self
145    }
146}
147
148impl<T> Sub for VolField<T>
149where
150    T: Sub<Output = T> + Clone + Default,
151{
152    type Output = VolField<T>;
153    fn sub(mut self, rhs: Self) -> Self::Output {
154        // See `Add` above: do not rebuild the name (avoids unbounded growth).
155        self.internal = self.internal - rhs.internal;
156        for (l, r) in self.boundary.iter_mut().zip(rhs.boundary) {
157            l.values = l.values.clone() - r.values;
158        }
159        self
160    }
161}
162
163impl<T> Neg for VolField<T>
164where
165    T: Neg<Output = T> + Clone,
166{
167    type Output = VolField<T>;
168    fn neg(mut self) -> Self::Output {
169        // See `Add` above: do not rebuild the name (avoids unbounded growth).
170        self.internal = -self.internal;
171        for p in self.boundary.iter_mut() {
172            p.values = -p.values.clone();
173        }
174        self
175    }
176}
177
178impl<T> Mul<f64> for VolField<T>
179where
180    T: Mul<f64, Output = T> + Clone,
181{
182    type Output = VolField<T>;
183    fn mul(mut self, s: f64) -> Self::Output {
184        self.internal = self.internal * s;
185        for p in self.boundary.iter_mut() {
186            p.values = p.values.clone() * s;
187        }
188        self
189    }
190}
191
192impl<T> Div<f64> for VolField<T>
193where
194    T: Mul<f64, Output = T> + Clone,
195{
196    type Output = VolField<T>;
197    fn div(self, s: f64) -> Self::Output {
198        self * (1.0 / s)
199    }
200}
201
202// Pointwise scalar * scalar
203impl Mul for VolScalarField {
204    type Output = Self;
205    fn mul(mut self, rhs: Self) -> Self {
206        self.internal = self.internal * rhs.internal;
207        for (l, r) in self.boundary.iter_mut().zip(rhs.boundary) {
208            l.values = l.values.clone() * r.values;
209        }
210        self
211    }
212}
213
214// Scalar field * vector field
215impl Mul<VolVectorField> for VolScalarField {
216    type Output = VolVectorField;
217    fn mul(self, mut rhs: VolVectorField) -> VolVectorField {
218        rhs.internal = rhs.internal * self.internal;
219        for (r, l) in rhs.boundary.iter_mut().zip(self.boundary) {
220            r.values = r.values.clone() * l.values;
221        }
222        rhs
223    }
224}
225
226// f64 * VolField<T>
227impl<T: Mul<f64, Output = T> + Clone> Mul<VolField<T>> for f64 {
228    type Output = VolField<T>;
229    fn mul(self, mut rhs: VolField<T>) -> VolField<T> {
230        rhs.internal = rhs.internal * self;
231        for p in rhs.boundary.iter_mut() {
232            p.values = p.values.clone() * self;
233        }
234        rhs
235    }
236}
237
238// ── Tests ─────────────────────────────────────────────────────────────────────
239
240#[cfg(test)]
241mod tests {
242    use super::*;
243    use crate::mesh::fv_mesh::{BoundaryPatch, FvMeshBuilder, PatchKind};
244
245    fn unit_mesh() -> Arc<FvMesh> {
246        Arc::new(
247            FvMeshBuilder::new()
248                .n_cells(2)
249                .n_internal_faces(1)
250                .owner(vec![0, 1, 0])
251                .neighbour(vec![1])
252                .patches(vec![
253                    BoundaryPatch::new("right", 1, 1, PatchKind::Wall),
254                    BoundaryPatch::new("left", 2, 1, PatchKind::Wall),
255                ])
256                .cell_volumes(vec![1.0, 1.0])
257                .cell_centres(vec![
258                    Vector3::new(0.25, 0.0, 0.0),
259                    Vector3::new(0.75, 0.0, 0.0),
260                ])
261                .face_area_vectors(vec![
262                    Vector3::new(1.0, 0.0, 0.0),
263                    Vector3::new(1.0, 0.0, 0.0),
264                    Vector3::new(-1.0, 0.0, 0.0),
265                ])
266                .face_centres(vec![
267                    Vector3::new(0.5, 0.0, 0.0),
268                    Vector3::new(1.0, 0.0, 0.0),
269                    Vector3::new(0.0, 0.0, 0.0),
270                ])
271                .build()
272                .unwrap(),
273        )
274    }
275
276    #[test]
277    fn uniform_vol_scalar_field() {
278        let m = unit_mesh();
279        let p = VolScalarField::uniform("p", m, 101325.0);
280        assert_eq!(p.internal.len(), 2);
281        assert!((p.internal[0] - 101325.0).abs() < 1e-8);
282        assert_eq!(p.boundary.len(), 2);
283    }
284
285    #[test]
286    fn add_vol_scalar_fields() {
287        let m = unit_mesh();
288        let a = VolScalarField::uniform("a", m.clone(), 1.0);
289        let b = VolScalarField::uniform("b", m, 2.0);
290        let c = a + b;
291        assert!((c.internal[0] - 3.0).abs() < 1e-15);
292        assert!((c.internal[1] - 3.0).abs() < 1e-15);
293    }
294
295    #[test]
296    fn mul_scalar_field_by_f64() {
297        let m = unit_mesh();
298        let a = VolScalarField::uniform("a", m, 3.0);
299        let b = a * 2.0;
300        assert!((b.internal[0] - 6.0).abs() < 1e-15);
301    }
302
303    #[test]
304    fn neg_vol_vector_field() {
305        let m = unit_mesh();
306        let u = VolVectorField::uniform("U", m, Vector3::new(1.0, 2.0, 3.0));
307        let neg_u = -u;
308        assert_eq!(neg_u.internal[0], Vector3::new(-1.0, -2.0, -3.0));
309    }
310}
outram_foam_basic_lib/fields/vol_field.rs

outram_foam_basic_lib/fields/
vol_field.rs
