outram_foam_basic_lib/fields/
surface_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::Vector3;
29
30/// A surface field: one value per *internal* face in the internal field, plus
31/// one `PatchField` per boundary patch.
32///
33/// Mirrors `Foam::surfaceScalarField` / `Foam::SurfaceField<Type>`.
34///
35/// ## Why `internal` has length `n_internal_faces`, not `n_faces`
36///
37/// In OpenFOAM, `surfaceScalarField.internalField()` only covers the internal
38/// faces; boundary-face values live in `boundaryField()[patch]`.  This matches
39/// the LDU matrix structure: `lower` and `upper` arrays have length
40/// `n_internal_faces`.
41#[derive(Debug, Clone)]
42pub struct SurfaceField<T: Clone> {
43    /// Field name (diagnostic label, e.g. `"phi"`).
44    pub name: String,
45    /// Mesh this field is defined on.
46    pub mesh: Arc<FvMesh>,
47    /// Face values for all internal faces; length == `mesh.n_internal_faces`.
48    pub internal: Field<T>,
49    /// One entry per boundary patch; `boundary[i].values` has length
50    /// `mesh.patches[i].size`.
51    pub boundary: Vec<PatchField<T>>,
52}
53
54// ── Type aliases ──────────────────────────────────────────────────────────────
55
56/// Scalar surface field: one `f64` per face (e.g. face flux `phi` `[m³/s]`).
57pub type SurfaceScalarField = SurfaceField<f64>;
58/// Vector surface field: one `Vector3` per face.
59pub type SurfaceVectorField = SurfaceField<Vector3>;
60
61// ── Construction ─────────────────────────────────────────────────────────────
62
63impl<T: Clone> SurfaceField<T> {
64    /// Assemble a surface field from its internal-face values and per-patch
65    /// boundary fields. In debug builds, asserts that `internal` has length
66    /// `mesh.n_internal_faces` and that `boundary` has one entry per patch.
67    pub fn new(
68        name: impl Into<String>,
69        mesh: Arc<FvMesh>,
70        internal: Field<T>,
71        boundary: Vec<PatchField<T>>,
72    ) -> Self {
73        debug_assert_eq!(
74            internal.len(),
75            mesh.n_internal_faces,
76            "SurfaceField internal field length must equal n_internal_faces"
77        );
78        debug_assert_eq!(
79            boundary.len(),
80            mesh.patches.len(),
81            "SurfaceField boundary length must equal number of patches"
82        );
83        Self {
84            name: name.into(),
85            mesh,
86            internal,
87            boundary,
88        }
89    }
90}
91
92impl SurfaceScalarField {
93    /// Scalar surface field with all internal-face values zero and zero-gradient
94    /// boundary patches.
95    pub fn zeros(name: impl Into<String>, mesh: Arc<FvMesh>) -> Self {
96        let n_int = mesh.n_internal_faces;
97        let boundary = mesh
98            .patches
99            .iter()
100            .map(|p| PatchField::zero_gradient(p.size))
101            .collect();
102        Self::new(name, mesh, Field::zeros(n_int), boundary)
103    }
104
105    /// Scalar surface field with all internal-face values set to `value` and
106    /// zero-gradient boundary patches.
107    pub fn uniform(name: impl Into<String>, mesh: Arc<FvMesh>, value: f64) -> Self {
108        let n_int = mesh.n_internal_faces;
109        let boundary = mesh
110            .patches
111            .iter()
112            .map(|p| PatchField::zero_gradient(p.size))
113            .collect();
114        Self::new(name, mesh, Field::uniform(n_int, value), boundary)
115    }
116}
117
118impl SurfaceVectorField {
119    /// Vector surface field with all internal-face values zero and zero-gradient
120    /// boundary patches.
121    pub fn zero(name: impl Into<String>, mesh: Arc<FvMesh>) -> Self {
122        let n_int = mesh.n_internal_faces;
123        let boundary = mesh
124            .patches
125            .iter()
126            .map(|p| PatchField::zero_gradient_vec(p.size))
127            .collect();
128        Self::new(name, mesh, Field::zero_vec(n_int), boundary)
129    }
130}
131
132// ── Value access across all faces ─────────────────────────────────────────────
133
134impl<T: Clone + Default> SurfaceField<T> {
135    /// Value at any face: internal face → from `internal`; boundary face →
136    /// from the appropriate patch's `values`.
137    pub fn face_value(&self, f: usize) -> T {
138        if self.mesh.is_internal_face(f) {
139            self.internal[f].clone()
140        } else {
141            let (pi, fi) = self
142                .mesh
143                .patch_for_face(f)
144                .expect("face index out of range for boundary");
145            self.boundary[pi].values[fi].clone()
146        }
147    }
148}
149
150// ── Arithmetic ────────────────────────────────────────────────────────────────
151
152impl<T> Add for SurfaceField<T>
153where
154    T: Add<Output = T> + Clone + Default,
155{
156    type Output = Self;
157    fn add(mut self, rhs: Self) -> Self {
158        // NB: deliberately keep `self.name` unchanged. Rebuilding a compound
159        // name here makes the `name` String grow exponentially when a field is
160        // reassigned from an expression containing itself in a solver loop
161        // (see the matching note in `vol_field.rs`). The name is only a label.
162        self.internal = self.internal + rhs.internal;
163        for (l, r) in self.boundary.iter_mut().zip(rhs.boundary) {
164            l.values = l.values.clone() + r.values;
165        }
166        self
167    }
168}
169
170impl<T> Sub for SurfaceField<T>
171where
172    T: Sub<Output = T> + Clone + Default,
173{
174    type Output = Self;
175    fn sub(mut self, rhs: Self) -> Self {
176        // See `Add` above: do not rebuild the name (avoids unbounded growth).
177        self.internal = self.internal - rhs.internal;
178        for (l, r) in self.boundary.iter_mut().zip(rhs.boundary) {
179            l.values = l.values.clone() - r.values;
180        }
181        self
182    }
183}
184
185impl<T> Neg for SurfaceField<T>
186where
187    T: Neg<Output = T> + Clone,
188{
189    type Output = Self;
190    fn neg(mut self) -> Self {
191        // See `Add` above: do not rebuild the name (avoids unbounded growth).
192        self.internal = -self.internal;
193        for p in self.boundary.iter_mut() {
194            p.values = -p.values.clone();
195        }
196        self
197    }
198}
199
200impl<T> Mul<f64> for SurfaceField<T>
201where
202    T: Mul<f64, Output = T> + Clone,
203{
204    type Output = Self;
205    fn mul(mut self, s: f64) -> Self {
206        self.internal = self.internal * s;
207        for p in self.boundary.iter_mut() {
208            p.values = p.values.clone() * s;
209        }
210        self
211    }
212}
213
214impl<T> Div<f64> for SurfaceField<T>
215where
216    T: Mul<f64, Output = T> + Clone,
217{
218    type Output = Self;
219    fn div(self, s: f64) -> Self {
220        self * (1.0 / s)
221    }
222}
223
224// Pointwise scalar * scalar
225impl Mul for SurfaceScalarField {
226    type Output = Self;
227    fn mul(mut self, rhs: Self) -> Self {
228        self.internal = self.internal * rhs.internal;
229        for (l, r) in self.boundary.iter_mut().zip(rhs.boundary) {
230            l.values = l.values.clone() * r.values;
231        }
232        self
233    }
234}
235
236// f64 * SurfaceField<T>
237impl<T: Mul<f64, Output = T> + Clone> Mul<SurfaceField<T>> for f64 {
238    type Output = SurfaceField<T>;
239    fn mul(self, mut rhs: SurfaceField<T>) -> SurfaceField<T> {
240        rhs.internal = rhs.internal * self;
241        for p in rhs.boundary.iter_mut() {
242            p.values = p.values.clone() * self;
243        }
244        rhs
245    }
246}
247
248// ── Tests ─────────────────────────────────────────────────────────────────────
249
250#[cfg(test)]
251mod tests {
252    use super::*;
253    use crate::mesh::fv_mesh::{BoundaryPatch, FvMeshBuilder, PatchKind};
254
255    fn unit_mesh() -> Arc<FvMesh> {
256        Arc::new(
257            FvMeshBuilder::new()
258                .n_cells(2)
259                .n_internal_faces(1)
260                .owner(vec![0, 1, 0])
261                .neighbour(vec![1])
262                .patches(vec![
263                    BoundaryPatch::new("right", 1, 1, PatchKind::Wall),
264                    BoundaryPatch::new("left", 2, 1, PatchKind::Wall),
265                ])
266                .cell_volumes(vec![1.0, 1.0])
267                .cell_centres(vec![
268                    Vector3::new(0.25, 0.0, 0.0),
269                    Vector3::new(0.75, 0.0, 0.0),
270                ])
271                .face_area_vectors(vec![
272                    Vector3::new(1.0, 0.0, 0.0),
273                    Vector3::new(1.0, 0.0, 0.0),
274                    Vector3::new(-1.0, 0.0, 0.0),
275                ])
276                .face_centres(vec![
277                    Vector3::new(0.5, 0.0, 0.0),
278                    Vector3::new(1.0, 0.0, 0.0),
279                    Vector3::new(0.0, 0.0, 0.0),
280                ])
281                .build()
282                .unwrap(),
283        )
284    }
285
286    #[test]
287    fn zeros_surface_scalar_field() {
288        let m = unit_mesh();
289        let phi = SurfaceScalarField::zeros("phi", m.clone());
290        // Internal: 1 face; boundary: 2 patches × 1 face each
291        assert_eq!(phi.internal.len(), 1);
292        assert_eq!(phi.boundary.len(), 2);
293        assert_eq!(phi.boundary[0].values.len(), 1);
294    }
295
296    #[test]
297    fn face_value_internal() {
298        let m = unit_mesh();
299        let phi = SurfaceScalarField::uniform("phi", m.clone(), 5.0);
300        assert!((phi.face_value(0) - 5.0).abs() < 1e-15); // internal face 0
301    }
302
303    #[test]
304    fn add_surface_fields() {
305        let m = unit_mesh();
306        let a = SurfaceScalarField::uniform("a", m.clone(), 2.0);
307        let b = SurfaceScalarField::uniform("b", m, 3.0);
308        let c = a + b;
309        assert!((c.internal[0] - 5.0).abs() < 1e-15);
310    }
311
312    #[test]
313    fn mul_surface_field_by_scalar() {
314        let m = unit_mesh();
315        let a = SurfaceScalarField::uniform("a", m, 4.0);
316        let b = a * 0.5;
317        assert!((b.internal[0] - 2.0).abs() < 1e-15);
318    }
319}
outram_foam_basic_lib/fields/surface_field.rs

outram_foam_basic_lib/fields/
surface_field.rs
