1#[allow(unused_imports)]
6use super::functions_2::*;
7use crate::{Error, Result};
8use oxiphysics_core::math::Vec3;
9use std::fs::File;
10use std::io::{BufWriter, Write};
11use std::path::Path;
12
13#[allow(unused_imports)]
14use super::functions::*;
15use super::functions::{VTK_HEX, VTK_QUAD, VTK_TET, VTK_TRIANGLE};
16
17#[allow(dead_code)]
20pub struct VtuMultiPiece;
21#[allow(dead_code)]
22impl VtuMultiPiece {
23 pub fn write(pieces: &[&VtuWriter]) -> String {
25 let mut s = String::from("<?xml version=\"1.0\"?>\n");
26 s.push_str(
27 "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">\n",
28 );
29 s.push_str(" <UnstructuredGrid>\n");
30 for w in pieces {
31 let xml = w.to_xml();
32 if let Some(start) = xml.find("<Piece")
33 && let Some(end) = xml.rfind("</Piece>")
34 {
35 s.push_str(" ");
36 s.push_str(&xml[start..end + 8]);
37 s.push('\n');
38 }
39 }
40 s.push_str(" </UnstructuredGrid>\n");
41 s.push_str("</VTKFile>\n");
42 s
43 }
44 pub fn total_points(pieces: &[&VtuWriter]) -> usize {
46 pieces.iter().map(|w| w.num_points()).sum()
47 }
48 pub fn total_cells(pieces: &[&VtuWriter]) -> usize {
50 pieces.iter().map(|w| w.num_cells()).sum()
51 }
52}
53#[allow(dead_code)]
55#[derive(Debug, Clone)]
56pub struct VtuAnnotation {
57 pub key: String,
59 pub value: String,
61}
62#[allow(dead_code)]
63pub(super) struct PointVector {
64 pub(super) name: String,
65 pub(super) vectors: Vec<[f64; 3]>,
66}
67#[allow(dead_code)]
71pub struct VtuWriterLegacy;
72#[allow(dead_code)]
73impl VtuWriterLegacy {
74 pub fn write(
78 path: &str,
79 positions: &[Vec3],
80 cells: &[[usize; 4]],
81 point_data: &[(&str, &[f64])],
82 cell_data: &[(&str, &[f64])],
83 ) -> Result<()> {
84 let file = File::create(Path::new(path))?;
85 let mut w = BufWriter::new(file);
86 let npoints = positions.len();
87 let ncells = cells.len();
88 writeln!(w, "<?xml version=\"1.0\"?>")?;
89 writeln!(
90 w,
91 "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">"
92 )?;
93 writeln!(w, " <UnstructuredGrid>")?;
94 writeln!(
95 w,
96 " <Piece NumberOfPoints=\"{}\" NumberOfCells=\"{}\">",
97 npoints, ncells
98 )?;
99 if !point_data.is_empty() {
100 writeln!(w, " <PointData>")?;
101 for (name, vals) in point_data {
102 writeln!(
103 w,
104 " <DataArray type=\"Float64\" Name=\"{}\" format=\"ascii\">",
105 name
106 )?;
107 write!(w, " ")?;
108 for (i, v) in vals.iter().enumerate() {
109 if i > 0 {
110 write!(w, " ")?;
111 }
112 write!(w, "{}", v)?;
113 }
114 writeln!(w)?;
115 writeln!(w, " </DataArray>")?;
116 }
117 writeln!(w, " </PointData>")?;
118 }
119 if !cell_data.is_empty() {
120 writeln!(w, " <CellData>")?;
121 for (name, vals) in cell_data {
122 writeln!(
123 w,
124 " <DataArray type=\"Float64\" Name=\"{}\" format=\"ascii\">",
125 name
126 )?;
127 write!(w, " ")?;
128 for (i, v) in vals.iter().enumerate() {
129 if i > 0 {
130 write!(w, " ")?;
131 }
132 write!(w, "{}", v)?;
133 }
134 writeln!(w)?;
135 writeln!(w, " </DataArray>")?;
136 }
137 writeln!(w, " </CellData>")?;
138 }
139 writeln!(w, " <Points>")?;
140 writeln!(
141 w,
142 " <DataArray type=\"Float64\" NumberOfComponents=\"3\" format=\"ascii\">"
143 )?;
144 for p in positions {
145 writeln!(w, " {} {} {}", p.x, p.y, p.z)?;
146 }
147 writeln!(w, " </DataArray>")?;
148 writeln!(w, " </Points>")?;
149 writeln!(w, " <Cells>")?;
150 writeln!(
151 w,
152 " <DataArray type=\"Int32\" Name=\"connectivity\" format=\"ascii\">"
153 )?;
154 for c in cells {
155 writeln!(w, " {} {} {} {}", c[0], c[1], c[2], c[3])?;
156 }
157 writeln!(w, " </DataArray>")?;
158 writeln!(
159 w,
160 " <DataArray type=\"Int32\" Name=\"offsets\" format=\"ascii\">"
161 )?;
162 write!(w, " ")?;
163 for i in 0..ncells {
164 if i > 0 {
165 write!(w, " ")?;
166 }
167 write!(w, "{}", (i + 1) * 4)?;
168 }
169 writeln!(w)?;
170 writeln!(w, " </DataArray>")?;
171 writeln!(
172 w,
173 " <DataArray type=\"UInt8\" Name=\"types\" format=\"ascii\">"
174 )?;
175 write!(w, " ")?;
176 for i in 0..ncells {
177 if i > 0 {
178 write!(w, " ")?;
179 }
180 write!(w, "10")?;
181 }
182 writeln!(w)?;
183 writeln!(w, " </DataArray>")?;
184 writeln!(w, " </Cells>")?;
185 writeln!(w, " </Piece>")?;
186 writeln!(w, " </UnstructuredGrid>")?;
187 writeln!(w, "</VTKFile>")?;
188 w.flush()?;
189 Ok(())
190 }
191}
192#[allow(dead_code)]
193pub(super) struct CellVector {
194 pub(super) name: String,
195 pub(super) vectors: Vec<[f64; 3]>,
196}
197#[allow(dead_code)]
202pub struct VtuWriter {
203 pub(super) points: Vec<[f64; 3]>,
204 pub(super) cells: Vec<Vec<usize>>,
205 pub(super) cell_types: Vec<u8>,
206 pub(super) point_scalars: Vec<PointScalar>,
207 pub(super) point_vectors: Vec<PointVector>,
208 pub(super) cell_scalars: Vec<CellScalar>,
209 pub(super) cell_vectors: Vec<CellVector>,
210}
211#[allow(dead_code)]
212impl VtuWriter {
213 pub fn new() -> Self {
215 Self {
216 points: Vec::new(),
217 cells: Vec::new(),
218 cell_types: Vec::new(),
219 point_scalars: Vec::new(),
220 point_vectors: Vec::new(),
221 cell_scalars: Vec::new(),
222 cell_vectors: Vec::new(),
223 }
224 }
225 pub fn add_point(&mut self, p: [f64; 3]) -> usize {
227 let idx = self.points.len();
228 self.points.push(p);
229 idx
230 }
231 pub fn add_points(&mut self, pts: &[[f64; 3]]) -> usize {
233 let first = self.points.len();
234 self.points.extend_from_slice(pts);
235 first
236 }
237 pub fn num_points(&self) -> usize {
239 self.points.len()
240 }
241 pub fn num_cells(&self) -> usize {
243 self.cells.len()
244 }
245 pub fn add_cell(&mut self, connectivity: Vec<usize>, cell_type_id: u8) {
247 self.cells.push(connectivity);
248 self.cell_types.push(cell_type_id);
249 }
250 pub fn add_triangle(&mut self, i0: usize, i1: usize, i2: usize) {
252 self.add_cell(vec![i0, i1, i2], VTK_TRIANGLE);
253 }
254 pub fn add_quad(&mut self, i0: usize, i1: usize, i2: usize, i3: usize) {
256 self.add_cell(vec![i0, i1, i2, i3], VTK_QUAD);
257 }
258 pub fn add_tet(&mut self, i0: usize, i1: usize, i2: usize, i3: usize) {
260 self.add_cell(vec![i0, i1, i2, i3], VTK_TET);
261 }
262 #[allow(clippy::too_many_arguments)]
264 pub fn add_hex(
265 &mut self,
266 i0: usize,
267 i1: usize,
268 i2: usize,
269 i3: usize,
270 i4: usize,
271 i5: usize,
272 i6: usize,
273 i7: usize,
274 ) {
275 self.add_cell(vec![i0, i1, i2, i3, i4, i5, i6, i7], VTK_HEX);
276 }
277 pub fn add_point_data_scalar(&mut self, name: &str, values: &[f64]) {
279 self.point_scalars.push(PointScalar {
280 name: name.to_string(),
281 values: values.to_vec(),
282 });
283 }
284 pub fn add_point_data_vector(&mut self, name: &str, vectors: &[[f64; 3]]) {
286 self.point_vectors.push(PointVector {
287 name: name.to_string(),
288 vectors: vectors.to_vec(),
289 });
290 }
291 #[allow(dead_code)]
296 pub fn write_cell_data_tensor(name: &str, tensors: &[[f64; 6]]) -> String {
297 let mut s = String::new();
298 s.push_str(
299 &format!(
300 " <DataArray type=\"Float64\" Name=\"{}\" NumberOfComponents=\"6\" format=\"ascii\">\n ",
301 name
302 ),
303 );
304 let vals: Vec<String> = tensors
305 .iter()
306 .flat_map(|t| t.iter())
307 .map(|v| format!("{}", v))
308 .collect();
309 s.push_str(&vals.join(" "));
310 s.push_str("\n </DataArray>\n");
311 s
312 }
313 #[allow(dead_code)]
318 pub fn write_pvtu_parallel(
319 piece_files: &[&str],
320 point_data_names: &[&str],
321 cell_data_names: &[&str],
322 ) -> String {
323 let mut s = String::new();
324 s.push_str("<?xml version=\"1.0\"?>\n");
325 s.push_str(
326 "<VTKFile type=\"PUnstructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">\n",
327 );
328 s.push_str(" <PUnstructuredGrid GhostLevel=\"0\">\n");
329 s.push_str(" <PPoints>\n");
330 s.push_str(
331 " <PDataArray type=\"Float64\" NumberOfComponents=\"3\" format=\"ascii\"/>\n",
332 );
333 s.push_str(" </PPoints>\n");
334 if !point_data_names.is_empty() {
335 s.push_str(" <PPointData>\n");
336 for name in point_data_names {
337 s.push_str(&format!(
338 " <PDataArray type=\"Float64\" Name=\"{}\" format=\"ascii\"/>\n",
339 name
340 ));
341 }
342 s.push_str(" </PPointData>\n");
343 }
344 if !cell_data_names.is_empty() {
345 s.push_str(" <PCellData>\n");
346 for name in cell_data_names {
347 s.push_str(&format!(
348 " <PDataArray type=\"Float64\" Name=\"{}\" format=\"ascii\"/>\n",
349 name
350 ));
351 }
352 s.push_str(" </PCellData>\n");
353 }
354 for file in piece_files {
355 s.push_str(&format!(" <Piece Source=\"{}\"/>\n", file));
356 }
357 s.push_str(" </PUnstructuredGrid>\n");
358 s.push_str("</VTKFile>\n");
359 s
360 }
361 pub fn add_cell_data_scalar(&mut self, name: &str, values: &[f64]) {
363 self.cell_scalars.push(CellScalar {
364 name: name.to_string(),
365 values: values.to_vec(),
366 });
367 }
368 pub fn add_cell_data_vector(&mut self, name: &str, vectors: &[[f64; 3]]) {
370 self.cell_vectors.push(CellVector {
371 name: name.to_string(),
372 vectors: vectors.to_vec(),
373 });
374 }
375 pub fn bounding_box(&self) -> Option<([f64; 3], [f64; 3])> {
378 if self.points.is_empty() {
379 return None;
380 }
381 let mut min = self.points[0];
382 let mut max = self.points[0];
383 for p in &self.points[1..] {
384 for d in 0..3 {
385 if p[d] < min[d] {
386 min[d] = p[d];
387 }
388 if p[d] > max[d] {
389 max[d] = p[d];
390 }
391 }
392 }
393 Some((min, max))
394 }
395 pub fn to_xml(&self) -> String {
397 let npoints = self.points.len();
398 let ncells = self.cells.len();
399 let mut s = String::new();
400 s.push_str("<?xml version=\"1.0\"?>\n");
401 s.push_str(
402 "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">\n",
403 );
404 s.push_str(" <UnstructuredGrid>\n");
405 s.push_str(&format!(
406 " <Piece NumberOfPoints=\"{}\" NumberOfCells=\"{}\">\n",
407 npoints, ncells
408 ));
409 let has_pdata = !self.point_scalars.is_empty() || !self.point_vectors.is_empty();
410 if has_pdata {
411 s.push_str(" <PointData>\n");
412 for ps in &self.point_scalars {
413 s.push_str(
414 &format!(
415 " <DataArray type=\"Float64\" Name=\"{}\" NumberOfComponents=\"1\" format=\"ascii\">\n ",
416 ps.name
417 ),
418 );
419 let vals: Vec<String> = ps.values.iter().map(|v| format!("{}", v)).collect();
420 s.push_str(&vals.join(" "));
421 s.push_str("\n </DataArray>\n");
422 }
423 for pv in &self.point_vectors {
424 s.push_str(
425 &format!(
426 " <DataArray type=\"Float64\" Name=\"{}\" NumberOfComponents=\"3\" format=\"ascii\">\n ",
427 pv.name
428 ),
429 );
430 let vals: Vec<String> = pv
431 .vectors
432 .iter()
433 .flat_map(|v| v.iter())
434 .map(|x| format!("{}", x))
435 .collect();
436 s.push_str(&vals.join(" "));
437 s.push_str("\n </DataArray>\n");
438 }
439 s.push_str(" </PointData>\n");
440 }
441 let has_cdata = !self.cell_scalars.is_empty() || !self.cell_vectors.is_empty();
442 if has_cdata {
443 s.push_str(" <CellData>\n");
444 for cs in &self.cell_scalars {
445 s.push_str(
446 &format!(
447 " <DataArray type=\"Float64\" Name=\"{}\" NumberOfComponents=\"1\" format=\"ascii\">\n ",
448 cs.name
449 ),
450 );
451 let vals: Vec<String> = cs.values.iter().map(|v| format!("{}", v)).collect();
452 s.push_str(&vals.join(" "));
453 s.push_str("\n </DataArray>\n");
454 }
455 for cv in &self.cell_vectors {
456 s.push_str(
457 &format!(
458 " <DataArray type=\"Float64\" Name=\"{}\" NumberOfComponents=\"3\" format=\"ascii\">\n ",
459 cv.name
460 ),
461 );
462 let vals: Vec<String> = cv
463 .vectors
464 .iter()
465 .flat_map(|v| v.iter())
466 .map(|x| format!("{}", x))
467 .collect();
468 s.push_str(&vals.join(" "));
469 s.push_str("\n </DataArray>\n");
470 }
471 s.push_str(" </CellData>\n");
472 }
473 s.push_str(" <Points>\n");
474 s.push_str(
475 " <DataArray type=\"Float64\" NumberOfComponents=\"3\" format=\"ascii\">\n",
476 );
477 for p in &self.points {
478 s.push_str(&format!(" {} {} {}\n", p[0], p[1], p[2]));
479 }
480 s.push_str(" </DataArray>\n");
481 s.push_str(" </Points>\n");
482 s.push_str(" <Cells>\n");
483 s.push_str(
484 " <DataArray type=\"Int32\" Name=\"connectivity\" format=\"ascii\">\n ",
485 );
486 let conn: Vec<String> = self
487 .cells
488 .iter()
489 .flat_map(|c| c.iter())
490 .map(|i| format!("{}", i))
491 .collect();
492 s.push_str(&conn.join(" "));
493 s.push_str("\n </DataArray>\n");
494 s.push_str(
495 " <DataArray type=\"Int32\" Name=\"offsets\" format=\"ascii\">\n ",
496 );
497 let mut offset = 0usize;
498 let offsets: Vec<String> = self
499 .cells
500 .iter()
501 .map(|c| {
502 offset += c.len();
503 format!("{}", offset)
504 })
505 .collect();
506 s.push_str(&offsets.join(" "));
507 s.push_str("\n </DataArray>\n");
508 s.push_str(
509 " <DataArray type=\"UInt8\" Name=\"types\" format=\"ascii\">\n ",
510 );
511 let type_strs: Vec<String> = self.cell_types.iter().map(|t| format!("{}", t)).collect();
512 s.push_str(&type_strs.join(" "));
513 s.push_str("\n </DataArray>\n");
514 s.push_str(" </Cells>\n");
515 s.push_str(" </Piece>\n");
516 s.push_str(" </UnstructuredGrid>\n");
517 s.push_str("</VTKFile>\n");
518 s
519 }
520 pub fn write_to_file(&self, path: &str) -> Result<()> {
522 let file = File::create(Path::new(path))?;
523 let mut w = BufWriter::new(file);
524 write!(w, "{}", self.to_xml())?;
525 w.flush()?;
526 Ok(())
527 }
528 pub fn validate(&self) -> Result<()> {
531 let np = self.points.len();
532 for (ci, cell) in self.cells.iter().enumerate() {
533 for &idx in cell {
534 if idx >= np {
535 return Err(Error::Parse(format!(
536 "cell {} references point index {} but only {} points exist",
537 ci, idx, np
538 )));
539 }
540 }
541 }
542 for ps in &self.point_scalars {
543 if ps.values.len() != np {
544 return Err(Error::Parse(format!(
545 "point scalar '{}' has {} values but {} points exist",
546 ps.name,
547 ps.values.len(),
548 np
549 )));
550 }
551 }
552 for pv in &self.point_vectors {
553 if pv.vectors.len() != np {
554 return Err(Error::Parse(format!(
555 "point vector '{}' has {} vectors but {} points exist",
556 pv.name,
557 pv.vectors.len(),
558 np
559 )));
560 }
561 }
562 let nc = self.cells.len();
563 for cs in &self.cell_scalars {
564 if cs.values.len() != nc {
565 return Err(Error::Parse(format!(
566 "cell scalar '{}' has {} values but {} cells exist",
567 cs.name,
568 cs.values.len(),
569 nc
570 )));
571 }
572 }
573 for cv in &self.cell_vectors {
574 if cv.vectors.len() != nc {
575 return Err(Error::Parse(format!(
576 "cell vector '{}' has {} vectors but {} cells exist",
577 cv.name,
578 cv.vectors.len(),
579 nc
580 )));
581 }
582 }
583 Ok(())
584 }
585 pub fn points(&self) -> &[[f64; 3]] {
587 &self.points
588 }
589 pub fn cell_types(&self) -> &[u8] {
591 &self.cell_types
592 }
593 pub fn cells_per_cell(&self) -> &[Vec<usize>] {
595 &self.cells
596 }
597 pub fn cells_flat(&self) -> Vec<usize> {
599 self.cells.iter().flat_map(|c| c.iter().cloned()).collect()
600 }
601}
602#[allow(dead_code)]
604pub struct VtuFieldOps;
605#[allow(dead_code)]
606impl VtuFieldOps {
607 pub fn add(a: &[f64], b: &[f64]) -> Vec<f64> {
609 a.iter().zip(b.iter()).map(|(x, y)| x + y).collect()
610 }
611 pub fn sub(a: &[f64], b: &[f64]) -> Vec<f64> {
613 a.iter().zip(b.iter()).map(|(x, y)| x - y).collect()
614 }
615 pub fn mul(a: &[f64], b: &[f64]) -> Vec<f64> {
617 a.iter().zip(b.iter()).map(|(x, y)| x * y).collect()
618 }
619 pub fn scale(a: &[f64], s: f64) -> Vec<f64> {
621 a.iter().map(|x| x * s).collect()
622 }
623 pub fn vector_magnitude(v: &[f64]) -> Vec<f64> {
625 v.chunks(3)
626 .map(|c| {
627 if c.len() == 3 {
628 (c[0] * c[0] + c[1] * c[1] + c[2] * c[2]).sqrt()
629 } else {
630 0.0
631 }
632 })
633 .collect()
634 }
635 pub fn normalize_vectors(v: &[f64]) -> Vec<f64> {
637 let mut out = Vec::with_capacity(v.len());
638 for c in v.chunks(3) {
639 if c.len() == 3 {
640 let mag = (c[0] * c[0] + c[1] * c[1] + c[2] * c[2]).sqrt();
641 if mag > 1e-30 {
642 out.push(c[0] / mag);
643 out.push(c[1] / mag);
644 out.push(c[2] / mag);
645 } else {
646 out.push(0.0);
647 out.push(0.0);
648 out.push(1.0);
649 }
650 }
651 }
652 out
653 }
654 pub fn clamp(a: &[f64], lo: f64, hi: f64) -> Vec<f64> {
656 a.iter().map(|&x| x.max(lo).min(hi)).collect()
657 }
658 pub fn dot_product(a: &[f64], b: &[f64]) -> Vec<f64> {
660 a.chunks(3)
661 .zip(b.chunks(3))
662 .map(|(u, v)| {
663 if u.len() == 3 && v.len() == 3 {
664 u[0] * v[0] + u[1] * v[1] + u[2] * v[2]
665 } else {
666 0.0
667 }
668 })
669 .collect()
670 }
671 pub fn cross_product(a: &[f64], b: &[f64]) -> Vec<f64> {
673 let mut out = Vec::with_capacity(a.len());
674 for (u, v) in a.chunks(3).zip(b.chunks(3)) {
675 if u.len() == 3 && v.len() == 3 {
676 out.push(u[1] * v[2] - u[2] * v[1]);
677 out.push(u[2] * v[0] - u[0] * v[2]);
678 out.push(u[0] * v[1] - u[1] * v[0]);
679 }
680 }
681 out
682 }
683 pub fn stats(a: &[f64]) -> (f64, f64, f64) {
685 if a.is_empty() {
686 return (0.0, 0.0, 0.0);
687 }
688 let mut min_v = a[0];
689 let mut max_v = a[0];
690 let mut sum = 0.0_f64;
691 for &v in a {
692 if v < min_v {
693 min_v = v;
694 }
695 if v > max_v {
696 max_v = v;
697 }
698 sum += v;
699 }
700 (min_v, max_v, sum / a.len() as f64)
701 }
702 pub fn threshold(a: &[f64], threshold: f64) -> Vec<f64> {
704 a.iter()
705 .map(|&v| if v > threshold { 1.0 } else { 0.0 })
706 .collect()
707 }
708}
709#[allow(dead_code)]
710pub(super) struct PointScalar {
711 pub(super) name: String,
712 pub(super) values: Vec<f64>,
713}
714#[allow(dead_code)]
715pub(super) struct CellScalar {
716 pub(super) name: String,
717 pub(super) values: Vec<f64>,
718}
719#[allow(dead_code)]
721pub struct VtuReader {
722 pub(super) points: Vec<[f64; 3]>,
723 pub(super) num_cells: usize,
724 pub(super) cell_types: Vec<u8>,
725 pub(super) cell_connectivity: Vec<Vec<usize>>,
726 pub(super) point_scalar_names: Vec<String>,
727 pub(super) cell_scalar_names: Vec<String>,
728}
729#[allow(dead_code)]
730impl VtuReader {
731 pub fn from_xml(data: &str) -> Result<Self> {
733 let mut points = Vec::new();
734 let mut in_points_array = false;
735 let mut cell_types = Vec::new();
736 let mut cell_connectivity = Vec::new();
737 let mut num_cells = 0usize;
738 let mut point_scalar_names = Vec::new();
739 let mut cell_scalar_names = Vec::new();
740 let mut in_cell_types = false;
741 let mut in_connectivity = false;
742 let mut in_point_data = false;
743 let mut in_cell_data = false;
744 for line in data.lines() {
745 let trimmed = line.trim();
746 if trimmed.starts_with("<PointData") {
747 in_point_data = true;
748 continue;
749 }
750 if trimmed.starts_with("</PointData>") {
751 in_point_data = false;
752 continue;
753 }
754 if trimmed.starts_with("<CellData") {
755 in_cell_data = true;
756 continue;
757 }
758 if trimmed.starts_with("</CellData>") {
759 in_cell_data = false;
760 continue;
761 }
762 if (in_point_data || in_cell_data)
763 && trimmed.contains("Name=\"")
764 && let Some(start) = trimmed.find("Name=\"")
765 {
766 let rest = &trimmed[start + 6..];
767 if let Some(end) = rest.find('"') {
768 let name = rest[..end].to_string();
769 if in_point_data {
770 point_scalar_names.push(name);
771 } else {
772 cell_scalar_names.push(name);
773 }
774 }
775 }
776 if trimmed.contains("NumberOfCells=")
777 && let Some(start) = trimmed.find("NumberOfCells=\"")
778 {
779 let rest = &trimmed[start + 15..];
780 if let Some(end) = rest.find('"') {
781 num_cells = rest[..end].parse::<usize>().unwrap_or(0);
782 }
783 }
784 if trimmed.contains("NumberOfComponents=\"3\"")
785 && trimmed.contains("Float64")
786 && !in_point_data
787 && !in_cell_data
788 {
789 in_points_array = true;
790 continue;
791 }
792 if in_points_array && trimmed.starts_with("</DataArray>") {
793 in_points_array = false;
794 continue;
795 }
796 if in_points_array && !trimmed.starts_with('<') && !trimmed.is_empty() {
797 let nums: Vec<&str> = trimmed.split_whitespace().collect();
798 if nums.len() >= 3 {
799 let mut i = 0;
800 while i + 2 < nums.len() {
801 let x = nums[i]
802 .parse::<f64>()
803 .map_err(|e| Error::Parse(e.to_string()))?;
804 let y = nums[i + 1]
805 .parse::<f64>()
806 .map_err(|e| Error::Parse(e.to_string()))?;
807 let z = nums[i + 2]
808 .parse::<f64>()
809 .map_err(|e| Error::Parse(e.to_string()))?;
810 points.push([x, y, z]);
811 i += 3;
812 }
813 }
814 }
815 if trimmed.contains("Name=\"types\"") {
816 in_cell_types = true;
817 continue;
818 }
819 if in_cell_types && trimmed.starts_with("</DataArray>") {
820 in_cell_types = false;
821 continue;
822 }
823 if in_cell_types && !trimmed.starts_with('<') && !trimmed.is_empty() {
824 for tok in trimmed.split_whitespace() {
825 if let Ok(t) = tok.parse::<u8>() {
826 cell_types.push(t);
827 }
828 }
829 }
830 if trimmed.contains("Name=\"connectivity\"") {
831 in_connectivity = true;
832 continue;
833 }
834 if in_connectivity && trimmed.starts_with("</DataArray>") {
835 in_connectivity = false;
836 continue;
837 }
838 if in_connectivity && !trimmed.starts_with('<') && !trimmed.is_empty() {
839 let indices: Vec<usize> = trimmed
840 .split_whitespace()
841 .filter_map(|s| s.parse::<usize>().ok())
842 .collect();
843 if !indices.is_empty() {
844 cell_connectivity.push(indices);
845 }
846 }
847 }
848 Ok(Self {
849 points,
850 num_cells,
851 cell_types,
852 cell_connectivity,
853 point_scalar_names,
854 cell_scalar_names,
855 })
856 }
857 pub fn points(&self) -> &[[f64; 3]] {
859 &self.points
860 }
861 pub fn num_cells(&self) -> usize {
863 self.num_cells
864 }
865 pub fn cell_types(&self) -> &[u8] {
867 &self.cell_types
868 }
869 pub fn point_data_names(&self) -> &[String] {
871 &self.point_scalar_names
872 }
873 pub fn cell_data_names(&self) -> &[String] {
875 &self.cell_scalar_names
876 }
877 #[allow(dead_code)]
882 pub fn read_cell_data(data: &str, array_name: &str) -> std::result::Result<Vec<f64>, String> {
883 let mut in_cell_data = false;
884 let mut target_found = false;
885 let mut collecting = false;
886 let mut values = Vec::new();
887 for line in data.lines() {
888 let trimmed = line.trim();
889 if trimmed.starts_with("<CellData") {
890 in_cell_data = true;
891 continue;
892 }
893 if trimmed.starts_with("</CellData>") {
894 in_cell_data = false;
895 target_found = false;
896 collecting = false;
897 continue;
898 }
899 if in_cell_data && trimmed.contains("Name=\"") {
900 if let Some(start) = trimmed.find("Name=\"") {
901 let rest = &trimmed[start + 6..];
902 if let Some(end) = rest.find('"') {
903 let name = &rest[..end];
904 target_found = name == array_name;
905 }
906 }
907 if target_found {
908 collecting = true;
909 }
910 continue;
911 }
912 if collecting && trimmed.starts_with("</DataArray>") {
913 collecting = false;
914 continue;
915 }
916 if collecting && !trimmed.starts_with('<') && !trimmed.is_empty() {
917 for tok in trimmed.split_whitespace() {
918 if let Ok(v) = tok.parse::<f64>() {
919 values.push(v);
920 }
921 }
922 }
923 }
924 if values.is_empty() && !data.contains(&format!("Name=\"{}\"", array_name)) {
925 Err(format!("cell data array '{}' not found", array_name))
926 } else {
927 Ok(values)
928 }
929 }
930 pub fn bounding_box(&self) -> Option<([f64; 3], [f64; 3])> {
932 if self.points.is_empty() {
933 return None;
934 }
935 let mut min = self.points[0];
936 let mut max = self.points[0];
937 for p in &self.points[1..] {
938 for d in 0..3 {
939 if p[d] < min[d] {
940 min[d] = p[d];
941 }
942 if p[d] > max[d] {
943 max[d] = p[d];
944 }
945 }
946 }
947 Some((min, max))
948 }
949}
950#[allow(dead_code)]
954pub struct PvtuWriter;
955#[allow(dead_code)]
956impl PvtuWriter {
957 pub fn write(
962 piece_files: &[&str],
963 point_data_names: &[&str],
964 cell_data_names: &[&str],
965 ) -> String {
966 let mut s = String::new();
967 s.push_str("<?xml version=\"1.0\"?>\n");
968 s.push_str(
969 "<VTKFile type=\"PUnstructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">\n",
970 );
971 s.push_str(" <PUnstructuredGrid>\n");
972 if !point_data_names.is_empty() {
973 s.push_str(" <PPointData>\n");
974 for name in point_data_names {
975 s.push_str(&format!(
976 " <PDataArray type=\"Float64\" Name=\"{}\" NumberOfComponents=\"1\"/>\n",
977 name
978 ));
979 }
980 s.push_str(" </PPointData>\n");
981 }
982 if !cell_data_names.is_empty() {
983 s.push_str(" <PCellData>\n");
984 for name in cell_data_names {
985 s.push_str(&format!(
986 " <PDataArray type=\"Float64\" Name=\"{}\" NumberOfComponents=\"1\"/>\n",
987 name
988 ));
989 }
990 s.push_str(" </PCellData>\n");
991 }
992 s.push_str(" <PPoints>\n");
993 s.push_str(" <PDataArray type=\"Float64\" NumberOfComponents=\"3\"/>\n");
994 s.push_str(" </PPoints>\n");
995 for path in piece_files {
996 s.push_str(&format!(" <Piece Source=\"{}\"/>\n", path));
997 }
998 s.push_str(" </PUnstructuredGrid>\n");
999 s.push_str("</VTKFile>\n");
1000 s
1001 }
1002}
1003#[allow(dead_code)]
1005pub struct VtuTimeSeries {
1006 pub(super) prefix: String,
1007 pub(super) directory: String,
1008 pub(super) timesteps: Vec<(f64, String)>,
1009}
1010#[allow(dead_code)]
1011impl VtuTimeSeries {
1012 pub fn new(directory: &str, prefix: &str) -> Self {
1017 Self {
1018 prefix: prefix.to_string(),
1019 directory: directory.to_string(),
1020 timesteps: Vec::new(),
1021 }
1022 }
1023 pub fn write_frame(&mut self, time: f64, writer: &VtuWriter) -> Result<String> {
1025 let frame_idx = self.timesteps.len();
1026 let filename = format!("{}_{:04}.vtu", self.prefix, frame_idx);
1027 let filepath = format!("{}/{}", self.directory, filename);
1028 writer.write_to_file(&filepath)?;
1029 self.timesteps.push((time, filename));
1030 Ok(filepath)
1031 }
1032 pub fn write_pvd(&self) -> Result<String> {
1034 let pvd_path = format!("{}/{}.pvd", self.directory, self.prefix);
1035 let entries: Vec<(f64, &str)> = self
1036 .timesteps
1037 .iter()
1038 .map(|(t, f)| (*t, f.as_str()))
1039 .collect();
1040 let pvd_xml = PvdWriter::write(&entries);
1041 let file = File::create(Path::new(&pvd_path))?;
1042 let mut w = BufWriter::new(file);
1043 write!(w, "{}", pvd_xml)?;
1044 w.flush()?;
1045 Ok(pvd_path)
1046 }
1047 pub fn num_frames(&self) -> usize {
1049 self.timesteps.len()
1050 }
1051 pub fn times(&self) -> Vec<f64> {
1053 self.timesteps.iter().map(|(t, _)| *t).collect()
1054 }
1055}
1056#[allow(dead_code)]
1058#[derive(Debug, Clone)]
1059pub struct VtuCellQuality {
1060 pub cell_index: usize,
1062 pub cell_type: u8,
1064 pub aspect_ratio: f64,
1066 pub min_edge: f64,
1068 pub max_edge: f64,
1070 pub scaled_jacobian: f64,
1072}
1073#[allow(dead_code)]
1078pub struct PvdWriter;
1079#[allow(dead_code)]
1080impl PvdWriter {
1081 pub fn write(timesteps: &[(f64, &str)]) -> String {
1085 let mut s = String::new();
1086 s.push_str("<?xml version=\"1.0\"?>\n");
1087 s.push_str("<VTKFile type=\"Collection\" version=\"0.1\">\n");
1088 s.push_str(" <Collection>\n");
1089 for (t, path) in timesteps {
1090 s.push_str(&format!(
1091 " <DataSet timestep=\"{}\" part=\"0\" file=\"{}\"/>\n",
1092 t, path
1093 ));
1094 }
1095 s.push_str(" </Collection>\n");
1096 s.push_str("</VTKFile>\n");
1097 s
1098 }
1099 pub fn write_parallel(timesteps: &[(f64, Vec<&str>)]) -> String {
1103 let mut s = String::new();
1104 s.push_str("<?xml version=\"1.0\"?>\n");
1105 s.push_str("<VTKFile type=\"Collection\" version=\"0.1\">\n");
1106 s.push_str(" <Collection>\n");
1107 for (t, paths) in timesteps {
1108 for (part, path) in paths.iter().enumerate() {
1109 s.push_str(&format!(
1110 " <DataSet timestep=\"{}\" part=\"{}\" file=\"{}\"/>\n",
1111 t, part, path
1112 ));
1113 }
1114 }
1115 s.push_str(" </Collection>\n");
1116 s.push_str("</VTKFile>\n");
1117 s
1118 }
1119}