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