outram_foam_basic_lib/io/
poly_mesh.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
22//! `constant/polyMesh` read/write.
23//!
24//! The crate's [`FvMesh`] stores only flat finite-volume geometry
25//! (cell/face centres, areas, owner/neighbour) with no point/face-vertex
26//! connectivity. OpenFOAM's on-disk `polyMesh`, by contrast, is defined by its
27//! **connectivity**: `points` (vertices), `faces` (vertex loops),
28//! `owner`/`neighbour` (cell adjacency), and `boundary` (patches). This module
29//! therefore defines a connectivity-carrying [`PolyMesh`] as the I/O
30//! representation and computes full FV geometry from it via
31//! [`PolyMesh::to_fv_mesh`] — the same divergence-theorem pyramid
32//! decomposition OpenFOAM's `primitiveMesh` uses (mirrored from
33//! `outram-foam-mesh`'s `block_mesh` / `poly_dual_mesh`).
34//!
35//! ## Files
36//!
37//! | file | class | contents |
38//! |---|---|---|
39//! | `points`    | `vectorField`      | vertex coordinates `[m]` |
40//! | `faces`     | `faceList`         | each face as a vertex-index loop |
41//! | `owner`     | `labelList`        | owner cell per face |
42//! | `neighbour` | `labelList`        | neighbour cell per internal face |
43//! | `boundary`  | `polyBoundaryMesh` | patches (`type`, `nFaces`, `startFace`) |
44//!
45//! Faces are ordered OpenFOAM-style: internal faces first
46//! (`[0, n_internal_faces)`), then boundary faces grouped by patch.
47
48use std::path::Path;
49
50use super::dict::{
51    fmt_scalar, tokenize, write_header, Token, BANNER, FOOTER, SEPARATOR,
52};
53use super::IoError;
54use crate::mesh::{BoundaryPatch, FvMesh, FvMeshBuilder, PatchKind};
55use crate::primitives::Vector3;
56
57/// A single mesh face: its point-index loop plus owner / neighbour cells.
58///
59/// `verts` is wound so the face normal points **from `owner` toward
60/// `neighbour`** (outward from the owner cell). Boundary faces have
61/// `neighbour == None`.
62#[derive(Debug, Clone, PartialEq)]
63pub struct MeshFace {
64    /// Ordered point indices (into [`PolyMesh::points`]) forming the face loop.
65    pub verts: Vec<usize>,
66    /// Owning cell index.
67    pub owner: usize,
68    /// Neighbour cell index (internal faces only).
69    pub neighbour: Option<usize>,
70}
71
72/// A connectivity-carrying poly-mesh — the on-disk `polyMesh` representation.
73///
74/// Faces are ordered internal-first (`[0, n_internal_faces)`), then boundary
75/// faces grouped by patch, matching OpenFOAM. Call [`PolyMesh::to_fv_mesh`]
76/// to obtain the geometry-carrying [`FvMesh`].
77#[derive(Debug, Clone, PartialEq)]
78pub struct PolyMesh {
79    /// Mesh points `[m]`.
80    pub points: Vec<Vector3>,
81    /// All faces, internal first then boundary.
82    pub faces: Vec<MeshFace>,
83    /// Number of internal faces (leading internal entries in `faces`).
84    pub n_internal_faces: usize,
85    /// Number of cells.
86    pub n_cells: usize,
87    /// Boundary patches, covering `[n_internal_faces, faces.len())`.
88    pub patches: Vec<BoundaryPatch>,
89}
90
91impl PolyMesh {
92    /// Number of points.
93    pub fn n_points(&self) -> usize {
94        self.points.len()
95    }
96    /// Total number of faces (internal + boundary).
97    pub fn n_faces(&self) -> usize {
98        self.faces.len()
99    }
100    /// Number of boundary faces.
101    pub fn n_boundary_faces(&self) -> usize {
102        self.faces.len() - self.n_internal_faces
103    }
104
105    /// Total mesh volume `[m^3]` — the sum of all cell volumes.
106    pub fn total_volume(&self) -> f64 {
107        self.cell_geometry().0.iter().sum()
108    }
109
110    /// Compute per-cell `(volumes, centres)` `[m^3]`, `[m]` via the OpenFOAM
111    /// pyramid decomposition (`primitiveMeshCellCentresAndVols`).
112    fn cell_geometry(&self) -> (Vec<f64>, Vec<Vector3>) {
113        let mut face_ctr = Vec::with_capacity(self.faces.len());
114        let mut face_sf = Vec::with_capacity(self.faces.len());
115        for f in &self.faces {
116            let (c, sf) = face_centre_and_area(&self.points, &f.verts);
117            face_ctr.push(c);
118            face_sf.push(sf);
119        }
120
121        // First pass: estimate each cell centre as the mean of its face centres.
122        let mut c_est = vec![Vector3::ZERO; self.n_cells];
123        let mut n_cell_faces = vec![0.0f64; self.n_cells];
124        for (fi, f) in self.faces.iter().enumerate() {
125            c_est[f.owner] += face_ctr[fi];
126            n_cell_faces[f.owner] += 1.0;
127            if let Some(nb) = f.neighbour {
128                c_est[nb] += face_ctr[fi];
129                n_cell_faces[nb] += 1.0;
130            }
131        }
132        for c in 0..self.n_cells {
133            if n_cell_faces[c] > 0.0 {
134                c_est[c] /= n_cell_faces[c];
135            }
136        }
137
138        // Second pass: sum pyramid volumes / centroids about the estimate.
139        let mut vol = vec![0.0f64; self.n_cells];
140        let mut ctr = vec![Vector3::ZERO; self.n_cells];
141        for (fi, f) in self.faces.iter().enumerate() {
142            let cf = face_ctr[fi];
143            {
144                let cell = f.owner;
145                let sf = face_sf[fi];
146                let pyr3 = sf.dot(cf - c_est[cell]);
147                let pyr_ctr = cf * 0.75 + c_est[cell] * 0.25;
148                vol[cell] += pyr3;
149                ctr[cell] += pyr_ctr * pyr3;
150            }
151            if let Some(cell) = f.neighbour {
152                let sf = face_sf[fi] * -1.0;
153                let pyr3 = sf.dot(cf - c_est[cell]);
154                let pyr_ctr = cf * 0.75 + c_est[cell] * 0.25;
155                vol[cell] += pyr3;
156                ctr[cell] += pyr_ctr * pyr3;
157            }
158        }
159        for c in 0..self.n_cells {
160            if vol[c].abs() > f64::EPSILON {
161                ctr[c] /= vol[c];
162            } else {
163                ctr[c] = c_est[c];
164            }
165            vol[c] /= 3.0;
166        }
167        (vol, ctr)
168    }
169
170    /// Convert to the geometry-carrying [`FvMesh`], computing cell
171    /// volumes/centres, face-area vectors, face areas, and face centres from
172    /// the point/face topology.
173    ///
174    /// # Errors
175    /// [`IoError::Mesh`] if the assembled mesh fails `FvMesh::validate`.
176    pub fn to_fv_mesh(&self) -> Result<FvMesh, IoError> {
177        let (cell_volumes, cell_centres) = self.cell_geometry();
178
179        let mut owner = Vec::with_capacity(self.faces.len());
180        let mut neighbour = Vec::with_capacity(self.n_internal_faces);
181        let mut face_area_vectors = Vec::with_capacity(self.faces.len());
182        let mut face_centres = Vec::with_capacity(self.faces.len());
183        for f in &self.faces {
184            let (c, sf) = face_centre_and_area(&self.points, &f.verts);
185            owner.push(f.owner);
186            if let Some(nb) = f.neighbour {
187                neighbour.push(nb);
188            }
189            face_area_vectors.push(sf);
190            face_centres.push(c);
191        }
192
193        FvMeshBuilder::new()
194            .n_cells(self.n_cells)
195            .n_internal_faces(self.n_internal_faces)
196            .owner(owner)
197            .neighbour(neighbour)
198            .patches(self.patches.clone())
199            .cell_volumes(cell_volumes)
200            .cell_centres(cell_centres)
201            .face_area_vectors(face_area_vectors)
202            .face_centres(face_centres)
203            .build()
204            .map_err(|e| IoError::Mesh(e.to_string()))
205    }
206
207    // ── Reading ─────────────────────────────────────────────────────────────
208
209    /// Read a `polyMesh` from a `constant/polyMesh` directory.
210    ///
211    /// Reads `points`, `faces`, `owner`, `neighbour`, and `boundary`.
212    pub fn read(dir: impl AsRef<Path>) -> Result<Self, IoError> {
213        let dir = dir.as_ref();
214        let points = read_points(&dir.join("points"))?;
215        let face_loops = read_faces(&dir.join("faces"))?;
216        let owner = read_labels(&dir.join("owner"))?;
217        let neighbour = read_labels(&dir.join("neighbour"))?;
218        let raw_patches = read_boundary(&dir.join("boundary"))?;
219
220        if owner.len() != face_loops.len() {
221            return Err(IoError::parse(
222                "polyMesh",
223                format!(
224                    "owner has {} entries but faces has {}",
225                    owner.len(),
226                    face_loops.len()
227                ),
228            ));
229        }
230        let n_internal_faces = neighbour.len();
231        let mut faces = Vec::with_capacity(face_loops.len());
232        for (fi, verts) in face_loops.into_iter().enumerate() {
233            let nb = if fi < n_internal_faces {
234                Some(neighbour[fi])
235            } else {
236                None
237            };
238            faces.push(MeshFace {
239                verts,
240                owner: owner[fi],
241                neighbour: nb,
242            });
243        }
244
245        // Cell count = highest referenced cell index + 1.
246        let mut n_cells = 0usize;
247        for f in &faces {
248            n_cells = n_cells.max(f.owner + 1);
249            if let Some(nb) = f.neighbour {
250                n_cells = n_cells.max(nb + 1);
251            }
252        }
253
254        let patches = raw_patches
255            .into_iter()
256            .map(|p| BoundaryPatch::new(p.name, p.start_face, p.n_faces, kind_from_str(&p.kind)))
257            .collect();
258
259        Ok(Self {
260            points,
261            faces,
262            n_internal_faces,
263            n_cells,
264            patches,
265        })
266    }
267
268    // ── Writing ─────────────────────────────────────────────────────────────
269
270    /// Write the `polyMesh` files into `dir`, creating it if necessary.
271    pub fn write(&self, dir: impl AsRef<Path>) -> Result<(), IoError> {
272        let dir = dir.as_ref();
273        std::fs::create_dir_all(dir).map_err(|e| IoError::io(dir.display().to_string(), e))?;
274        write_file(&dir.join("points"), &self.points_string())?;
275        write_file(&dir.join("faces"), &self.faces_string())?;
276        write_file(&dir.join("owner"), &self.owner_string())?;
277        write_file(&dir.join("neighbour"), &self.neighbour_string())?;
278        write_file(&dir.join("boundary"), &self.boundary_string())?;
279        Ok(())
280    }
281
282    fn points_string(&self) -> String {
283        let mut s = header_str("vectorField", "points");
284        s.push_str(&format!("{}\n(\n", self.points.len()));
285        for p in &self.points {
286            s.push_str(&format!(
287                "({} {} {})\n",
288                fmt_scalar(p.x),
289                fmt_scalar(p.y),
290                fmt_scalar(p.z)
291            ));
292        }
293        s.push_str(")\n\n");
294        s.push_str(FOOTER);
295        s
296    }
297
298    fn faces_string(&self) -> String {
299        let mut s = header_str("faceList", "faces");
300        s.push_str(&format!("{}\n(\n", self.faces.len()));
301        for f in &self.faces {
302            s.push_str(&format!("{}(", f.verts.len()));
303            for (i, v) in f.verts.iter().enumerate() {
304                if i > 0 {
305                    s.push(' ');
306                }
307                s.push_str(&v.to_string());
308            }
309            s.push_str(")\n");
310        }
311        s.push_str(")\n\n");
312        s.push_str(FOOTER);
313        s
314    }
315
316    fn owner_string(&self) -> String {
317        let mut s = header_str("labelList", "owner");
318        s.push_str(&format!("{}\n(\n", self.faces.len()));
319        for f in &self.faces {
320            s.push_str(&f.owner.to_string());
321            s.push('\n');
322        }
323        s.push_str(")\n\n");
324        s.push_str(FOOTER);
325        s
326    }
327
328    fn neighbour_string(&self) -> String {
329        let mut s = header_str("labelList", "neighbour");
330        s.push_str(&format!("{}\n(\n", self.n_internal_faces));
331        for f in self.faces.iter().take(self.n_internal_faces) {
332            s.push_str(&f.neighbour.unwrap_or(0).to_string());
333            s.push('\n');
334        }
335        s.push_str(")\n\n");
336        s.push_str(FOOTER);
337        s
338    }
339
340    fn boundary_string(&self) -> String {
341        let mut s = header_str("polyBoundaryMesh", "boundary");
342        s.push_str(&format!("{}\n(\n", self.patches.len()));
343        for p in &self.patches {
344            s.push_str(&format!("    {}\n    {{\n", p.name));
345            s.push_str(&format!("        type            {};\n", kind_to_str(p.kind)));
346            s.push_str(&format!("        nFaces          {};\n", p.size));
347            s.push_str(&format!("        startFace       {};\n", p.start));
348            s.push_str("    }\n");
349        }
350        s.push_str(")\n\n");
351        s.push_str(FOOTER);
352        s
353    }
354}
355
356// ── Geometry (mirrors outram-foam-mesh::block_mesh) ─────────────────────────
357
358/// Face centre `[m]` and area vector `Sf` `[m^2]` from a vertex loop, using the
359/// OpenFOAM decomposition of the polygon into triangles about its centroid.
360fn face_centre_and_area(points: &[Vector3], verts: &[usize]) -> (Vector3, Vector3) {
361    let n = verts.len();
362    debug_assert!(n >= 3);
363    if n == 3 {
364        let a = points[verts[0]];
365        let b = points[verts[1]];
366        let c = points[verts[2]];
367        let centre = (a + b + c) / 3.0;
368        let area = (b - a).cross(c - a) * 0.5;
369        return (centre, area);
370    }
371    let mut c_est = Vector3::ZERO;
372    for &v in verts {
373        c_est += points[v];
374    }
375    c_est /= n as f64;
376
377    let mut sum_a = 0.0;
378    let mut sum_n = Vector3::ZERO;
379    let mut sum_ac = Vector3::ZERO;
380    for i in 0..n {
381        let p1 = points[verts[i]];
382        let p2 = points[verts[(i + 1) % n]];
383        let mid = (p1 + p2 + c_est) / 3.0;
384        let n_tri = (p2 - p1).cross(c_est - p1);
385        let a_tri = n_tri.mag();
386        sum_a += a_tri;
387        sum_n += n_tri;
388        sum_ac += mid * a_tri;
389    }
390    let centre = if sum_a > f64::EPSILON {
391        sum_ac / sum_a
392    } else {
393        c_est
394    };
395    (centre, sum_n * 0.5)
396}
397
398// ── PatchKind <-> OpenFOAM type strings ─────────────────────────────────────
399
400fn kind_from_str(s: &str) -> PatchKind {
401    match s {
402        "wall" => PatchKind::Wall,
403        "symmetry" | "symmetryPlane" => PatchKind::Symmetry,
404        "empty" => PatchKind::Empty,
405        "wedge" => PatchKind::Wedge,
406        "cyclic" => PatchKind::Cyclic,
407        "processor" => PatchKind::Processor,
408        _ => PatchKind::Patch,
409    }
410}
411
412fn kind_to_str(k: PatchKind) -> &'static str {
413    match k {
414        PatchKind::Patch => "patch",
415        PatchKind::Wall => "wall",
416        PatchKind::Symmetry => "symmetry",
417        PatchKind::Empty => "empty",
418        PatchKind::Wedge => "wedge",
419        PatchKind::Cyclic => "cyclic",
420        PatchKind::Processor => "processor",
421    }
422}
423
424// ── Low-level token cursor over a polyMesh list file ────────────────────────
425
426/// Owns its tokens (no lifetime parameters, per crate design rules).
427struct ListCursor {
428    tokens: Vec<Token>,
429    pos: usize,
430    context: String,
431}
432
433impl ListCursor {
434    fn new(text: &str, context: impl Into<String>) -> Self {
435        let mut c = Self {
436            tokens: tokenize(text),
437            pos: 0,
438            context: context.into(),
439        };
440        c.skip_header();
441        c
442    }
443
444    fn err(&self, m: impl Into<String>) -> IoError {
445        IoError::parse(self.context.clone(), m.into())
446    }
447
448    fn peek(&self) -> Option<&str> {
449        self.tokens.get(self.pos).map(|t| t.text.as_str())
450    }
451
452    fn next(&mut self) -> Option<String> {
453        let t = self.tokens.get(self.pos).map(|t| t.text.clone());
454        if t.is_some() {
455            self.pos += 1;
456        }
457        t
458    }
459
460    fn expect(&mut self, tok: &str) -> Result<(), IoError> {
461        match self.next() {
462            Some(t) if t == tok => Ok(()),
463            other => Err(self.err(format!("expected `{tok}`, found {other:?}"))),
464        }
465    }
466
467    /// Skip a leading `FoamFile { … }` header block if present.
468    fn skip_header(&mut self) {
469        if self.peek() == Some("FoamFile") {
470            self.next(); // FoamFile
471            // consume up to and including the matching '}'
472            let mut depth = 0usize;
473            while let Some(t) = self.next() {
474                match t.as_str() {
475                    "{" => depth += 1,
476                    "}" => {
477                        depth -= 1;
478                        if depth == 0 {
479                            break;
480                        }
481                    }
482                    _ => {}
483                }
484            }
485        }
486    }
487
488    fn parse_usize(&mut self) -> Result<usize, IoError> {
489        let t = self
490            .next()
491            .ok_or_else(|| self.err("expected an integer, found end of input"))?;
492        t.parse::<usize>()
493            .map_err(|_| self.err(format!("expected an integer, found `{t}`")))
494    }
495
496    fn parse_f64(&mut self) -> Result<f64, IoError> {
497        let t = self
498            .next()
499            .ok_or_else(|| self.err("expected a number, found end of input"))?;
500        t.parse::<f64>()
501            .map_err(|_| self.err(format!("expected a number, found `{t}`")))
502    }
503}
504
505fn read_points(path: &Path) -> Result<Vec<Vector3>, IoError> {
506    let text = read_to_string(path)?;
507    let mut c = ListCursor::new(&text, path.display().to_string());
508    let count = c.parse_usize()?;
509    c.expect("(")?;
510    let mut pts = Vec::with_capacity(count);
511    for _ in 0..count {
512        c.expect("(")?;
513        let x = c.parse_f64()?;
514        let y = c.parse_f64()?;
515        let z = c.parse_f64()?;
516        c.expect(")")?;
517        pts.push(Vector3::new(x, y, z));
518    }
519    c.expect(")")?;
520    Ok(pts)
521}
522
523fn read_faces(path: &Path) -> Result<Vec<Vec<usize>>, IoError> {
524    let text = read_to_string(path)?;
525    let mut c = ListCursor::new(&text, path.display().to_string());
526    let count = c.parse_usize()?;
527    c.expect("(")?;
528    let mut faces = Vec::with_capacity(count);
529    for _ in 0..count {
530        let nverts = c.parse_usize()?;
531        c.expect("(")?;
532        let mut verts = Vec::with_capacity(nverts);
533        for _ in 0..nverts {
534            verts.push(c.parse_usize()?);
535        }
536        c.expect(")")?;
537        faces.push(verts);
538    }
539    c.expect(")")?;
540    Ok(faces)
541}
542
543fn read_labels(path: &Path) -> Result<Vec<usize>, IoError> {
544    let text = read_to_string(path)?;
545    let mut c = ListCursor::new(&text, path.display().to_string());
546    let count = c.parse_usize()?;
547    c.expect("(")?;
548    let mut out = Vec::with_capacity(count);
549    for _ in 0..count {
550        out.push(c.parse_usize()?);
551    }
552    c.expect(")")?;
553    Ok(out)
554}
555
556struct RawPatch {
557    name: String,
558    kind: String,
559    n_faces: usize,
560    start_face: usize,
561}
562
563fn read_boundary(path: &Path) -> Result<Vec<RawPatch>, IoError> {
564    let text = read_to_string(path)?;
565    let mut c = ListCursor::new(&text, path.display().to_string());
566    let count = c.parse_usize()?;
567    c.expect("(")?;
568    let mut patches = Vec::with_capacity(count);
569    for _ in 0..count {
570        let name = c
571            .next()
572            .ok_or_else(|| c.err("expected a patch name"))?
573            .to_string();
574        c.expect("{")?;
575        let mut kind = String::from("patch");
576        let mut n_faces = 0usize;
577        let mut start_face = 0usize;
578        // Parse `key value ;` entries; recognise type/nFaces/startFace, skip
579        // the rest (inGroups, physicalType, …), tolerating nested `( … )`.
580        while c.peek() != Some("}") {
581            let key = c
582                .next()
583                .ok_or_else(|| c.err("unexpected end of input in patch entry"))?
584                .to_string();
585            match key.as_str() {
586                "type" => {
587                    kind = c
588                        .next()
589                        .ok_or_else(|| c.err("expected a patch type"))?
590                        .to_string();
591                    c.expect(";")?;
592                }
593                "nFaces" => {
594                    n_faces = c.parse_usize()?;
595                    c.expect(";")?;
596                }
597                "startFace" => {
598                    start_face = c.parse_usize()?;
599                    c.expect(";")?;
600                }
601                _ => {
602                    // Skip an arbitrary value up to ';', balancing parens.
603                    let mut depth = 0usize;
604                    loop {
605                        match c.next().as_deref() {
606                            Some("(") => depth += 1,
607                            Some(")") => depth = depth.saturating_sub(1),
608                            Some(";") if depth == 0 => break,
609                            Some(_) => {}
610                            None => return Err(c.err("unterminated patch entry")),
611                        }
612                    }
613                }
614            }
615        }
616        c.expect("}")?;
617        patches.push(RawPatch {
618            name,
619            kind,
620            n_faces,
621            start_face,
622        });
623    }
624    c.expect(")")?;
625    Ok(patches)
626}
627
628// ── small file helpers ──────────────────────────────────────────────────────
629
630fn read_to_string(path: &Path) -> Result<String, IoError> {
631    std::fs::read_to_string(path).map_err(|e| IoError::io(path.display().to_string(), e))
632}
633
634fn write_file(path: &Path, contents: &str) -> Result<(), IoError> {
635    std::fs::write(path, contents).map_err(|e| IoError::io(path.display().to_string(), e))
636}
637
638/// Build a `constant/polyMesh` file header string (banner + FoamFile + separator).
639fn header_str(class: &str, object: &str) -> String {
640    let mut s = String::new();
641    s.push_str(BANNER);
642    let header = super::dict::FoamHeader::standard_with_location(class, "constant/polyMesh", object);
643    write_header(&mut s, &header);
644    s.push_str(SEPARATOR);
645    s.push_str("\n\n");
646    s
647}
outram_foam_basic_lib/io/poly_mesh.rs

outram_foam_basic_lib/io/
poly_mesh.rs
