Skip to main content

cmtool_core/ensight_gold/
geo.rs

1// SPDX-License-Identifier: GPL-3.0-or-later
2
3use crate::coordinates::*;
4use crate::{
5    ensight_gold::{Reader, reader::EnsightGoldReader, types::ElementsType},
6    utils,
7};
8use std::{
9    io::{Error, ErrorKind},
10    path::Path,
11    str::FromStr,
12};
13
14const MAXIMAL_NUMBER_OF_MESH_ELEMENT_TYPE: usize = 20;
15const MAXIMAL_NUMBER_OF_PART: usize = 20;
16
17///Marks "the section holds no further item", as opposed to a read failure. Sections are only
18///delimited by what the next line holds, so the readers have to report the end as an error and
19///the callers must not confuse it with a truncated or unreadable file.
20#[derive(Debug)]
21struct SectionEnd;
22
23impl std::fmt::Display for SectionEnd {
24    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
25        write!(f, "end of section")
26    }
27}
28
29impl std::error::Error for SectionEnd {}
30
31fn section_end() -> Error {
32    Error::other(SectionEnd)
33}
34
35fn is_section_end(error: &Error) -> bool {
36    error
37        .get_ref()
38        .is_some_and(|inner| inner.is::<SectionEnd>())
39}
40
41///The first line of a section is the only place where the file is allowed to end
42fn read_section_header(reader: &mut EnsightGoldReader) -> std::io::Result<String> {
43    match reader.get_line_string() {
44        Ok(line) => Ok(line),
45        Err(error) if error.kind() == ErrorKind::UnexpectedEof => Err(section_end()),
46        Err(error) => Err(error),
47    }
48}
49
50#[derive(Debug, Default, Clone)]
51pub(crate) struct MeshElementType {
52    pub n_nodes: usize,
53    pub(crate) n_elements: usize,
54    pub etype: ElementsType,
55    pub vertices: Vec<usize>,
56}
57
58pub struct Part {
59    id: u32,
60    name: String,
61    vertex_coordinates: Vec<f64>,
62    pub n_vertex: usize,
63    pub(crate) elements: Vec<MeshElementType>,
64}
65
66impl std::fmt::Display for Part {
67    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
68        writeln!(f, "Part #{}: {}", self.id, self.name)?;
69        writeln!(f, "  Number of vertices: {}", self.n_vertex)?;
70        writeln!(f, "  Vertex coordinates: {}", self.vertex_coordinates.len())?;
71        writeln!(f, "  Number of elements: {}", self.elements.len())?;
72
73        Ok(())
74    }
75}
76
77#[derive(Debug)]
78pub struct Geometry {
79    pub(crate) parts: Vec<Part>,
80}
81
82impl std::fmt::Display for Geometry {
83    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
84        writeln!(f, "Geometry with {} parts", self.parts.len())?;
85        for p in &self.parts {
86            writeln!(f, "{}", p)?;
87        }
88
89        Ok(())
90    }
91}
92
93impl MeshElementType {
94    pub fn get_vertex(&self, i_element: usize, i_vertex: usize) -> usize {
95        self.vertices[i_element * self.n_nodes + i_vertex]
96    }
97
98    pub fn read(reader: &mut EnsightGoldReader, ignore_element_id: bool) -> std::io::Result<Self> {
99        let element_type = read_section_header(reader)?;
100
101        //The next part starts here, this one has no more elements
102        if element_type.contains("part") {
103            reader.rollback()?;
104            return Err(section_end());
105        }
106
107        let etype = ElementsType::from_str(&element_type)
108            .map_err(|_| Error::new(ErrorKind::Unsupported, "Missing 'ElementsType' in header"))?;
109
110        let n_elements = reader.read_i32()? as usize;
111
112        if ignore_element_id {
113            for _ in 0..n_elements {
114                reader.ignore_line()?;
115            }
116        }
117        //Regular
118        if etype != ElementsType::Nfaced {
119            // let n_vertex = etype.number_vertex();
120            let n_nodes = etype.node_count() as usize;
121            let mut max = 0;
122            let mut vertices = Vec::with_capacity(n_elements * n_nodes);
123            for _ in 0..n_elements * n_nodes {
124                let cn = reader.read_i32()?;
125                vertices.push(cn as usize);
126                if cn > max {
127                    max = cn;
128                }
129            }
130            Ok(Self {
131                n_nodes,
132                n_elements,
133                etype,
134                vertices,
135            })
136        } else {
137            let n_nodes = etype.node_count() as usize;
138            for _ in 0..n_nodes * n_elements {
139                println!("{}", reader.read_i32()? as usize);
140            }
141            todo!()
142        }
143    }
144}
145
146impl std::fmt::Debug for Part {
147    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
148        write!(f, "Part {} with  {} elements", self.id, self.elements.len())
149    }
150}
151
152impl Part {
153    pub fn get_vertex_coordinates(&self, k_vertex: usize, k_xyz: usize) -> f64 {
154        self.vertex_coordinates[utils::linear_index_coordinates_matrix(k_vertex, k_xyz)]
155    }
156
157    pub fn get_vertex_coordinates_vec(&self, k_vertex: usize) -> Coords3 {
158        // [
159        //     self.vertex_coordinates[utils::linear_index_coordinates_matrix(k_vertex, 0)],
160        //     self.vertex_coordinates[utils::linear_index_coordinates_matrix(k_vertex, 1)],
161        //     self.vertex_coordinates[utils::linear_index_coordinates_matrix(k_vertex, 2)],
162        // ]
163        let offset = k_vertex * 3;
164        [
165            self.vertex_coordinates[offset],
166            self.vertex_coordinates[offset + 1],
167            self.vertex_coordinates[offset + 2],
168        ]
169    }
170
171    pub fn get_vertex_coordinates_slice(&self, k_vertex: usize) -> &Coords3 {
172        let offset = k_vertex * 3;
173        self.vertex_coordinates[offset..offset + 3]
174            .try_into()
175            .expect("Part slice vertex")
176    }
177    pub fn set_vertex_coordinates(&mut self, k_vertex: usize, k_xyz: usize, value: f64) {
178        self.vertex_coordinates[utils::linear_index_coordinates_matrix(k_vertex, k_xyz)] = value;
179    }
180
181    pub fn read(
182        reader: &mut EnsightGoldReader,
183        ignore_node_id: bool,
184        ignore_element_id: bool,
185    ) -> std::io::Result<Self> {
186        let buf = read_section_header(reader)?;
187
188        //Not a part header: the geometry holds no further part
189        if !buf.contains("part") {
190            return Err(section_end());
191        }
192
193        let id = reader.read_i32()? as u32;
194        reader.ignore_line()?;
195
196        reader.check_lines_contains("coordinates")?;
197
198        let n_nodes = reader.read_i32()? as usize;
199
200        if ignore_node_id {
201            for _ in 0..n_nodes {
202                reader.ignore_line()?;
203            }
204        }
205
206        let mut current_part = Self {
207            id,
208            name: "part".to_string(),
209            vertex_coordinates: vec![0.; n_nodes * 3],
210            elements: Vec::with_capacity(MAXIMAL_NUMBER_OF_MESH_ELEMENT_TYPE),
211            n_vertex: n_nodes,
212        };
213        for i_xyz in 0..3 {
214            for i_vertex in 0..n_nodes {
215                current_part.vertex_coordinates
216                    [utils::linear_index_coordinates_matrix(i_vertex, i_xyz)] =
217                    reader.read_f32()? as f64;
218            }
219        }
220
221        loop {
222            match MeshElementType::read(reader, ignore_element_id) {
223                Ok(element) => current_part.elements.push(element),
224                Err(error) if is_section_end(&error) => break,
225                Err(error) => return Err(error),
226            }
227        }
228
229        Ok(current_part)
230    }
231}
232
233impl Geometry {
234    pub fn new(path: &Path) -> std::io::Result<Self> {
235        let mut reader = Reader::new(path)?;
236        Self::read(&mut reader)
237    }
238
239    pub fn number_of_part(&self) -> usize {
240        self.parts.len()
241    }
242
243    pub fn get_part_by_id(&self, id: u32) -> Option<&Part> {
244        self.parts.iter().find(|&p| p.id == id)
245    }
246
247    // pub fn get_part_position(&self, part:&Part) -> Option<usize> {
248    //     self.parts.iter().position(|p| p.id == part.id)
249    // }
250
251    pub fn read(reader: &mut EnsightGoldReader) -> std::io::Result<Self> {
252        //SKIP
253        for _ in 0..3 {
254            reader.ignore_line()?;
255        }
256
257        let (ignore_node_id, ignore_element_id) = {
258            let _node_id_choice = reader.get_line_string()?;
259
260            let ignore_node_id = false; //TODO find in node_id_choice: assign
261
262            let _element_id_choice = reader.get_line_string()?;
263            let ignore_element_id = false; //TODO find in element_id_choice: assign
264            (ignore_node_id, ignore_element_id)
265        };
266
267        // println!("{} {}", node_id_choice, element_id_choice);
268
269        let buf = reader.get_line_string()?;
270        if buf.contains("extents") {
271            println!("IGNORE EXTENTS");
272            reader.ignore_bytes(6)?;
273        } else {
274            reader.rollback()?;
275        }
276
277        let mut parts: Vec<Part> = Vec::with_capacity(MAXIMAL_NUMBER_OF_PART);
278
279        loop {
280            match Part::read(reader, ignore_node_id, ignore_element_id) {
281                Ok(part) => parts.push(part),
282                Err(error) if is_section_end(&error) => break,
283                Err(error) => return Err(error),
284            }
285        }
286
287        if reader.check_eof()? {
288            Ok(Self { parts })
289        } else {
290            Err(Error::new(
291                ErrorKind::InvalidData,
292                format!(
293                    "{}: trailing data after the last part, {} part(s) read",
294                    reader.path().display(),
295                    parts.len()
296                ),
297            ))
298        }
299    }
300}
301
302#[cfg(test)]
303mod test {
304    use super::*;
305
306    const LINE_SIZE: usize = 80;
307    const N_NODES: i32 = 4;
308
309    fn line(content: &str) -> Vec<u8> {
310        let mut buffer = vec![0u8; LINE_SIZE];
311        buffer[..content.len()].copy_from_slice(content.as_bytes());
312        buffer
313    }
314
315    ///One part of four nodes holding a single tetra, the smallest geometry the reader accepts.
316    ///Returns the file and the offsets of the coordinate and vertex blocks.
317    fn minimal_geometry() -> (Vec<u8>, usize, usize) {
318        let mut bytes = Vec::new();
319
320        for _ in 0..3 {
321            bytes.extend(line("description"));
322        }
323        bytes.extend(line("node id off"));
324        bytes.extend(line("element id off"));
325
326        //Read once while looking for "extents", then rolled back and read as the part header
327        bytes.extend(line("part"));
328        bytes.extend(1i32.to_le_bytes());
329        bytes.extend(line("a part"));
330        bytes.extend(line("coordinates"));
331        bytes.extend(N_NODES.to_le_bytes());
332
333        let coordinates_offset = bytes.len();
334        for value in 0..3 * N_NODES {
335            bytes.extend((value as f32).to_le_bytes());
336        }
337
338        bytes.extend(line("tetra4"));
339        bytes.extend(1i32.to_le_bytes());
340
341        let vertices_offset = bytes.len();
342        for vertex in 0..N_NODES {
343            bytes.extend(vertex.to_le_bytes());
344        }
345
346        (bytes, coordinates_offset, vertices_offset)
347    }
348
349    fn read_geometry(name: &str, bytes: &[u8]) -> std::io::Result<Geometry> {
350        let path = std::path::PathBuf::from("/tmp").join(name);
351        std::fs::write(&path, bytes).unwrap();
352        let geometry = Geometry::new(&path);
353        let _ = std::fs::remove_file(&path);
354        geometry
355    }
356
357    #[test]
358    fn test_read_minimal_geometry() {
359        let (bytes, _, _) = minimal_geometry();
360
361        let geometry = read_geometry("test_geo_minimal.geo", &bytes).expect("geometry");
362
363        assert_eq!(geometry.number_of_part(), 1);
364        assert_eq!(geometry.parts[0].n_vertex, N_NODES as usize);
365        assert_eq!(geometry.parts[0].elements.len(), 1);
366        assert_eq!(geometry.parts[0].elements[0].n_elements, 1);
367    }
368
369    ///A file cut inside the coordinates used to be reported as a geometry without any part
370    #[test]
371    fn test_truncated_coordinates_is_an_error() {
372        let (bytes, coordinates_offset, _) = minimal_geometry();
373
374        let result = read_geometry(
375            "test_geo_truncated_coordinates.geo",
376            &bytes[..coordinates_offset + 8],
377        );
378
379        assert!(result.is_err(), "truncated coordinates read as a geometry");
380    }
381
382    ///A file cut inside the element vertices used to yield a part without its elements
383    #[test]
384    fn test_truncated_elements_is_an_error() {
385        let (bytes, _, vertices_offset) = minimal_geometry();
386
387        let result = read_geometry(
388            "test_geo_truncated_elements.geo",
389            &bytes[..vertices_offset + 4],
390        );
391
392        assert!(result.is_err(), "truncated elements read as a geometry");
393    }
394}