outram_foam_basic_lib/mesh/region_interface.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::sync::Arc;
23
24use crate::fields::boundary::bc::{BoundaryCondition, PatchField};
25use crate::fields::field::Field;
26use crate::fields::vol_field::VolScalarField;
27use crate::mesh::fv_mesh::FvMesh;
28use crate::primitives::Vector3;
29
30/// Face-to-face mapping between two mesh patches at a shared interface.
31///
32/// Used by `chtMultiRegionFoam`-style solvers where a fluid region and a
33/// solid region share an interface. Each side has a patch (identified by
34/// mesh + patch index); the `face_map` gives the paired face index on side B
35/// for each face on side A.
36///
37/// For matching meshes (same layout, same face count) `face_map[i] = i`.
38/// For non-matching meshes (different refinements) the map is built by
39/// nearest-face-centre search (see `from_face_centres`).
40#[derive(Debug, Clone)]
41pub struct RegionInterface {
42 pub mesh_a: Arc<FvMesh>,
43 pub patch_a: usize,
44 pub mesh_b: Arc<FvMesh>,
45 pub patch_b: usize,
46 /// `face_map[fi_a]` = `fi_b` on the B-side patch.
47 pub face_map: Vec<usize>,
48}
49
50impl RegionInterface {
51 /// Construct a matching interface: face `i` on A is coupled to face `i` on B.
52 ///
53 /// Panics if the two patches have different sizes.
54 pub fn matching(
55 mesh_a: Arc<FvMesh>,
56 patch_a: usize,
57 mesh_b: Arc<FvMesh>,
58 patch_b: usize,
59 ) -> Self {
60 let n_a = mesh_a.patches[patch_a].size;
61 let n_b = mesh_b.patches[patch_b].size;
62 assert_eq!(n_a, n_b, "matching interface requires equal patch sizes");
63 Self {
64 mesh_a,
65 patch_a,
66 mesh_b,
67 patch_b,
68 face_map: (0..n_a).collect(),
69 }
70 }
71
72 /// Construct a non-matching interface via nearest-face-centre search.
73 ///
74 /// For each face on patch A, the closest face centre on patch B is found
75 /// and stored in `face_map`. O(n_A × n_B) — acceptable for boundary
76 /// patches which are typically small.
77 pub fn from_face_centres(
78 mesh_a: Arc<FvMesh>,
79 patch_a: usize,
80 mesh_b: Arc<FvMesh>,
81 patch_b: usize,
82 ) -> Self {
83 let patch_info_a = &mesh_a.patches[patch_a];
84 let patch_info_b = &mesh_b.patches[patch_b];
85
86 let centres_b: Vec<Vector3> = (0..patch_info_b.size)
87 .map(|fi| mesh_b.face_centres[patch_info_b.start + fi])
88 .collect();
89
90 let face_map: Vec<usize> = (0..patch_info_a.size)
91 .map(|fi_a| {
92 let ca = mesh_a.face_centres[patch_info_a.start + fi_a];
93 centres_b
94 .iter()
95 .enumerate()
96 .min_by(|(_, cb1), (_, cb2)| {
97 (ca - **cb1)
98 .mag_sqr()
99 .partial_cmp(&(ca - **cb2).mag_sqr())
100 .unwrap()
101 })
102 .map(|(i, _)| i)
103 .unwrap_or(0)
104 })
105 .collect();
106
107 Self {
108 mesh_a,
109 patch_a,
110 mesh_b,
111 patch_b,
112 face_map,
113 }
114 }
115
116 /// Exchange temperature boundary values at the interface.
117 ///
118 /// Returns `(bc_a_new, bc_b_new)` where each is a `PatchField` of
119 /// `FixedValue` temperatures drawn from the cell-centre values of the
120 /// opposite region's cells adjacent to the interface.
121 ///
122 /// This implements the simplest coupling: Dirichlet T on each side equals
123 /// the adjacent cell T of the other side. A more accurate coupling would
124 /// use the heat-flux equality condition; that requires knowing `kappa` on
125 /// both sides (available from `SolidThermo::kappa()`).
126 pub fn exchange_temperature(
127 &self,
128 t_a: &VolScalarField,
129 t_b: &VolScalarField,
130 ) -> (PatchField<f64>, PatchField<f64>) {
131 let pa = &self.mesh_a.patches[self.patch_a];
132 let pb = &self.mesh_b.patches[self.patch_b];
133
134 // Values on A-side face come from adjacent B-side cells
135 let vals_a: Vec<f64> = (0..pa.size)
136 .map(|fi_a| {
137 let fi_b = self.face_map[fi_a];
138 let gf_b = pb.start + fi_b;
139 let cell_b = self.mesh_b.owner[gf_b];
140 t_b.internal[cell_b]
141 })
142 .collect();
143
144 // Values on B-side face come from adjacent A-side cells
145 let n_b = pb.size;
146 let mut vals_b = vec![0.0; n_b];
147 for fi_a in 0..pa.size {
148 let fi_b = self.face_map[fi_a];
149 let gf_a = pa.start + fi_a;
150 let cell_a = self.mesh_a.owner[gf_a];
151 vals_b[fi_b] = t_a.internal[cell_a];
152 }
153
154 let bc_a = PatchField {
155 bc: BoundaryCondition::FixedValue(0.0), // sentinel; values override
156 values: Field::new(vals_a),
157 };
158 let bc_b = PatchField {
159 bc: BoundaryCondition::FixedValue(0.0),
160 values: Field::new(vals_b),
161 };
162 (bc_a, bc_b)
163 }
164}
165
166#[cfg(test)]
167mod tests {
168 use super::*;
169 use crate::fields::vol_field::VolScalarField;
170 use crate::mesh::fv_mesh::{BoundaryPatch, FvMeshBuilder, PatchKind};
171
172 fn simple_mesh(_hot: bool) -> Arc<FvMesh> {
173 Arc::new(
174 FvMeshBuilder::new()
175 .n_cells(2)
176 .n_internal_faces(1)
177 .owner(vec![0, 1, 0])
178 .neighbour(vec![1])
179 .patches(vec![
180 BoundaryPatch::new("interface", 1, 1, PatchKind::Wall),
181 BoundaryPatch::new("outer", 2, 1, PatchKind::Wall),
182 ])
183 .cell_volumes(vec![1.0, 1.0])
184 .cell_centres(vec![
185 Vector3::new(0.25, 0.0, 0.0),
186 Vector3::new(0.75, 0.0, 0.0),
187 ])
188 .face_area_vectors(vec![
189 Vector3::new(1.0, 0.0, 0.0),
190 Vector3::new(1.0, 0.0, 0.0),
191 Vector3::new(-1.0, 0.0, 0.0),
192 ])
193 .face_centres(vec![
194 Vector3::new(0.5, 0.0, 0.0),
195 Vector3::new(1.0, 0.0, 0.0),
196 Vector3::new(0.0, 0.0, 0.0),
197 ])
198 .build()
199 .unwrap(),
200 )
201 }
202
203 #[test]
204 fn exchange_temperature_matching() {
205 let m_a = simple_mesh(true);
206 let m_b = simple_mesh(false);
207 let iface = RegionInterface::matching(m_a.clone(), 0, m_b.clone(), 0);
208
209 let t_a = VolScalarField::uniform("T_a", m_a.clone(), 500.0);
210 let t_b = VolScalarField::uniform("T_b", m_b.clone(), 300.0);
211 let (bc_a, bc_b) = iface.exchange_temperature(&t_a, &t_b);
212
213 // A-side BC gets B's adjacent cell temperature (300)
214 assert!((bc_a.values[0] - 300.0).abs() < 1e-10);
215 // B-side BC gets A's adjacent cell temperature (500)
216 assert!((bc_b.values[0] - 500.0).abs() < 1e-10);
217 }
218
219 #[test]
220 fn non_matching_interface_nearest_centre() {
221 let m_a = simple_mesh(true);
222 let m_b = simple_mesh(false);
223 let iface = RegionInterface::from_face_centres(m_a.clone(), 0, m_b.clone(), 0);
224 // Should map face 0 → face 0 (same centre)
225 assert_eq!(iface.face_map[0], 0);
226 }
227}