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

outram_foam_basic_lib/fields/
vol_field.rs
