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

outram_foam_basic_lib/fields/
surface_field.rs
