1use crate::Result;
6use oxiphysics_core::math::Vec3;
7use std::fs::File;
8use std::io::{BufWriter, Write};
9use std::path::Path;
10
11pub struct VtuGrid {
17 pub points: Vec<[f64; 3]>,
19 pub cells: Vec<Vec<usize>>,
21 pub cell_types: Vec<VtkCellType>,
23 pub point_data: Vec<VtkDataArray>,
25 pub cell_data: Vec<VtkDataArray>,
27}
28impl VtuGrid {
29 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 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 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 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 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 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 pub fn n_points(&self) -> usize {
73 self.points.len()
74 }
75 pub fn n_cells(&self) -> usize {
77 self.cells.len()
78 }
79 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 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 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 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 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#[derive(Debug, Clone)]
267pub struct VtkFieldRecord {
268 pub name: String,
270 pub values: Vec<f64>,
272}
273impl VtkFieldRecord {
274 pub fn new(name: impl Into<String>, values: Vec<f64>) -> Self {
276 Self {
277 name: name.into(),
278 values,
279 }
280 }
281 pub fn len(&self) -> usize {
283 self.values.len()
284 }
285 pub fn is_empty(&self) -> bool {
287 self.values.is_empty()
288 }
289}
290#[derive(Debug, Clone, Copy)]
292pub enum VtkCellType {
293 Vertex = 1,
295 Line = 3,
297 Triangle = 5,
299 Quad = 9,
301 Tetra = 10,
303 Hexahedron = 12,
305 Wedge = 13,
307 Pyramid = 14,
309}
310#[derive(Debug, Clone)]
312pub struct VtkLegacyData {
313 pub title: String,
315 pub dataset_type: String,
317 pub points: Vec<[f64; 3]>,
319 pub point_scalars: Vec<(String, Vec<f64>)>,
321}
322impl VtkLegacyData {
323 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#[derive(Debug, Clone, Default)]
338pub struct VtkFieldData {
339 pub records: Vec<VtkFieldRecord>,
341}
342impl VtkFieldData {
343 pub fn new() -> Self {
345 Self {
346 records: Vec::new(),
347 }
348 }
349 pub fn add(&mut self, name: impl Into<String>, values: Vec<f64>) {
351 self.records.push(VtkFieldRecord::new(name, values));
352 }
353 pub fn add_scalar(&mut self, name: impl Into<String>, value: f64) {
355 self.add(name, vec![value]);
356 }
357 pub fn get(&self, name: &str) -> Option<&VtkFieldRecord> {
359 self.records.iter().find(|r| r.name == name)
360 }
361 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}
377pub struct VtkWriter;
379impl VtkWriter {
380 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 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 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}
466pub struct VtkTimeSeries {
468 pub times: Vec<f64>,
470 pub grids: Vec<VtuGrid>,
472 pub base_name: String,
474}
475impl VtkTimeSeries {
476 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 pub fn push(&mut self, t: f64, grid: VtuGrid) {
486 self.times.push(t);
487 self.grids.push(grid);
488 }
489 pub fn n_steps(&self) -> usize {
491 self.times.len()
492 }
493 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 pub fn vtu_string(&self, i: usize) -> Option<String> {
507 self.grids.get(i).map(|g| g.to_vtu_string())
508 }
509 pub fn estimated_size_bytes(&self) -> usize {
511 self.grids.iter().map(|g| g.n_points() * 30).sum()
512 }
513}
514#[derive(Debug, Clone)]
516pub struct PvdEntry {
517 pub time: f64,
519 pub filename: String,
521 pub part: u32,
523 pub group: String,
525}
526impl PvdEntry {
527 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}
537pub struct VtkMultiBlock {
542 pub blocks: Vec<VtkBlock>,
544 pub title: String,
546}
547impl VtkMultiBlock {
548 pub fn new(title: impl Into<String>) -> Self {
550 Self {
551 blocks: Vec::new(),
552 title: title.into(),
553 }
554 }
555 pub fn add_block(&mut self, name: impl Into<String>, grid: VtuGrid) {
557 self.blocks.push(VtkBlock::new(name, grid));
558 }
559 pub fn n_blocks(&self) -> usize {
561 self.blocks.len()
562 }
563 pub fn total_points(&self) -> usize {
565 self.blocks.iter().map(|b| b.grid.n_points()).sum()
566 }
567 pub fn total_cells(&self) -> usize {
569 self.blocks.iter().map(|b| b.grid.n_cells()).sum()
570 }
571 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 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}
612pub struct VtkPolyDataGrid {
614 pub points: Vec<[f64; 3]>,
616 pub lines: Vec<[usize; 2]>,
618 pub triangles: Vec<[usize; 3]>,
620 pub point_data: Vec<VtkDataArray>,
622 pub cell_data: Vec<VtkDataArray>,
624}
625impl VtkPolyDataGrid {
626 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 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 pub fn add_triangle(&mut self, a: usize, b: usize, c: usize) {
644 self.triangles.push([a, b, c]);
645 }
646 pub fn add_line(&mut self, a: usize, b: usize) {
648 self.lines.push([a, b]);
649 }
650 pub fn n_points(&self) -> usize {
652 self.points.len()
653 }
654 pub fn n_triangles(&self) -> usize {
656 self.triangles.len()
657 }
658 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 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#[derive(Debug, Clone)]
699pub enum VtkDataArray {
700 Scalar {
702 name: String,
704 values: Vec<f64>,
706 },
707 Vector3 {
709 name: String,
711 values: Vec<[f64; 3]>,
713 },
714 Integer {
716 name: String,
718 values: Vec<i64>,
720 },
721}
722impl VtkDataArray {
723 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 pub fn n_components(&self) -> usize {
733 match self {
734 Self::Scalar { .. } => 1,
735 Self::Vector3 { .. } => 3,
736 Self::Integer { .. } => 1,
737 }
738 }
739 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 pub fn is_empty(&self) -> bool {
749 self.len() == 0
750 }
751}
752pub struct VtkBlock {
754 pub name: String,
756 pub grid: VtuGrid,
758}
759impl VtkBlock {
760 pub fn new(name: impl Into<String>, grid: VtuGrid) -> Self {
762 Self {
763 name: name.into(),
764 grid,
765 }
766 }
767}
768pub struct VtkRectilinearGrid {
770 pub x_coords: Vec<f64>,
772 pub y_coords: Vec<f64>,
774 pub z_coords: Vec<f64>,
776 pub point_data: Vec<VtkDataArray>,
778}
779impl VtkRectilinearGrid {
780 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 pub fn n_points(&self) -> usize {
791 self.x_coords.len() * self.y_coords.len() * self.z_coords.len()
792 }
793 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 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 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}
872pub struct VtkStructuredGrid {
877 pub dims: [usize; 3],
879 pub points: Vec<[f64; 3]>,
881 pub point_data: Vec<VtkDataArray>,
883}
884impl VtkStructuredGrid {
885 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 pub fn n_points(&self) -> usize {
895 self.dims[0] * self.dims[1] * self.dims[2]
896 }
897 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 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 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}