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

outram_foam_basic_lib/io/
poly_mesh.rs
