Skip to main content

oxiphysics_io/vtk/
types.rs

1//! Auto-generated module
2//!
3//! 🤖 Generated with [SplitRS](https://github.com/cool-japan/splitrs)
4
5use crate::Result;
6use oxiphysics_core::math::Vec3;
7use std::fs::File;
8use std::io::{BufWriter, Write};
9use std::path::Path;
10
11/// VTK XML Unstructured Grid.
12///
13/// Build up a grid with [`VtuGrid::add_point`] / [`VtuGrid::add_cell`], attach
14/// field data via the `add_point_*` / `add_cell_*` helpers, then call
15/// [`VtuGrid::to_vtu_string`] to obtain a ready-to-write VTU XML string.
16pub struct VtuGrid {
17    /// 3-D coordinates of every point.
18    pub points: Vec<[f64; 3]>,
19    /// Connectivity list: each entry is the ordered point indices of one cell.
20    pub cells: Vec<Vec<usize>>,
21    /// VTK cell type for each cell (parallel to [`VtuGrid::cells`]).
22    pub cell_types: Vec<VtkCellType>,
23    /// Per-point data arrays.
24    pub point_data: Vec<VtkDataArray>,
25    /// Per-cell data arrays.
26    pub cell_data: Vec<VtkDataArray>,
27}
28impl VtuGrid {
29    /// Create an empty grid.
30    pub fn new() -> Self {
31        Self {
32            points: Vec::new(),
33            cells: Vec::new(),
34            cell_types: Vec::new(),
35            point_data: Vec::new(),
36            cell_data: Vec::new(),
37        }
38    }
39    /// Append a point and return its zero-based index.
40    pub fn add_point(&mut self, p: [f64; 3]) -> usize {
41        let idx = self.points.len();
42        self.points.push(p);
43        idx
44    }
45    /// Append a cell defined by `connectivity` (point indices) and a `cell_type`.
46    pub fn add_cell(&mut self, connectivity: Vec<usize>, cell_type: VtkCellType) {
47        self.cells.push(connectivity);
48        self.cell_types.push(cell_type);
49    }
50    /// Attach a per-point scalar field.
51    pub fn add_point_scalar(&mut self, name: &str, values: Vec<f64>) {
52        self.point_data.push(VtkDataArray::Scalar {
53            name: name.to_owned(),
54            values,
55        });
56    }
57    /// Attach a per-point 3-component vector field.
58    pub fn add_point_vector(&mut self, name: &str, values: Vec<[f64; 3]>) {
59        self.point_data.push(VtkDataArray::Vector3 {
60            name: name.to_owned(),
61            values,
62        });
63    }
64    /// Attach a per-cell scalar field.
65    pub fn add_cell_scalar(&mut self, name: &str, values: Vec<f64>) {
66        self.cell_data.push(VtkDataArray::Scalar {
67            name: name.to_owned(),
68            values,
69        });
70    }
71    /// Number of points in the grid.
72    pub fn n_points(&self) -> usize {
73        self.points.len()
74    }
75    /// Number of cells in the grid.
76    pub fn n_cells(&self) -> usize {
77        self.cells.len()
78    }
79    /// Serialise the grid to a VTU XML string.
80    ///
81    /// The output is ASCII-encoded and compatible with ParaView / VisIt.
82    pub fn to_vtu_string(&self) -> String {
83        let mut s = String::new();
84        s.push_str("<?xml version=\"1.0\"?>\n");
85        s.push_str(
86            "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">\n",
87        );
88        s.push_str("  <UnstructuredGrid>\n");
89        s.push_str(&format!(
90            "    <Piece NumberOfPoints=\"{}\" NumberOfCells=\"{}\">\n",
91            self.n_points(),
92            self.n_cells()
93        ));
94        s.push_str("      <Points>\n");
95        s.push_str(
96            "        <DataArray type=\"Float64\" NumberOfComponents=\"3\" format=\"ascii\">\n",
97        );
98        for p in &self.points {
99            s.push_str(&format!("          {} {} {}\n", p[0], p[1], p[2]));
100        }
101        s.push_str("        </DataArray>\n");
102        s.push_str("      </Points>\n");
103        s.push_str("      <Cells>\n");
104        s.push_str("        <DataArray type=\"Int64\" Name=\"connectivity\" format=\"ascii\">\n");
105        s.push_str("          ");
106        let mut first = true;
107        for conn in &self.cells {
108            for &idx in conn {
109                if !first {
110                    s.push(' ');
111                }
112                s.push_str(&idx.to_string());
113                first = false;
114            }
115        }
116        s.push('\n');
117        s.push_str("        </DataArray>\n");
118        s.push_str("        <DataArray type=\"Int64\" Name=\"offsets\" format=\"ascii\">\n");
119        s.push_str("          ");
120        let mut offset: usize = 0;
121        for (i, conn) in self.cells.iter().enumerate() {
122            if i > 0 {
123                s.push(' ');
124            }
125            offset += conn.len();
126            s.push_str(&offset.to_string());
127        }
128        s.push('\n');
129        s.push_str("        </DataArray>\n");
130        s.push_str("        <DataArray type=\"UInt8\" Name=\"types\" format=\"ascii\">\n");
131        s.push_str("          ");
132        for (i, ct) in self.cell_types.iter().enumerate() {
133            if i > 0 {
134                s.push(' ');
135            }
136            s.push_str(&(*ct as u8).to_string());
137        }
138        s.push('\n');
139        s.push_str("        </DataArray>\n");
140        s.push_str("      </Cells>\n");
141        if !self.point_data.is_empty() {
142            s.push_str("      <PointData>\n");
143            for arr in &self.point_data {
144                s.push_str(&Self::data_array_xml(arr));
145            }
146            s.push_str("      </PointData>\n");
147        }
148        if !self.cell_data.is_empty() {
149            s.push_str("      <CellData>\n");
150            for arr in &self.cell_data {
151                s.push_str(&Self::data_array_xml(arr));
152            }
153            s.push_str("      </CellData>\n");
154        }
155        s.push_str("    </Piece>\n");
156        s.push_str("  </UnstructuredGrid>\n");
157        s.push_str("</VTKFile>\n");
158        s
159    }
160    /// Render a [`VtkDataArray`] as a ``DataArray` XML element.
161    fn data_array_xml(arr: &VtkDataArray) -> String {
162        let mut s = String::new();
163        match arr {
164            VtkDataArray::Scalar { name, values } => {
165                s.push_str(
166                    &format!(
167                        "        <DataArray type=\"Float64\" Name=\"{}\" NumberOfComponents=\"1\" format=\"ascii\">\n",
168                        name
169                    ),
170                );
171                s.push_str("          ");
172                for (i, v) in values.iter().enumerate() {
173                    if i > 0 {
174                        s.push(' ');
175                    }
176                    s.push_str(&v.to_string());
177                }
178                s.push('\n');
179                s.push_str("        </DataArray>\n");
180            }
181            VtkDataArray::Vector3 { name, values } => {
182                s.push_str(
183                    &format!(
184                        "        <DataArray type=\"Float64\" Name=\"{}\" NumberOfComponents=\"3\" format=\"ascii\">\n",
185                        name
186                    ),
187                );
188                for v in values {
189                    s.push_str(&format!("          {} {} {}\n", v[0], v[1], v[2]));
190                }
191                s.push_str("        </DataArray>\n");
192            }
193            VtkDataArray::Integer { name, values } => {
194                s.push_str(
195                    &format!(
196                        "        <DataArray type=\"Int64\" Name=\"{}\" NumberOfComponents=\"1\" format=\"ascii\">\n",
197                        name
198                    ),
199                );
200                s.push_str("          ");
201                for (i, v) in values.iter().enumerate() {
202                    if i > 0 {
203                        s.push(' ');
204                    }
205                    s.push_str(&v.to_string());
206                }
207                s.push('\n');
208                s.push_str("        </DataArray>\n");
209            }
210        }
211        s
212    }
213    /// Write a PVD collection file that references a series of VTU snapshots.
214    ///
215    /// Returns the PVD XML as a `String`. The caller is responsible for writing
216    /// it to disk and ensuring that the referenced VTU files are co-located.
217    ///
218    /// # Arguments
219    /// * `base_name` – e.g. `"simulation"` (currently used as the collection title).
220    /// * `time_steps` – slice of `(time, vtu_filename)` pairs.
221    pub fn write_pvd_collection(base_name: &str, time_steps: &[(f64, String)]) -> String {
222        let mut s = String::new();
223        s.push_str("<?xml version=\"1.0\"?>\n");
224        s.push_str(
225            &format!(
226                "<VTKFile type=\"Collection\" version=\"0.1\" byte_order=\"LittleEndian\">\n  <!-- {} -->\n",
227                base_name
228            ),
229        );
230        s.push_str("  <Collection>\n");
231        for (time, filename) in time_steps {
232            s.push_str(&format!(
233                "    <DataSet timestep=\"{}\" group=\"\" part=\"0\" file=\"{}\"/>\n",
234                time, filename
235            ));
236        }
237        s.push_str("  </Collection>\n");
238        s.push_str("</VTKFile>\n");
239        s
240    }
241    /// Create a point-cloud grid: each position becomes a `Vertex` cell.
242    pub fn from_points(positions: &[[f64; 3]]) -> Self {
243        let mut grid = Self::new();
244        for &p in positions {
245            let idx = grid.add_point(p);
246            grid.add_cell(vec![idx], VtkCellType::Vertex);
247        }
248        grid
249    }
250    /// Create a grid from a tetrahedral mesh.
251    ///
252    /// `nodes` are the 3-D coordinates; `elements` contains 4-node connectivity
253    /// arrays, one per tetrahedron.
254    pub fn from_tet_mesh(nodes: &[[f64; 3]], elements: &[[usize; 4]]) -> Self {
255        let mut grid = Self::new();
256        for &n in nodes {
257            grid.add_point(n);
258        }
259        for &[a, b, c, d] in elements {
260            grid.add_cell(vec![a, b, c, d], VtkCellType::Tetra);
261        }
262        grid
263    }
264}
265/// A key-value pair in a VTK field data section.
266#[derive(Debug, Clone)]
267pub struct VtkFieldRecord {
268    /// Field name.
269    pub name: String,
270    /// Scalar values.
271    pub values: Vec<f64>,
272}
273impl VtkFieldRecord {
274    /// Create a field record.
275    pub fn new(name: impl Into<String>, values: Vec<f64>) -> Self {
276        Self {
277            name: name.into(),
278            values,
279        }
280    }
281    /// Number of values.
282    pub fn len(&self) -> usize {
283        self.values.len()
284    }
285    /// Returns true if the record has no values.
286    pub fn is_empty(&self) -> bool {
287        self.values.is_empty()
288    }
289}
290/// VTK cell types for use with [`VtuGrid`].
291#[derive(Debug, Clone, Copy)]
292pub enum VtkCellType {
293    /// VTK_VERTEX (1)
294    Vertex = 1,
295    /// VTK_LINE (3)
296    Line = 3,
297    /// VTK_TRIANGLE (5)
298    Triangle = 5,
299    /// VTK_QUAD (9)
300    Quad = 9,
301    /// VTK_TETRA (10)
302    Tetra = 10,
303    /// VTK_HEXAHEDRON (12)
304    Hexahedron = 12,
305    /// VTK_WEDGE (13)
306    Wedge = 13,
307    /// VTK_PYRAMID (14)
308    Pyramid = 14,
309}
310/// A parsed VTK legacy ASCII dataset.
311#[derive(Debug, Clone)]
312pub struct VtkLegacyData {
313    /// Dataset title line.
314    pub title: String,
315    /// Dataset type (e.g. "POLYDATA", "UNSTRUCTURED_GRID").
316    pub dataset_type: String,
317    /// Point coordinates (x, y, z).
318    pub points: Vec<[f64; 3]>,
319    /// Named scalar point data arrays.
320    pub point_scalars: Vec<(String, Vec<f64>)>,
321}
322impl VtkLegacyData {
323    /// Create an empty VTK legacy data container.
324    pub fn empty() -> Self {
325        Self {
326            title: String::new(),
327            dataset_type: String::new(),
328            points: Vec::new(),
329            point_scalars: Vec::new(),
330        }
331    }
332}
333/// A collection of VTK field data records.
334///
335/// Field data is global (not attached to points or cells) and is useful for
336/// simulation metadata such as time, timestep index, and solver parameters.
337#[derive(Debug, Clone, Default)]
338pub struct VtkFieldData {
339    /// The field records in this collection.
340    pub records: Vec<VtkFieldRecord>,
341}
342impl VtkFieldData {
343    /// Create an empty field data collection.
344    pub fn new() -> Self {
345        Self {
346            records: Vec::new(),
347        }
348    }
349    /// Add a named field.
350    pub fn add(&mut self, name: impl Into<String>, values: Vec<f64>) {
351        self.records.push(VtkFieldRecord::new(name, values));
352    }
353    /// Add a single scalar value as a 1-element field.
354    pub fn add_scalar(&mut self, name: impl Into<String>, value: f64) {
355        self.add(name, vec![value]);
356    }
357    /// Look up a field by name.
358    pub fn get(&self, name: &str) -> Option<&VtkFieldRecord> {
359        self.records.iter().find(|r| r.name == name)
360    }
361    /// Serialize as a VTK ASCII `FIELD` section.
362    pub fn to_vtk_field_string(&self) -> String {
363        if self.records.is_empty() {
364            return String::new();
365        }
366        let mut s = String::new();
367        s.push_str(&format!("FIELD FieldData {}\n", self.records.len()));
368        for rec in &self.records {
369            s.push_str(&format!("{} {} 1 float\n", rec.name, rec.values.len()));
370            let vals: Vec<String> = rec.values.iter().map(|v| v.to_string()).collect();
371            s.push_str(&vals.join(" "));
372            s.push('\n');
373        }
374        s
375    }
376}
377/// Writer for legacy VTK (.vtk) files.
378pub struct VtkWriter;
379impl VtkWriter {
380    /// Write a point cloud to a legacy VTK file.
381    pub fn write_points(path: &str, positions: &[Vec3]) -> Result<()> {
382        let file = File::create(Path::new(path))?;
383        let mut w = BufWriter::new(file);
384        writeln!(w, "# vtk DataFile Version 3.0")?;
385        writeln!(w, "OxiPhysics point cloud")?;
386        writeln!(w, "ASCII")?;
387        writeln!(w, "DATASET POLYDATA")?;
388        writeln!(w, "POINTS {} float", positions.len())?;
389        for p in positions {
390            writeln!(w, "{} {} {}", p.x, p.y, p.z)?;
391        }
392        w.flush()?;
393        Ok(())
394    }
395    /// Write an unstructured grid with tetrahedral cells to a legacy VTK file.
396    ///
397    /// Optionally includes scalar and vector point data.
398    pub fn write_unstructured_grid(
399        path: &str,
400        positions: &[Vec3],
401        cells: &[[usize; 4]],
402        scalars: Option<(&str, &[f64])>,
403        vectors: Option<(&str, &[Vec3])>,
404    ) -> Result<()> {
405        let file = File::create(Path::new(path))?;
406        let mut w = BufWriter::new(file);
407        writeln!(w, "# vtk DataFile Version 3.0")?;
408        writeln!(w, "OxiPhysics unstructured grid")?;
409        writeln!(w, "ASCII")?;
410        writeln!(w, "DATASET UNSTRUCTURED_GRID")?;
411        writeln!(w, "POINTS {} float", positions.len())?;
412        for p in positions {
413            writeln!(w, "{} {} {}", p.x, p.y, p.z)?;
414        }
415        let ncells = cells.len();
416        let cell_size = ncells * 5;
417        writeln!(w, "CELLS {} {}", ncells, cell_size)?;
418        for c in cells {
419            writeln!(w, "4 {} {} {} {}", c[0], c[1], c[2], c[3])?;
420        }
421        writeln!(w, "CELL_TYPES {}", ncells)?;
422        for _ in 0..ncells {
423            writeln!(w, "10")?;
424        }
425        let has_data = scalars.is_some() || vectors.is_some();
426        if has_data {
427            writeln!(w, "POINT_DATA {}", positions.len())?;
428        }
429        if let Some((name, vals)) = scalars {
430            writeln!(w, "SCALARS {} float 1", name)?;
431            writeln!(w, "LOOKUP_TABLE default")?;
432            for v in vals {
433                writeln!(w, "{}", v)?;
434            }
435        }
436        if let Some((name, vecs)) = vectors {
437            writeln!(w, "VECTORS {} float", name)?;
438            for v in vecs {
439                writeln!(w, "{} {} {}", v.x, v.y, v.z)?;
440            }
441        }
442        w.flush()?;
443        Ok(())
444    }
445    /// Write polygon (triangle) surface data to a legacy VTK file.
446    pub fn write_polydata(path: &str, positions: &[Vec3], triangles: &[[usize; 3]]) -> Result<()> {
447        let file = File::create(Path::new(path))?;
448        let mut w = BufWriter::new(file);
449        writeln!(w, "# vtk DataFile Version 3.0")?;
450        writeln!(w, "OxiPhysics polydata")?;
451        writeln!(w, "ASCII")?;
452        writeln!(w, "DATASET POLYDATA")?;
453        writeln!(w, "POINTS {} float", positions.len())?;
454        for p in positions {
455            writeln!(w, "{} {} {}", p.x, p.y, p.z)?;
456        }
457        let ntri = triangles.len();
458        writeln!(w, "POLYGONS {} {}", ntri, ntri * 4)?;
459        for t in triangles {
460            writeln!(w, "3 {} {} {}", t[0], t[1], t[2])?;
461        }
462        w.flush()?;
463        Ok(())
464    }
465}
466/// A time series of `VtuGrid` snapshots.
467pub struct VtkTimeSeries {
468    /// The simulation time of each snapshot.
469    pub times: Vec<f64>,
470    /// The grid snapshot at each time.
471    pub grids: Vec<VtuGrid>,
472    /// Base file name (without extension) for output files.
473    pub base_name: String,
474}
475impl VtkTimeSeries {
476    /// Create an empty time series.
477    pub fn new(base_name: &str) -> Self {
478        Self {
479            times: Vec::new(),
480            grids: Vec::new(),
481            base_name: base_name.to_owned(),
482        }
483    }
484    /// Append a snapshot at time `t`.
485    pub fn push(&mut self, t: f64, grid: VtuGrid) {
486        self.times.push(t);
487        self.grids.push(grid);
488    }
489    /// Number of snapshots.
490    pub fn n_steps(&self) -> usize {
491        self.times.len()
492    }
493    /// Generate a PVD collection file referencing all snapshots.
494    ///
495    /// File names are `{base_name}_{step:05}.vtu`.
496    pub fn to_pvd_string(&self) -> String {
497        let entries: Vec<(f64, String)> = self
498            .times
499            .iter()
500            .enumerate()
501            .map(|(i, &t)| (t, format!("{}_{:05}.vtu", self.base_name, i)))
502            .collect();
503        VtuGrid::write_pvd_collection(&self.base_name, &entries)
504    }
505    /// Get the VTU XML string for snapshot `i`.
506    pub fn vtu_string(&self, i: usize) -> Option<String> {
507        self.grids.get(i).map(|g| g.to_vtu_string())
508    }
509    /// Return estimated total data size in bytes (rough: 30 chars per point).
510    pub fn estimated_size_bytes(&self) -> usize {
511        self.grids.iter().map(|g| g.n_points() * 30).sum()
512    }
513}
514/// A single time step entry in a PVD collection.
515#[derive(Debug, Clone)]
516pub struct PvdEntry {
517    /// Simulation time.
518    pub time: f64,
519    /// Path to the VTU file for this time step.
520    pub filename: String,
521    /// Part name (optional, default "0").
522    pub part: u32,
523    /// Group name (optional).
524    pub group: String,
525}
526impl PvdEntry {
527    /// Create a new PVD entry.
528    pub fn new(time: f64, filename: impl Into<String>) -> Self {
529        Self {
530            time,
531            filename: filename.into(),
532            part: 0,
533            group: String::new(),
534        }
535    }
536}
537/// A ParaView-compatible multi-block dataset.
538///
539/// Wraps several [`VtuGrid`] blocks, each with a name.  Serializes to a
540/// `<vtkMultiBlockDataSet>` XML document.
541pub struct VtkMultiBlock {
542    /// The named blocks comprising this dataset.
543    pub blocks: Vec<VtkBlock>,
544    /// Optional dataset title.
545    pub title: String,
546}
547impl VtkMultiBlock {
548    /// Create an empty multi-block dataset.
549    pub fn new(title: impl Into<String>) -> Self {
550        Self {
551            blocks: Vec::new(),
552            title: title.into(),
553        }
554    }
555    /// Add a block.
556    pub fn add_block(&mut self, name: impl Into<String>, grid: VtuGrid) {
557        self.blocks.push(VtkBlock::new(name, grid));
558    }
559    /// Number of blocks.
560    pub fn n_blocks(&self) -> usize {
561        self.blocks.len()
562    }
563    /// Total number of points across all blocks.
564    pub fn total_points(&self) -> usize {
565        self.blocks.iter().map(|b| b.grid.n_points()).sum()
566    }
567    /// Total number of cells across all blocks.
568    pub fn total_cells(&self) -> usize {
569        self.blocks.iter().map(|b| b.grid.n_cells()).sum()
570    }
571    /// Serialize to a VTK multi-block VTM XML string.
572    ///
573    /// Each block references an external VTU file `{block_name}.vtu`.
574    pub fn to_vtm_string(&self) -> String {
575        let mut s = String::new();
576        s.push_str("<?xml version=\"1.0\"?>\n");
577        s.push_str(
578            "<VTKFile type=\"vtkMultiBlockDataSet\" version=\"1.0\" byte_order=\"LittleEndian\">\n",
579        );
580        s.push_str("  <vtkMultiBlockDataSet>\n");
581        for (i, block) in self.blocks.iter().enumerate() {
582            s.push_str(&format!(
583                "    <DataSet index=\"{}\" name=\"{}\" file=\"{}.vtu\"/>\n",
584                i, block.name, block.name
585            ));
586        }
587        s.push_str("  </vtkMultiBlockDataSet>\n");
588        s.push_str("</VTKFile>\n");
589        s
590    }
591    /// Write all VTU block files and a VTM index to the given directory.
592    ///
593    /// Returns the list of written file paths.
594    pub fn write_to_dir(&self, dir: &str) -> crate::Result<Vec<String>> {
595        use std::io::Write;
596        let mut written = Vec::new();
597        for block in &self.blocks {
598            let path = format!("{}/{}.vtu", dir, block.name);
599            let xml = block.grid.to_vtu_string();
600            let mut f = std::fs::File::create(&path)?;
601            f.write_all(xml.as_bytes())?;
602            written.push(path);
603        }
604        let vtm_path = format!("{}/{}.vtm", dir, self.title);
605        let vtm = self.to_vtm_string();
606        let mut f = std::fs::File::create(&vtm_path)?;
607        f.write_all(vtm.as_bytes())?;
608        written.push(vtm_path);
609        Ok(written)
610    }
611}
612/// A VTK polydata dataset storing points, lines, and polygons.
613pub struct VtkPolyDataGrid {
614    /// 3-D point coordinates.
615    pub points: Vec<[f64; 3]>,
616    /// Line connectivity.
617    pub lines: Vec<[usize; 2]>,
618    /// Triangle connectivity.
619    pub triangles: Vec<[usize; 3]>,
620    /// Per-point data arrays.
621    pub point_data: Vec<VtkDataArray>,
622    /// Per-cell (poly) data arrays.
623    pub cell_data: Vec<VtkDataArray>,
624}
625impl VtkPolyDataGrid {
626    /// Create an empty polydata grid.
627    pub fn new() -> Self {
628        Self {
629            points: Vec::new(),
630            lines: Vec::new(),
631            triangles: Vec::new(),
632            point_data: Vec::new(),
633            cell_data: Vec::new(),
634        }
635    }
636    /// Add a point and return its index.
637    pub fn add_point(&mut self, p: [f64; 3]) -> usize {
638        let idx = self.points.len();
639        self.points.push(p);
640        idx
641    }
642    /// Add a triangle.
643    pub fn add_triangle(&mut self, a: usize, b: usize, c: usize) {
644        self.triangles.push([a, b, c]);
645    }
646    /// Add a line.
647    pub fn add_line(&mut self, a: usize, b: usize) {
648        self.lines.push([a, b]);
649    }
650    /// Number of points.
651    pub fn n_points(&self) -> usize {
652        self.points.len()
653    }
654    /// Number of triangles.
655    pub fn n_triangles(&self) -> usize {
656        self.triangles.len()
657    }
658    /// Compute normals for each triangle (cross product, normalised).
659    pub fn compute_triangle_normals(&self) -> Vec<[f64; 3]> {
660        self.triangles
661            .iter()
662            .map(|&[a, b, c]| {
663                let pa = self.points[a];
664                let pb = self.points[b];
665                let pc = self.points[c];
666                let ab = [pb[0] - pa[0], pb[1] - pa[1], pb[2] - pa[2]];
667                let ac = [pc[0] - pa[0], pc[1] - pa[1], pc[2] - pa[2]];
668                let n = [
669                    ab[1] * ac[2] - ab[2] * ac[1],
670                    ab[2] * ac[0] - ab[0] * ac[2],
671                    ab[0] * ac[1] - ab[1] * ac[0],
672                ];
673                let len = (n[0] * n[0] + n[1] * n[1] + n[2] * n[2]).sqrt();
674                if len < 1e-12 {
675                    [0.0, 0.0, 1.0]
676                } else {
677                    [n[0] / len, n[1] / len, n[2] / len]
678                }
679            })
680            .collect()
681    }
682    /// Convert to a VTU unstructured grid (triangles become Triangle cells).
683    pub fn to_vtu_grid(&self) -> VtuGrid {
684        let mut g = VtuGrid::new();
685        for &p in &self.points {
686            g.add_point(p);
687        }
688        for &[a, b, c] in &self.triangles {
689            g.add_cell(vec![a, b, c], VtkCellType::Triangle);
690        }
691        for &[a, b] in &self.lines {
692            g.add_cell(vec![a, b], VtkCellType::Line);
693        }
694        g
695    }
696}
697/// A scalar, vector, or integer data array attached to points or cells.
698#[derive(Debug, Clone)]
699pub enum VtkDataArray {
700    /// Named scalar (1-component) field.
701    Scalar {
702        /// Field name.
703        name: String,
704        /// One value per point or cell.
705        values: Vec<f64>,
706    },
707    /// Named 3-component vector field.
708    Vector3 {
709        /// Field name.
710        name: String,
711        /// One `\[x, y, z\]` tuple per point or cell.
712        values: Vec<[f64; 3]>,
713    },
714    /// Named integer field.
715    Integer {
716        /// Field name.
717        name: String,
718        /// One value per point or cell.
719        values: Vec<i64>,
720    },
721}
722impl VtkDataArray {
723    /// Returns the name of the data array.
724    pub fn name(&self) -> &str {
725        match self {
726            Self::Scalar { name, .. } => name,
727            Self::Vector3 { name, .. } => name,
728            Self::Integer { name, .. } => name,
729        }
730    }
731    /// Returns the number of components per tuple (1 for scalar/integer, 3 for vector).
732    pub fn n_components(&self) -> usize {
733        match self {
734            Self::Scalar { .. } => 1,
735            Self::Vector3 { .. } => 3,
736            Self::Integer { .. } => 1,
737        }
738    }
739    /// Returns the number of tuples in the array.
740    pub fn len(&self) -> usize {
741        match self {
742            Self::Scalar { values, .. } => values.len(),
743            Self::Vector3 { values, .. } => values.len(),
744            Self::Integer { values, .. } => values.len(),
745        }
746    }
747    /// Returns `true` if the array contains no tuples.
748    pub fn is_empty(&self) -> bool {
749        self.len() == 0
750    }
751}
752/// A named block in a multi-block dataset.
753pub struct VtkBlock {
754    /// Block name.
755    pub name: String,
756    /// The grid for this block.
757    pub grid: VtuGrid,
758}
759impl VtkBlock {
760    /// Create a named block from a grid.
761    pub fn new(name: impl Into<String>, grid: VtuGrid) -> Self {
762        Self {
763            name: name.into(),
764            grid,
765        }
766    }
767}
768/// A VTK rectilinear grid (`.vtr`): axes are independently sampled.
769pub struct VtkRectilinearGrid {
770    /// Coordinate values along the X axis.
771    pub x_coords: Vec<f64>,
772    /// Coordinate values along the Y axis.
773    pub y_coords: Vec<f64>,
774    /// Coordinate values along the Z axis.
775    pub z_coords: Vec<f64>,
776    /// Per-point data arrays (ordered i,j,k like structured grid).
777    pub point_data: Vec<VtkDataArray>,
778}
779impl VtkRectilinearGrid {
780    /// Create from independent coordinate arrays.
781    pub fn new(x_coords: Vec<f64>, y_coords: Vec<f64>, z_coords: Vec<f64>) -> Self {
782        Self {
783            x_coords,
784            y_coords,
785            z_coords,
786            point_data: Vec::new(),
787        }
788    }
789    /// Total number of grid points.
790    pub fn n_points(&self) -> usize {
791        self.x_coords.len() * self.y_coords.len() * self.z_coords.len()
792    }
793    /// Grid dimensions `\[ni, nj, nk\]`.
794    pub fn dims(&self) -> [usize; 3] {
795        [
796            self.x_coords.len(),
797            self.y_coords.len(),
798            self.z_coords.len(),
799        ]
800    }
801    /// Append a per-point scalar field.
802    pub fn add_point_scalar(&mut self, name: &str, values: Vec<f64>) {
803        self.point_data.push(VtkDataArray::Scalar {
804            name: name.to_owned(),
805            values,
806        });
807    }
808    /// Serialise to VTK XML `.vtr` string.
809    pub fn to_vtr_string(&self) -> String {
810        let [ni, nj, nk] = self.dims();
811        let mut s = String::new();
812        s.push_str("<?xml version=\"1.0\"?>\n");
813        s.push_str(
814            "<VTKFile type=\"RectilinearGrid\" version=\"0.1\" byte_order=\"LittleEndian\">\n",
815        );
816        s.push_str(&format!(
817            "  <RectilinearGrid WholeExtent=\"0 {} 0 {} 0 {}\">\n",
818            ni.saturating_sub(1),
819            nj.saturating_sub(1),
820            nk.saturating_sub(1)
821        ));
822        s.push_str(&format!(
823            "    <Piece Extent=\"0 {} 0 {} 0 {}\">\n",
824            ni.saturating_sub(1),
825            nj.saturating_sub(1),
826            nk.saturating_sub(1)
827        ));
828        s.push_str("      <Coordinates>\n");
829        for (label, coords) in [
830            ("x", &self.x_coords),
831            ("y", &self.y_coords),
832            ("z", &self.z_coords),
833        ] {
834            s.push_str(&format!(
835                "        <DataArray type=\"Float64\" Name=\"{}\" format=\"ascii\">\n          ",
836                label
837            ));
838            for (i, v) in coords.iter().enumerate() {
839                if i > 0 {
840                    s.push(' ');
841                }
842                s.push_str(&v.to_string());
843            }
844            s.push_str("\n        </DataArray>\n");
845        }
846        s.push_str("      </Coordinates>\n");
847        if !self.point_data.is_empty() {
848            s.push_str("      <PointData>\n");
849            for arr in &self.point_data {
850                if let VtkDataArray::Scalar { name, values } = arr {
851                    s.push_str(
852                        &format!(
853                            "        <DataArray type=\"Float64\" Name=\"{}\" format=\"ascii\">\n          ",
854                            name
855                        ),
856                    );
857                    for (i, v) in values.iter().enumerate() {
858                        if i > 0 {
859                            s.push(' ');
860                        }
861                        s.push_str(&v.to_string());
862                    }
863                    s.push_str("\n        </DataArray>\n");
864                }
865            }
866            s.push_str("      </PointData>\n");
867        }
868        s.push_str("    </Piece>\n  </RectilinearGrid>\n</VTKFile>\n");
869        s
870    }
871}
872/// A VTK structured grid (`.vts`) with explicit point coordinates.
873///
874/// Dimensions are `(ni, nj, nk)` in index space.  Points are ordered
875/// with i varying fastest, then j, then k.
876pub struct VtkStructuredGrid {
877    /// Dimensions `\[ni, nj, nk\]`.
878    pub dims: [usize; 3],
879    /// 3-D coordinates of every point (ordered i,j,k).
880    pub points: Vec<[f64; 3]>,
881    /// Per-point data arrays.
882    pub point_data: Vec<VtkDataArray>,
883}
884impl VtkStructuredGrid {
885    /// Create an empty structured grid with the given dimensions.
886    pub fn new(ni: usize, nj: usize, nk: usize) -> Self {
887        Self {
888            dims: [ni, nj, nk],
889            points: Vec::with_capacity(ni * nj * nk),
890            point_data: Vec::new(),
891        }
892    }
893    /// Number of points (ni × nj × nk).
894    pub fn n_points(&self) -> usize {
895        self.dims[0] * self.dims[1] * self.dims[2]
896    }
897    /// Append a per-point scalar field.
898    pub fn add_point_scalar(&mut self, name: &str, values: Vec<f64>) {
899        self.point_data.push(VtkDataArray::Scalar {
900            name: name.to_owned(),
901            values,
902        });
903    }
904    /// Serialise to VTK XML `.vts` format string.
905    pub fn to_vts_string(&self) -> String {
906        let mut s = String::new();
907        s.push_str("<?xml version=\"1.0\"?>\n");
908        s.push_str(
909            "<VTKFile type=\"StructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">\n",
910        );
911        s.push_str(&format!(
912            "  <StructuredGrid WholeExtent=\"0 {} 0 {} 0 {}\">\n",
913            self.dims[0].saturating_sub(1),
914            self.dims[1].saturating_sub(1),
915            self.dims[2].saturating_sub(1)
916        ));
917        s.push_str(&format!(
918            "    <Piece Extent=\"0 {} 0 {} 0 {}\">\n",
919            self.dims[0].saturating_sub(1),
920            self.dims[1].saturating_sub(1),
921            self.dims[2].saturating_sub(1)
922        ));
923        s.push_str("      <Points>\n");
924        s.push_str(
925            "        <DataArray type=\"Float64\" NumberOfComponents=\"3\" format=\"ascii\">\n",
926        );
927        for p in &self.points {
928            s.push_str(&format!("          {} {} {}\n", p[0], p[1], p[2]));
929        }
930        s.push_str("        </DataArray>\n      </Points>\n");
931        if !self.point_data.is_empty() {
932            s.push_str("      <PointData>\n");
933            for arr in &self.point_data {
934                if let VtkDataArray::Scalar { name, values } = arr {
935                    s.push_str(
936                        &format!(
937                            "        <DataArray type=\"Float64\" Name=\"{}\" NumberOfComponents=\"1\" format=\"ascii\">\n          ",
938                            name
939                        ),
940                    );
941                    for (i, v) in values.iter().enumerate() {
942                        if i > 0 {
943                            s.push(' ');
944                        }
945                        s.push_str(&v.to_string());
946                    }
947                    s.push_str("\n        </DataArray>\n");
948                }
949            }
950            s.push_str("      </PointData>\n");
951        }
952        s.push_str("    </Piece>\n  </StructuredGrid>\n</VTKFile>\n");
953        s
954    }
955    /// Build a uniform Cartesian grid covering `\[x0,x1\] x \[y0,y1\] x \[z0,z1\]`.
956    pub fn uniform(
957        x0: f64,
958        x1: f64,
959        ni: usize,
960        y0: f64,
961        y1: f64,
962        nj: usize,
963        z0: f64,
964        z1: f64,
965        nk: usize,
966    ) -> Self {
967        let mut g = Self::new(ni, nj, nk);
968        let dx = if ni > 1 {
969            (x1 - x0) / (ni - 1) as f64
970        } else {
971            0.0
972        };
973        let dy = if nj > 1 {
974            (y1 - y0) / (nj - 1) as f64
975        } else {
976            0.0
977        };
978        let dz = if nk > 1 {
979            (z1 - z0) / (nk - 1) as f64
980        } else {
981            0.0
982        };
983        for k in 0..nk {
984            for j in 0..nj {
985                for i in 0..ni {
986                    g.points
987                        .push([x0 + i as f64 * dx, y0 + j as f64 * dy, z0 + k as f64 * dz]);
988                }
989            }
990        }
991        g
992    }
993}