mesh_sieve/data/
bc.rs

1//! Boundary condition helpers using label queries.
2
3use crate::data::constrained_section::ConstrainedSection;
4use crate::data::section::Section;
5use crate::data::storage::Storage;
6use crate::mesh_error::MeshSieveError;
7use crate::physics::fvm::{
8    BoundaryBranchError, BoundaryCondition, FvBoundaryBranch, boundary_branch_for_face_checked,
9};
10use crate::topology::cache::InvalidateCache;
11use crate::topology::coastal::CoastalLabelQueries;
12use crate::topology::labels::LabelSet;
13use crate::topology::point::PointId;
14use std::collections::{HashMap, HashSet};
15
16#[derive(Clone, Debug, PartialEq, Eq)]
17pub enum CoastalBoundaryAssemblyError {
18    InvalidBoundaryFaceRole(BoundaryBranchError),
19}
20
21/// Label query selector for boundary condition application.
22#[derive(Clone, Debug, PartialEq, Eq)]
23pub struct LabelQuery {
24    name: String,
25    value: i32,
26}
27
28impl LabelQuery {
29    /// Create a new label query.
30    pub fn new(name: impl Into<String>, value: i32) -> Self {
31        Self {
32            name: name.into(),
33            value,
34        }
35    }
36
37    /// Label name.
38    pub fn name(&self) -> &str {
39        &self.name
40    }
41
42    /// Label value.
43    pub fn value(&self) -> i32 {
44        self.value
45    }
46}
47
48/// Per-field DOF indices for packed multi-field sections.
49#[derive(Clone, Debug, PartialEq, Eq)]
50pub struct FieldDofIndices {
51    /// Field index in the packed layout.
52    pub field: usize,
53    /// DOF indices within the field slice.
54    pub dof_indices: Vec<usize>,
55}
56
57impl FieldDofIndices {
58    /// Create a new field DOF index set.
59    pub fn new(field: usize, dof_indices: impl Into<Vec<usize>>) -> Self {
60        Self {
61            field,
62            dof_indices: dof_indices.into(),
63        }
64    }
65}
66
67fn label_points(labels: &LabelSet, query: &LabelQuery) -> Vec<PointId> {
68    labels.stratum_points(query.name(), query.value())
69}
70
71/// Boundary face sets derived from canonical coastal labels.
72#[derive(Clone, Debug, Default, PartialEq, Eq)]
73pub struct CoastalBoundaryFaceSets {
74    pub free_surface: Vec<PointId>,
75    pub bed: Vec<PointId>,
76    pub open: Vec<PointId>,
77    pub inflow: Vec<PointId>,
78    pub outflow: Vec<PointId>,
79    pub tidal: Vec<PointId>,
80}
81
82/// Build boundary-face sets by intersecting supplied boundary faces with coastal labels.
83pub fn coastal_boundary_face_sets(
84    labels: &LabelSet,
85    boundary_faces: impl IntoIterator<Item = PointId>,
86) -> CoastalBoundaryFaceSets {
87    let faces: HashSet<_> = boundary_faces.into_iter().collect();
88    let filter = |pts: Vec<PointId>| -> Vec<PointId> {
89        let mut v: Vec<_> = pts.into_iter().filter(|p| faces.contains(p)).collect();
90        v.sort_unstable();
91        v
92    };
93    CoastalBoundaryFaceSets {
94        free_surface: filter(labels.free_surface_points()),
95        bed: filter(labels.bed_points()),
96        open: filter(labels.open_boundary_points()),
97        inflow: filter(labels.inflow_points()),
98        outflow: filter(labels.outflow_points()),
99        tidal: filter(labels.tidal_points()),
100    }
101}
102
103/// Assign per-face boundary conditions using coastal boundary class/role labels.
104pub fn map_coastal_boundary_conditions(
105    labels: &LabelSet,
106    boundary_faces: impl IntoIterator<Item = PointId>,
107    branch_closure: impl Fn(FvBoundaryBranch, PointId) -> BoundaryCondition,
108) -> HashMap<PointId, BoundaryCondition> {
109    let mut out = HashMap::new();
110    for face in boundary_faces {
111        if let Ok(branch) = boundary_branch_for_face_checked(labels, face) {
112            out.insert(face, branch_closure(branch, face));
113        }
114    }
115    out
116}
117
118/// Resolved coastal boundary data ready for FV flux assembly.
119#[derive(Clone, Debug, Default, PartialEq)]
120pub struct CoastalBoundaryAssembly {
121    pub face_sets: CoastalBoundaryFaceSets,
122    pub boundary_conditions: HashMap<PointId, BoundaryCondition>,
123    pub branches: HashMap<PointId, FvBoundaryBranch>,
124}
125
126/// High-level coastal API: resolve boundary face groups and BC closures in one pass.
127pub fn resolve_coastal_boundary_assembly(
128    labels: &LabelSet,
129    boundary_faces: impl IntoIterator<Item = PointId>,
130    branch_closure: impl Fn(FvBoundaryBranch, PointId) -> BoundaryCondition,
131) -> Result<CoastalBoundaryAssembly, CoastalBoundaryAssemblyError> {
132    let faces: Vec<_> = boundary_faces.into_iter().collect();
133    let face_sets = coastal_boundary_face_sets(labels, faces.iter().copied());
134    let mut boundary_conditions = HashMap::new();
135    let mut branches = HashMap::new();
136    for face in faces {
137        let branch = boundary_branch_for_face_checked(labels, face)
138            .map_err(CoastalBoundaryAssemblyError::InvalidBoundaryFaceRole)?;
139        boundary_conditions.insert(face, branch_closure(branch, face));
140        branches.insert(face, branch);
141    }
142    Ok(CoastalBoundaryAssembly {
143        face_sets,
144        boundary_conditions,
145        branches,
146    })
147}
148
149fn field_offsets(field_dofs: &[usize]) -> Vec<usize> {
150    let mut offsets = Vec::with_capacity(field_dofs.len());
151    let mut total = 0usize;
152    for dof in field_dofs {
153        offsets.push(total);
154        total += *dof;
155    }
156    offsets
157}
158
159/// Apply Dirichlet constraints directly to a section for all points matching a label query.
160pub fn apply_dirichlet_to_section<V, S, F>(
161    section: &mut Section<V, S>,
162    labels: &LabelSet,
163    query: &LabelQuery,
164    dof_indices: &[usize],
165    mut value: F,
166) -> Result<(), MeshSieveError>
167where
168    V: Clone,
169    S: Storage<V>,
170    F: FnMut(PointId, usize) -> V,
171{
172    for point in label_points(labels, query) {
173        let slice = section.try_restrict_mut(point)?;
174        let len = slice.len();
175        for &index in dof_indices {
176            if index >= len {
177                return Err(MeshSieveError::ConstraintIndexOutOfBounds { point, index, len });
178            }
179            slice[index] = value(point, index);
180        }
181    }
182    section.invalidate_cache();
183    Ok(())
184}
185
186/// Apply Dirichlet constraints directly to a section using packed per-field DOF indices.
187pub fn apply_dirichlet_to_section_fields<V, S, F>(
188    section: &mut Section<V, S>,
189    labels: &LabelSet,
190    query: &LabelQuery,
191    field_dofs: &[usize],
192    field_indices: &[FieldDofIndices],
193    mut value: F,
194) -> Result<(), MeshSieveError>
195where
196    V: Clone,
197    S: Storage<V>,
198    F: FnMut(PointId, usize, usize) -> V,
199{
200    let offsets = field_offsets(field_dofs);
201    for point in label_points(labels, query) {
202        let slice = section.try_restrict_mut(point)?;
203        let len = slice.len();
204        for field_spec in field_indices {
205            let field = field_spec.field;
206            let field_len =
207                *field_dofs
208                    .get(field)
209                    .ok_or_else(|| MeshSieveError::SectionAccess {
210                        point,
211                        source: Box::new(std::io::Error::new(
212                            std::io::ErrorKind::InvalidInput,
213                            "field index out of bounds",
214                        )),
215                    })?;
216            let base = offsets[field];
217            for &dof in &field_spec.dof_indices {
218                if dof >= field_len {
219                    return Err(MeshSieveError::ConstraintIndexOutOfBounds {
220                        point,
221                        index: dof,
222                        len: field_len,
223                    });
224                }
225                let index = base + dof;
226                if index >= len {
227                    return Err(MeshSieveError::ConstraintIndexOutOfBounds { point, index, len });
228                }
229                slice[index] = value(point, field, dof);
230            }
231        }
232    }
233    section.invalidate_cache();
234    Ok(())
235}
236
237/// Apply Dirichlet constraints to a constrained section for all points matching a label query.
238pub fn apply_dirichlet_to_constrained_section<V, S, F>(
239    section: &mut ConstrainedSection<V, S>,
240    labels: &LabelSet,
241    query: &LabelQuery,
242    dof_indices: &[usize],
243    mut value: F,
244) -> Result<(), MeshSieveError>
245where
246    V: Clone,
247    S: Storage<V>,
248    F: FnMut(PointId, usize) -> V,
249{
250    for point in label_points(labels, query) {
251        for &index in dof_indices {
252            let val = value(point, index);
253            section.insert_constraint(point, index, val)?;
254        }
255    }
256    section.apply_constraints()
257}
258
259/// Apply Dirichlet constraints to a constrained section using packed per-field DOF indices.
260pub fn apply_dirichlet_to_constrained_section_fields<V, S, F>(
261    section: &mut ConstrainedSection<V, S>,
262    labels: &LabelSet,
263    query: &LabelQuery,
264    field_dofs: &[usize],
265    field_indices: &[FieldDofIndices],
266    mut value: F,
267) -> Result<(), MeshSieveError>
268where
269    V: Clone,
270    S: Storage<V>,
271    F: FnMut(PointId, usize, usize) -> V,
272{
273    let offsets = field_offsets(field_dofs);
274    for point in label_points(labels, query) {
275        for field_spec in field_indices {
276            let field = field_spec.field;
277            let field_len =
278                *field_dofs
279                    .get(field)
280                    .ok_or_else(|| MeshSieveError::SectionAccess {
281                        point,
282                        source: Box::new(std::io::Error::new(
283                            std::io::ErrorKind::InvalidInput,
284                            "field index out of bounds",
285                        )),
286                    })?;
287            let base = offsets[field];
288            for &dof in &field_spec.dof_indices {
289                if dof >= field_len {
290                    return Err(MeshSieveError::ConstraintIndexOutOfBounds {
291                        point,
292                        index: dof,
293                        len: field_len,
294                    });
295                }
296                let index = base + dof;
297                let val = value(point, field, dof);
298                section.insert_constraint(point, index, val)?;
299            }
300        }
301    }
302    section.apply_constraints()
303}
304
305#[cfg(test)]
306mod tests {
307    use super::*;
308    use crate::discretization::runtime::FluxStencil;
309    use crate::physics::fvm::{FvmInputs, flux_activity_mask_from_wet_dry};
310    use crate::topology::coastal::{
311        BOUNDARY_CLASS_LABEL, BOUNDARY_ROLE_LABEL, BoundaryClass, OpenBoundaryRole,
312        WET_DRY_MASK_LABEL, WetDryMask,
313    };
314
315    fn p(id: u64) -> PointId {
316        PointId::new(id).unwrap()
317    }
318
319    #[test]
320    fn resolves_all_boundary_branches_to_bc_closures() {
321        let mut labels = LabelSet::new();
322        labels.set_label(
323            p(10),
324            BOUNDARY_CLASS_LABEL,
325            BoundaryClass::FreeSurface.code(),
326        );
327        labels.set_label(p(11), BOUNDARY_CLASS_LABEL, BoundaryClass::Bed.code());
328        labels.set_label(p(12), BOUNDARY_CLASS_LABEL, BoundaryClass::Open.code());
329        labels.set_label(p(12), BOUNDARY_ROLE_LABEL, OpenBoundaryRole::Inflow.code());
330        labels.set_label(p(13), BOUNDARY_CLASS_LABEL, BoundaryClass::Open.code());
331        labels.set_label(p(13), BOUNDARY_ROLE_LABEL, OpenBoundaryRole::Outflow.code());
332        labels.set_label(p(14), BOUNDARY_CLASS_LABEL, BoundaryClass::Open.code());
333        labels.set_label(p(14), BOUNDARY_ROLE_LABEL, OpenBoundaryRole::Tidal.code());
334
335        let resolved = resolve_coastal_boundary_assembly(
336            &labels,
337            [p(10), p(11), p(12), p(13), p(14)],
338            |branch, _| match branch {
339                FvBoundaryBranch::FreeSurface => BoundaryCondition::Neumann { gradient: 1.0 },
340                FvBoundaryBranch::Bed => BoundaryCondition::Neumann { gradient: -1.0 },
341                FvBoundaryBranch::Inflow => BoundaryCondition::Dirichlet { value: 10.0 },
342                FvBoundaryBranch::Outflow => BoundaryCondition::Robin {
343                    alpha: 1.0,
344                    beta: 0.5,
345                    gamma: 0.2,
346                },
347                FvBoundaryBranch::Tidal => BoundaryCondition::Dirichlet { value: 7.0 },
348                FvBoundaryBranch::Open => unreachable!(),
349            },
350        )
351        .unwrap();
352        assert_eq!(resolved.face_sets.free_surface, vec![p(10)]);
353        assert_eq!(resolved.face_sets.bed, vec![p(11)]);
354        assert_eq!(resolved.face_sets.inflow, vec![p(12)]);
355        assert_eq!(resolved.face_sets.outflow, vec![p(13)]);
356        assert_eq!(resolved.face_sets.tidal, vec![p(14)]);
357        assert_eq!(resolved.boundary_conditions.len(), 5);
358    }
359
360    #[test]
361    fn errors_on_missing_open_role_in_assembly_faces() {
362        let mut labels = LabelSet::new();
363        labels.set_label(p(12), BOUNDARY_CLASS_LABEL, BoundaryClass::Open.code());
364        let err = resolve_coastal_boundary_assembly(&labels, [p(12)], |_, _| {
365            BoundaryCondition::Dirichlet { value: 0.0 }
366        })
367        .unwrap_err();
368        assert!(matches!(
369            err,
370            CoastalBoundaryAssemblyError::InvalidBoundaryFaceRole(BoundaryBranchError::OpenBoundaryMissingRole { face }) if face == p(12)
371        ));
372    }
373
374    #[test]
375    fn wet_dry_mask_permutations_and_mixed_boundary_sets() {
376        let c0 = p(1);
377        let c1 = p(2);
378        let c2 = p(3);
379        let b0 = p(10);
380        let b1 = p(11);
381        let i0 = p(20);
382        let inputs = FvmInputs::new(
383            [
384                FluxStencil {
385                    face: b0,
386                    left: c0,
387                    right: None,
388                },
389                FluxStencil {
390                    face: b1,
391                    left: c1,
392                    right: None,
393                },
394                FluxStencil {
395                    face: i0,
396                    left: c1,
397                    right: Some(c2),
398                },
399            ],
400            vec![],
401            vec![],
402        );
403        let mut labels = LabelSet::new();
404        labels.set_label(b0, BOUNDARY_CLASS_LABEL, BoundaryClass::Open.code());
405        labels.set_label(b0, BOUNDARY_ROLE_LABEL, OpenBoundaryRole::Inflow.code());
406        labels.set_label(b1, BOUNDARY_CLASS_LABEL, BoundaryClass::Bed.code());
407        labels.set_label(i0, BOUNDARY_CLASS_LABEL, BoundaryClass::Open.code());
408        labels.set_label(c1, WET_DRY_MASK_LABEL, WetDryMask::Dry.code());
409        labels.set_label(c2, WET_DRY_MASK_LABEL, WetDryMask::Wet.code());
410        let mask = flux_activity_mask_from_wet_dry(&inputs, &labels);
411        assert_eq!(mask.boundary_faces_active.get(&b0), Some(&true));
412        assert_eq!(mask.boundary_faces_active.get(&b1), Some(&false));
413        assert_eq!(mask.near_boundary_faces_active.get(&i0), Some(&false));
414    }
415}
mesh_sieve/data/bc.rs

mesh_sieve/data/
bc.rs
