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

ndgrid/io/
gmsh.rs
