ndgrid/io/
gmsh.rs

1//! Gmsh I/O
2
3use crate::traits::{Builder, Entity, Geometry, GmshExport, GmshImport, Grid, Point};
4use itertools::izip;
5use ndelement::types::ReferenceCellType;
6use num::Zero;
7use num::traits::FromBytes;
8use std::collections::HashMap;
9use std::fs::File;
10use std::io::{BufRead, BufReader, Read};
11use std::str::FromStr;
12
13fn get_permutation_to_gmsh(cell_type: ReferenceCellType, degree: usize) -> Vec<usize> {
14    match cell_type {
15        ReferenceCellType::Triangle => match degree {
16            1 => vec![0, 1, 2],
17            2 => vec![0, 1, 2, 5, 3, 4],
18            3 => vec![0, 1, 2, 7, 8, 3, 4, 6, 5, 9],
19            4 => vec![0, 1, 2, 9, 10, 11, 3, 4, 5, 8, 7, 6, 12, 13, 14],
20            5 => vec![
21                0, 1, 2, 11, 12, 13, 14, 3, 4, 5, 6, 10, 9, 8, 7, 15, 16, 17, 18, 19, 20,
22            ],
23            _ => {
24                panic!("Unsupported degree");
25            }
26        },
27        ReferenceCellType::Quadrilateral => match degree {
28            1 => vec![0, 1, 3, 2],
29            2 => vec![0, 1, 3, 2, 4, 6, 7, 5, 8],
30            _ => {
31                panic!("Unsupported degree");
32            }
33        },
34        ReferenceCellType::Tetrahedron => match degree {
35            1 => vec![0, 1, 2, 3],
36            _ => {
37                panic!("Unsupported degree");
38            }
39        },
40        ReferenceCellType::Hexahedron => match degree {
41            1 => vec![0, 1, 3, 2, 4, 5, 7, 6],
42            _ => {
43                panic!("Unsupported degree");
44            }
45        },
46        _ => {
47            panic!("Unsupported cell type.");
48        }
49    }
50}
51
52fn get_gmsh_cell(cell_type: ReferenceCellType, degree: usize) -> usize {
53    match cell_type {
54        ReferenceCellType::Triangle => match degree {
55            1 => 2,
56            2 => 9,
57            3 => 21,
58            4 => 23,
59            5 => 25,
60            _ => {
61                panic!("Unsupported degree");
62            }
63        },
64        ReferenceCellType::Quadrilateral => match degree {
65            1 => 3,
66            2 => 10,
67            _ => {
68                panic!("Unsupported degree");
69            }
70        },
71        ReferenceCellType::Tetrahedron => match degree {
72            1 => 4,
73            _ => {
74                panic!("Unsupported degree");
75            }
76        },
77        ReferenceCellType::Hexahedron => match degree {
78            1 => 5,
79            _ => {
80                panic!("Unsupported degree");
81            }
82        },
83        _ => {
84            panic!("Unsupported cell type.");
85        }
86    }
87}
88
89fn interpret_gmsh_cell(gmsh_cell: usize) -> (ReferenceCellType, usize) {
90    match gmsh_cell {
91        2 => (ReferenceCellType::Triangle, 1),
92        9 => (ReferenceCellType::Triangle, 2),
93        21 => (ReferenceCellType::Triangle, 3),
94        23 => (ReferenceCellType::Triangle, 4),
95        25 => (ReferenceCellType::Triangle, 5),
96        3 => (ReferenceCellType::Quadrilateral, 1),
97        10 => (ReferenceCellType::Quadrilateral, 2),
98        4 => (ReferenceCellType::Tetrahedron, 1),
99        5 => (ReferenceCellType::Hexahedron, 1),
100        _ => {
101            panic!("Unsupported cell type.");
102        }
103    }
104}
105
106/// Convert gmsh element node ordering to ndgrid ordering
107fn permute_gmsh_to_ndgrid(
108    gmsh_cell_type: usize,
109    node_tags: &[usize],
110) -> (Vec<usize>, ReferenceCellType, usize) {
111    let (cell_type, degree) = interpret_gmsh_cell(gmsh_cell_type);
112    let gmsh_perm = get_permutation_to_gmsh(cell_type, degree);
113    let mut cell = vec![0; node_tags.len()];
114    for (i, j) in gmsh_perm.iter().enumerate() {
115        cell[*j] = node_tags[i];
116    }
117    (cell, cell_type, degree)
118}
119
120/// Get the expected number of nodes for a gmsh cell type
121fn gmsh_node_count(gmsh_cell_type: usize) -> usize {
122    let (cell_type, degree) = interpret_gmsh_cell(gmsh_cell_type);
123    get_permutation_to_gmsh(cell_type, degree).len()
124}
125
126/// Parse coordinate strings into typed values
127fn parse_coordinates<T: FromStr>(values: &[&str]) -> Vec<T> {
128    values
129        .iter()
130        .map(|c| T::from_str(c).unwrap_or_else(|_| panic!("Could not parse coordinate: {c}")))
131        .collect()
132}
133
134/// Parse gmsh binary data
135fn parse<T: FromBytes>(data: &[u8], gmsh_data_size: usize, is_le: bool) -> T
136where
137    <T as FromBytes>::Bytes: Sized,
138{
139    let mut buf: T::Bytes = unsafe { std::mem::zeroed() };
140    let target_data_size = std::mem::size_of::<T::Bytes>();
141    let buf_bytes: &mut [u8] = unsafe {
142        std::slice::from_raw_parts_mut((&mut buf as *mut _) as *mut u8, target_data_size)
143    };
144
145    if is_le {
146        buf_bytes[..gmsh_data_size].copy_from_slice(data);
147        T::from_le_bytes(&buf)
148    } else {
149        buf_bytes[target_data_size - gmsh_data_size..].copy_from_slice(data);
150        T::from_be_bytes(&buf)
151    }
152}
153
154impl<G: Grid<EntityDescriptor = ReferenceCellType>> GmshExport for G {
155    fn to_gmsh_string(&self) -> String {
156        let tdim = self.topology_dim();
157        let gdim = self.geometry_dim();
158
159        let mut points = HashMap::new();
160        for cell_type in self.cell_types() {
161            for cell in self.entity_iter(*cell_type) {
162                for point in cell.geometry().points() {
163                    let mut p = vec![G::T::zero(); gdim];
164                    point.coords(&mut p);
165                    points.insert(point.index(), p);
166                }
167            }
168        }
169        let mut points = points.iter().collect::<Vec<_>>();
170        points.sort_by(|i, j| i.0.cmp(j.0));
171        let node_count = points.len();
172
173        let mut gmsh_s = String::from("");
174        gmsh_s.push_str("$MeshFormat\n");
175        gmsh_s.push_str("4.1 0 8\n");
176        gmsh_s.push_str("$EndMeshFormat\n");
177        gmsh_s.push_str("$Nodes\n");
178        gmsh_s.push_str(&format!("1 {node_count} 1 {node_count}\n"));
179        gmsh_s.push_str(&format!("{tdim} 1 0 {node_count}\n"));
180        for (i, _) in &points {
181            gmsh_s.push_str(&format!("{}\n", *i + 1));
182        }
183        for (_, coords) in &points {
184            for n in 0..3 {
185                if n != 0 {
186                    gmsh_s.push(' ');
187                }
188                gmsh_s.push_str(&format!(
189                    "{}",
190                    if let Some(c) = coords.get(n) {
191                        *c
192                    } else {
193                        G::T::zero()
194                    }
195                ));
196            }
197            gmsh_s.push('\n');
198        }
199        gmsh_s.push_str("$EndNodes\n");
200        gmsh_s.push_str("$Elements\n");
201
202        let cell_count = self
203            .entity_types(tdim)
204            .iter()
205            .map(|t| self.entity_count(*t))
206            .sum::<usize>();
207
208        let mut elements = vec![];
209        let mut cells_by_element = vec![];
210        for t in self.cell_types() {
211            for (index, cell) in self.entity_iter(*t).enumerate() {
212                let element = (cell.entity_type(), cell.geometry().degree());
213                if !elements.contains(&element) {
214                    elements.push(element);
215                    cells_by_element.push(vec![]);
216                }
217                cells_by_element[elements.iter().position(|i| *i == element).unwrap()].push(index);
218            }
219        }
220
221        gmsh_s.push_str(&format!("{} {cell_count} 1 {cell_count}\n", elements.len()));
222
223        let mut next_index = 1;
224        for ((cell_type, degree), cells) in izip!(elements, cells_by_element) {
225            let gmsh_perm = get_permutation_to_gmsh(cell_type, degree);
226
227            gmsh_s.push_str(&format!(
228                "{tdim} 1 {} {}\n",
229                get_gmsh_cell(cell_type, degree),
230                cells.len()
231            ));
232
233            for index in cells.iter() {
234                gmsh_s.push_str(&format!("{}", {
235                    let idx = next_index;
236                    next_index += 1;
237                    idx
238                },));
239                let entity = self.entity(cell_type, *index).unwrap();
240                let point_indices = entity
241                    .geometry()
242                    .points()
243                    .map(|i| i.index())
244                    .collect::<Vec<_>>();
245                for j in &gmsh_perm {
246                    gmsh_s.push_str(&format!(" {}", point_indices[*j] + 1));
247                }
248                gmsh_s.push('\n');
249            }
250        }
251
252        gmsh_s.push_str("$EndElements\n");
253
254        gmsh_s
255    }
256}
257
258/// Get a section from a gmsh string
259fn gmsh_section(s: &str, section: &str) -> String {
260    let a = s.split(&format!("${section}\n")).collect::<Vec<_>>();
261    if a.len() <= 1 {
262        panic!("Section not found: {section}");
263    }
264    String::from(a[1].split(&format!("\n$End{section}")).collect::<Vec<_>>()[0])
265}
266
267/// Import nodes from Gmsh v1/v2 ASCII format (shared between v1 and v2)
268fn import_nodes_v1_v2<T, B>(builder: &mut B, s: &str)
269where
270    T: FromStr,
271    B: Builder<T = T>,
272{
273    let nodes = gmsh_section(s, "Nodes");
274    let mut nodes = nodes.lines();
275
276    let Some(num_nodes) = nodes.next() else {
277        panic!("Could not read num nodes");
278    };
279    let num_nodes = num_nodes
280        .parse::<usize>()
281        .expect("Could not parse num nodes");
282
283    for _ in 0..num_nodes {
284        let Some(line) = nodes.next() else {
285            panic!("Unable to read node line");
286        };
287        let line = line.split(" ").collect::<Vec<&str>>();
288
289        let (tag, coord_strs) = line.split_at(1);
290        let tag = tag[0].parse::<usize>().expect("Could not parse node tag");
291        let coords: Vec<T> = parse_coordinates(coord_strs);
292
293        builder.add_point(tag, &coords[..builder.gdim()]);
294    }
295}
296
297impl<T, B> GmshImport for B
298where
299    T: FromStr + FromBytes + std::fmt::Debug,
300    <T as FromBytes>::Bytes: Sized,
301    B: Builder<T = T, EntityDescriptor = ReferenceCellType>,
302{
303    fn import_from_gmsh_v1(&mut self, s: String) {
304        import_nodes_v1_v2(self, &s);
305
306        // Load elements
307        let elements = gmsh_section(&s, "Elements");
308        let mut elements = elements.lines();
309
310        let Some(num_elements) = elements.next() else {
311            panic!("Could not read num nodes");
312        };
313        let num_elements = num_elements
314            .parse::<usize>()
315            .expect("Could not parse num nodes");
316
317        for _ in 0..num_elements {
318            let Some(line) = elements.next() else {
319                panic!("Unable to read element line");
320            };
321            let line = line.split(" ").collect::<Vec<&str>>();
322
323            let (a, node_tags) = line.split_at(5);
324            let [tag, element_type, ..] = a
325                .iter()
326                .map(|i| {
327                    i.parse::<usize>()
328                        .expect("Could not parse element tag and type")
329                })
330                .collect::<Vec<_>>()[..]
331            else {
332                panic!("Unrecognised format");
333            };
334
335            let node_tags = node_tags
336                .iter()
337                .map(|i| {
338                    i.parse::<usize>()
339                        .expect("Could not parse nodes for element")
340                })
341                .collect::<Vec<_>>();
342            let (cell, cell_type, degree) = permute_gmsh_to_ndgrid(element_type, &node_tags);
343            self.add_cell_from_nodes_and_type(tag, &cell, cell_type, degree);
344        }
345    }
346
347    fn import_from_gmsh_string_v2(&mut self, s: String) {
348        import_nodes_v1_v2(self, &s);
349
350        // Load elements
351        let elements = gmsh_section(&s, "Elements");
352        let mut elements = elements.lines();
353
354        let Some(num_elements) = elements.next() else {
355            panic!("Could not read num nodes");
356        };
357        let num_elements = num_elements
358            .parse::<usize>()
359            .expect("Could not parse num nodes");
360
361        for _ in 0..num_elements {
362            let Some(line) = elements.next() else {
363                panic!("Unable to read element line");
364            };
365            let line = line.split(" ").collect::<Vec<&str>>();
366
367            let (a, rem_line) = line.split_at(3);
368            let [tag, element_type, num_tags] = a
369                .iter()
370                .map(|i| {
371                    i.parse::<usize>()
372                        .expect("Could not parse element tag and type")
373                })
374                .collect::<Vec<_>>()[..]
375            else {
376                panic!("Unrecognised format");
377            };
378
379            let (_, node_tags) = rem_line.split_at(num_tags);
380            let node_tags = node_tags
381                .iter()
382                .map(|i| {
383                    i.parse::<usize>()
384                        .expect("Could not parse nodes for element")
385                })
386                .collect::<Vec<_>>();
387            let (cell, cell_type, degree) = permute_gmsh_to_ndgrid(element_type, &node_tags);
388            self.add_cell_from_nodes_and_type(tag, &cell, cell_type, degree);
389        }
390    }
391
392    fn import_from_gmsh_binary_v2(
393        &mut self,
394        mut reader: BufReader<File>,
395        data_size: usize,
396        is_le: bool,
397    ) {
398        let mut line = String::new();
399        let mut buf = Vec::new();
400
401        macro_rules! read_exact {
402            ($size: expr, $msg: expr) => {{
403                buf.resize($size, 0);
404                reader.read_exact(&mut buf).expect($msg);
405            }};
406        }
407
408        const GMSH_INT_SIZE: usize = 4;
409
410        loop {
411            let Ok(num_bytes) = reader.read_line(&mut line) else {
412                continue;
413            };
414
415            // EOF reached.
416            if num_bytes == 0 {
417                break;
418            }
419
420            match line.trim() {
421                // Load all nodes.
422                "$Nodes" => {
423                    line.clear();
424                    let Ok(_) = reader.read_line(&mut line) else {
425                        panic!("Unable to read num nodes");
426                    };
427                    let num_nodes = line
428                        .trim()
429                        .parse::<usize>()
430                        .expect("Could not parse num nodes");
431
432                    for _ in 0..num_nodes {
433                        read_exact!(GMSH_INT_SIZE, "Unable to read node tag");
434                        let tag = parse::<usize>(&buf, GMSH_INT_SIZE, is_le);
435
436                        read_exact!(3 * data_size, "Unable to read node coords");
437                        let coords = buf
438                            .chunks(data_size)
439                            .map(|i| parse::<T>(i, data_size, is_le))
440                            .collect::<Vec<_>>();
441
442                        self.add_point(tag, &coords[..self.gdim()]);
443                    }
444
445                    line.clear();
446                }
447                // Load all elements.
448                "$Elements" => {
449                    line.clear();
450                    let Ok(_) = reader.read_line(&mut line) else {
451                        panic!("Unable to read num elements")
452                    };
453                    let num_elements = line
454                        .trim()
455                        .parse::<usize>()
456                        .expect("Could not parse num elements");
457
458                    for _ in 0..num_elements {
459                        read_exact!(3 * GMSH_INT_SIZE, "Unable to element tag and type");
460                        let [elm_type, _num_elm_follow, num_tags] = buf
461                            .chunks(GMSH_INT_SIZE)
462                            .map(|i| parse::<usize>(i, GMSH_INT_SIZE, is_le))
463                            .collect::<Vec<_>>()[..]
464                        else {
465                            panic!("Could not parse element tag and type")
466                        };
467
468                        read_exact!(GMSH_INT_SIZE, "Unable to read num tags");
469                        let tag = parse::<usize>(&buf, 4, is_le);
470
471                        // Skip tags
472                        read_exact!(num_tags * GMSH_INT_SIZE, "Unable to read element tags");
473
474                        let (cell_type, degree) = interpret_gmsh_cell(elm_type);
475                        let gmsh_perm = get_permutation_to_gmsh(cell_type, degree);
476
477                        read_exact!(
478                            gmsh_perm.len() * GMSH_INT_SIZE,
479                            "Unable to read element node tags"
480                        );
481                        let node_tags = buf
482                            .chunks(GMSH_INT_SIZE)
483                            .map(|i| parse::<usize>(i, GMSH_INT_SIZE, is_le))
484                            .collect::<Vec<_>>();
485
486                        let (cell, cell_type, degree) =
487                            permute_gmsh_to_ndgrid(elm_type, &node_tags);
488                        self.add_cell_from_nodes_and_type(tag, &cell, cell_type, degree);
489                    }
490                    line.clear();
491                }
492                _ => {
493                    line.clear();
494                }
495            }
496        }
497    }
498
499    fn import_from_gmsh_string_v4(&mut self, s: String) {
500        // Load nodes
501        let nodes = gmsh_section(&s, "Nodes");
502        let nodes = nodes.lines().collect::<Vec<_>>();
503
504        let [num_entity_blocks, _num_nodes, _min_node_tag, _max_node_tag] = nodes[0]
505            .trim()
506            .split(" ")
507            .map(|i| i.parse::<usize>().unwrap())
508            .collect::<Vec<_>>()[..]
509        else {
510            panic!("Unrecognised gmsh format for node blocks");
511        };
512
513        let mut line_n = 1;
514        for _ in 0..num_entity_blocks {
515            let [entity_dim, _entity_tag, parametric, num_nodes_in_block] = nodes[line_n]
516                .trim()
517                .split(" ")
518                .map(|i| i.parse::<usize>().unwrap())
519                .collect::<Vec<_>>()[..]
520            else {
521                panic!("Unrecognised gmsh format for nodes");
522            };
523            line_n += 1;
524            let tags = &nodes[line_n..line_n + num_nodes_in_block];
525            let coords = &nodes[line_n + num_nodes_in_block..line_n + 2 * num_nodes_in_block];
526            for (t, c) in izip!(tags, coords) {
527                let coord_strs: Vec<&str> = c.trim().split(" ").collect();
528                let pt: Vec<T> = parse_coordinates(&coord_strs);
529                let tag = t.parse::<usize>().unwrap();
530                self.add_point(tag, &pt[..self.gdim()]);
531
532                if parametric == 1 && entity_dim >= 1 {
533                    let param_count = entity_dim.min(3);
534                    if pt.len() >= 3 + param_count {
535                        self.add_point_parametric_coords(tag, entity_dim, &pt[3..3 + param_count]);
536                    }
537                }
538            }
539            line_n += num_nodes_in_block * 2;
540        }
541
542        // Load elements
543        let elements = gmsh_section(&s, "Elements");
544        let elements = elements.lines().collect::<Vec<_>>();
545
546        let [
547            num_entity_blocks,
548            _num_elements,
549            _min_element_tag,
550            _max_element_tag,
551        ] = elements[0]
552            .trim()
553            .split(" ")
554            .map(|i| i.parse::<usize>().unwrap())
555            .collect::<Vec<_>>()[..]
556        else {
557            panic!("Unrecognised gmsh format");
558        };
559
560        let mut line_n = 1;
561        for _ in 0..num_entity_blocks {
562            let [
563                _entity_dim,
564                _entity_tag,
565                element_type,
566                num_elements_in_block,
567            ] = elements[line_n]
568                .trim()
569                .split(" ")
570                .map(|i| i.parse::<usize>().unwrap())
571                .collect::<Vec<_>>()[..]
572            else {
573                panic!("Unrecognised gmsh format");
574            };
575
576            line_n += 1;
577            for line in &elements[line_n..line_n + num_elements_in_block] {
578                let line = line
579                    .trim()
580                    .split(" ")
581                    .map(|i| i.parse::<usize>().unwrap())
582                    .collect::<Vec<_>>();
583                let node_tags = &line[1..];
584                let (cell, cell_type, degree) = permute_gmsh_to_ndgrid(element_type, node_tags);
585                self.add_cell_from_nodes_and_type(line[0], &cell, cell_type, degree);
586            }
587
588            line_n += num_elements_in_block;
589        }
590    }
591
592    fn import_from_gmsh_binary_v4(
593        &mut self,
594        mut reader: BufReader<File>,
595        data_size: usize,
596        is_le: bool,
597    ) {
598        let mut line = String::new();
599        let mut buf = Vec::new();
600
601        macro_rules! read_exact {
602            ($size: expr, $msg: expr) => {{
603                buf.resize($size, 0);
604                reader.read_exact(&mut buf).expect($msg);
605            }};
606        }
607
608        const GMSH_INT_SIZE: usize = 4;
609
610        loop {
611            let Ok(num_bytes) = reader.read_line(&mut line) else {
612                continue;
613            };
614
615            // EOF reached.
616            if num_bytes == 0 {
617                break;
618            }
619
620            match line.trim() {
621                // Load all nodes.
622                "$Nodes" => {
623                    read_exact!(4 * data_size, "Unable to read node section info");
624                    let [num_entity_blocks, _num_nodes, _min_node_tag, _max_node_tag] = buf
625                        .chunks(data_size)
626                        .map(|i| parse::<usize>(i, data_size, is_le))
627                        .collect::<Vec<_>>()[..]
628                    else {
629                        panic!("Could not parse node section info")
630                    };
631
632                    for _ in 0..num_entity_blocks {
633                        read_exact!(3 * GMSH_INT_SIZE, "Unable to read node entity block info");
634                        let [entity_dim, _entity_tag, parametric] = buf
635                            .chunks(GMSH_INT_SIZE)
636                            .map(|i| parse::<usize>(i, GMSH_INT_SIZE, is_le))
637                            .collect::<Vec<_>>()[..]
638                        else {
639                            panic!("Could not parse node entity block info")
640                        };
641
642                        read_exact!(data_size, "Unable to read num nodes in block");
643                        let num_nodes_in_block = parse::<usize>(&buf, data_size, is_le);
644
645                        read_exact!(num_nodes_in_block * data_size, "Unable to read node tags");
646                        let tags = buf
647                            .chunks(data_size)
648                            .map(|i| parse::<usize>(i, data_size, is_le))
649                            .collect::<Vec<_>>();
650
651                        let param_count = if parametric == 1 && entity_dim >= 1 {
652                            entity_dim.min(3)
653                        } else {
654                            0
655                        };
656                        let total_coords_per_node = 3 + param_count;
657
658                        read_exact!(
659                            total_coords_per_node * num_nodes_in_block * data_size,
660                            "Unable to read node coords"
661                        );
662                        let all_coords = buf
663                            .chunks(total_coords_per_node * data_size)
664                            .map(|i| {
665                                i.chunks(data_size)
666                                    .map(|j| parse::<T>(j, data_size, is_le))
667                                    .collect::<Vec<_>>()
668                            })
669                            .collect::<Vec<_>>();
670
671                        for (t, c) in izip!(tags, all_coords) {
672                            self.add_point(t, &c[..self.gdim()]);
673
674                            if parametric == 1 && entity_dim >= 1 && c.len() > 3 {
675                                self.add_point_parametric_coords(
676                                    t,
677                                    entity_dim,
678                                    &c[3..3 + param_count],
679                                );
680                            }
681                        }
682                    }
683
684                    line.clear();
685                }
686                // Load all elements.
687                "$Elements" => {
688                    read_exact!(4 * data_size, "Unable to read element section info");
689                    let [
690                        num_entity_blocks,
691                        _num_elements,
692                        _min_element_tag,
693                        _max_element_tag,
694                    ] = buf
695                        .chunks(data_size)
696                        .map(|i| parse::<usize>(i, data_size, is_le))
697                        .collect::<Vec<_>>()[..]
698                    else {
699                        panic!("Could not parse element section info")
700                    };
701
702                    for _ in 0..num_entity_blocks {
703                        read_exact!(
704                            3 * GMSH_INT_SIZE,
705                            "Unable to read element entity block info"
706                        );
707                        let [_entity_dim, _entity_tag, element_type] = buf
708                            .chunks(GMSH_INT_SIZE)
709                            .map(|i| parse::<usize>(i, GMSH_INT_SIZE, is_le))
710                            .collect::<Vec<_>>()[..]
711                        else {
712                            panic!("Could not parse element entity block info")
713                        };
714
715                        read_exact!(data_size, "Unable to read num elements in block");
716                        let num_elements_in_block = parse::<usize>(&buf, data_size, is_le);
717
718                        let node_count = gmsh_node_count(element_type);
719
720                        for _ in 0..num_elements_in_block {
721                            read_exact!(data_size, "Unable to read element tag");
722                            let tag = parse::<usize>(&buf, data_size, is_le);
723
724                            read_exact!(data_size * node_count, "Unable to read node tags");
725                            let node_tags = buf
726                                .chunks(data_size)
727                                .map(|i| parse::<usize>(i, data_size, is_le))
728                                .collect::<Vec<_>>();
729
730                            let (cell, cell_type, degree) =
731                                permute_gmsh_to_ndgrid(element_type, &node_tags);
732                            self.add_cell_from_nodes_and_type(tag, &cell, cell_type, degree);
733                        }
734                    }
735                    line.clear();
736                }
737                _ => {
738                    line.clear();
739                }
740            }
741        }
742    }
743}
744
745#[cfg(test)]
746mod test {
747    use super::*;
748    use crate::{
749        MixedGridBuilder, SingleElementGridBuilder,
750        shapes::regular_sphere,
751        traits::{Builder, Topology},
752    };
753    use approx::*;
754
755    #[test]
756    fn test_regular_sphere_gmsh_io() {
757        let g = regular_sphere::<f64>(2, ReferenceCellType::Triangle);
758        g.export_as_gmsh("_test_io_sphere.msh");
759
760        let g = regular_sphere::<f64>(2, ReferenceCellType::Quadrilateral);
761        g.export_as_gmsh("_test_io_sphere_quads.msh");
762    }
763
764    #[test]
765    fn test_export_quads() {
766        let mut b = SingleElementGridBuilder::<f64>::new(3, (ReferenceCellType::Quadrilateral, 1));
767        b.add_point(0, &[0.0, 0.0, 0.0]);
768        b.add_point(1, &[1.0, 0.0, 0.0]);
769        b.add_point(2, &[0.0, 1.0, 0.0]);
770        b.add_point(3, &[1.0, 1.0, 0.0]);
771        b.add_point(4, &[0.0, 0.0, 1.0]);
772        b.add_point(5, &[1.0, 0.0, 1.0]);
773        b.add_point(6, &[0.0, 1.0, 1.0]);
774        b.add_point(7, &[1.0, 1.0, 1.0]);
775        b.add_cell(0, &[0, 2, 1, 3]);
776        b.add_cell(1, &[0, 1, 4, 5]);
777        b.add_cell(2, &[0, 4, 2, 6]);
778        b.add_cell(3, &[1, 3, 5, 7]);
779        b.add_cell(4, &[2, 6, 3, 7]);
780        b.add_cell(5, &[4, 5, 6, 7]);
781        let g = b.create_grid();
782        g.export_as_gmsh("_test_io_cube.msh");
783    }
784
785    #[test]
786    fn test_export_tetrahedra() {
787        let mut b = SingleElementGridBuilder::<f64>::new(3, (ReferenceCellType::Tetrahedron, 1));
788        b.add_point(0, &[0.0, 0.0, 0.0]);
789        b.add_point(1, &[1.0, 0.0, 0.0]);
790        b.add_point(2, &[0.0, 1.0, 0.0]);
791        b.add_point(3, &[1.0, 1.0, 0.0]);
792        b.add_point(4, &[0.0, 0.0, 1.0]);
793        b.add_point(5, &[1.0, 0.0, 1.0]);
794        b.add_point(6, &[0.0, 1.0, 1.0]);
795        b.add_point(7, &[1.0, 1.0, 1.0]);
796        b.add_cell(0, &[0, 1, 5, 7]);
797        b.add_cell(1, &[0, 2, 6, 7]);
798        b.add_cell(2, &[0, 4, 5, 7]);
799        b.add_cell(3, &[0, 1, 3, 7]);
800        b.add_cell(4, &[0, 2, 3, 7]);
801        b.add_cell(5, &[0, 4, 6, 7]);
802        let g = b.create_grid();
803        g.export_as_gmsh("_test_io_tetrahedra.msh");
804    }
805
806    #[test]
807    fn test_hexahedra() {
808        let mut b = SingleElementGridBuilder::<f64>::new(3, (ReferenceCellType::Hexahedron, 1));
809        b.add_point(0, &[0.0, 0.0, 0.0]);
810        b.add_point(1, &[1.0, 0.0, 0.0]);
811        b.add_point(2, &[0.0, 1.0, 0.0]);
812        b.add_point(3, &[1.0, 1.0, 0.0]);
813        b.add_point(4, &[0.0, 0.0, 1.0]);
814        b.add_point(5, &[1.0, 0.0, 1.0]);
815        b.add_point(6, &[0.0, 1.0, 1.0]);
816        b.add_point(7, &[1.0, 1.0, 1.0]);
817        b.add_point(8, &[0.0, 0.0, 2.0]);
818        b.add_point(9, &[1.0, 0.0, 2.0]);
819        b.add_point(10, &[0.0, 1.0, 2.0]);
820        b.add_point(11, &[1.0, 1.0, 1.5]);
821        b.add_cell(1, &[0, 1, 2, 3, 4, 5, 6, 7]);
822        b.add_cell(2, &[4, 5, 6, 7, 8, 9, 10, 11]);
823        let g = b.create_grid();
824        g.export_as_gmsh("_test_io_hexahedra.msh");
825
826        let mut b = SingleElementGridBuilder::<f64>::new(3, (ReferenceCellType::Hexahedron, 1));
827        b.import_from_gmsh("_test_io_hexahedra.msh");
828        let g2 = b.create_grid();
829
830        let mut p1 = [0.0; 3];
831        let mut p2 = [0.0; 3];
832        for (v1, v2) in izip!(
833            g.entity_iter(ReferenceCellType::Point),
834            g2.entity_iter(ReferenceCellType::Point)
835        ) {
836            for (pt1, pt2) in izip!(v1.geometry().points(), v2.geometry().points()) {
837                pt1.coords(&mut p1);
838                pt2.coords(&mut p2);
839                for (c1, c2) in izip!(&p1, &p2) {
840                    assert_relative_eq!(c1, c2, epsilon = 1e-10);
841                }
842            }
843        }
844        for (h1, h2) in izip!(
845            g.entity_iter(ReferenceCellType::Hexahedron),
846            g2.entity_iter(ReferenceCellType::Hexahedron)
847        ) {
848            for (v1, v2) in izip!(
849                h1.topology().sub_entity_iter(ReferenceCellType::Point),
850                h2.topology().sub_entity_iter(ReferenceCellType::Point)
851            ) {
852                assert_eq!(v1, v2);
853            }
854        }
855    }
856
857    #[test]
858    fn test_high_order_triangles() {
859        let mut b = SingleElementGridBuilder::<f64>::new(3, (ReferenceCellType::Triangle, 5));
860        b.add_point(0, &[0.0, 0.0, 0.0]);
861        b.add_point(1, &[5.0, 0.0, 0.2]);
862        b.add_point(2, &[0.0, 5.0, 0.4]);
863        b.add_point(3, &[4.0, 1.0, -0.1]);
864        b.add_point(4, &[3.0, 2.0, 0.2]);
865        b.add_point(5, &[2.0, 3.0, -0.3]);
866        b.add_point(6, &[1.0, 4.0, -0.5]);
867        b.add_point(7, &[0.0, 1.0, 0.6]);
868        b.add_point(8, &[0.0, 2.0, 0.2]);
869        b.add_point(9, &[0.0, 3.0, 0.1]);
870        b.add_point(10, &[0.0, 4.0, -0.2]);
871        b.add_point(11, &[1.0, 0.0, -0.3]);
872        b.add_point(12, &[2.0, 0.0, -0.4]);
873        b.add_point(13, &[3.0, 0.0, -0.5]);
874        b.add_point(14, &[4.0, 0.0, -0.2]);
875        b.add_point(15, &[1.0, 1.0, 0.1]);
876        b.add_point(16, &[2.0, 1.0, 0.1]);
877        b.add_point(17, &[3.0, 1.0, 0.1]);
878        b.add_point(18, &[2.0, 1.0, 0.2]);
879        b.add_point(19, &[2.0, 2.0, 0.1]);
880        b.add_point(20, &[3.0, 1.0, 0.1]);
881        b.add_cell(
882            1,
883            &[
884                0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
885            ],
886        );
887        let g = b.create_grid();
888        g.export_as_gmsh("_test_io_high_order_triangle.msh");
889
890        let mut b = SingleElementGridBuilder::<f64>::new(3, (ReferenceCellType::Triangle, 5));
891        b.import_from_gmsh("_test_io_high_order_triangle.msh");
892        let g2 = b.create_grid();
893
894        let mut p1 = [0.0; 3];
895        let mut p2 = [0.0; 3];
896        for (v1, v2) in izip!(
897            g.entity_iter(ReferenceCellType::Point),
898            g2.entity_iter(ReferenceCellType::Point)
899        ) {
900            v1.geometry().points().next().unwrap().coords(&mut p1);
901            v2.geometry().points().next().unwrap().coords(&mut p2);
902        }
903        for (v1, v2) in izip!(
904            g.entity_iter(ReferenceCellType::Point),
905            g2.entity_iter(ReferenceCellType::Point)
906        ) {
907            v1.geometry().points().next().unwrap().coords(&mut p1);
908            v2.geometry().points().next().unwrap().coords(&mut p2);
909            for (c1, c2) in izip!(&p1, &p2) {
910                assert_relative_eq!(c1, c2, epsilon = 1e-10);
911            }
912        }
913        for (h1, h2) in izip!(
914            g.entity_iter(ReferenceCellType::Triangle),
915            g2.entity_iter(ReferenceCellType::Triangle)
916        ) {
917            for (v1, v2) in izip!(
918                h1.topology().sub_entity_iter(ReferenceCellType::Point),
919                h2.topology().sub_entity_iter(ReferenceCellType::Point)
920            ) {
921                assert_eq!(v1, v2);
922            }
923        }
924    }
925
926    #[test]
927    fn test_export_mixed() {
928        let mut b = MixedGridBuilder::<f64>::new(2);
929        b.add_point(0, &[0.0, 0.0]);
930        b.add_point(1, &[1.0, 0.0]);
931        b.add_point(2, &[0.0, 1.0]);
932        b.add_point(3, &[1.0, 1.0]);
933        b.add_point(4, &[2.0, 0.5]);
934        b.add_point(5, &[1.6, 0.9]);
935        b.add_point(6, &[1.0, 0.5]);
936        b.add_point(7, &[1.6, 0.1]);
937
938        b.add_cell(0, (ReferenceCellType::Quadrilateral, 1, &[0, 1, 2, 3]));
939        b.add_cell(1, (ReferenceCellType::Triangle, 2, &[1, 4, 3, 5, 6, 7]));
940        let g = b.create_grid();
941        g.export_as_gmsh("_test_io_mixed.msh");
942
943        let mut b = MixedGridBuilder::<f64>::new(2);
944        b.import_from_gmsh("_test_io_mixed.msh");
945        let _g = b.create_grid();
946    }
947
948    #[test]
949    fn test_import_mixed() {
950        let mut b = MixedGridBuilder::<f64>::new(2);
951        b.import_from_gmsh("meshes/mixed.msh");
952        let _g = b.create_grid();
953    }
954    #[test]
955    fn test_import_mixed_triangle() {
956        let mut b = MixedGridBuilder::<f64>::new(3);
957        b.import_from_gmsh("meshes/sphere_triangle.msh");
958        let _g = b.create_grid();
959    }
960
961    #[test]
962    fn test_import_triangle() {
963        let mut b = SingleElementGridBuilder::<f64>::new(3, (ReferenceCellType::Triangle, 1));
964        b.import_from_gmsh("meshes/sphere_triangle.msh");
965        let _g = b.create_grid();
966    }
967
968    #[test]
969    fn test_import_quadrilateral() {
970        let mut b = SingleElementGridBuilder::<f64>::new(3, (ReferenceCellType::Quadrilateral, 1));
971        b.import_from_gmsh("meshes/cube_quadrilateral.msh");
972        let _g = b.create_grid();
973    }
974
975    #[test]
976    fn test_import_tetrahedron() {
977        let mut b = SingleElementGridBuilder::<f64>::new(3, (ReferenceCellType::Tetrahedron, 1));
978        b.import_from_gmsh("meshes/cube_tetrahedron.msh");
979        let _g = b.create_grid();
980    }
981
982    #[test]
983    #[should_panic]
984    fn test_import_wrong_cell() {
985        let mut b = SingleElementGridBuilder::<f64>::new(3, (ReferenceCellType::Quadrilateral, 1));
986        b.import_from_gmsh("meshes/cube_tetrahedron.msh");
987        let _g = b.create_grid();
988    }
989
990    #[test]
991    #[should_panic]
992    fn test_import_wrong_degree() {
993        let mut b = SingleElementGridBuilder::<f64>::new(3, (ReferenceCellType::Triangle, 2));
994        b.import_from_gmsh("meshes/sphere_triangle.msh");
995        let _g = b.create_grid();
996    }
997
998    #[test]
999    fn test_import_triangle_bin() {
1000        let mut b = SingleElementGridBuilder::<f64>::new(3, (ReferenceCellType::Triangle, 1));
1001        b.import_from_gmsh("meshes/sphere_triangle_bin.msh");
1002        let _g = b.create_grid();
1003    }
1004
1005    #[test]
1006    fn test_import_quadrilateral_bin() {
1007        let mut b = SingleElementGridBuilder::<f64>::new(3, (ReferenceCellType::Quadrilateral, 1));
1008        b.import_from_gmsh("meshes/cube_quadrilateral_bin.msh");
1009        let _g = b.create_grid();
1010    }
1011
1012    #[test]
1013    fn test_import_tetrahedron_bin() {
1014        let mut b = SingleElementGridBuilder::<f64>::new(3, (ReferenceCellType::Tetrahedron, 1));
1015        b.import_from_gmsh("meshes/cube_tetrahedron_bin.msh");
1016        let _g = b.create_grid();
1017    }
1018
1019    #[test]
1020    #[should_panic]
1021    fn test_import_wrong_cell_bin() {
1022        let mut b = SingleElementGridBuilder::<f64>::new(3, (ReferenceCellType::Quadrilateral, 1));
1023        b.import_from_gmsh("meshes/cube_tetrahedron_bin.msh");
1024        let _g = b.create_grid();
1025    }
1026
1027    #[test]
1028    #[should_panic]
1029    fn test_import_wrong_degree_bin() {
1030        let mut b = SingleElementGridBuilder::<f64>::new(3, (ReferenceCellType::Triangle, 2));
1031        b.import_from_gmsh("meshes/sphere_triangle_bin.msh");
1032        let _g = b.create_grid();
1033    }
1034}
ndgrid/io/gmsh.rs

ndgrid/io/
gmsh.rs
