Skip to main content

oxiphysics_io/vtk/
functions.rs

1//! Auto-generated module
2//!
3//! 🤖 Generated with [SplitRS](https://github.com/cool-japan/splitrs)
4
5use std::io::Write;
6
7use super::types::{PvdEntry, VtkCellType, VtkLegacyData, VtuGrid};
8
9#[cfg(test)]
10mod tests {
11
12    use crate::VtkWriter;
13
14    use oxiphysics_core::math::Vec3;
15    use std::path::Path;
16    #[test]
17    fn test_vtk_point_cloud_write() {
18        let path = "/tmp/oxiphy_test_points.vtk";
19        let pts = vec![
20            Vec3::new(0.0, 0.0, 0.0),
21            Vec3::new(1.0, 0.0, 0.0),
22            Vec3::new(0.0, 1.0, 0.0),
23        ];
24        VtkWriter::write_points(path, &pts).unwrap();
25        assert!(Path::new(path).exists());
26        let content = std::fs::read_to_string(path).unwrap();
27        assert!(content.contains("POINTS 3 float"));
28        std::fs::remove_file(path).ok();
29    }
30    #[test]
31    fn test_vtk_unstructured_grid_write() {
32        let path = "/tmp/oxiphy_test_ugrid.vtk";
33        let pts = vec![
34            Vec3::new(0.0, 0.0, 0.0),
35            Vec3::new(1.0, 0.0, 0.0),
36            Vec3::new(0.0, 1.0, 0.0),
37            Vec3::new(0.0, 0.0, 1.0),
38        ];
39        let cells = vec![[0, 1, 2, 3]];
40        let scalars = vec![1.0, 2.0, 3.0, 4.0];
41        VtkWriter::write_unstructured_grid(path, &pts, &cells, Some(("pressure", &scalars)), None)
42            .unwrap();
43        let content = std::fs::read_to_string(path).unwrap();
44        assert!(content.contains("UNSTRUCTURED_GRID"));
45        assert!(content.contains("CELL_TYPES 1"));
46        std::fs::remove_file(path).ok();
47    }
48    #[test]
49    fn test_vtk_header_format() {
50        let path = "/tmp/oxiphy_test_header.vtk";
51        let pts = vec![Vec3::new(0.0, 0.0, 0.0)];
52        VtkWriter::write_points(path, &pts).unwrap();
53        let content = std::fs::read_to_string(path).unwrap();
54        let lines: Vec<&str> = content.lines().collect();
55        assert_eq!(lines[0], "# vtk DataFile Version 3.0");
56        assert_eq!(lines[2], "ASCII");
57        std::fs::remove_file(path).ok();
58    }
59    #[test]
60    fn test_vtk_polydata_write() {
61        let path = "/tmp/oxiphy_test_polydata.vtk";
62        let pts = vec![
63            Vec3::new(0.0, 0.0, 0.0),
64            Vec3::new(1.0, 0.0, 0.0),
65            Vec3::new(0.0, 1.0, 0.0),
66        ];
67        let tris = vec![[0, 1, 2]];
68        VtkWriter::write_polydata(path, &pts, &tris).unwrap();
69        let content = std::fs::read_to_string(path).unwrap();
70        assert!(content.contains("POLYGONS 1 4"));
71        std::fs::remove_file(path).ok();
72    }
73    #[test]
74    fn test_vtk_write_read_roundtrip() {
75        let path = "/tmp/oxiphy_test_roundtrip.vtk";
76        let positions = vec![
77            Vec3::new(1.0, 2.0, 3.0),
78            Vec3::new(4.0, 5.0, 6.0),
79            Vec3::new(7.0, 8.0, 9.0),
80            Vec3::new(0.5, 1.5, 2.5),
81        ];
82        let scalars = vec![10.0, 20.0, 30.0, 40.0];
83        let cells: Vec<[usize; 4]> = vec![];
84        VtkWriter::write_unstructured_grid(
85            path,
86            &positions,
87            &cells,
88            Some(("density", &scalars)),
89            None,
90        )
91        .unwrap();
92        let content = std::fs::read_to_string(path).unwrap();
93        assert!(
94            content.contains("POINTS 4 float"),
95            "expected POINTS 4 float header"
96        );
97        assert!(
98            content.contains("SCALARS density float 1"),
99            "expected scalars section"
100        );
101        let mut parsed: Vec<[f64; 3]> = Vec::new();
102        let mut in_points = false;
103        for line in content.lines() {
104            if line.starts_with("POINTS") {
105                in_points = true;
106                continue;
107            }
108            if in_points {
109                let parts: Vec<&str> = line.split_whitespace().collect();
110                if parts.len() == 3 {
111                    if let (Ok(x), Ok(y), Ok(z)) = (
112                        parts[0].parse::<f64>(),
113                        parts[1].parse::<f64>(),
114                        parts[2].parse::<f64>(),
115                    ) {
116                        parsed.push([x, y, z]);
117                    }
118                } else {
119                    in_points = false;
120                }
121            }
122        }
123        assert_eq!(parsed.len(), 4, "expected 4 parsed points");
124        let tol = 1e-4;
125        assert!((parsed[0][0] - 1.0).abs() < tol);
126        assert!((parsed[0][1] - 2.0).abs() < tol);
127        assert!((parsed[0][2] - 3.0).abs() < tol);
128        assert!((parsed[1][0] - 4.0).abs() < tol);
129        assert!((parsed[2][1] - 8.0).abs() < tol);
130        assert!((parsed[3][2] - 2.5).abs() < tol);
131        std::fs::remove_file(path).ok();
132    }
133}
134/// Validate a `VtuGrid` for consistency.
135///
136/// Checks that all cell connectivity indices are valid and all point data
137/// arrays have the correct length.
138#[allow(dead_code)]
139pub fn validate_vtu_grid(grid: &VtuGrid) -> Vec<String> {
140    let mut errors = Vec::new();
141    let n_pts = grid.n_points();
142    for (ci, conn) in grid.cells.iter().enumerate() {
143        for &idx in conn {
144            if idx >= n_pts {
145                errors.push(format!("Cell {ci}: point index {idx} >= n_points {n_pts}"));
146            }
147        }
148    }
149    for arr in &grid.point_data {
150        if arr.len() != n_pts {
151            errors.push(format!(
152                "Point data '{}': len {} != n_points {}",
153                arr.name(),
154                arr.len(),
155                n_pts
156            ));
157        }
158    }
159    let n_cells = grid.n_cells();
160    for arr in &grid.cell_data {
161        if arr.len() != n_cells {
162            errors.push(format!(
163                "Cell data '{}': len {} != n_cells {}",
164                arr.name(),
165                arr.len(),
166                n_cells
167            ));
168        }
169    }
170    errors
171}
172#[cfg(test)]
173mod vtu_grid_tests {
174    use super::*;
175
176    use crate::vtk::VtkPolyDataGrid;
177    use crate::vtk::VtkRectilinearGrid;
178    use crate::vtk::VtkStructuredGrid;
179    use crate::vtk::VtkTimeSeries;
180    use crate::vtk::types::*;
181
182    #[test]
183    fn test_vtu_empty_grid() {
184        let grid = VtuGrid::new();
185        let xml = grid.to_vtu_string();
186        assert!(
187            xml.starts_with("<?xml version=\"1.0\"?>"),
188            "missing XML header"
189        );
190        assert!(xml.contains("<VTKFile"), "missing VTKFile element");
191        assert!(xml.contains("UnstructuredGrid"), "missing UnstructuredGrid");
192    }
193    #[test]
194    fn test_vtu_single_point() {
195        let mut grid = VtuGrid::new();
196        grid.add_point([1.0, 2.0, 3.0]);
197        let xml = grid.to_vtu_string();
198        assert!(
199            xml.contains("NumberOfPoints=\"1\""),
200            "expected NumberOfPoints=\"1\""
201        );
202    }
203    #[test]
204    fn test_vtu_tetrahedron() {
205        let mut grid = VtuGrid::new();
206        grid.add_point([0.0, 0.0, 0.0]);
207        grid.add_point([1.0, 0.0, 0.0]);
208        grid.add_point([0.0, 1.0, 0.0]);
209        grid.add_point([0.0, 0.0, 1.0]);
210        grid.add_cell(vec![0, 1, 2, 3], VtkCellType::Tetra);
211        let xml = grid.to_vtu_string();
212        assert!(
213            xml.contains(">10<") || xml.contains("          10"),
214            "expected cell type 10 in output; got:\n{}",
215            xml
216        );
217    }
218    #[test]
219    fn test_vtu_point_data() {
220        let mut grid = VtuGrid::new();
221        grid.add_point([0.0, 0.0, 0.0]);
222        grid.add_point([1.0, 0.0, 0.0]);
223        grid.add_point([0.0, 1.0, 0.0]);
224        grid.add_cell(vec![0, 1, 2], VtkCellType::Triangle);
225        grid.add_point_scalar("pressure", vec![1.0, 2.0, 3.0]);
226        let xml = grid.to_vtu_string();
227        assert!(
228            xml.contains("Name=\"pressure\""),
229            "pressure field missing from PointData"
230        );
231        assert!(xml.contains("<PointData>"), "PointData section missing");
232    }
233    #[test]
234    fn test_vtu_from_points() {
235        let pts = [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]];
236        let grid = VtuGrid::from_points(&pts);
237        assert_eq!(grid.n_points(), 3);
238        assert_eq!(grid.n_cells(), 3);
239        let xml = grid.to_vtu_string();
240        assert!(xml.contains("NumberOfCells=\"3\""));
241    }
242    #[test]
243    fn test_vtu_from_tet_mesh() {
244        let nodes = [
245            [0.0, 0.0, 0.0],
246            [1.0, 0.0, 0.0],
247            [0.0, 1.0, 0.0],
248            [0.0, 0.0, 1.0],
249        ];
250        let elems = [[0, 1, 2, 3]];
251        let grid = VtuGrid::from_tet_mesh(&nodes, &elems);
252        assert_eq!(grid.n_points(), 4);
253        assert_eq!(grid.n_cells(), 1);
254        let xml = grid.to_vtu_string();
255        assert!(xml.contains("NumberOfCells=\"1\""));
256        assert!(xml.contains("10"), "expected cell type 10 for tet");
257    }
258    #[test]
259    fn test_pvd_collection() {
260        let steps = vec![
261            (0.0_f64, "sim_0000.vtu".to_owned()),
262            (1.0_f64, "sim_0001.vtu".to_owned()),
263        ];
264        let pvd = VtuGrid::write_pvd_collection("simulation", &steps);
265        assert!(
266            pvd.contains("timestep=\"0\"")
267                || pvd.contains("timestep=\"0.0\"")
268                || pvd.contains("timestep=\"0\""),
269            "missing first timestep entry"
270        );
271        let count = pvd.matches("<DataSet").count();
272        assert_eq!(count, 2, "expected 2 DataSet entries, got {}", count);
273    }
274    #[test]
275    fn test_vtu_offsets() {
276        let mut grid = VtuGrid::new();
277        for i in 0..8 {
278            grid.add_point([i as f64, 0.0, 0.0]);
279        }
280        grid.add_cell(vec![0, 1, 2, 3], VtkCellType::Tetra);
281        grid.add_cell(vec![4, 5, 6, 7], VtkCellType::Tetra);
282        let xml = grid.to_vtu_string();
283        assert!(
284            xml.contains("4 8"),
285            "expected offsets '4 8' in output;\n{}",
286            xml
287        );
288    }
289    #[test]
290    fn test_structured_grid_uniform_point_count() {
291        let g = VtkStructuredGrid::uniform(0.0, 1.0, 3, 0.0, 1.0, 3, 0.0, 1.0, 3);
292        assert_eq!(g.n_points(), 27);
293        assert_eq!(g.points.len(), 27);
294    }
295    #[test]
296    fn test_structured_grid_single_point() {
297        let g = VtkStructuredGrid::uniform(0.0, 0.0, 1, 0.0, 0.0, 1, 0.0, 0.0, 1);
298        assert_eq!(g.n_points(), 1);
299    }
300    #[test]
301    fn test_structured_grid_to_vts_contains_header() {
302        let g = VtkStructuredGrid::uniform(0.0, 1.0, 2, 0.0, 1.0, 2, 0.0, 1.0, 2);
303        let s = g.to_vts_string();
304        assert!(s.contains("<VTKFile type=\"StructuredGrid\""));
305        assert!(s.contains("<Piece"));
306    }
307    #[test]
308    fn test_structured_grid_with_scalar() {
309        let mut g = VtkStructuredGrid::uniform(0.0, 1.0, 2, 0.0, 1.0, 2, 0.0, 1.0, 2);
310        let vals: Vec<f64> = (0..g.n_points()).map(|i| i as f64).collect();
311        g.add_point_scalar("index", vals);
312        let s = g.to_vts_string();
313        assert!(s.contains("Name=\"index\""));
314    }
315    #[test]
316    fn test_rectilinear_grid_n_points() {
317        let g = VtkRectilinearGrid::new(vec![0.0, 0.5, 1.0], vec![0.0, 1.0], vec![0.0]);
318        assert_eq!(g.n_points(), 6);
319    }
320    #[test]
321    fn test_rectilinear_grid_dims() {
322        let g = VtkRectilinearGrid::new(vec![0.0, 1.0], vec![0.0, 1.0, 2.0], vec![0.0, 0.5, 1.0]);
323        assert_eq!(g.dims(), [2, 3, 3]);
324    }
325    #[test]
326    fn test_rectilinear_grid_vtr_string() {
327        let g = VtkRectilinearGrid::new(vec![0.0, 1.0], vec![0.0, 1.0], vec![0.0, 1.0]);
328        let s = g.to_vtr_string();
329        assert!(s.contains("<VTKFile type=\"RectilinearGrid\""));
330        assert!(s.contains("<Coordinates>"));
331        assert!(s.contains("Name=\"x\""));
332        assert!(s.contains("Name=\"y\""));
333        assert!(s.contains("Name=\"z\""));
334    }
335    #[test]
336    fn test_polydata_add_triangle() {
337        let mut g = VtkPolyDataGrid::new();
338        g.add_point([0.0, 0.0, 0.0]);
339        g.add_point([1.0, 0.0, 0.0]);
340        g.add_point([0.0, 1.0, 0.0]);
341        g.add_triangle(0, 1, 2);
342        assert_eq!(g.n_triangles(), 1);
343        assert_eq!(g.n_points(), 3);
344    }
345    #[test]
346    fn test_polydata_triangle_normals() {
347        let mut g = VtkPolyDataGrid::new();
348        g.add_point([0.0, 0.0, 0.0]);
349        g.add_point([1.0, 0.0, 0.0]);
350        g.add_point([0.0, 1.0, 0.0]);
351        g.add_triangle(0, 1, 2);
352        let normals = g.compute_triangle_normals();
353        assert_eq!(normals.len(), 1);
354        assert!(
355            (normals[0][2]).abs() > 0.9,
356            "Normal z-component should be ~1"
357        );
358    }
359    #[test]
360    fn test_polydata_to_vtu_grid() {
361        let mut g = VtkPolyDataGrid::new();
362        g.add_point([0.0, 0.0, 0.0]);
363        g.add_point([1.0, 0.0, 0.0]);
364        g.add_point([0.0, 1.0, 0.0]);
365        g.add_triangle(0, 1, 2);
366        let vtu = g.to_vtu_grid();
367        assert_eq!(vtu.n_points(), 3);
368        assert_eq!(vtu.n_cells(), 1);
369    }
370    #[test]
371    fn test_time_series_push_and_count() {
372        let mut ts = VtkTimeSeries::new("sim");
373        let g1 = VtuGrid::from_points(&[[0.0, 0.0, 0.0]]);
374        let g2 = VtuGrid::from_points(&[[1.0, 0.0, 0.0]]);
375        ts.push(0.0, g1);
376        ts.push(0.1, g2);
377        assert_eq!(ts.n_steps(), 2);
378    }
379    #[test]
380    fn test_time_series_pvd_string() {
381        let mut ts = VtkTimeSeries::new("run");
382        ts.push(0.0, VtuGrid::new());
383        ts.push(1.0, VtuGrid::new());
384        let pvd = ts.to_pvd_string();
385        assert!(
386            pvd.contains("<DataSet"),
387            "PVD should contain DataSet entries"
388        );
389        let count = pvd.matches("<DataSet").count();
390        assert_eq!(count, 2);
391    }
392    #[test]
393    fn test_time_series_vtu_string() {
394        let mut ts = VtkTimeSeries::new("s");
395        let mut g = VtuGrid::new();
396        g.add_point([1.0, 2.0, 3.0]);
397        ts.push(0.5, g);
398        let vtu = ts.vtu_string(0).expect("should have step 0");
399        assert!(
400            vtu.contains("1 2 3") || vtu.contains("1.0 2.0 3.0") || vtu.contains("1 2 3"),
401            "VTU should contain coordinates"
402        );
403    }
404    #[test]
405    fn test_time_series_none_out_of_range() {
406        let ts = VtkTimeSeries::new("s");
407        assert!(ts.vtu_string(0).is_none());
408    }
409    #[test]
410    fn test_validate_vtu_grid_valid() {
411        let mut grid = VtuGrid::new();
412        grid.add_point([0.0, 0.0, 0.0]);
413        grid.add_point([1.0, 0.0, 0.0]);
414        grid.add_point([0.0, 1.0, 0.0]);
415        grid.add_cell(vec![0, 1, 2], VtkCellType::Triangle);
416        grid.add_point_scalar("p", vec![1.0, 2.0, 3.0]);
417        let errors = validate_vtu_grid(&grid);
418        assert!(
419            errors.is_empty(),
420            "Valid grid should have no errors: {:?}",
421            errors
422        );
423    }
424    #[test]
425    fn test_validate_vtu_grid_bad_connectivity() {
426        let mut grid = VtuGrid::new();
427        grid.add_point([0.0, 0.0, 0.0]);
428        grid.cells.push(vec![0, 5]);
429        grid.cell_types.push(VtkCellType::Line);
430        let errors = validate_vtu_grid(&grid);
431        assert!(
432            !errors.is_empty(),
433            "Out-of-range index should produce an error"
434        );
435    }
436    #[test]
437    fn test_validate_vtu_grid_mismatched_point_data() {
438        let mut grid = VtuGrid::new();
439        grid.add_point([0.0, 0.0, 0.0]);
440        grid.add_point([1.0, 0.0, 0.0]);
441        grid.add_point_scalar("bad", vec![1.0]);
442        let errors = validate_vtu_grid(&grid);
443        assert!(
444            !errors.is_empty(),
445            "Mismatched point data should produce error"
446        );
447    }
448}
449/// Write a legacy VTK unstructured grid in binary (big-endian) format.
450///
451/// VTK legacy binary files use big-endian byte order.
452/// Supports arbitrary cell types and float64 point coordinates.
453#[allow(dead_code)]
454pub fn write_vtk_binary_unstructured(
455    path: &str,
456    points: &[[f64; 3]],
457    cells: &[Vec<usize>],
458    cell_types: &[u8],
459) -> crate::Result<()> {
460    let file = std::fs::File::create(path)?;
461    let mut w = std::io::BufWriter::new(file);
462    writeln!(w, "# vtk DataFile Version 3.0")?;
463    writeln!(w, "OxiPhysics binary unstructured grid")?;
464    writeln!(w, "BINARY")?;
465    writeln!(w, "DATASET UNSTRUCTURED_GRID")?;
466    writeln!(w, "POINTS {} double", points.len())?;
467    for p in points {
468        for &coord in p {
469            w.write_all(&coord.to_be_bytes())?;
470        }
471    }
472    writeln!(w)?;
473    let total_entries: usize = cells.iter().map(|c| 1 + c.len()).sum();
474    writeln!(w, "CELLS {} {}", cells.len(), total_entries)?;
475    for cell in cells {
476        let n = cell.len() as i32;
477        w.write_all(&n.to_be_bytes())?;
478        for &idx in cell {
479            let i = idx as i32;
480            w.write_all(&i.to_be_bytes())?;
481        }
482    }
483    writeln!(w)?;
484    writeln!(w, "CELL_TYPES {}", cell_types.len())?;
485    for &ct in cell_types {
486        let ct_i32 = ct as i32;
487        w.write_all(&ct_i32.to_be_bytes())?;
488    }
489    writeln!(w)?;
490    w.flush()?;
491    Ok(())
492}
493/// Parse a legacy VTK ASCII file, extracting points and point scalars.
494///
495/// This is a simplified reader that handles POINTS and POINT_DATA/SCALARS sections.
496#[allow(dead_code)]
497pub fn read_vtk_legacy_ascii(content: &str) -> std::result::Result<VtkLegacyData, String> {
498    let mut lines = content.lines();
499    let magic = lines.next().ok_or("Missing magic line")?;
500    if !magic.starts_with("# vtk DataFile") {
501        return Err(format!("Not a VTK file: {magic}"));
502    }
503    let title = lines.next().ok_or("Missing title")?.to_string();
504    let format = lines.next().ok_or("Missing format line")?;
505    if format.trim() != "ASCII" {
506        return Err(format!("Only ASCII format supported, got: {format}"));
507    }
508    let ds_line = lines.next().ok_or("Missing DATASET line")?;
509    let dataset_type = ds_line
510        .strip_prefix("DATASET ")
511        .ok_or("Expected DATASET line")?
512        .trim()
513        .to_string();
514    let mut result = VtkLegacyData {
515        title,
516        dataset_type,
517        points: Vec::new(),
518        point_scalars: Vec::new(),
519    };
520    let remaining: Vec<&str> = lines.collect();
521    let mut i = 0;
522    while i < remaining.len() {
523        let line = remaining[i].trim();
524        if line.starts_with("POINTS ") {
525            let parts: Vec<&str> = line.split_whitespace().collect();
526            let n: usize = parts.get(1).and_then(|s| s.parse().ok()).unwrap_or(0);
527            i += 1;
528            let mut pts = Vec::with_capacity(n);
529            let mut coords_read = 0;
530            let mut current = [0.0f64; 3];
531            while coords_read < n * 3 && i < remaining.len() {
532                for tok in remaining[i].split_whitespace() {
533                    if let Ok(v) = tok.parse::<f64>() {
534                        let idx = coords_read % 3;
535                        current[idx] = v;
536                        if idx == 2 {
537                            pts.push(current);
538                            current = [0.0; 3];
539                        }
540                        coords_read += 1;
541                    }
542                }
543                i += 1;
544            }
545            result.points = pts;
546            continue;
547        }
548        if line.starts_with("SCALARS ") {
549            let parts: Vec<&str> = line.split_whitespace().collect();
550            let name = parts.get(1).unwrap_or(&"scalar").to_string();
551            i += 1;
552            if i < remaining.len() && remaining[i].trim().starts_with("LOOKUP_TABLE") {
553                i += 1;
554            }
555            let n = result.points.len();
556            let mut vals = Vec::with_capacity(n);
557            while vals.len() < n && i < remaining.len() {
558                for tok in remaining[i].split_whitespace() {
559                    if let Ok(v) = tok.parse::<f64>() {
560                        vals.push(v);
561                    }
562                }
563                i += 1;
564            }
565            result.point_scalars.push((name, vals));
566            continue;
567        }
568        i += 1;
569    }
570    Ok(result)
571}
572/// Write a ParaView Data (PVD) collection file pointing to a sequence of VTU files.
573#[allow(dead_code)]
574pub fn write_pvd_file(path: &str, entries: &[PvdEntry]) -> crate::Result<()> {
575    let file = std::fs::File::create(path)?;
576    let mut w = std::io::BufWriter::new(file);
577    writeln!(w, r#"<?xml version="1.0"?>"#)?;
578    writeln!(
579        w,
580        r#"<VTKFile type="Collection" version="0.1" byte_order="LittleEndian">"#
581    )?;
582    writeln!(w, r#"  <Collection>"#)?;
583    for entry in entries {
584        writeln!(
585            w,
586            r#"    <DataSet timestep="{:.6e}" group="{}" part="{}" file="{}"/>"#,
587            entry.time, entry.group, entry.part, entry.filename
588        )?;
589    }
590    writeln!(w, r#"  </Collection>"#)?;
591    writeln!(w, r#"</VTKFile>"#)?;
592    w.flush()?;
593    Ok(())
594}
595/// Generate a sequence of PVD entries for evenly spaced time steps.
596#[allow(dead_code)]
597pub fn pvd_entries_uniform(
598    base_name: &str,
599    extension: &str,
600    t_start: f64,
601    dt: f64,
602    n_steps: usize,
603) -> Vec<PvdEntry> {
604    (0..n_steps)
605        .map(|i| {
606            let t = t_start + i as f64 * dt;
607            let fname = format!("{base_name}_{i:06}.{extension}");
608            PvdEntry::new(t, fname)
609        })
610        .collect()
611}
612/// Write a triangle surface mesh as a VTK POLYDATA file (ASCII legacy).
613#[allow(dead_code)]
614pub fn write_vtk_polydata(
615    path: &str,
616    vertices: &[[f64; 3]],
617    triangles: &[[usize; 3]],
618    normals: Option<&[[f64; 3]]>,
619) -> crate::Result<()> {
620    let file = std::fs::File::create(path)?;
621    let mut w = std::io::BufWriter::new(file);
622    writeln!(w, "# vtk DataFile Version 3.0")?;
623    writeln!(w, "OxiPhysics surface mesh")?;
624    writeln!(w, "ASCII")?;
625    writeln!(w, "DATASET POLYDATA")?;
626    writeln!(w, "POINTS {} double", vertices.len())?;
627    for v in vertices {
628        writeln!(w, "{:.10e} {:.10e} {:.10e}", v[0], v[1], v[2])?;
629    }
630    let n_tris = triangles.len();
631    writeln!(w, "POLYGONS {} {}", n_tris, n_tris * 4)?;
632    for tri in triangles {
633        writeln!(w, "3 {} {} {}", tri[0], tri[1], tri[2])?;
634    }
635    if let Some(nrms) = normals {
636        writeln!(w, "POINT_DATA {}", vertices.len())?;
637        writeln!(w, "NORMALS normals double")?;
638        for n in nrms {
639            writeln!(w, "{:.10e} {:.10e} {:.10e}", n[0], n[1], n[2])?;
640        }
641    }
642    w.flush()?;
643    Ok(())
644}
645/// Write VTK CELL_DATA with a symmetric 3×3 stress tensor per cell.
646///
647/// VTK represents tensors as 9-component arrays (row-major).
648#[allow(dead_code)]
649pub fn write_vtk_cell_stress(
650    path: &str,
651    points: &[[f64; 3]],
652    cells: &[Vec<usize>],
653    cell_types: &[u8],
654    stress_tensors: &[[[f64; 3]; 3]],
655) -> crate::Result<()> {
656    let file = std::fs::File::create(path)?;
657    let mut w = std::io::BufWriter::new(file);
658    writeln!(w, "# vtk DataFile Version 3.0")?;
659    writeln!(w, "OxiPhysics cell stress data")?;
660    writeln!(w, "ASCII")?;
661    writeln!(w, "DATASET UNSTRUCTURED_GRID")?;
662    writeln!(w, "POINTS {} double", points.len())?;
663    for p in points {
664        writeln!(w, "{:.10e} {:.10e} {:.10e}", p[0], p[1], p[2])?;
665    }
666    let total_entries: usize = cells.iter().map(|c| 1 + c.len()).sum();
667    writeln!(w, "CELLS {} {}", cells.len(), total_entries)?;
668    for cell in cells {
669        let indices: Vec<String> = cell.iter().map(|i| i.to_string()).collect();
670        writeln!(w, "{} {}", cell.len(), indices.join(" "))?;
671    }
672    writeln!(w, "CELL_TYPES {}", cell_types.len())?;
673    for &ct in cell_types {
674        writeln!(w, "{ct}")?;
675    }
676    if !stress_tensors.is_empty() {
677        writeln!(w, "CELL_DATA {}", stress_tensors.len())?;
678        writeln!(w, "TENSORS stress double")?;
679        for sigma in stress_tensors {
680            for row in sigma {
681                writeln!(w, "{:.10e} {:.10e} {:.10e}", row[0], row[1], row[2])?;
682            }
683            writeln!(w)?;
684        }
685    }
686    w.flush()?;
687    Ok(())
688}
689/// Write VTK POINT_DATA with velocity vectors for each point.
690#[allow(dead_code)]
691pub fn write_vtk_point_velocity(
692    path: &str,
693    points: &[[f64; 3]],
694    velocities: &[[f64; 3]],
695) -> crate::Result<()> {
696    let file = std::fs::File::create(path)?;
697    let mut w = std::io::BufWriter::new(file);
698    writeln!(w, "# vtk DataFile Version 3.0")?;
699    writeln!(w, "OxiPhysics velocity field")?;
700    writeln!(w, "ASCII")?;
701    writeln!(w, "DATASET POLYDATA")?;
702    writeln!(w, "POINTS {} double", points.len())?;
703    for p in points {
704        writeln!(w, "{:.10e} {:.10e} {:.10e}", p[0], p[1], p[2])?;
705    }
706    writeln!(w, "POINT_DATA {}", points.len())?;
707    writeln!(w, "VECTORS velocity double")?;
708    for v in velocities {
709        writeln!(w, "{:.10e} {:.10e} {:.10e}", v[0], v[1], v[2])?;
710    }
711    w.flush()?;
712    Ok(())
713}
714#[cfg(test)]
715mod tests_vtk_extended {
716    use super::*;
717
718    use crate::vtk::types::*;
719
720    #[test]
721    fn test_vtk_legacy_data_empty() {
722        let d = VtkLegacyData::empty();
723        assert!(d.points.is_empty());
724        assert!(d.point_scalars.is_empty());
725    }
726    #[test]
727    fn test_read_vtk_legacy_ascii_points() {
728        let content = "# vtk DataFile Version 3.0\nTest\nASCII\nDATASET POLYDATA\n\
729            POINTS 3 float\n0.0 0.0 0.0\n1.0 0.0 0.0\n0.0 1.0 0.0\n";
730        let data = read_vtk_legacy_ascii(content).unwrap();
731        assert_eq!(data.points.len(), 3);
732        assert_eq!(data.title, "Test");
733        assert_eq!(data.dataset_type, "POLYDATA");
734    }
735    #[test]
736    fn test_read_vtk_legacy_ascii_point_values() {
737        let content = "# vtk DataFile Version 3.0\nTest\nASCII\nDATASET POLYDATA\n\
738            POINTS 2 float\n0.0 0.0 0.0\n1.0 0.0 0.0\n";
739        let data = read_vtk_legacy_ascii(content).unwrap();
740        assert!((data.points[1][0] - 1.0).abs() < 1e-10);
741    }
742    #[test]
743    fn test_read_vtk_legacy_ascii_scalars() {
744        let content = "# vtk DataFile Version 3.0\nTest\nASCII\nDATASET POLYDATA\n\
745            POINTS 2 float\n0.0 0.0 0.0\n1.0 0.0 0.0\n\
746            POINT_DATA 2\nSCALARS pressure float 1\nLOOKUP_TABLE default\n0.5 1.5\n";
747        let data = read_vtk_legacy_ascii(content).unwrap();
748        assert_eq!(data.point_scalars.len(), 1);
749        assert_eq!(data.point_scalars[0].0, "pressure");
750        assert!((data.point_scalars[0].1[0] - 0.5).abs() < 1e-10);
751        assert!((data.point_scalars[0].1[1] - 1.5).abs() < 1e-10);
752    }
753    #[test]
754    fn test_read_vtk_legacy_ascii_bad_magic() {
755        let content = "NOTAVTK\nTitle\nASCII\nDATASET POLYDATA\n";
756        assert!(read_vtk_legacy_ascii(content).is_err());
757    }
758    #[test]
759    fn test_read_vtk_legacy_ascii_binary_unsupported() {
760        let content = "# vtk DataFile Version 3.0\nTitle\nBINARY\nDATASET POLYDATA\n";
761        assert!(read_vtk_legacy_ascii(content).is_err());
762    }
763    #[test]
764    fn test_pvd_entry_new() {
765        let e = PvdEntry::new(0.5, "step_000000.vtu");
766        assert!((e.time - 0.5).abs() < 1e-12);
767        assert_eq!(e.filename, "step_000000.vtu");
768        assert_eq!(e.part, 0);
769    }
770    #[test]
771    fn test_pvd_entries_uniform_count() {
772        let entries = pvd_entries_uniform("sim", "vtu", 0.0, 0.01, 5);
773        assert_eq!(entries.len(), 5);
774    }
775    #[test]
776    fn test_pvd_entries_uniform_times() {
777        let entries = pvd_entries_uniform("sim", "vtu", 1.0, 0.1, 3);
778        assert!((entries[0].time - 1.0).abs() < 1e-12);
779        assert!((entries[1].time - 1.1).abs() < 1e-12);
780        assert!((entries[2].time - 1.2).abs() < 1e-12);
781    }
782    #[test]
783    fn test_pvd_entries_uniform_filenames() {
784        let entries = pvd_entries_uniform("flow", "vtu", 0.0, 1.0, 2);
785        assert_eq!(entries[0].filename, "flow_000000.vtu");
786        assert_eq!(entries[1].filename, "flow_000001.vtu");
787    }
788    #[test]
789    fn test_write_pvd_file() {
790        let entries = pvd_entries_uniform("sim", "vtu", 0.0, 0.01, 3);
791        let path = "/tmp/test_oxiphysics_pvd.pvd";
792        write_pvd_file(path, &entries).unwrap();
793        let content = std::fs::read_to_string(path).unwrap();
794        assert!(content.contains("Collection"));
795        assert!(content.contains("timestep"));
796        assert!(content.contains("sim_000000.vtu"));
797        std::fs::remove_file(path).ok();
798    }
799    #[test]
800    fn test_write_vtk_polydata_no_normals() {
801        let verts = vec![[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]];
802        let tris = vec![[0, 1, 2]];
803        let path = "/tmp/test_oxiphysics_polydata.vtk";
804        write_vtk_polydata(path, &verts, &tris, None).unwrap();
805        let content = std::fs::read_to_string(path).unwrap();
806        assert!(content.contains("POLYDATA"));
807        assert!(content.contains("POINTS 3"));
808        assert!(content.contains("POLYGONS 1"));
809        std::fs::remove_file(path).ok();
810    }
811    #[test]
812    fn test_write_vtk_polydata_with_normals() {
813        let verts = vec![[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]];
814        let tris = vec![[0, 1, 2]];
815        let normals = vec![[0.0, 0.0, 1.0]; 3];
816        let path = "/tmp/test_oxiphysics_polydata_normals.vtk";
817        write_vtk_polydata(path, &verts, &tris, Some(&normals)).unwrap();
818        let content = std::fs::read_to_string(path).unwrap();
819        assert!(content.contains("NORMALS"));
820        std::fs::remove_file(path).ok();
821    }
822    #[test]
823    fn test_write_vtk_cell_stress() {
824        let points = vec![[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.5, 1.0, 0.0]];
825        let cells = vec![vec![0_usize, 1, 2]];
826        let cell_types = vec![5_u8];
827        let stress = vec![[[1.0, 0.5, 0.0], [0.5, 2.0, 0.0], [0.0, 0.0, 0.5]]];
828        let path = "/tmp/test_oxiphysics_stress.vtk";
829        write_vtk_cell_stress(path, &points, &cells, &cell_types, &stress).unwrap();
830        let content = std::fs::read_to_string(path).unwrap();
831        assert!(content.contains("TENSORS stress"));
832        std::fs::remove_file(path).ok();
833    }
834    #[test]
835    fn test_write_vtk_point_velocity() {
836        let points = vec![[0.0, 0.0, 0.0], [1.0, 0.0, 0.0]];
837        let velocities = vec![[0.1, 0.2, 0.3], [0.4, 0.5, 0.6]];
838        let path = "/tmp/test_oxiphysics_velocity.vtk";
839        write_vtk_point_velocity(path, &points, &velocities).unwrap();
840        let content = std::fs::read_to_string(path).unwrap();
841        assert!(content.contains("VECTORS velocity"));
842        std::fs::remove_file(path).ok();
843    }
844    #[test]
845    fn test_write_vtk_binary_unstructured() {
846        let points = vec![[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.5, 1.0, 0.0]];
847        let cells = vec![vec![0_usize, 1, 2]];
848        let cell_types = vec![5_u8];
849        let path = "/tmp/test_oxiphysics_binary.vtk";
850        write_vtk_binary_unstructured(path, &points, &cells, &cell_types).unwrap();
851        let meta = std::fs::metadata(path).unwrap();
852        assert!(meta.len() > 0);
853        std::fs::remove_file(path).ok();
854    }
855    #[test]
856    fn test_write_vtk_binary_header_ascii() {
857        let points = vec![[0.5, 0.5, 0.5]];
858        let cells = vec![vec![0_usize]];
859        let cell_types = vec![1_u8];
860        let path = "/tmp/test_oxiphysics_binary2.vtk";
861        write_vtk_binary_unstructured(path, &points, &cells, &cell_types).unwrap();
862        let bytes = std::fs::read(path).unwrap();
863        let header_str = std::str::from_utf8(&bytes[..5]).unwrap_or("");
864        assert_eq!(header_str, "# vtk");
865        std::fs::remove_file(path).ok();
866    }
867}
868/// Returns the expected number of nodes for a given [`VtkCellType`].
869///
870/// Returns `0` for unknown types.
871#[allow(dead_code)]
872pub fn cell_type_node_count(ct: VtkCellType) -> usize {
873    match ct {
874        VtkCellType::Vertex => 1,
875        VtkCellType::Line => 2,
876        VtkCellType::Triangle => 3,
877        VtkCellType::Quad => 4,
878        VtkCellType::Tetra => 4,
879        VtkCellType::Hexahedron => 8,
880        VtkCellType::Wedge => 6,
881        VtkCellType::Pyramid => 5,
882    }
883}
884/// Return a human-readable name for a [`VtkCellType`].
885#[allow(dead_code)]
886pub fn cell_type_name(ct: VtkCellType) -> &'static str {
887    match ct {
888        VtkCellType::Vertex => "Vertex",
889        VtkCellType::Line => "Line",
890        VtkCellType::Triangle => "Triangle",
891        VtkCellType::Quad => "Quad",
892        VtkCellType::Tetra => "Tetra",
893        VtkCellType::Hexahedron => "Hexahedron",
894        VtkCellType::Wedge => "Wedge",
895        VtkCellType::Pyramid => "Pyramid",
896    }
897}
898/// Compute axis-aligned bounding box of a [`VtuGrid`]'s points.
899///
900/// Returns `([min_x, min_y, min_z], [max_x, max_y, max_z])`, or `None` if
901/// the grid has no points.
902#[allow(dead_code)]
903pub fn grid_bounding_box(grid: &VtuGrid) -> Option<([f64; 3], [f64; 3])> {
904    if grid.points.is_empty() {
905        return None;
906    }
907    let mut lo = [f64::INFINITY; 3];
908    let mut hi = [f64::NEG_INFINITY; 3];
909    for p in &grid.points {
910        for d in 0..3 {
911            if p[d] < lo[d] {
912                lo[d] = p[d];
913            }
914            if p[d] > hi[d] {
915                hi[d] = p[d];
916            }
917        }
918    }
919    Some((lo, hi))
920}
921/// Compute centroid of all points in a [`VtuGrid`].
922#[allow(dead_code)]
923pub fn grid_centroid(grid: &VtuGrid) -> Option<[f64; 3]> {
924    let n = grid.points.len();
925    if n == 0 {
926        return None;
927    }
928    let sum = grid.points.iter().fold([0.0_f64; 3], |acc, p| {
929        [acc[0] + p[0], acc[1] + p[1], acc[2] + p[2]]
930    });
931    Some([sum[0] / n as f64, sum[1] / n as f64, sum[2] / n as f64])
932}
933/// Merge two [`VtuGrid`]s into one (points and cells are concatenated).
934///
935/// Point indices in `b`'s cells are offset by `a.n_points()`.
936#[allow(dead_code)]
937pub fn merge_vtu_grids(a: &VtuGrid, b: &VtuGrid) -> VtuGrid {
938    let mut out = VtuGrid::new();
939    for &p in &a.points {
940        out.add_point(p);
941    }
942    for &p in &b.points {
943        out.add_point(p);
944    }
945    for (conn, &ct) in a.cells.iter().zip(a.cell_types.iter()) {
946        out.add_cell(conn.clone(), ct);
947    }
948    let offset = a.n_points();
949    for (conn, &ct) in b.cells.iter().zip(b.cell_types.iter()) {
950        let shifted: Vec<usize> = conn.iter().map(|&i| i + offset).collect();
951        out.add_cell(shifted, ct);
952    }
953    out
954}
955/// Compute per-cell volume for tetrahedral cells in a [`VtuGrid`].
956///
957/// Non-tetrahedral cells get volume `0.0`.
958#[allow(dead_code)]
959pub fn compute_cell_volumes(grid: &VtuGrid) -> Vec<f64> {
960    grid.cells
961        .iter()
962        .zip(grid.cell_types.iter())
963        .map(|(conn, &ct)| {
964            if matches!(ct, VtkCellType::Tetra) && conn.len() == 4 {
965                let p = |i: usize| grid.points[conn[i]];
966                let [a, b, c, d] = [p(0), p(1), p(2), p(3)];
967                let e1 = [b[0] - a[0], b[1] - a[1], b[2] - a[2]];
968                let e2 = [c[0] - a[0], c[1] - a[1], c[2] - a[2]];
969                let e3 = [d[0] - a[0], d[1] - a[1], d[2] - a[2]];
970                let cross = [
971                    e2[1] * e3[2] - e2[2] * e3[1],
972                    e2[2] * e3[0] - e2[0] * e3[2],
973                    e2[0] * e3[1] - e2[1] * e3[0],
974                ];
975                (e1[0] * cross[0] + e1[1] * cross[1] + e1[2] * cross[2]).abs() / 6.0
976            } else {
977                0.0
978            }
979        })
980        .collect()
981}
982/// Compute the Euclidean distance between two 3-D points.
983#[allow(dead_code)]
984pub fn point_distance(a: [f64; 3], b: [f64; 3]) -> f64 {
985    let dx = b[0] - a[0];
986    let dy = b[1] - a[1];
987    let dz = b[2] - a[2];
988    (dx * dx + dy * dy + dz * dz).sqrt()
989}
990#[cfg(test)]
991mod vtk_extended_tests {
992    use super::*;
993    use crate::VtkDataArray;
994
995    use crate::vtk::VtkFieldData;
996    use crate::vtk::VtkMultiBlock;
997
998    use crate::vtk::VtkTimeSeries;
999    use crate::vtk::types::*;
1000
1001    #[test]
1002    fn test_field_data_add_and_get() {
1003        let mut fd = VtkFieldData::new();
1004        fd.add("time", vec![1.5, 2.5]);
1005        let rec = fd.get("time").expect("should find 'time'");
1006        assert_eq!(rec.values, vec![1.5, 2.5]);
1007    }
1008    #[test]
1009    fn test_field_data_add_scalar() {
1010        let mut fd = VtkFieldData::new();
1011        fd.add_scalar("step", 42.0);
1012        let rec = fd.get("step").expect("should find 'step'");
1013        assert_eq!(rec.values.len(), 1);
1014        assert!((rec.values[0] - 42.0).abs() < 1e-12);
1015    }
1016    #[test]
1017    fn test_field_data_get_missing_returns_none() {
1018        let fd = VtkFieldData::new();
1019        assert!(fd.get("missing").is_none());
1020    }
1021    #[test]
1022    fn test_field_data_to_vtk_string_contains_name() {
1023        let mut fd = VtkFieldData::new();
1024        fd.add("pressure", vec![101325.0]);
1025        let s = fd.to_vtk_field_string();
1026        assert!(
1027            s.contains("pressure"),
1028            "field string should contain field name"
1029        );
1030        assert!(
1031            s.contains("FIELD"),
1032            "field string should contain FIELD keyword"
1033        );
1034    }
1035    #[test]
1036    fn test_field_data_empty_string_empty() {
1037        let fd = VtkFieldData::new();
1038        assert!(fd.to_vtk_field_string().is_empty());
1039    }
1040    #[test]
1041    fn test_field_data_multiple_records() {
1042        let mut fd = VtkFieldData::new();
1043        fd.add_scalar("dt", 0.001);
1044        fd.add_scalar("iter", 100.0);
1045        let s = fd.to_vtk_field_string();
1046        assert!(s.contains("dt"));
1047        assert!(s.contains("iter"));
1048        assert!(s.contains("FieldData 2"));
1049    }
1050    #[test]
1051    fn test_multi_block_empty() {
1052        let mb = VtkMultiBlock::new("empty");
1053        assert_eq!(mb.n_blocks(), 0);
1054        assert_eq!(mb.total_points(), 0);
1055    }
1056    #[test]
1057    fn test_multi_block_add_block() {
1058        let mut mb = VtkMultiBlock::new("test");
1059        let g = VtuGrid::from_points(&[[0.0, 0.0, 0.0], [1.0, 0.0, 0.0]]);
1060        mb.add_block("domain_a", g);
1061        assert_eq!(mb.n_blocks(), 1);
1062        assert_eq!(mb.total_points(), 2);
1063    }
1064    #[test]
1065    fn test_multi_block_total_cells() {
1066        let mut mb = VtkMultiBlock::new("m");
1067        let mut g1 = VtuGrid::new();
1068        g1.add_point([0.0, 0.0, 0.0]);
1069        g1.add_point([1.0, 0.0, 0.0]);
1070        g1.add_cell(vec![0, 1], VtkCellType::Line);
1071        let g2 = VtuGrid::from_points(&[[5.0, 0.0, 0.0]]);
1072        mb.add_block("a", g1);
1073        mb.add_block("b", g2);
1074        assert_eq!(mb.total_cells(), 2);
1075    }
1076    #[test]
1077    fn test_multi_block_vtm_string_contains_blocks() {
1078        let mut mb = VtkMultiBlock::new("sim");
1079        mb.add_block("fluid", VtuGrid::new());
1080        mb.add_block("solid", VtuGrid::new());
1081        let vtm = mb.to_vtm_string();
1082        assert!(vtm.contains("fluid"));
1083        assert!(vtm.contains("solid"));
1084        assert!(vtm.contains("vtkMultiBlockDataSet"));
1085    }
1086    #[test]
1087    fn test_multi_block_write_to_dir() {
1088        let mut mb = VtkMultiBlock::new("test_mb");
1089        mb.add_block("part0", VtuGrid::from_points(&[[0.0, 0.0, 0.0]]));
1090        let written = mb.write_to_dir("/tmp").unwrap();
1091        assert!(!written.is_empty());
1092        for f in &written {
1093            assert!(std::path::Path::new(f).exists(), "file should exist: {f}");
1094            std::fs::remove_file(f).ok();
1095        }
1096    }
1097    #[test]
1098    fn test_cell_type_node_counts() {
1099        assert_eq!(cell_type_node_count(VtkCellType::Vertex), 1);
1100        assert_eq!(cell_type_node_count(VtkCellType::Line), 2);
1101        assert_eq!(cell_type_node_count(VtkCellType::Triangle), 3);
1102        assert_eq!(cell_type_node_count(VtkCellType::Quad), 4);
1103        assert_eq!(cell_type_node_count(VtkCellType::Tetra), 4);
1104        assert_eq!(cell_type_node_count(VtkCellType::Hexahedron), 8);
1105        assert_eq!(cell_type_node_count(VtkCellType::Wedge), 6);
1106        assert_eq!(cell_type_node_count(VtkCellType::Pyramid), 5);
1107    }
1108    #[test]
1109    fn test_cell_type_names() {
1110        assert_eq!(cell_type_name(VtkCellType::Triangle), "Triangle");
1111        assert_eq!(cell_type_name(VtkCellType::Tetra), "Tetra");
1112        assert_eq!(cell_type_name(VtkCellType::Hexahedron), "Hexahedron");
1113    }
1114    #[test]
1115    fn test_grid_bounding_box_empty_is_none() {
1116        let grid = VtuGrid::new();
1117        assert!(grid_bounding_box(&grid).is_none());
1118    }
1119    #[test]
1120    fn test_grid_bounding_box_values() {
1121        let mut g = VtuGrid::new();
1122        g.add_point([1.0, 2.0, 3.0]);
1123        g.add_point([4.0, 5.0, 6.0]);
1124        let (lo, hi) = grid_bounding_box(&g).unwrap();
1125        assert!((lo[0] - 1.0).abs() < 1e-12);
1126        assert!((hi[0] - 4.0).abs() < 1e-12);
1127    }
1128    #[test]
1129    fn test_grid_centroid_single_point() {
1130        let mut g = VtuGrid::new();
1131        g.add_point([3.0, 4.0, 5.0]);
1132        let c = grid_centroid(&g).unwrap();
1133        assert!((c[0] - 3.0).abs() < 1e-12);
1134    }
1135    #[test]
1136    fn test_grid_centroid_two_points() {
1137        let mut g = VtuGrid::new();
1138        g.add_point([0.0, 0.0, 0.0]);
1139        g.add_point([2.0, 4.0, 6.0]);
1140        let c = grid_centroid(&g).unwrap();
1141        assert!((c[0] - 1.0).abs() < 1e-12);
1142        assert!((c[1] - 2.0).abs() < 1e-12);
1143    }
1144    #[test]
1145    fn test_grid_centroid_empty_is_none() {
1146        let g = VtuGrid::new();
1147        assert!(grid_centroid(&g).is_none());
1148    }
1149    #[test]
1150    fn test_merge_vtu_grids_point_count() {
1151        let a = VtuGrid::from_points(&[[0.0, 0.0, 0.0], [1.0, 0.0, 0.0]]);
1152        let b = VtuGrid::from_points(&[[2.0, 0.0, 0.0]]);
1153        let merged = merge_vtu_grids(&a, &b);
1154        assert_eq!(merged.n_points(), 3);
1155    }
1156    #[test]
1157    fn test_merge_vtu_grids_cell_count() {
1158        let a = VtuGrid::from_points(&[[0.0, 0.0, 0.0]]);
1159        let b = VtuGrid::from_points(&[[1.0, 0.0, 0.0], [2.0, 0.0, 0.0]]);
1160        let merged = merge_vtu_grids(&a, &b);
1161        assert_eq!(merged.n_cells(), 3);
1162    }
1163    #[test]
1164    fn test_merge_vtu_grids_indices_offset() {
1165        let mut a = VtuGrid::new();
1166        a.add_point([0.0, 0.0, 0.0]);
1167        a.add_point([1.0, 0.0, 0.0]);
1168        a.add_cell(vec![0, 1], VtkCellType::Line);
1169        let mut b = VtuGrid::new();
1170        b.add_point([2.0, 0.0, 0.0]);
1171        b.add_point([3.0, 0.0, 0.0]);
1172        b.add_cell(vec![0, 1], VtkCellType::Line);
1173        let merged = merge_vtu_grids(&a, &b);
1174        assert_eq!(merged.cells[1], vec![2, 3]);
1175    }
1176    #[test]
1177    fn test_compute_cell_volumes_tet() {
1178        let nodes = [
1179            [0.0, 0.0, 0.0],
1180            [1.0, 0.0, 0.0],
1181            [0.0, 1.0, 0.0],
1182            [0.0, 0.0, 1.0],
1183        ];
1184        let grid = VtuGrid::from_tet_mesh(&nodes, &[[0, 1, 2, 3]]);
1185        let vols = compute_cell_volumes(&grid);
1186        assert_eq!(vols.len(), 1);
1187        assert!(
1188            (vols[0] - 1.0 / 6.0).abs() < 1e-10,
1189            "tet volume should be 1/6, got {}",
1190            vols[0]
1191        );
1192    }
1193    #[test]
1194    fn test_compute_cell_volumes_vertex_zero() {
1195        let grid = VtuGrid::from_points(&[[0.0, 0.0, 0.0]]);
1196        let vols = compute_cell_volumes(&grid);
1197        assert_eq!(vols[0], 0.0);
1198    }
1199    #[test]
1200    fn test_compute_cell_volumes_multiple_tets() {
1201        let nodes = [
1202            [0.0, 0.0, 0.0],
1203            [1.0, 0.0, 0.0],
1204            [0.0, 1.0, 0.0],
1205            [0.0, 0.0, 1.0],
1206            [2.0, 0.0, 0.0],
1207            [3.0, 0.0, 0.0],
1208            [2.0, 1.0, 0.0],
1209            [2.0, 0.0, 1.0],
1210        ];
1211        let elements = [[0, 1, 2, 3], [4, 5, 6, 7]];
1212        let grid = VtuGrid::from_tet_mesh(&nodes, &elements);
1213        let vols = compute_cell_volumes(&grid);
1214        assert_eq!(vols.len(), 2);
1215        for v in &vols {
1216            assert!(
1217                (v - 1.0 / 6.0).abs() < 1e-10,
1218                "each tet should have volume 1/6"
1219            );
1220        }
1221    }
1222    #[test]
1223    fn test_point_distance_known() {
1224        let d = point_distance([0.0, 0.0, 0.0], [3.0, 4.0, 0.0]);
1225        assert!((d - 5.0).abs() < 1e-12);
1226    }
1227    #[test]
1228    fn test_point_distance_same_point_is_zero() {
1229        let d = point_distance([1.0, 2.0, 3.0], [1.0, 2.0, 3.0]);
1230        assert!(d.abs() < 1e-12);
1231    }
1232    #[test]
1233    fn test_vtk_data_array_scalar_len() {
1234        let arr = VtkDataArray::Scalar {
1235            name: "p".into(),
1236            values: vec![1.0, 2.0, 3.0],
1237        };
1238        assert_eq!(arr.len(), 3);
1239        assert_eq!(arr.n_components(), 1);
1240        assert!(!arr.is_empty());
1241    }
1242    #[test]
1243    fn test_vtk_data_array_vector3_n_components() {
1244        let arr = VtkDataArray::Vector3 {
1245            name: "v".into(),
1246            values: vec![[1.0, 0.0, 0.0], [0.0, 1.0, 0.0]],
1247        };
1248        assert_eq!(arr.n_components(), 3);
1249        assert_eq!(arr.len(), 2);
1250    }
1251    #[test]
1252    fn test_vtk_data_array_integer() {
1253        let arr = VtkDataArray::Integer {
1254            name: "id".into(),
1255            values: vec![1, 2, 3, 4],
1256        };
1257        assert_eq!(arr.len(), 4);
1258        assert_eq!(arr.name(), "id");
1259    }
1260    #[test]
1261    fn test_vtk_data_array_empty() {
1262        let arr = VtkDataArray::Scalar {
1263            name: "empty".into(),
1264            values: vec![],
1265        };
1266        assert!(arr.is_empty());
1267    }
1268    #[test]
1269    fn test_vtu_cell_scalar_data_in_xml() {
1270        let mut grid = VtuGrid::new();
1271        for i in 0..4 {
1272            grid.add_point([i as f64, 0.0, 0.0]);
1273        }
1274        grid.add_cell(vec![0, 1, 2, 3], VtkCellType::Quad);
1275        grid.add_cell_scalar("material_id", vec![7.0]);
1276        let xml = grid.to_vtu_string();
1277        assert!(xml.contains("<CellData>"), "CellData section missing");
1278        assert!(xml.contains("material_id"), "material_id field missing");
1279    }
1280    #[test]
1281    fn test_vtu_vector_point_data_in_xml() {
1282        let mut grid = VtuGrid::new();
1283        grid.add_point([0.0, 0.0, 0.0]);
1284        grid.add_point([1.0, 0.0, 0.0]);
1285        grid.add_point([0.0, 1.0, 0.0]);
1286        grid.add_cell(vec![0, 1, 2], VtkCellType::Triangle);
1287        grid.add_point_vector("velocity", vec![[1.0, 0.0, 0.0]; 3]);
1288        let xml = grid.to_vtu_string();
1289        assert!(xml.contains("velocity"), "velocity field should be in XML");
1290        assert!(xml.contains("NumberOfComponents=\"3\""));
1291    }
1292    #[test]
1293    fn test_time_series_estimated_size() {
1294        let mut ts = VtkTimeSeries::new("s");
1295        let mut g = VtuGrid::new();
1296        for i in 0..10 {
1297            g.add_point([i as f64, 0.0, 0.0]);
1298        }
1299        ts.push(0.0, g);
1300        let sz = ts.estimated_size_bytes();
1301        assert!(sz > 0);
1302    }
1303    #[test]
1304    fn test_time_series_increasing_times() {
1305        let mut ts = VtkTimeSeries::new("flow");
1306        ts.push(0.0, VtuGrid::new());
1307        ts.push(1.0, VtuGrid::new());
1308        ts.push(2.5, VtuGrid::new());
1309        assert_eq!(ts.n_steps(), 3);
1310        assert!((ts.times[2] - 2.5).abs() < 1e-12);
1311    }
1312    #[test]
1313    fn test_validate_vtu_grid_cell_data_mismatch() {
1314        let mut grid = VtuGrid::new();
1315        grid.add_point([0.0, 0.0, 0.0]);
1316        grid.add_point([1.0, 0.0, 0.0]);
1317        grid.add_point([0.0, 1.0, 0.0]);
1318        grid.add_cell(vec![0, 1, 2], VtkCellType::Triangle);
1319        grid.cell_data.push(VtkDataArray::Scalar {
1320            name: "bad_cell".into(),
1321            values: vec![1.0, 2.0],
1322        });
1323        let errors = validate_vtu_grid(&grid);
1324        assert!(
1325            !errors.is_empty(),
1326            "mismatched cell data should produce error"
1327        );
1328    }
1329    #[test]
1330    fn test_validate_vtu_grid_empty_valid() {
1331        let grid = VtuGrid::new();
1332        let errors = validate_vtu_grid(&grid);
1333        assert!(errors.is_empty(), "empty grid should be valid");
1334    }
1335}