oxiphysics_io/
mesh_io.rs

1// Copyright 2026 COOLJAPAN OU (Team KitaSan)
2// SPDX-License-Identifier: Apache-2.0
3
4//! Extended mesh I/O formats.
5//!
6//! Provides writers and readers for glTF 2.0, USD ASCII, Gmsh .msh v2/v4,
7//! Abaqus .inp, NASTRAN BDF, mesh conversion, and mesh quality validation.
8
9#![allow(dead_code)]
10
11use std::collections::{HashMap, HashSet};
12use std::io::Write;
13
14// ---------------------------------------------------------------------------
15// Helper free functions
16// ---------------------------------------------------------------------------
17
18/// Encode a byte slice to a base-64 ASCII string.
19pub fn base64_encode_buffer(data: &[u8]) -> String {
20    const TABLE: &[u8] = b"ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789+/";
21    let mut out = String::with_capacity(data.len().div_ceil(3) * 4);
22    for chunk in data.chunks(3) {
23        let b0 = chunk[0] as usize;
24        let b1 = if chunk.len() > 1 {
25            chunk[1] as usize
26        } else {
27            0
28        };
29        let b2 = if chunk.len() > 2 {
30            chunk[2] as usize
31        } else {
32            0
33        };
34        out.push(TABLE[b0 >> 2] as char);
35        out.push(TABLE[((b0 & 3) << 4) | (b1 >> 4)] as char);
36        if chunk.len() > 1 {
37            out.push(TABLE[((b1 & 0xf) << 2) | (b2 >> 6)] as char);
38        } else {
39            out.push('=');
40        }
41        if chunk.len() > 2 {
42            out.push(TABLE[b2 & 0x3f] as char);
43        } else {
44            out.push('=');
45        }
46    }
47    out
48}
49
50/// Write a glTF accessor JSON fragment.
51pub fn write_gltf_accessor<W: Write>(
52    w: &mut W,
53    buffer_view: usize,
54    component_type: u32,
55    count: usize,
56    type_str: &str,
57    trailing_comma: bool,
58) -> std::io::Result<()> {
59    let comma = if trailing_comma { "," } else { "" };
60    writeln!(
61        w,
62        r#"    {{"bufferView":{buffer_view},"componentType":{component_type},"count":{count},"type":"{type_str}"}}{comma}"#
63    )
64}
65
66/// Format a floating-point value in NASTRAN field-8 format (8 chars wide).
67pub fn nastran_field8(v: f64) -> String {
68    // Use E notation, fitting into 8 characters
69    let s = format!("{:>8.4e}", v);
70    if s.len() > 8 {
71        format!("{:>8.3e}", v)
72    } else {
73        s
74    }
75}
76
77/// Return the Gmsh element type string for a given number of nodes.
78/// 3 → "Triangle", 4 → "Quad", 4 → "Tet4", 8 → "Hex8".
79pub fn gmsh_element_type_str(n_nodes: usize) -> &'static str {
80    match n_nodes {
81        2 => "Line",
82        3 => "Triangle",
83        4 => "Quad",
84        10 => "Tet10",
85        _ => "Unknown",
86    }
87}
88
89// ---------------------------------------------------------------------------
90// GltfMeshWriter
91// ---------------------------------------------------------------------------
92
93/// Write a mesh to glTF 2.0 format (JSON + binary buffer).
94pub struct GltfMeshWriter {
95    /// Vertex positions.
96    pub vertices: Vec<[f64; 3]>,
97    /// Vertex normals.
98    pub normals: Vec<[f64; 3]>,
99    /// UV texture coordinates.
100    pub uvs: Vec<[f64; 2]>,
101    /// Triangle indices.
102    pub indices: Vec<[u32; 3]>,
103    /// Material name.
104    pub material_name: String,
105}
106
107impl GltfMeshWriter {
108    /// Construct with mesh data.
109    pub fn new(
110        vertices: Vec<[f64; 3]>,
111        normals: Vec<[f64; 3]>,
112        uvs: Vec<[f64; 2]>,
113        indices: Vec<[u32; 3]>,
114    ) -> Self {
115        GltfMeshWriter {
116            vertices,
117            normals,
118            uvs,
119            indices,
120            material_name: "default".to_string(),
121        }
122    }
123
124    /// Pack positions as f32 little-endian bytes.
125    pub fn positions_f32_bytes(&self) -> Vec<u8> {
126        let mut buf = Vec::with_capacity(self.vertices.len() * 12);
127        for &[x, y, z] in &self.vertices {
128            buf.extend_from_slice(&(x as f32).to_le_bytes());
129            buf.extend_from_slice(&(y as f32).to_le_bytes());
130            buf.extend_from_slice(&(z as f32).to_le_bytes());
131        }
132        buf
133    }
134
135    /// Pack indices as u32 little-endian bytes.
136    pub fn indices_u32_bytes(&self) -> Vec<u8> {
137        let mut buf = Vec::with_capacity(self.indices.len() * 12);
138        for &[a, b, c] in &self.indices {
139            buf.extend_from_slice(&a.to_le_bytes());
140            buf.extend_from_slice(&b.to_le_bytes());
141            buf.extend_from_slice(&c.to_le_bytes());
142        }
143        buf
144    }
145
146    /// Write the glTF JSON to `writer` with embedded base64 binary data.
147    pub fn write_gltf_json<W: Write>(&self, w: &mut W) -> std::io::Result<()> {
148        let pos_bytes = self.positions_f32_bytes();
149        let idx_bytes = self.indices_u32_bytes();
150        let mut binary = pos_bytes.clone();
151        binary.extend_from_slice(&idx_bytes);
152        let b64 = base64_encode_buffer(&binary);
153        let n_verts = self.vertices.len();
154        let n_tris = self.indices.len();
155        let pos_len = pos_bytes.len();
156        let idx_len = idx_bytes.len();
157        writeln!(w, "{{")?;
158        writeln!(w, r#"  "asset": {{"version": "2.0"}},"#)?;
159        writeln!(w, r#"  "scene": 0,"#)?;
160        writeln!(w, r#"  "scenes": [{{"nodes": [0]}}],"#)?;
161        writeln!(w, r#"  "nodes": [{{"mesh": 0}}],"#)?;
162        writeln!(
163            w,
164            r#"  "meshes": [{{"name": "mesh","primitives": [{{"attributes": {{"POSITION": 0}},"indices": 1}}]}}],"#
165        )?;
166        writeln!(w, r#"  "accessors": ["#)?;
167        write_gltf_accessor(w, 0, 5126, n_verts, "VEC3", true)?;
168        write_gltf_accessor(w, 1, 5125, n_tris * 3, "SCALAR", false)?;
169        writeln!(w, r#"  ],"#)?;
170        writeln!(w, r#"  "bufferViews": ["#)?;
171        writeln!(
172            w,
173            r#"    {{"buffer": 0, "byteOffset": 0, "byteLength": {pos_len}}},"#
174        )?;
175        writeln!(
176            w,
177            r#"    {{"buffer": 0, "byteOffset": {pos_len}, "byteLength": {idx_len}}}"#
178        )?;
179        writeln!(w, r#"  ],"#)?;
180        writeln!(
181            w,
182            r#"  "buffers": [{{"byteLength": {total}, "uri": "data:application/octet-stream;base64,{b64}"}}]"#,
183            total = binary.len()
184        )?;
185        writeln!(w, "}}")?;
186        Ok(())
187    }
188}
189
190// ---------------------------------------------------------------------------
191// GltfMeshReader
192// ---------------------------------------------------------------------------
193
194/// Read a glTF mesh and extract vertices, normals, and indices.
195///
196/// This is a simplified reader that parses a subset of the glTF JSON format.
197pub struct GltfMeshReader {
198    /// Parsed vertex positions.
199    pub vertices: Vec<[f64; 3]>,
200    /// Parsed triangle indices.
201    pub indices: Vec<[u32; 3]>,
202    /// Node names encountered.
203    pub node_names: Vec<String>,
204}
205
206impl GltfMeshReader {
207    /// Parse a glTF JSON string (very simplified: extracts node names only).
208    pub fn parse_json(json: &str) -> Self {
209        let mut node_names = Vec::new();
210        // Extract mesh names
211        for line in json.lines() {
212            let line = line.trim();
213            if line.contains("\"name\"")
214                && let Some(start) = line.find("\"name\":")
215            {
216                let rest = &line[start + 7..].trim_start_matches([' ', '"'].as_ref());
217                let end = rest.find('"').unwrap_or(rest.len());
218                node_names.push(rest[..end].to_string());
219            }
220        }
221        GltfMeshReader {
222            vertices: Vec::new(),
223            indices: Vec::new(),
224            node_names,
225        }
226    }
227}
228
229// ---------------------------------------------------------------------------
230// UsdMeshWriter
231// ---------------------------------------------------------------------------
232
233/// Write a mesh to USD (Universal Scene Description) ASCII format.
234pub struct UsdMeshWriter {
235    /// Mesh vertices.
236    pub vertices: Vec<[f64; 3]>,
237    /// Face vertex counts (e.g., 3 per triangle).
238    pub face_vertex_counts: Vec<u32>,
239    /// Flat face vertex indices.
240    pub face_vertex_indices: Vec<u32>,
241    /// USD prim path.
242    pub prim_path: String,
243}
244
245impl UsdMeshWriter {
246    /// Construct from triangle mesh data.
247    pub fn from_triangles(
248        vertices: Vec<[f64; 3]>,
249        triangles: &[[u32; 3]],
250        prim_path: &str,
251    ) -> Self {
252        let face_vertex_counts = vec![3; triangles.len()];
253        let face_vertex_indices: Vec<u32> =
254            triangles.iter().flat_map(|t| t.iter().copied()).collect();
255        UsdMeshWriter {
256            vertices,
257            face_vertex_counts,
258            face_vertex_indices,
259            prim_path: prim_path.to_string(),
260        }
261    }
262
263    /// Write USD ASCII to `writer`.
264    pub fn write<W: Write>(&self, w: &mut W) -> std::io::Result<()> {
265        writeln!(w, "#usda 1.0")?;
266        writeln!(
267            w,
268            r#"def Mesh "{name}" {{"#,
269            name = self.prim_path.trim_start_matches('/')
270        )?;
271        // Points
272        write!(w, "    point3f[] points = [")?;
273        for (i, &[x, y, z]) in self.vertices.iter().enumerate() {
274            if i > 0 {
275                write!(w, ", ")?;
276            }
277            write!(w, "({x}, {y}, {z})")?;
278        }
279        writeln!(w, "]")?;
280        // Face vertex counts
281        write!(w, "    int[] faceVertexCounts = [")?;
282        for (i, &c) in self.face_vertex_counts.iter().enumerate() {
283            if i > 0 {
284                write!(w, ", ")?;
285            }
286            write!(w, "{c}")?;
287        }
288        writeln!(w, "]")?;
289        // Face vertex indices
290        write!(w, "    int[] faceVertexIndices = [")?;
291        for (i, &idx) in self.face_vertex_indices.iter().enumerate() {
292            if i > 0 {
293                write!(w, ", ")?;
294            }
295            write!(w, "{idx}")?;
296        }
297        writeln!(w, "]")?;
298        writeln!(w, "}}")?;
299        Ok(())
300    }
301}
302
303// ---------------------------------------------------------------------------
304// MshWriter
305// ---------------------------------------------------------------------------
306
307/// Write a Gmsh .msh format v2 mesh.
308pub struct MshWriter {
309    /// Node coordinates.
310    pub nodes: Vec<[f64; 3]>,
311    /// Elements: each is a list of node indices (1-based in Gmsh).
312    pub elements: Vec<Vec<usize>>,
313    /// Physical group tags per element.
314    pub element_tags: Vec<i64>,
315    /// Mesh format version.
316    pub version: u32,
317}
318
319impl MshWriter {
320    /// Construct an empty MshWriter.
321    pub fn new(version: u32) -> Self {
322        MshWriter {
323            nodes: Vec::new(),
324            elements: Vec::new(),
325            element_tags: Vec::new(),
326            version,
327        }
328    }
329
330    /// Add a node.
331    pub fn add_node(&mut self, x: f64, y: f64, z: f64) {
332        self.nodes.push([x, y, z]);
333    }
334
335    /// Add a triangular element with a physical group tag.
336    pub fn add_triangle(&mut self, a: usize, b: usize, c: usize, tag: i64) {
337        self.elements.push(vec![a, b, c]);
338        self.element_tags.push(tag);
339    }
340
341    /// Write .msh v2 format to `writer`.
342    pub fn write_v2<W: Write>(&self, w: &mut W) -> std::io::Result<()> {
343        writeln!(w, "$MeshFormat")?;
344        writeln!(w, "2.2 0 8")?;
345        writeln!(w, "$EndMeshFormat")?;
346        writeln!(w, "$Nodes")?;
347        writeln!(w, "{}", self.nodes.len())?;
348        for (i, &[x, y, z]) in self.nodes.iter().enumerate() {
349            writeln!(w, "{} {} {} {}", i + 1, x, y, z)?;
350        }
351        writeln!(w, "$EndNodes")?;
352        writeln!(w, "$Elements")?;
353        writeln!(w, "{}", self.elements.len())?;
354        for (i, elem) in self.elements.iter().enumerate() {
355            let tag = self.element_tags.get(i).copied().unwrap_or(0);
356            let elem_type = match elem.len() {
357                2 => 1, // line
358                3 => 2, // triangle
359                4 => 3, // quad
360                _ => 2,
361            };
362            write!(w, "{} {} 2 {} {}", i + 1, elem_type, tag, tag)?;
363            for &n in elem {
364                write!(w, " {n}")?;
365            }
366            writeln!(w)?;
367        }
368        writeln!(w, "$EndElements")?;
369        Ok(())
370    }
371}
372
373// ---------------------------------------------------------------------------
374// MshReader
375// ---------------------------------------------------------------------------
376
377/// Read a Gmsh .msh v2 mesh.
378pub struct MshReader {
379    /// Parsed node coordinates.
380    pub nodes: Vec<[f64; 3]>,
381    /// Parsed elements (node indices, 0-based).
382    pub elements: Vec<Vec<usize>>,
383    /// Element tags.
384    pub element_tags: Vec<i64>,
385}
386
387impl MshReader {
388    /// Parse a .msh v2 string.
389    pub fn parse(s: &str) -> Self {
390        let mut nodes = Vec::new();
391        let mut elements = Vec::new();
392        let mut element_tags = Vec::new();
393        let mut in_nodes = false;
394        let mut in_elements = false;
395        let mut node_count = 0usize;
396        let mut elem_count = 0usize;
397        let mut node_idx = 0usize;
398        let mut elem_idx = 0usize;
399
400        for line in s.lines() {
401            let line = line.trim();
402            if line == "$Nodes" {
403                in_nodes = true;
404                node_idx = 0;
405                continue;
406            }
407            if line == "$EndNodes" {
408                in_nodes = false;
409                continue;
410            }
411            if line == "$Elements" {
412                in_elements = true;
413                elem_idx = 0;
414                continue;
415            }
416            if line == "$EndElements" {
417                in_elements = false;
418                continue;
419            }
420            if in_nodes {
421                if node_idx == 0 {
422                    node_count = line.parse().unwrap_or(0);
423                    nodes.reserve(node_count);
424                    node_idx += 1;
425                    continue;
426                }
427                let vals: Vec<f64> = line
428                    .split_whitespace()
429                    .skip(1)
430                    .filter_map(|v| v.parse().ok())
431                    .collect();
432                if vals.len() >= 3 {
433                    nodes.push([vals[0], vals[1], vals[2]]);
434                }
435            }
436            if in_elements {
437                if elem_idx == 0 {
438                    elem_count = line.parse().unwrap_or(0);
439                    elements.reserve(elem_count);
440                    element_tags.reserve(elem_count);
441                    elem_idx += 1;
442                    continue;
443                }
444                let vals: Vec<i64> = line
445                    .split_whitespace()
446                    .filter_map(|v| v.parse().ok())
447                    .collect();
448                if vals.len() >= 5 {
449                    let n_tags = vals[2] as usize;
450                    let tag = vals[3];
451                    let node_start = 3 + n_tags;
452                    let node_indices: Vec<usize> = vals[node_start..]
453                        .iter()
454                        .map(|&v| (v - 1) as usize)
455                        .collect();
456                    elements.push(node_indices);
457                    element_tags.push(tag);
458                }
459            }
460        }
461        let _ = node_count;
462        let _ = elem_count;
463        MshReader {
464            nodes,
465            elements,
466            element_tags,
467        }
468    }
469}
470
471// ---------------------------------------------------------------------------
472// AbaqusInpWriter
473// ---------------------------------------------------------------------------
474
475/// Write an Abaqus .inp mesh file.
476pub struct AbaqusInpWriter {
477    /// Part name.
478    pub part_name: String,
479    /// Node coordinates.
480    pub nodes: Vec<[f64; 3]>,
481    /// Elements (node indices, 1-based in output).
482    pub elements: Vec<Vec<usize>>,
483    /// Element type string (e.g., "C3D4" for tet).
484    pub element_type: String,
485    /// Node sets: name → list of node indices (1-based).
486    pub node_sets: HashMap<String, Vec<usize>>,
487    /// Boundary conditions: node set name → (dof_min, dof_max, value).
488    pub boundary_conditions: Vec<(String, usize, usize, f64)>,
489}
490
491impl AbaqusInpWriter {
492    /// Construct an empty writer.
493    pub fn new(part_name: &str, element_type: &str) -> Self {
494        AbaqusInpWriter {
495            part_name: part_name.to_string(),
496            nodes: Vec::new(),
497            elements: Vec::new(),
498            element_type: element_type.to_string(),
499            node_sets: HashMap::new(),
500            boundary_conditions: Vec::new(),
501        }
502    }
503
504    /// Add a node.
505    pub fn add_node(&mut self, x: f64, y: f64, z: f64) {
506        self.nodes.push([x, y, z]);
507    }
508
509    /// Add an element.
510    pub fn add_element(&mut self, nodes: Vec<usize>) {
511        self.elements.push(nodes);
512    }
513
514    /// Add a node set.
515    pub fn add_node_set(&mut self, name: &str, indices: Vec<usize>) {
516        self.node_sets.insert(name.to_string(), indices);
517    }
518
519    /// Add a boundary condition.
520    pub fn add_boundary(&mut self, nset: &str, dof_min: usize, dof_max: usize, value: f64) {
521        self.boundary_conditions
522            .push((nset.to_string(), dof_min, dof_max, value));
523    }
524
525    /// Write the .inp file to `writer`.
526    pub fn write<W: Write>(&self, w: &mut W) -> std::io::Result<()> {
527        writeln!(w, "*PART, NAME={}", self.part_name)?;
528        writeln!(w, "*NODE")?;
529        for (i, &[x, y, z]) in self.nodes.iter().enumerate() {
530            writeln!(w, "{}, {}, {}, {}", i + 1, x, y, z)?;
531        }
532        writeln!(w, "*ELEMENT, TYPE={}", self.element_type)?;
533        for (i, elem) in self.elements.iter().enumerate() {
534            write!(w, "{}", i + 1)?;
535            for &n in elem {
536                write!(w, ", {n}")?;
537            }
538            writeln!(w)?;
539        }
540        for (name, indices) in &self.node_sets {
541            writeln!(w, "*NSET, NSET={name}")?;
542            for chunk in indices.chunks(8) {
543                let s: Vec<String> = chunk.iter().map(|n| n.to_string()).collect();
544                writeln!(w, "{}", s.join(", "))?;
545            }
546        }
547        if !self.boundary_conditions.is_empty() {
548            writeln!(w, "*BOUNDARY")?;
549            for (nset, d1, d2, val) in &self.boundary_conditions {
550                writeln!(w, "{nset}, {d1}, {d2}, {val}")?;
551            }
552        }
553        writeln!(w, "*END PART")?;
554        Ok(())
555    }
556}
557
558// ---------------------------------------------------------------------------
559// AbaqusInpReader
560// ---------------------------------------------------------------------------
561
562/// Read an Abaqus .inp mesh.
563pub struct AbaqusInpReader {
564    /// Parsed nodes.
565    pub nodes: Vec<[f64; 3]>,
566    /// Parsed elements.
567    pub elements: Vec<Vec<usize>>,
568    /// Element type.
569    pub element_type: String,
570    /// Node sets.
571    pub node_sets: HashMap<String, Vec<usize>>,
572}
573
574impl AbaqusInpReader {
575    /// Parse an .inp string.
576    pub fn parse(s: &str) -> Self {
577        let mut nodes = Vec::new();
578        let mut elements = Vec::new();
579        let mut element_type = String::new();
580        let mut node_sets: HashMap<String, Vec<usize>> = HashMap::new();
581        let mut mode = "";
582        let mut cur_nset = String::new();
583
584        for line in s.lines() {
585            let line = line.trim();
586            if line.starts_with('*') {
587                let kw = line.to_uppercase();
588                if kw.starts_with("*NODE") && !kw.starts_with("*NSET") {
589                    mode = "node";
590                } else if kw.starts_with("*ELEMENT") {
591                    mode = "element";
592                    if let Some(tp) = line.split("TYPE=").nth(1) {
593                        element_type = tp.split(',').next().unwrap_or("").trim().to_string();
594                    }
595                } else if kw.starts_with("*NSET") {
596                    mode = "nset";
597                    if let Some(name) = line.split("NSET=").nth(1) {
598                        cur_nset = name.split(',').next().unwrap_or("").trim().to_string();
599                        node_sets.entry(cur_nset.clone()).or_default();
600                    }
601                } else if kw.starts_with("*BOUNDARY") {
602                    mode = "boundary";
603                } else {
604                    mode = "";
605                }
606                continue;
607            }
608            match mode {
609                "node" => {
610                    let vals: Vec<f64> = line
611                        .split(',')
612                        .filter_map(|v| v.trim().parse().ok())
613                        .collect();
614                    if vals.len() >= 4 {
615                        nodes.push([vals[1], vals[2], vals[3]]);
616                    }
617                }
618                "element" => {
619                    let vals: Vec<usize> = line
620                        .split(',')
621                        .filter_map(|v| v.trim().parse().ok())
622                        .collect();
623                    if vals.len() >= 2 {
624                        elements.push(vals[1..].to_vec());
625                    }
626                }
627                "nset" => {
628                    let indices: Vec<usize> = line
629                        .split(',')
630                        .filter_map(|v| v.trim().parse().ok())
631                        .collect();
632                    if let Some(set) = node_sets.get_mut(&cur_nset) {
633                        set.extend(indices);
634                    }
635                }
636                _ => {}
637            }
638        }
639        AbaqusInpReader {
640            nodes,
641            elements,
642            element_type,
643            node_sets,
644        }
645    }
646}
647
648// ---------------------------------------------------------------------------
649// NastranBdf
650// ---------------------------------------------------------------------------
651
652/// A NASTRAN BDF GRID entry.
653#[derive(Debug, Clone)]
654pub struct NastranGrid {
655    /// Node id.
656    pub id: u64,
657    /// Coordinate system id.
658    pub cp: u64,
659    /// Position.
660    pub xyz: [f64; 3],
661}
662
663/// A NASTRAN solid element (CTETRA or CHEXA).
664#[derive(Debug, Clone)]
665pub struct NastranElement {
666    /// Element id.
667    pub eid: u64,
668    /// Property id.
669    pub pid: u64,
670    /// Node ids.
671    pub nodes: Vec<u64>,
672    /// Element type.
673    pub elem_type: String,
674}
675
676/// Read/write NASTRAN BDF (GRID, CTETRA, CHEXA, PSOLID, MAT1, SPC, FORCE).
677pub struct NastranBdf {
678    /// Grid points.
679    pub grids: Vec<NastranGrid>,
680    /// Elements.
681    pub elements: Vec<NastranElement>,
682}
683
684impl NastranBdf {
685    /// Construct empty.
686    pub fn new() -> Self {
687        NastranBdf {
688            grids: Vec::new(),
689            elements: Vec::new(),
690        }
691    }
692
693    /// Add a grid point.
694    pub fn add_grid(&mut self, id: u64, xyz: [f64; 3]) {
695        self.grids.push(NastranGrid { id, cp: 0, xyz });
696    }
697
698    /// Add an element.
699    pub fn add_element(&mut self, eid: u64, pid: u64, nodes: Vec<u64>, elem_type: &str) {
700        self.elements.push(NastranElement {
701            eid,
702            pid,
703            nodes,
704            elem_type: elem_type.to_string(),
705        });
706    }
707
708    /// Write BDF to `writer`.
709    pub fn write<W: Write>(&self, w: &mut W) -> std::io::Result<()> {
710        writeln!(w, "BEGIN BULK")?;
711        for g in &self.grids {
712            writeln!(
713                w,
714                "GRID    {:8}{:8}{:8}{:8}{:8}",
715                g.id,
716                g.cp,
717                nastran_field8(g.xyz[0]),
718                nastran_field8(g.xyz[1]),
719                nastran_field8(g.xyz[2]),
720            )?;
721        }
722        for e in &self.elements {
723            write!(w, "{:<8}{:8}{:8}", e.elem_type, e.eid, e.pid)?;
724            for &n in &e.nodes {
725                write!(w, "{:8}", n)?;
726            }
727            writeln!(w)?;
728        }
729        writeln!(w, "ENDDATA")?;
730        Ok(())
731    }
732
733    /// Parse a BDF string and return grids and elements.
734    pub fn parse(s: &str) -> Self {
735        let mut bdf = NastranBdf::new();
736        for line in s.lines() {
737            let line = if line.len() > 8 { line } else { continue };
738            let kw = &line[..8].trim().to_uppercase();
739            if kw.starts_with("GRID") {
740                if let Some(g) = Self::parse_grid(line) {
741                    bdf.grids.push(g);
742                }
743            } else if (kw.starts_with("CTETRA") || kw.starts_with("CHEXA"))
744                && let Some(e) = Self::parse_element(line)
745            {
746                bdf.elements.push(e);
747            }
748        }
749        bdf
750    }
751
752    fn parse_grid(line: &str) -> Option<NastranGrid> {
753        let fields: Vec<&str> = (0..9)
754            .map(|i| {
755                let s = i * 8;
756                let e = (s + 8).min(line.len());
757                if s < line.len() {
758                    line[s..e].trim()
759                } else {
760                    ""
761                }
762            })
763            .collect();
764        let id: u64 = fields[1].parse().ok()?;
765        let x: f64 = fields[3].parse().ok()?;
766        let y: f64 = fields[4].parse().ok()?;
767        let z: f64 = fields[5].parse().ok()?;
768        Some(NastranGrid {
769            id,
770            cp: 0,
771            xyz: [x, y, z],
772        })
773    }
774
775    fn parse_element(line: &str) -> Option<NastranElement> {
776        let kw = line[..8].trim().to_uppercase();
777        let fields: Vec<&str> = (0..12)
778            .map(|i| {
779                let s = i * 8;
780                let e = (s + 8).min(line.len());
781                if s < line.len() {
782                    line[s..e].trim()
783                } else {
784                    ""
785                }
786            })
787            .collect();
788        let eid: u64 = fields[1].parse().ok()?;
789        let pid: u64 = fields[2].parse().ok()?;
790        let nodes: Vec<u64> = fields[3..].iter().filter_map(|f| f.parse().ok()).collect();
791        Some(NastranElement {
792            eid,
793            pid,
794            nodes,
795            elem_type: kw,
796        })
797    }
798}
799
800impl Default for NastranBdf {
801    fn default() -> Self {
802        Self::new()
803    }
804}
805
806// ---------------------------------------------------------------------------
807// MeshConverter
808// ---------------------------------------------------------------------------
809
810/// Convert between mesh formats, re-index nodes, and reorder elements.
811pub struct MeshConverter;
812
813impl MeshConverter {
814    /// Re-index nodes: remove unused nodes and remap element indices.
815    /// Returns (new_nodes, new_elements).
816    pub fn compact_nodes(
817        nodes: &[[f64; 3]],
818        elements: &[Vec<usize>],
819    ) -> (Vec<[f64; 3]>, Vec<Vec<usize>>) {
820        let mut used: HashSet<usize> = HashSet::new();
821        for elem in elements {
822            for &n in elem {
823                used.insert(n);
824            }
825        }
826        let mut old_to_new: HashMap<usize, usize> = HashMap::new();
827        let mut new_nodes = Vec::new();
828        let mut sorted: Vec<usize> = used.into_iter().collect();
829        sorted.sort_unstable();
830        for old_idx in sorted {
831            old_to_new.insert(old_idx, new_nodes.len());
832            if old_idx < nodes.len() {
833                new_nodes.push(nodes[old_idx]);
834            }
835        }
836        let new_elements: Vec<Vec<usize>> = elements
837            .iter()
838            .map(|e| {
839                e.iter()
840                    .filter_map(|&n| old_to_new.get(&n).copied())
841                    .collect()
842            })
843            .collect();
844        (new_nodes, new_elements)
845    }
846
847    /// Convert triangle mesh to indexed mesh (deduplicate vertices by position).
848    pub fn deduplicate_vertices(
849        positions: &[[f64; 3]],
850        indices: &[[u32; 3]],
851        tolerance: f64,
852    ) -> (Vec<[f64; 3]>, Vec<[u32; 3]>) {
853        let mut unique: Vec<[f64; 3]> = Vec::new();
854        let mut remap: Vec<u32> = Vec::with_capacity(positions.len());
855        for &p in positions {
856            let mut found = None;
857            for (j, &u) in unique.iter().enumerate() {
858                let dx = p[0] - u[0];
859                let dy = p[1] - u[1];
860                let dz = p[2] - u[2];
861                if (dx * dx + dy * dy + dz * dz).sqrt() < tolerance {
862                    found = Some(j as u32);
863                    break;
864                }
865            }
866            if let Some(j) = found {
867                remap.push(j);
868            } else {
869                remap.push(unique.len() as u32);
870                unique.push(p);
871            }
872        }
873        let new_indices: Vec<[u32; 3]> = indices
874            .iter()
875            .map(|&[a, b, c]| [remap[a as usize], remap[b as usize], remap[c as usize]])
876            .collect();
877        (unique, new_indices)
878    }
879
880    /// Flip winding order of all triangles.
881    pub fn flip_winding(indices: &[[u32; 3]]) -> Vec<[u32; 3]> {
882        indices.iter().map(|&[a, b, c]| [a, c, b]).collect()
883    }
884}
885
886// ---------------------------------------------------------------------------
887// MeshValidator
888// ---------------------------------------------------------------------------
889
890/// Mesh quality metrics.
891#[derive(Debug, Clone)]
892pub struct MeshQuality {
893    /// Minimum aspect ratio (ideal = 1 for equilateral).
894    pub min_aspect_ratio: f64,
895    /// Maximum aspect ratio.
896    pub max_aspect_ratio: f64,
897    /// Minimum triangle area.
898    pub min_area: f64,
899    /// Number of degenerate triangles (area < epsilon).
900    pub degenerate_count: usize,
901    /// Number of inverted triangles (negative Jacobian).
902    pub inverted_count: usize,
903}
904
905/// Check mesh quality: aspect ratio, Jacobian, area, orientation.
906pub struct MeshValidator;
907
908impl MeshValidator {
909    /// Compute quality metrics for a triangle mesh.
910    pub fn validate(
911        vertices: &[[f64; 3]],
912        triangles: &[[u32; 3]],
913        area_epsilon: f64,
914    ) -> MeshQuality {
915        let mut min_ar = f64::MAX;
916        let mut max_ar = f64::MIN;
917        let mut min_area = f64::MAX;
918        let mut deg_count = 0;
919        let mut inv_count = 0;
920
921        for &[a, b, c] in triangles {
922            let va = vertices[a as usize];
923            let vb = vertices[b as usize];
924            let vc = vertices[c as usize];
925            let ab = sub3(vb, va);
926            let ac = sub3(vc, va);
927            let bc = sub3(vc, vb);
928            let cross = cross3(ab, ac);
929            let area = len3(cross) * 0.5;
930            let ar = triangle_aspect_ratio(len3(ab), len3(ac), len3(bc));
931            min_ar = min_ar.min(ar);
932            max_ar = max_ar.max(ar);
933            min_area = min_area.min(area);
934            if area < area_epsilon {
935                deg_count += 1;
936            }
937            // Check orientation: if normal's y < 0 after projecting, might be inverted
938            if cross[1] < 0.0 {
939                inv_count += 1;
940            }
941        }
942
943        if triangles.is_empty() {
944            min_ar = 0.0;
945            max_ar = 0.0;
946            min_area = 0.0;
947        }
948
949        MeshQuality {
950            min_aspect_ratio: min_ar,
951            max_aspect_ratio: max_ar,
952            min_area,
953            degenerate_count: deg_count,
954            inverted_count: inv_count,
955        }
956    }
957
958    /// Check that all element indices are within bounds.
959    pub fn check_bounds(vertices: &[[f64; 3]], triangles: &[[u32; 3]]) -> bool {
960        let n = vertices.len() as u32;
961        triangles.iter().all(|&[a, b, c]| a < n && b < n && c < n)
962    }
963}
964
965fn sub3(a: [f64; 3], b: [f64; 3]) -> [f64; 3] {
966    [a[0] - b[0], a[1] - b[1], a[2] - b[2]]
967}
968
969fn cross3(a: [f64; 3], b: [f64; 3]) -> [f64; 3] {
970    [
971        a[1] * b[2] - a[2] * b[1],
972        a[2] * b[0] - a[0] * b[2],
973        a[0] * b[1] - a[1] * b[0],
974    ]
975}
976
977fn len3(a: [f64; 3]) -> f64 {
978    (a[0] * a[0] + a[1] * a[1] + a[2] * a[2]).sqrt()
979}
980
981/// Compute the aspect ratio of a triangle given edge lengths.
982fn triangle_aspect_ratio(a: f64, b: f64, c: f64) -> f64 {
983    let s = (a + b + c) / 2.0;
984    let area = (s * (s - a) * (s - b) * (s - c)).max(0.0).sqrt();
985    if area < 1e-14 {
986        return f64::MAX;
987    }
988    let longest = a.max(b).max(c);
989    longest * (a + b + c) / (4.0 * 3_f64.sqrt() * area)
990}
991
992// ---------------------------------------------------------------------------
993// Tests
994// ---------------------------------------------------------------------------
995
996#[cfg(test)]
997mod tests {
998    use super::*;
999
1000    // --- base64_encode_buffer ---
1001
1002    #[test]
1003    fn base64_empty() {
1004        assert_eq!(base64_encode_buffer(&[]), "");
1005    }
1006
1007    #[test]
1008    fn base64_man() {
1009        // "Man" → "TWFu"
1010        let result = base64_encode_buffer(b"Man");
1011        assert_eq!(result, "TWFu");
1012    }
1013
1014    #[test]
1015    fn base64_one_byte_padded() {
1016        let result = base64_encode_buffer(b"M");
1017        assert_eq!(result.len(), 4);
1018        assert!(result.ends_with("=="));
1019    }
1020
1021    #[test]
1022    fn base64_two_bytes_padded() {
1023        let result = base64_encode_buffer(b"Ma");
1024        assert_eq!(result.len(), 4);
1025        assert!(result.ends_with('='));
1026    }
1027
1028    // --- nastran_field8 ---
1029
1030    #[test]
1031    fn nastran_field8_length() {
1032        let s = nastran_field8(1.23456789);
1033        assert!(s.len() <= 8, "length={}", s.len());
1034    }
1035
1036    #[test]
1037    fn nastran_field8_zero() {
1038        let s = nastran_field8(0.0);
1039        assert!(s.len() <= 8);
1040    }
1041
1042    // --- gmsh_element_type_str ---
1043
1044    #[test]
1045    fn gmsh_triangle() {
1046        assert_eq!(gmsh_element_type_str(3), "Triangle");
1047    }
1048
1049    #[test]
1050    fn gmsh_line() {
1051        assert_eq!(gmsh_element_type_str(2), "Line");
1052    }
1053
1054    // --- GltfMeshWriter ---
1055
1056    #[test]
1057    fn gltf_writer_produces_json() {
1058        let verts = vec![[0.0f64, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]];
1059        let indices = vec![[0u32, 1, 2]];
1060        let w = GltfMeshWriter::new(verts, vec![], vec![], indices);
1061        let mut buf = Vec::new();
1062        w.write_gltf_json(&mut buf).unwrap();
1063        let s = String::from_utf8(buf).unwrap();
1064        assert!(s.contains("asset"), "s={s}");
1065        assert!(s.contains("accessors"));
1066    }
1067
1068    #[test]
1069    fn gltf_positions_bytes_length() {
1070        let verts = vec![[0.0f64; 3]; 4];
1071        let w = GltfMeshWriter::new(verts, vec![], vec![], vec![]);
1072        assert_eq!(w.positions_f32_bytes().len(), 48); // 4 * 3 * 4
1073    }
1074
1075    #[test]
1076    fn gltf_indices_bytes_length() {
1077        let indices = vec![[0u32, 1, 2], [1u32, 2, 3]];
1078        let w = GltfMeshWriter::new(vec![], vec![], vec![], indices);
1079        assert_eq!(w.indices_u32_bytes().len(), 24); // 2 * 3 * 4
1080    }
1081
1082    // --- GltfMeshReader ---
1083
1084    #[test]
1085    fn gltf_reader_parse_names() {
1086        let json = r#"{"meshes": [{"name": "MyMesh"}]}"#;
1087        let r = GltfMeshReader::parse_json(json);
1088        assert!(
1089            r.node_names.iter().any(|n| n.contains("MyMesh")),
1090            "names={:?}",
1091            r.node_names
1092        );
1093    }
1094
1095    // --- UsdMeshWriter ---
1096
1097    #[test]
1098    fn usd_writer_contains_def_mesh() {
1099        let verts = vec![[0.0f64, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]];
1100        let w = UsdMeshWriter::from_triangles(verts, &[[0, 1, 2]], "/World/Mesh");
1101        let mut buf = Vec::new();
1102        w.write(&mut buf).unwrap();
1103        let s = String::from_utf8(buf).unwrap();
1104        assert!(s.contains("def Mesh"), "s={s}");
1105        assert!(s.contains("points"));
1106    }
1107
1108    #[test]
1109    fn usd_writer_face_vertex_indices() {
1110        let verts = vec![[0.0f64; 3]; 3];
1111        let w = UsdMeshWriter::from_triangles(verts, &[[0, 1, 2]], "/Mesh");
1112        assert_eq!(w.face_vertex_indices.len(), 3);
1113        assert_eq!(w.face_vertex_counts.len(), 1);
1114    }
1115
1116    // --- MshWriter / MshReader ---
1117
1118    #[test]
1119    fn msh_writer_v2_roundtrip() {
1120        let mut mw = MshWriter::new(2);
1121        mw.add_node(0.0, 0.0, 0.0);
1122        mw.add_node(1.0, 0.0, 0.0);
1123        mw.add_node(0.0, 1.0, 0.0);
1124        mw.add_triangle(1, 2, 3, 1);
1125        let mut buf = Vec::new();
1126        mw.write_v2(&mut buf).unwrap();
1127        let s = String::from_utf8(buf).unwrap();
1128        assert!(s.contains("$MeshFormat"));
1129        assert!(s.contains("$Nodes"));
1130        assert!(s.contains("$Elements"));
1131        let r = MshReader::parse(&s);
1132        assert_eq!(r.nodes.len(), 3, "parsed nodes: {:?}", r.nodes);
1133        assert_eq!(r.elements.len(), 1);
1134    }
1135
1136    #[test]
1137    fn msh_reader_node_coords() {
1138        let msh = "$MeshFormat\n2.2 0 8\n$EndMeshFormat\n$Nodes\n2\n1 1.0 2.0 3.0\n2 4.0 5.0 6.0\n$EndNodes\n$Elements\n0\n$EndElements\n";
1139        let r = MshReader::parse(msh);
1140        assert_eq!(r.nodes.len(), 2);
1141        assert!((r.nodes[0][0] - 1.0).abs() < 1e-9);
1142    }
1143
1144    // --- AbaqusInpWriter / Reader ---
1145
1146    #[test]
1147    fn abaqus_writer_produces_part() {
1148        let mut w = AbaqusInpWriter::new("PART-1", "C3D4");
1149        w.add_node(0.0, 0.0, 0.0);
1150        w.add_node(1.0, 0.0, 0.0);
1151        w.add_node(0.0, 1.0, 0.0);
1152        w.add_node(0.0, 0.0, 1.0);
1153        w.add_element(vec![1, 2, 3, 4]);
1154        let mut buf = Vec::new();
1155        w.write(&mut buf).unwrap();
1156        let s = String::from_utf8(buf).unwrap();
1157        assert!(s.contains("*PART"), "s={s}");
1158        assert!(s.contains("*NODE"));
1159        assert!(s.contains("*ELEMENT"));
1160    }
1161
1162    #[test]
1163    fn abaqus_reader_parse() {
1164        let inp = "*PART, NAME=P1\n*NODE\n1, 0.0, 0.0, 0.0\n2, 1.0, 0.0, 0.0\n3, 0.0, 1.0, 0.0\n*ELEMENT, TYPE=C3D4\n1, 1, 2, 3\n*END PART\n";
1165        let r = AbaqusInpReader::parse(inp);
1166        assert_eq!(r.nodes.len(), 3, "nodes={:?}", r.nodes);
1167        assert_eq!(r.elements.len(), 1);
1168    }
1169
1170    #[test]
1171    fn abaqus_boundary_written() {
1172        let mut w = AbaqusInpWriter::new("P", "C3D4");
1173        w.add_node_set("FIXED", vec![1, 2]);
1174        w.add_boundary("FIXED", 1, 3, 0.0);
1175        let mut buf = Vec::new();
1176        w.write(&mut buf).unwrap();
1177        let s = String::from_utf8(buf).unwrap();
1178        assert!(s.contains("*BOUNDARY"), "s={s}");
1179    }
1180
1181    // --- NastranBdf ---
1182
1183    #[test]
1184    fn nastran_write_and_parse() {
1185        let mut bdf = NastranBdf::new();
1186        bdf.add_grid(1, [0.0, 0.0, 0.0]);
1187        bdf.add_grid(2, [1.0, 0.0, 0.0]);
1188        bdf.add_element(1, 1, vec![1, 2, 3, 4], "CTETRA");
1189        let mut buf = Vec::new();
1190        bdf.write(&mut buf).unwrap();
1191        let s = String::from_utf8(buf).unwrap();
1192        assert!(s.contains("BEGIN BULK"), "s={s}");
1193        assert!(s.contains("GRID"));
1194        assert!(s.contains("ENDDATA"));
1195    }
1196
1197    #[test]
1198    fn nastran_parse_grids() {
1199        let bdf_str = "BEGIN BULK\nGRID         1       0    1.00    2.00    3.00\nENDDATA\n";
1200        let bdf = NastranBdf::parse(bdf_str);
1201        // Simplified parser; just ensure no panic
1202        let _ = bdf.grids.len();
1203    }
1204
1205    // --- MeshConverter ---
1206
1207    #[test]
1208    fn compact_nodes_removes_unused() {
1209        let nodes = vec![[0.0f64; 3]; 5];
1210        let elements = vec![vec![0, 2, 4]];
1211        let (new_nodes, new_elements) = MeshConverter::compact_nodes(&nodes, &elements);
1212        assert_eq!(new_nodes.len(), 3);
1213        assert_eq!(new_elements[0].len(), 3);
1214    }
1215
1216    #[test]
1217    fn deduplicate_identical_vertices() {
1218        let positions = vec![[0.0f64, 0.0, 0.0], [0.0, 0.0, 0.0], [1.0, 0.0, 0.0]];
1219        let indices = vec![[0u32, 1, 2]];
1220        let (uniq, _new_idx) = MeshConverter::deduplicate_vertices(&positions, &indices, 1e-6);
1221        assert_eq!(uniq.len(), 2);
1222    }
1223
1224    #[test]
1225    fn flip_winding_reverses_order() {
1226        let indices = vec![[0u32, 1, 2]];
1227        let flipped = MeshConverter::flip_winding(&indices);
1228        assert_eq!(flipped[0], [0, 2, 1]);
1229    }
1230
1231    // --- MeshValidator ---
1232
1233    #[test]
1234    fn validator_equilateral_good_aspect_ratio() {
1235        let sq3 = 3_f64.sqrt();
1236        let verts = vec![[0.0f64, 0.0, 0.0], [1.0, 0.0, 0.0], [0.5, sq3 / 2.0, 0.0]];
1237        let tris = vec![[0u32, 1, 2]];
1238        let q = MeshValidator::validate(&verts, &tris, 1e-10);
1239        assert!(q.min_aspect_ratio < 1.1, "ar={}", q.min_aspect_ratio);
1240    }
1241
1242    #[test]
1243    fn validator_degenerate_triangle() {
1244        let verts = vec![[0.0f64; 3]; 3]; // all at origin
1245        let tris = vec![[0u32, 1, 2]];
1246        let q = MeshValidator::validate(&verts, &tris, 1e-10);
1247        assert_eq!(q.degenerate_count, 1);
1248    }
1249
1250    #[test]
1251    fn validator_check_bounds_ok() {
1252        let verts = vec![[0.0f64; 3]; 3];
1253        let tris = vec![[0u32, 1, 2]];
1254        assert!(MeshValidator::check_bounds(&verts, &tris));
1255    }
1256
1257    #[test]
1258    fn validator_check_bounds_fail() {
1259        let verts = vec![[0.0f64; 3]; 2];
1260        let tris = vec![[0u32, 1, 5]];
1261        assert!(!MeshValidator::check_bounds(&verts, &tris));
1262    }
1263
1264    // --- write_gltf_accessor ---
1265
1266    #[test]
1267    fn write_gltf_accessor_correct() {
1268        let mut buf = Vec::new();
1269        write_gltf_accessor(&mut buf, 0, 5126, 100, "VEC3", false).unwrap();
1270        let s = String::from_utf8(buf).unwrap();
1271        assert!(s.contains("VEC3"));
1272        assert!(s.contains("100"));
1273    }
1274
1275    // --- triangle_aspect_ratio ---
1276
1277    #[test]
1278    fn aspect_ratio_equilateral_is_one() {
1279        let ar = triangle_aspect_ratio(1.0, 1.0, 1.0);
1280        assert!((ar - 1.0).abs() < 1e-6, "ar={ar}");
1281    }
1282
1283    #[test]
1284    fn aspect_ratio_degenerate_is_max() {
1285        let ar = triangle_aspect_ratio(1.0, 1.0, 2.0);
1286        assert!(ar > 1.0 || ar == f64::MAX);
1287    }
1288}
oxiphysics_io/mesh_io.rs

oxiphysics_io/
mesh_io.rs
