1#![allow(clippy::manual_div_ceil, clippy::too_many_arguments)]
6#[allow(unused_imports)]
7use super::functions::*;
8#[allow(unused_imports)]
9use super::functions_2::*;
10use std::io::Write;
11
12#[allow(dead_code)]
14#[derive(Debug, Clone, Copy, PartialEq)]
15pub enum XdmfTopologyType {
16 Triangle,
18 Tetrahedron,
20 Hexahedron,
22 Quadrilateral,
24 Mixed,
26}
27impl XdmfTopologyType {
28 pub fn xdmf_name(self) -> &'static str {
30 match self {
31 XdmfTopologyType::Triangle => "Triangle",
32 XdmfTopologyType::Tetrahedron => "Tetrahedron",
33 XdmfTopologyType::Hexahedron => "Hexahedron",
34 XdmfTopologyType::Quadrilateral => "Quadrilateral",
35 XdmfTopologyType::Mixed => "Mixed",
36 }
37 }
38 pub fn nodes_per_element(self) -> usize {
40 match self {
41 XdmfTopologyType::Triangle => 3,
42 XdmfTopologyType::Tetrahedron => 4,
43 XdmfTopologyType::Hexahedron => 8,
44 XdmfTopologyType::Quadrilateral => 4,
45 XdmfTopologyType::Mixed => 0,
46 }
47 }
48}
49#[derive(Debug, Clone)]
51#[allow(dead_code)]
52pub struct XdmfStep {
53 pub time: f64,
55 pub n_points: usize,
57 pub position_data: Vec<[f64; 3]>,
59 pub scalar_fields: Vec<(String, Vec<f64>)>,
61}
62#[derive(Debug, Clone)]
64#[allow(dead_code)]
65pub struct XdmfUniformGrid {
66 pub name: String,
68 pub dimensions: [usize; 3],
70 pub origin: [f64; 3],
72 pub spacing: [f64; 3],
74}
75#[derive(Debug, Clone, Default)]
78#[allow(dead_code)]
79pub struct XdmfMultiBlock {
80 pub(super) blocks: Vec<(String, String)>,
83}
84#[allow(dead_code)]
85impl XdmfMultiBlock {
86 pub fn new() -> Self {
88 Self::default()
89 }
90 pub fn add_block(&mut self, name: &str, grid_xml: &str) {
92 self.blocks.push((name.to_string(), grid_xml.to_string()));
93 }
94 pub fn len(&self) -> usize {
96 self.blocks.len()
97 }
98 pub fn is_empty(&self) -> bool {
100 self.blocks.is_empty()
101 }
102 pub fn to_xml(&self) -> String {
104 let mut s = String::new();
105 s.push_str("<?xml version=\"1.0\"?>\n");
106 s.push_str("<Xdmf Version=\"3.0\">\n");
107 s.push_str(" <Domain>\n");
108 s.push_str(
109 " <Grid Name=\"MultiBlock\" GridType=\"Collection\" CollectionType=\"Spatial\">\n",
110 );
111 for (_name, grid) in &self.blocks {
112 for line in grid.lines() {
113 s.push_str(" ");
114 s.push_str(line);
115 s.push('\n');
116 }
117 }
118 s.push_str(" </Grid>\n");
119 s.push_str(" </Domain>\n");
120 s.push_str("</Xdmf>\n");
121 s
122 }
123}
124#[derive(Debug, Clone)]
126#[allow(dead_code)]
127pub struct XdmfTimeSeriesHdf5 {
128 pub timesteps: Vec<f64>,
130 pub hdf5_paths: Vec<String>,
132 pub attribute_names: Vec<String>,
134}
135#[allow(dead_code)]
136impl XdmfTimeSeriesHdf5 {
137 pub fn new() -> Self {
139 Self {
140 timesteps: Vec::new(),
141 hdf5_paths: Vec::new(),
142 attribute_names: Vec::new(),
143 }
144 }
145 pub fn write_collection(
150 &self,
151 path: &str,
152 n_nodes: usize,
153 n_elements: usize,
154 topology: &str,
155 ) -> std::io::Result<()> {
156 let mut f = std::fs::File::create(path)?;
157 writeln!(f, "<?xml version=\"1.0\"?>")?;
158 writeln!(f, "<Xdmf Version=\"3.0\">")?;
159 writeln!(f, " <Domain>")?;
160 writeln!(
161 f,
162 " <Grid Name=\"TimeSeries\" GridType=\"Collection\" CollectionType=\"Temporal\">"
163 )?;
164 for (step_idx, &t) in self.timesteps.iter().enumerate() {
165 let h5path = self
166 .hdf5_paths
167 .get(step_idx)
168 .map(|s| s.as_str())
169 .unwrap_or("data.h5");
170 writeln!(
171 f,
172 " <Grid Name=\"step{}\" GridType=\"Uniform\">",
173 step_idx
174 )?;
175 writeln!(f, " <Time Value=\"{}\"/>", t)?;
176 writeln!(
177 f,
178 " <Topology TopologyType=\"{}\" NumberOfElements=\"{}\"/>",
179 topology, n_elements
180 )?;
181 writeln!(f, " <Geometry GeometryType=\"XYZ\">")?;
182 writeln!(
183 f,
184 " <DataItem Format=\"HDF\" Dimensions=\"{} 3\">{}:/coordinates</DataItem>",
185 n_nodes, h5path
186 )?;
187 writeln!(f, " </Geometry>")?;
188 for attr in &self.attribute_names {
189 writeln!(
190 f,
191 " <Attribute Name=\"{}\" AttributeType=\"Scalar\" Center=\"Node\">",
192 attr
193 )?;
194 writeln!(
195 f,
196 " <DataItem Format=\"HDF\" Dimensions=\"{}\">{}/{}</DataItem>",
197 n_nodes, h5path, attr
198 )?;
199 writeln!(f, " </Attribute>")?;
200 }
201 writeln!(f, " </Grid>")?;
202 }
203 writeln!(f, " </Grid>")?;
204 writeln!(f, " </Domain>")?;
205 writeln!(f, "</Xdmf>")?;
206 Ok(())
207 }
208}
209#[derive(Debug, Clone)]
213#[allow(dead_code)]
214pub struct XdmfStructuredGrid {
215 pub name: String,
217 pub ni: usize,
219 pub nj: usize,
221 pub nk: usize,
223 pub origin: [f64; 3],
225 pub dx: f64,
227 pub dy: f64,
229 pub dz: f64,
231 pub node_scalars: Vec<(String, Vec<f64>)>,
233 pub node_vectors: Vec<(String, Vec<[f64; 3]>)>,
235 pub cell_scalars: Vec<(String, Vec<f64>)>,
237}
238#[allow(dead_code)]
239impl XdmfStructuredGrid {
240 pub fn new(
242 name: &str,
243 ni: usize,
244 nj: usize,
245 nk: usize,
246 origin: [f64; 3],
247 dx: f64,
248 dy: f64,
249 dz: f64,
250 ) -> Self {
251 Self {
252 name: name.to_string(),
253 ni,
254 nj,
255 nk,
256 origin,
257 dx,
258 dy,
259 dz,
260 node_scalars: Vec::new(),
261 node_vectors: Vec::new(),
262 cell_scalars: Vec::new(),
263 }
264 }
265 pub fn n_nodes(&self) -> usize {
267 self.ni * self.nj * self.nk
268 }
269 pub fn n_cells(&self) -> usize {
271 (self.ni.saturating_sub(1)) * (self.nj.saturating_sub(1)) * (self.nk.saturating_sub(1))
272 }
273 pub fn add_node_scalar(&mut self, name: &str, data: Vec<f64>) {
275 self.node_scalars.push((name.to_string(), data));
276 }
277 pub fn add_node_vector(&mut self, name: &str, data: Vec<[f64; 3]>) {
279 self.node_vectors.push((name.to_string(), data));
280 }
281 pub fn add_cell_scalar(&mut self, name: &str, data: Vec<f64>) {
283 self.cell_scalars.push((name.to_string(), data));
284 }
285 pub fn to_xml(&self) -> String {
287 let mut s = String::new();
288 s.push_str("<?xml version=\"1.0\"?>\n");
289 s.push_str("<Xdmf Version=\"3.0\">\n");
290 s.push_str(" <Domain>\n");
291 s.push_str(&format!(
292 " <Grid Name=\"{}\" GridType=\"Uniform\">\n",
293 self.name
294 ));
295 s.push_str(&format!(
296 " <Topology TopologyType=\"3DCoRectMesh\" Dimensions=\"{} {} {}\"/>\n",
297 self.nk, self.nj, self.ni
298 ));
299 s.push_str(" <Geometry GeometryType=\"ORIGIN_DXDYDZ\">\n");
300 s.push_str(&format!(
301 " <DataItem Format=\"XML\" Dimensions=\"3\">{} {} {}</DataItem>\n",
302 self.origin[2], self.origin[1], self.origin[0]
303 ));
304 s.push_str(&format!(
305 " <DataItem Format=\"XML\" Dimensions=\"3\">{} {} {}</DataItem>\n",
306 self.dz, self.dy, self.dx
307 ));
308 s.push_str(" </Geometry>\n");
309 for (name, data) in &self.node_scalars {
310 s.push_str(&format!(
311 " <Attribute Name=\"{}\" AttributeType=\"Scalar\" Center=\"Node\">\n",
312 name
313 ));
314 s.push_str(&format!(
315 " <DataItem Format=\"XML\" Dimensions=\"{}\">\n ",
316 data.len()
317 ));
318 for (i, v) in data.iter().enumerate() {
319 if i > 0 {
320 s.push(' ');
321 }
322 s.push_str(&format!("{}", v));
323 }
324 s.push_str("\n </DataItem>\n </Attribute>\n");
325 }
326 for (name, data) in &self.node_vectors {
327 s.push_str(&format!(
328 " <Attribute Name=\"{}\" AttributeType=\"Vector\" Center=\"Node\">\n",
329 name
330 ));
331 s.push_str(&format!(
332 " <DataItem Format=\"XML\" Dimensions=\"{} 3\">\n",
333 data.len()
334 ));
335 for v in data {
336 s.push_str(&format!(" {} {} {}\n", v[0], v[1], v[2]));
337 }
338 s.push_str(" </DataItem>\n </Attribute>\n");
339 }
340 for (name, data) in &self.cell_scalars {
341 s.push_str(&format!(
342 " <Attribute Name=\"{}\" AttributeType=\"Scalar\" Center=\"Cell\">\n",
343 name
344 ));
345 s.push_str(&format!(
346 " <DataItem Format=\"XML\" Dimensions=\"{}\">\n ",
347 data.len()
348 ));
349 for (i, v) in data.iter().enumerate() {
350 if i > 0 {
351 s.push(' ');
352 }
353 s.push_str(&format!("{}", v));
354 }
355 s.push_str("\n </DataItem>\n </Attribute>\n");
356 }
357 s.push_str(" </Grid>\n");
358 s.push_str(" </Domain>\n");
359 s.push_str("</Xdmf>\n");
360 s
361 }
362}
363#[allow(dead_code)]
368pub struct XdmfWriter {
369 pub(super) buf: String,
370 pub(super) indent: usize,
371}
372#[allow(dead_code)]
373impl XdmfWriter {
374 pub fn new() -> Self {
376 let mut w = XdmfWriter {
377 buf: String::new(),
378 indent: 0,
379 };
380 w.buf.push_str("<?xml version=\"1.0\"?>\n");
381 w.buf.push_str("<Xdmf Version=\"3.0\">\n");
382 w.indent = 1;
383 w
384 }
385 fn push_line(&mut self, s: &str) {
386 let pad = " ".repeat(self.indent);
387 self.buf.push_str(&pad);
388 self.buf.push_str(s);
389 self.buf.push('\n');
390 }
391 pub fn open_domain(&mut self) {
393 self.push_line("<Domain>");
394 self.indent += 1;
395 }
396 pub fn close_domain(&mut self) {
398 self.indent = self.indent.saturating_sub(1);
399 self.push_line("</Domain>");
400 }
401 pub fn open_grid(&mut self, name: &str, grid_type: &str) {
403 self.push_line(&format!(
404 "<Grid Name=\"{}\" GridType=\"{}\">",
405 name, grid_type
406 ));
407 self.indent += 1;
408 }
409 pub fn open_temporal_collection(&mut self, name: &str) {
411 self.push_line(&format!(
412 "<Grid Name=\"{}\" GridType=\"Collection\" CollectionType=\"Temporal\">",
413 name
414 ));
415 self.indent += 1;
416 }
417 pub fn write_time(&mut self, t: f64) {
419 self.push_line(&format!("<Time Value=\"{}\"/>", t));
420 }
421 pub fn close_grid(&mut self) {
423 self.indent = self.indent.saturating_sub(1);
424 self.push_line("</Grid>");
425 }
426 pub fn write_polyvertex_topology(&mut self, n: usize) {
428 self.push_line(&format!(
429 "<Topology TopologyType=\"Polyvertex\" NumberOfElements=\"{}\"/>",
430 n
431 ));
432 }
433 pub fn write_unstructured_topology(
435 &mut self,
436 topo_type: &str,
437 n_elements: usize,
438 connectivity: &[usize],
439 npe: usize,
440 ) {
441 self.push_line(&format!(
442 "<Topology TopologyType=\"{}\" NumberOfElements=\"{}\">",
443 topo_type, n_elements
444 ));
445 self.indent += 1;
446 self.push_line(&format!(
447 "<DataItem Format=\"XML\" Dimensions=\"{} {}\">\n{} ",
448 n_elements,
449 npe,
450 " ".repeat(self.indent)
451 ));
452 let row_strings: Vec<String> = connectivity
453 .chunks(npe)
454 .map(|chunk| {
455 chunk
456 .iter()
457 .map(|i| i.to_string())
458 .collect::<Vec<_>>()
459 .join(" ")
460 })
461 .collect();
462 let pad = " ".repeat(self.indent);
463 for r in &row_strings {
464 self.buf.push_str(&pad);
465 self.buf.push_str(r);
466 self.buf.push('\n');
467 }
468 self.push_line("</DataItem>");
469 self.indent = self.indent.saturating_sub(1);
470 self.push_line("</Topology>");
471 }
472 pub fn write_xyz_geometry(&mut self, nodes: &[[f64; 3]]) {
474 self.push_line("<Geometry GeometryType=\"XYZ\">");
475 self.indent += 1;
476 self.push_line(&format!(
477 "<DataItem Format=\"XML\" Dimensions=\"{} 3\">",
478 nodes.len()
479 ));
480 self.indent += 1;
481 let pad = " ".repeat(self.indent);
482 for p in nodes {
483 self.buf.push_str(&pad);
484 self.buf.push_str(&format!("{} {} {}\n", p[0], p[1], p[2]));
485 }
486 self.indent = self.indent.saturating_sub(1);
487 self.push_line("</DataItem>");
488 self.indent = self.indent.saturating_sub(1);
489 self.push_line("</Geometry>");
490 }
491 pub fn write_scalar_attribute(&mut self, name: &str, center: &str, values: &[f64]) {
493 self.push_line(&format!(
494 "<Attribute Name=\"{}\" AttributeType=\"Scalar\" Center=\"{}\">",
495 name, center
496 ));
497 self.indent += 1;
498 self.push_line(&format!(
499 "<DataItem Format=\"XML\" Dimensions=\"{}\">",
500 values.len()
501 ));
502 self.indent += 1;
503 let pad = " ".repeat(self.indent);
504 self.buf.push_str(&pad);
505 for (i, v) in values.iter().enumerate() {
506 if i > 0 {
507 self.buf.push(' ');
508 }
509 self.buf.push_str(&format!("{}", v));
510 }
511 self.buf.push('\n');
512 self.indent = self.indent.saturating_sub(1);
513 self.push_line("</DataItem>");
514 self.indent = self.indent.saturating_sub(1);
515 self.push_line("</Attribute>");
516 }
517 pub fn write_vector_attribute(&mut self, name: &str, center: &str, vectors: &[[f64; 3]]) {
519 self.push_line(&format!(
520 "<Attribute Name=\"{}\" AttributeType=\"Vector\" Center=\"{}\">",
521 name, center
522 ));
523 self.indent += 1;
524 self.push_line(&format!(
525 "<DataItem Format=\"XML\" Dimensions=\"{} 3\">",
526 vectors.len()
527 ));
528 self.indent += 1;
529 let pad = " ".repeat(self.indent);
530 for v in vectors {
531 self.buf.push_str(&pad);
532 self.buf.push_str(&format!("{} {} {}\n", v[0], v[1], v[2]));
533 }
534 self.indent = self.indent.saturating_sub(1);
535 self.push_line("</DataItem>");
536 self.indent = self.indent.saturating_sub(1);
537 self.push_line("</Attribute>");
538 }
539 pub fn write_hdf5_attribute(
541 &mut self,
542 name: &str,
543 center: &str,
544 attr_type: &str,
545 dims: &str,
546 hdf5_ref: &str,
547 ) {
548 self.push_line(&format!(
549 "<Attribute Name=\"{}\" AttributeType=\"{}\" Center=\"{}\">",
550 name, attr_type, center
551 ));
552 self.indent += 1;
553 self.push_line(&format!(
554 "<DataItem Format=\"HDF\" Dimensions=\"{}\">{}</DataItem>",
555 dims, hdf5_ref
556 ));
557 self.indent = self.indent.saturating_sub(1);
558 self.push_line("</Attribute>");
559 }
560 pub fn finish(mut self) -> String {
562 self.buf.push_str("</Xdmf>\n");
563 self.buf
564 }
565 pub fn peek(&self) -> &str {
567 &self.buf
568 }
569}
570#[allow(dead_code)]
574pub struct Hdf5DataItemBuilder {
575 pub(super) filename: String,
576 pub(super) group: String,
577}
578#[allow(dead_code)]
579impl Hdf5DataItemBuilder {
580 pub fn new(filename: &str) -> Self {
582 Self {
583 filename: filename.to_owned(),
584 group: "/".to_owned(),
585 }
586 }
587 pub fn group(mut self, group: &str) -> Self {
589 self.group = group.to_owned();
590 self
591 }
592 pub fn build(&self, dataset_name: &str, dimensions: &str) -> String {
595 let path = if self.group == "/" {
596 format!("/{}", dataset_name)
597 } else {
598 format!("{}/{}", self.group, dataset_name)
599 };
600 format!(
601 "<DataItem Format=\"HDF\" Dimensions=\"{}\">\n {}:{}\n</DataItem>",
602 dimensions, self.filename, path
603 )
604 }
605}
606#[derive(Debug, Clone)]
608#[allow(dead_code)]
609pub struct XdmfGridEntry {
610 pub name: String,
612 pub time: f64,
614 pub nodes: Vec<[f64; 3]>,
616 pub connectivity: Vec<usize>,
618 pub topology_type: String,
620}
621impl XdmfGridEntry {
622 #[allow(dead_code)]
624 pub fn node_count(&self) -> usize {
625 self.nodes.len()
626 }
627 #[allow(dead_code)]
630 pub fn bounding_box(&self) -> ([f64; 3], [f64; 3]) {
631 let mut lo = [f64::INFINITY; 3];
632 let mut hi = [f64::NEG_INFINITY; 3];
633 for n in &self.nodes {
634 for k in 0..3 {
635 if n[k] < lo[k] {
636 lo[k] = n[k];
637 }
638 if n[k] > hi[k] {
639 hi[k] = n[k];
640 }
641 }
642 }
643 (lo, hi)
644 }
645 #[allow(dead_code)]
647 pub fn centroid(&self) -> [f64; 3] {
648 if self.nodes.is_empty() {
649 return [0.0; 3];
650 }
651 let mut s = [0.0_f64; 3];
652 for n in &self.nodes {
653 s[0] += n[0];
654 s[1] += n[1];
655 s[2] += n[2];
656 }
657 let inv = 1.0 / self.nodes.len() as f64;
658 [s[0] * inv, s[1] * inv, s[2] * inv]
659 }
660}
661#[derive(Debug, Clone)]
663#[allow(dead_code)]
664pub struct XdmfAttribute {
665 pub name: String,
667 pub center: String,
669 pub n_components: usize,
671 pub hdf5_path: String,
673}
674#[allow(dead_code)]
676#[derive(Debug, Clone)]
677pub struct XdmfMeshStep {
678 pub time: f64,
680 pub nodes: Vec<[f64; 3]>,
682 pub connectivity: Vec<usize>,
684 pub topology: XdmfTopologyType,
686 pub node_scalars: Vec<(String, Vec<f64>)>,
688 pub node_vectors: Vec<(String, Vec<[f64; 3]>)>,
690}
691impl XdmfMeshStep {
692 pub fn n_elements(&self) -> usize {
694 let npe = self.topology.nodes_per_element();
695 self.connectivity.len().checked_div(npe).unwrap_or(0)
696 }
697}
698#[derive(Debug, Clone, Default)]
700#[allow(dead_code)]
701pub struct XdmfDomainCollection {
702 pub grids: Vec<XdmfGridEntry>,
704}
705impl XdmfDomainCollection {
706 #[allow(dead_code)]
708 pub fn new() -> Self {
709 XdmfDomainCollection { grids: Vec::new() }
710 }
711 #[allow(dead_code)]
713 pub fn add_grid(&mut self, grid: XdmfGridEntry) {
714 self.grids.push(grid);
715 }
716 #[allow(dead_code)]
718 pub fn total_node_count(&self) -> usize {
719 self.grids.iter().map(|g| g.node_count()).sum()
720 }
721 #[allow(dead_code)]
723 pub fn find_grid(&self, name: &str) -> Option<&XdmfGridEntry> {
724 self.grids.iter().find(|g| g.name == name)
725 }
726 #[allow(dead_code)]
728 pub fn write_xml<W: std::io::Write>(&self, writer: &mut W) -> std::io::Result<()> {
729 writeln!(writer, "<?xml version=\"1.0\"?>")?;
730 writeln!(writer, "<Xdmf Version=\"3.0\">")?;
731 writeln!(writer, " <Domain>")?;
732 writeln!(
733 writer,
734 " <Grid Name=\"collection\" GridType=\"Collection\" CollectionType=\"Temporal\">"
735 )?;
736 for grid in &self.grids {
737 let nn = grid.node_count();
738 writeln!(
739 writer,
740 " <Grid Name=\"{}\" GridType=\"Uniform\">",
741 grid.name
742 )?;
743 writeln!(writer, " <Time Value=\"{}\"/>", grid.time)?;
744 writeln!(
745 writer,
746 " <Topology TopologyType=\"{}\" NumberOfElements=\"{}\"/>",
747 grid.topology_type,
748 grid.connectivity.len()
749 )?;
750 writeln!(writer, " <Geometry GeometryType=\"XYZ\">")?;
751 writeln!(
752 writer,
753 " <DataItem Format=\"XML\" Dimensions=\"{nn} 3\">"
754 )?;
755 for n in &grid.nodes {
756 writeln!(writer, " {} {} {}", n[0], n[1], n[2])?;
757 }
758 writeln!(writer, " </DataItem>")?;
759 writeln!(writer, " </Geometry>")?;
760 writeln!(writer, " </Grid>")?;
761 }
762 writeln!(writer, " </Grid>")?;
763 writeln!(writer, " </Domain>")?;
764 writeln!(writer, "</Xdmf>")?;
765 Ok(())
766 }
767 #[allow(dead_code)]
769 pub fn time_range(&self) -> (f64, f64) {
770 if self.grids.is_empty() {
771 return (0.0, 0.0);
772 }
773 let t_min = self
774 .grids
775 .iter()
776 .map(|g| g.time)
777 .fold(f64::INFINITY, f64::min);
778 let t_max = self
779 .grids
780 .iter()
781 .map(|g| g.time)
782 .fold(f64::NEG_INFINITY, f64::max);
783 (t_min, t_max)
784 }
785}
786#[derive(Debug, Clone)]
788#[allow(dead_code)]
789pub struct XdmfFieldDescriptor {
790 pub name: String,
792 pub attribute_type: String,
794 pub center: String,
796 pub data: Vec<f64>,
798 pub n_components: usize,
800}
801impl XdmfFieldDescriptor {
802 #[allow(dead_code)]
804 pub fn scalar(name: &str, data: Vec<f64>) -> Self {
805 XdmfFieldDescriptor {
806 name: name.to_string(),
807 attribute_type: "Scalar".to_string(),
808 center: "Node".to_string(),
809 data,
810 n_components: 1,
811 }
812 }
813 #[allow(dead_code)]
815 pub fn vector(name: &str, data: Vec<f64>) -> Self {
816 XdmfFieldDescriptor {
817 name: name.to_string(),
818 attribute_type: "Vector".to_string(),
819 center: "Node".to_string(),
820 data,
821 n_components: 3,
822 }
823 }
824 #[allow(dead_code)]
826 pub fn entry_count(&self) -> usize {
827 self.data.len().checked_div(self.n_components).unwrap_or(0)
828 }
829 #[allow(dead_code)]
831 pub fn data_lp_norm(&self, p: f64) -> f64 {
832 if p <= 0.0 || self.data.is_empty() {
833 return 0.0;
834 }
835 let sum: f64 = self.data.iter().map(|v| v.abs().powf(p)).sum();
836 sum.powf(1.0 / p)
837 }
838 #[allow(dead_code)]
840 pub fn max_abs(&self) -> f64 {
841 self.data
842 .iter()
843 .cloned()
844 .fold(0.0_f64, |acc, v| acc.max(v.abs()))
845 }
846 #[allow(dead_code)]
848 pub fn min_value(&self) -> f64 {
849 self.data.iter().cloned().fold(f64::INFINITY, f64::min)
850 }
851 #[allow(dead_code)]
853 pub fn max_value(&self) -> f64 {
854 self.data.iter().cloned().fold(f64::NEG_INFINITY, f64::max)
855 }
856}
857#[allow(dead_code)]
860#[derive(Debug, Clone, Default)]
861pub struct XdmfMeshTimeSeries {
862 pub steps: Vec<XdmfMeshStep>,
864}
865#[allow(dead_code)]
866impl XdmfMeshTimeSeries {
867 pub fn new() -> Self {
869 Self::default()
870 }
871 pub fn add_step(&mut self, step: XdmfMeshStep) {
873 self.steps.push(step);
874 }
875 pub fn len(&self) -> usize {
877 self.steps.len()
878 }
879 pub fn is_empty(&self) -> bool {
881 self.steps.is_empty()
882 }
883 pub fn times(&self) -> Vec<f64> {
885 self.steps.iter().map(|s| s.time).collect()
886 }
887 pub fn to_xml(&self) -> String {
889 let mut s = String::new();
890 s.push_str("<?xml version=\"1.0\"?>\n");
891 s.push_str("<Xdmf Version=\"3.0\">\n");
892 s.push_str(" <Domain>\n");
893 s.push_str(
894 " <Grid Name=\"MeshTimeSeries\" GridType=\"Collection\" CollectionType=\"Temporal\">\n",
895 );
896 for step in &self.steps {
897 let n_nodes = step.nodes.len();
898 let n_elements = step.n_elements();
899 let npe = step.topology.nodes_per_element();
900 s.push_str(" <Grid Name=\"mesh\" GridType=\"Uniform\">\n");
901 s.push_str(&format!(" <Time Value=\"{}\"/>\n", step.time));
902 s.push_str(&format!(
903 " <Topology TopologyType=\"{}\" NumberOfElements=\"{}\">\n",
904 step.topology.xdmf_name(),
905 n_elements
906 ));
907 s.push_str(&format!(
908 " <DataItem Format=\"XML\" Dimensions=\"{} {}\">\n",
909 n_elements, npe
910 ));
911 for chunk in step.connectivity.chunks(npe) {
912 let row: Vec<String> = chunk.iter().map(|&i| i.to_string()).collect();
913 s.push_str(&format!(" {}\n", row.join(" ")));
914 }
915 s.push_str(" </DataItem>\n");
916 s.push_str(" </Topology>\n");
917 s.push_str(" <Geometry GeometryType=\"XYZ\">\n");
918 s.push_str(&format!(
919 " <DataItem Format=\"XML\" Dimensions=\"{} 3\">\n",
920 n_nodes
921 ));
922 for p in &step.nodes {
923 s.push_str(&format!(" {} {} {}\n", p[0], p[1], p[2]));
924 }
925 s.push_str(" </DataItem>\n");
926 s.push_str(" </Geometry>\n");
927 for (name, values) in &step.node_scalars {
928 s.push_str(&format!(
929 " <Attribute Name=\"{}\" AttributeType=\"Scalar\" Center=\"Node\">\n",
930 name
931 ));
932 s.push_str(&format!(
933 " <DataItem Format=\"XML\" Dimensions=\"{}\">\n",
934 values.len()
935 ));
936 s.push_str(" ");
937 for (i, v) in values.iter().enumerate() {
938 if i > 0 {
939 s.push(' ');
940 }
941 s.push_str(&format!("{}", v));
942 }
943 s.push('\n');
944 s.push_str(" </DataItem>\n");
945 s.push_str(" </Attribute>\n");
946 }
947 for (name, vdata) in &step.node_vectors {
948 s.push_str(&format!(
949 " <Attribute Name=\"{}\" AttributeType=\"Vector\" Center=\"Node\">\n",
950 name
951 ));
952 s.push_str(&format!(
953 " <DataItem Format=\"XML\" Dimensions=\"{} 3\">\n",
954 vdata.len()
955 ));
956 for v in vdata {
957 s.push_str(&format!(" {} {} {}\n", v[0], v[1], v[2]));
958 }
959 s.push_str(" </DataItem>\n");
960 s.push_str(" </Attribute>\n");
961 }
962 s.push_str(" </Grid>\n");
963 }
964 s.push_str(" </Grid>\n");
965 s.push_str(" </Domain>\n");
966 s.push_str("</Xdmf>\n");
967 s
968 }
969}
970#[derive(Debug, Clone)]
972#[allow(dead_code)]
973pub struct XdmfTimeSeries {
974 pub steps: Vec<XdmfStep>,
976}
977#[allow(dead_code)]
978impl XdmfTimeSeries {
979 pub fn new() -> Self {
981 Self { steps: Vec::new() }
982 }
983 pub fn add_step(
985 &mut self,
986 time: f64,
987 positions: Vec<[f64; 3]>,
988 scalars: Vec<(String, Vec<f64>)>,
989 ) {
990 let n_points = positions.len();
991 self.steps.push(XdmfStep {
992 time,
993 n_points,
994 position_data: positions,
995 scalar_fields: scalars,
996 });
997 }
998 pub fn to_xml(&self) -> String {
1000 let mut s = String::new();
1001 s.push_str("<?xml version=\"1.0\"?>\n");
1002 s.push_str("<Xdmf Version=\"3.0\">\n");
1003 s.push_str(" <Domain>\n");
1004 s.push_str(
1005 " <Grid Name=\"TimeSeries\" GridType=\"Collection\" CollectionType=\"Temporal\">\n",
1006 );
1007 for step in &self.steps {
1008 s.push_str(" <Grid Name=\"particles\" GridType=\"Uniform\">\n");
1009 s.push_str(&format!(" <Time Value=\"{}\"/>\n", step.time));
1010 s.push_str(&format!(
1011 " <Topology TopologyType=\"Polyvertex\" NumberOfElements=\"{}\"/>\n",
1012 step.n_points
1013 ));
1014 s.push_str(" <Geometry GeometryType=\"XYZ\">\n");
1015 s.push_str(&format!(
1016 " <DataItem Format=\"XML\" Dimensions=\"{} 3\">\n",
1017 step.n_points
1018 ));
1019 for p in &step.position_data {
1020 s.push_str(&format!(" {} {} {}\n", p[0], p[1], p[2]));
1021 }
1022 s.push_str(" </DataItem>\n");
1023 s.push_str(" </Geometry>\n");
1024 for (name, values) in &step.scalar_fields {
1025 s.push_str(&format!(
1026 " <Attribute Name=\"{}\" AttributeType=\"Scalar\" Center=\"Node\">\n",
1027 name
1028 ));
1029 s.push_str(&format!(
1030 " <DataItem Format=\"XML\" Dimensions=\"{}\">\n",
1031 values.len()
1032 ));
1033 s.push_str(" ");
1034 for (i, v) in values.iter().enumerate() {
1035 if i > 0 {
1036 s.push(' ');
1037 }
1038 s.push_str(&format!("{}", v));
1039 }
1040 s.push('\n');
1041 s.push_str(" </DataItem>\n");
1042 s.push_str(" </Attribute>\n");
1043 }
1044 s.push_str(" </Grid>\n");
1045 }
1046 s.push_str(" </Grid>\n");
1047 s.push_str(" </Domain>\n");
1048 s.push_str("</Xdmf>\n");
1049 s
1050 }
1051}
1052impl XdmfTimeSeries {
1053 #[allow(dead_code)]
1058 pub fn add_step_with_vectors(
1059 &mut self,
1060 time: f64,
1061 positions: Vec<[f64; 3]>,
1062 scalars: Vec<(String, Vec<f64>)>,
1063 vectors: Vec<(String, Vec<[f64; 3]>)>,
1064 ) {
1065 let n_points = positions.len();
1066 let mut all_scalars = scalars;
1067 for (vname, vdata) in vectors {
1068 let mut xs = Vec::with_capacity(n_points);
1069 let mut ys = Vec::with_capacity(n_points);
1070 let mut zs = Vec::with_capacity(n_points);
1071 for v in &vdata {
1072 xs.push(v[0]);
1073 ys.push(v[1]);
1074 zs.push(v[2]);
1075 }
1076 all_scalars.push((format!("{}_x", vname), xs));
1077 all_scalars.push((format!("{}_y", vname), ys));
1078 all_scalars.push((format!("{}_z", vname), zs));
1079 }
1080 self.steps.push(XdmfStep {
1081 time,
1082 n_points,
1083 position_data: positions,
1084 scalar_fields: all_scalars,
1085 });
1086 }
1087 #[allow(dead_code)]
1089 pub fn times(&self) -> Vec<f64> {
1090 self.steps.iter().map(|s| s.time).collect()
1091 }
1092 #[allow(dead_code)]
1094 pub fn total_particle_count(&self) -> usize {
1095 self.steps.iter().map(|s| s.n_points).sum()
1096 }
1097}
1098impl XdmfTimeSeries {
1099 #[allow(dead_code)]
1104 pub fn add_frame(&mut self, time: f64, positions: Vec<[f64; 3]>) {
1105 self.add_step(time, positions, vec![]);
1106 }
1107 #[allow(dead_code)]
1109 pub fn add_frame_with_scalars(
1110 &mut self,
1111 time: f64,
1112 positions: Vec<[f64; 3]>,
1113 scalars: Vec<(String, Vec<f64>)>,
1114 ) {
1115 self.add_step(time, positions, scalars);
1116 }
1117 #[allow(dead_code)]
1122 pub fn write_xml<W: Write>(&self, writer: &mut W) -> std::io::Result<()> {
1123 let xml = self.to_xml();
1124 writer.write_all(xml.as_bytes())
1125 }
1126 #[allow(dead_code)]
1128 pub fn write_xml_to_file(&self, path: &str) -> std::io::Result<()> {
1129 let mut f = std::fs::File::create(path)?;
1130 self.write_xml(&mut f)
1131 }
1132}
1133#[derive(Debug, Clone)]
1135#[allow(dead_code)]
1136pub struct XdmfMeshPatch {
1137 pub name: String,
1139 pub element_ids: Vec<usize>,
1141 pub tag: Option<i32>,
1143}
1144impl XdmfMeshPatch {
1145 #[allow(dead_code)]
1147 pub fn new(name: &str, element_ids: Vec<usize>) -> Self {
1148 XdmfMeshPatch {
1149 name: name.to_string(),
1150 element_ids,
1151 tag: None,
1152 }
1153 }
1154 #[allow(dead_code)]
1156 pub fn element_count(&self) -> usize {
1157 self.element_ids.len()
1158 }
1159 #[allow(dead_code)]
1161 pub fn contains_element(&self, idx: usize) -> bool {
1162 self.element_ids.contains(&idx)
1163 }
1164 #[allow(dead_code)]
1166 pub fn merge(&mut self, other: &XdmfMeshPatch) {
1167 for &id in &other.element_ids {
1168 if !self.element_ids.contains(&id) {
1169 self.element_ids.push(id);
1170 }
1171 }
1172 }
1173 #[allow(dead_code)]
1175 pub fn to_debug_string(&self) -> String {
1176 format!(
1177 "patch \"{}\" [{} elements] tag={:?}",
1178 self.name,
1179 self.element_ids.len(),
1180 self.tag
1181 )
1182 }
1183}
1184#[allow(dead_code)]
1186pub struct XdmfReader;
1187#[allow(dead_code)]
1188impl XdmfReader {
1189 pub fn from_xml(data: &str) -> Result<Vec<XdmfStep>, std::io::Error> {
1194 let mut steps = Vec::new();
1195 let mut current_time: Option<f64> = None;
1196 let mut current_n: Option<usize> = None;
1197 for line in data.lines() {
1198 let trimmed = line.trim();
1199 if trimmed.starts_with("<Time")
1200 && let Some(start) = trimmed.find("Value=\"")
1201 {
1202 let rest = &trimmed[start + 7..];
1203 if let Some(end) = rest.find('"')
1204 && let Ok(t) = rest[..end].parse::<f64>()
1205 {
1206 current_time = Some(t);
1207 }
1208 }
1209 if trimmed.contains("NumberOfElements=\"")
1210 && let Some(start) = trimmed.find("NumberOfElements=\"")
1211 {
1212 let rest = &trimmed[start + 18..];
1213 if let Some(end) = rest.find('"')
1214 && let Ok(n) = rest[..end].parse::<usize>()
1215 {
1216 current_n = Some(n);
1217 }
1218 }
1219 if trimmed == "</Grid>"
1220 && let (Some(t), Some(n)) = (current_time.take(), current_n.take())
1221 {
1222 steps.push(XdmfStep {
1223 time: t,
1224 n_points: n,
1225 position_data: Vec::new(),
1226 scalar_fields: Vec::new(),
1227 });
1228 }
1229 }
1230 Ok(steps)
1231 }
1232}
1233#[allow(dead_code)]
1235#[derive(Debug, Clone)]
1236pub struct XdmfSchema {
1237 pub topology_type: String,
1239 pub required_attributes: Vec<String>,
1241}
1242impl XdmfSchema {
1243 pub fn new(topology_type: &str, required_attributes: Vec<String>) -> Self {
1245 Self {
1246 topology_type: topology_type.to_owned(),
1247 required_attributes,
1248 }
1249 }
1250 pub fn validate(&self, xml: &str) -> Vec<String> {
1254 let mut errors = Vec::new();
1255 if !xml.contains(&format!("TopologyType=\"{}\"", self.topology_type)) {
1256 errors.push(format!("missing TopologyType=\"{}\"", self.topology_type));
1257 }
1258 for attr in &self.required_attributes {
1259 if !xml.contains(&format!("Name=\"{}\"", attr)) {
1260 errors.push(format!("missing Attribute Name=\"{}\"", attr));
1261 }
1262 }
1263 errors
1264 }
1265}