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            self.add_point(tag, &coords[..self.gdim()]);
277        }
278
279        // Load elements
280        let elements = gmsh_section(&s, "Elements");
281        let mut elements = elements.lines();
282
283        let Some(num_elements) = elements.next() else {
284            panic!("Could not read num nodes");
285        };
286        let num_elements = num_elements
287            .parse::<usize>()
288            .expect("Could not parse num nodes");
289
290        for _ in 0..num_elements {
291            let Some(line) = elements.next() else {
292                panic!("Unable to read element line");
293            };
294            let line = line.split(" ").collect::<Vec<&str>>();
295
296            let (a, node_tags) = line.split_at(5);
297            let [tag, element_type, ..] = a
298                .iter()
299                .map(|i| {
300                    i.parse::<usize>()
301                        .expect("Could not parse element tag and type")
302                })
303                .collect::<Vec<_>>()[..]
304            else {
305                panic!("Unrecognised format");
306            };
307
308            let node_tags = node_tags
309                .iter()
310                .map(|i| {
311                    i.parse::<usize>()
312                        .expect("Could not parse nodes for element")
313                })
314                .collect::<Vec<_>>();
315            let (cell_type, degree) = interpret_gmsh_cell(element_type);
316            let gmsh_perm = get_permutation_to_gmsh(cell_type, degree);
317
318            let mut cell = vec![0; node_tags.len()];
319            for (i, j) in gmsh_perm.iter().enumerate() {
320                cell[*j] = node_tags[i];
321            }
322
323            self.add_cell_from_nodes_and_type(tag, &cell, cell_type, degree);
324        }
325    }
326
327    fn import_from_gmsh_string_v2(&mut self, s: String) {
328        // Load nodes
329        let nodes = gmsh_section(&s, "Nodes");
330        let mut nodes = nodes.lines();
331
332        let Some(num_nodes) = nodes.next() else {
333            panic!("Could not read num nodes");
334        };
335        let num_nodes = num_nodes
336            .parse::<usize>()
337            .expect("Could not parse num nodes");
338
339        for _ in 0..num_nodes {
340            let Some(line) = nodes.next() else {
341                panic!("Unable to read node line");
342            };
343            let line = line.split(" ").collect::<Vec<&str>>();
344
345            let (tag, coords) = line.split_at(1);
346            let tag = tag[0].parse::<usize>().expect("Could not parse node tag");
347            let coords = coords
348                .iter()
349                .map(|c| {
350                    if let Ok(d) = T::from_str(c) {
351                        d
352                    } else {
353                        panic!("Could not parse coords")
354                    }
355                })
356                .collect::<Vec<_>>();
357
358            self.add_point(tag, &coords[..self.gdim()]);
359        }
360
361        // Load elements
362        let elements = gmsh_section(&s, "Elements");
363        let mut elements = elements.lines();
364
365        let Some(num_elements) = elements.next() else {
366            panic!("Could not read num nodes");
367        };
368        let num_elements = num_elements
369            .parse::<usize>()
370            .expect("Could not parse num nodes");
371
372        for _ in 0..num_elements {
373            let Some(line) = elements.next() else {
374                panic!("Unable to read element line");
375            };
376            let line = line.split(" ").collect::<Vec<&str>>();
377
378            let (a, rem_line) = line.split_at(3);
379            let [tag, element_type, num_tags] = a
380                .iter()
381                .map(|i| {
382                    i.parse::<usize>()
383                        .expect("Could not parse element tag and type")
384                })
385                .collect::<Vec<_>>()[..]
386            else {
387                panic!("Unrecognised format");
388            };
389
390            let (_, node_tags) = rem_line.split_at(num_tags);
391            let node_tags = node_tags
392                .iter()
393                .map(|i| {
394                    i.parse::<usize>()
395                        .expect("Could not parse nodes for element")
396                })
397                .collect::<Vec<_>>();
398            let (cell_type, degree) = interpret_gmsh_cell(element_type);
399            let gmsh_perm = get_permutation_to_gmsh(cell_type, degree);
400
401            let mut cell = vec![0; node_tags.len()];
402            for (i, j) in gmsh_perm.iter().enumerate() {
403                cell[*j] = node_tags[i];
404            }
405
406            self.add_cell_from_nodes_and_type(tag, &cell, cell_type, degree);
407        }
408    }
409
410    fn import_from_gmsh_binary_v2(
411        &mut self,
412        mut reader: BufReader<File>,
413        data_size: usize,
414        is_le: bool,
415    ) {
416        let mut line = String::new();
417        let mut buf = Vec::new();
418
419        macro_rules! read_exact {
420            ($size: expr, $msg: expr) => {{
421                buf.resize($size, 0);
422                reader.read_exact(&mut buf).expect($msg);
423            }};
424        }
425
426        const GMSH_INT_SIZE: usize = 4;
427
428        loop {
429            let Ok(num_bytes) = reader.read_line(&mut line) else {
430                continue;
431            };
432
433            // EOF reached.
434            if num_bytes == 0 {
435                break;
436            }
437
438            match line.trim() {
439                // Load all nodes.
440                "$Nodes" => {
441                    line.clear();
442                    let Ok(_) = reader.read_line(&mut line) else {
443                        panic!("Unable to read num nodes");
444                    };
445                    let num_nodes = line
446                        .trim()
447                        .parse::<usize>()
448                        .expect("Could not parse num nodes");
449
450                    for _ in 0..num_nodes {
451                        read_exact!(GMSH_INT_SIZE, "Unable to read node tag");
452                        let tag = parse::<usize>(&buf, GMSH_INT_SIZE, is_le);
453
454                        read_exact!(3 * data_size, "Unable to read node coords");
455                        let coords = buf
456                            .chunks(data_size)
457                            .map(|i| parse::<T>(i, data_size, is_le))
458                            .collect::<Vec<_>>();
459
460                        self.add_point(tag, &coords[..self.gdim()]);
461                    }
462
463                    line.clear();
464                }
465                // Load all elements.
466                "$Elements" => {
467                    line.clear();
468                    let Ok(_) = reader.read_line(&mut line) else {
469                        panic!("Unable to read num elements")
470                    };
471                    let num_elements = line
472                        .trim()
473                        .parse::<usize>()
474                        .expect("Could not parse num elements");
475
476                    for _ in 0..num_elements {
477                        read_exact!(3 * GMSH_INT_SIZE, "Unable to element tag and type");
478                        let [elm_type, _num_elm_follow, num_tags] = buf
479                            .chunks(GMSH_INT_SIZE)
480                            .map(|i| parse::<usize>(i, GMSH_INT_SIZE, is_le))
481                            .collect::<Vec<_>>()[..]
482                        else {
483                            panic!("Could not parse element tag and type")
484                        };
485
486                        read_exact!(GMSH_INT_SIZE, "Unable to read num tags");
487                        let tag = parse::<usize>(&buf, 4, is_le);
488
489                        // Skip tags
490                        read_exact!(num_tags * GMSH_INT_SIZE, "Unable to read element tags");
491
492                        let (cell_type, degree) = interpret_gmsh_cell(elm_type);
493                        let gmsh_perm = get_permutation_to_gmsh(cell_type, degree);
494
495                        read_exact!(
496                            gmsh_perm.len() * GMSH_INT_SIZE,
497                            "Unable to read element node tags"
498                        );
499                        let node_tags = buf
500                            .chunks(GMSH_INT_SIZE)
501                            .map(|i| parse::<usize>(i, GMSH_INT_SIZE, is_le))
502                            .collect::<Vec<_>>();
503
504                        let mut cell = vec![0; node_tags.len()];
505                        for (i, j) in gmsh_perm.iter().enumerate() {
506                            cell[*j] = node_tags[i];
507                        }
508
509                        self.add_cell_from_nodes_and_type(tag, &cell, cell_type, degree);
510                    }
511                    line.clear();
512                }
513                _ => {
514                    line.clear();
515                }
516            }
517        }
518    }
519
520    fn import_from_gmsh_string_v4(&mut self, s: String) {
521        // Load nodes
522        let nodes = gmsh_section(&s, "Nodes");
523        let nodes = nodes.lines().collect::<Vec<_>>();
524
525        let [num_entity_blocks, _num_nodes, _min_node_tag, _max_node_tag] = nodes[0]
526            .trim()
527            .split(" ")
528            .map(|i| i.parse::<usize>().unwrap())
529            .collect::<Vec<_>>()[..]
530        else {
531            panic!("Unrecognised gmsh format for node blocks");
532        };
533
534        let mut line_n = 1;
535        for _ in 0..num_entity_blocks {
536            let [_entity_dim, _entity_tag, parametric, num_nodes_in_block] = nodes[line_n]
537                .trim()
538                .split(" ")
539                .map(|i| i.parse::<usize>().unwrap())
540                .collect::<Vec<_>>()[..]
541            else {
542                panic!("Unrecognised gmsh format for nodes");
543            };
544            if parametric == 1 {
545                unimplemented!("Parametric nodes currently not supported")
546            }
547            line_n += 1;
548            let tags = &nodes[line_n..line_n + num_nodes_in_block];
549            let coords = &nodes[line_n + num_nodes_in_block..line_n + 2 * num_nodes_in_block];
550            for (t, c) in izip!(tags, coords) {
551                let pt = c
552                    .trim()
553                    .split(" ")
554                    .map(|i| {
555                        if let Ok(j) = T::from_str(i) {
556                            j
557                        } else {
558                            panic!("Could not parse coordinate");
559                        }
560                    })
561                    .collect::<Vec<_>>();
562                self.add_point(t.parse::<usize>().unwrap(), &pt[..self.gdim()]);
563            }
564            line_n += num_nodes_in_block * 2;
565        }
566
567        // Load elements
568        let elements = gmsh_section(&s, "Elements");
569        let elements = elements.lines().collect::<Vec<_>>();
570
571        let [
572            num_entity_blocks,
573            _num_elements,
574            _min_element_tag,
575            _max_element_tag,
576        ] = elements[0]
577            .trim()
578            .split(" ")
579            .map(|i| i.parse::<usize>().unwrap())
580            .collect::<Vec<_>>()[..]
581        else {
582            panic!("Unrecognised gmsh format");
583        };
584
585        let mut line_n = 1;
586        for _ in 0..num_entity_blocks {
587            let [
588                _entity_dim,
589                _entity_tag,
590                element_type,
591                num_elements_in_block,
592            ] = elements[line_n]
593                .trim()
594                .split(" ")
595                .map(|i| i.parse::<usize>().unwrap())
596                .collect::<Vec<_>>()[..]
597            else {
598                panic!("Unrecognised gmsh format");
599            };
600            let (cell_type, degree) = interpret_gmsh_cell(element_type);
601            let gmsh_perm = get_permutation_to_gmsh(cell_type, degree);
602
603            line_n += 1;
604            for line in &elements[line_n..line_n + num_elements_in_block] {
605                let line = line
606                    .trim()
607                    .split(" ")
608                    .map(|i| i.parse::<usize>().unwrap())
609                    .collect::<Vec<_>>();
610                let mut cell = vec![0; line.len() - 1];
611                for (i, j) in gmsh_perm.iter().enumerate() {
612                    cell[*j] = line[i + 1];
613                }
614                self.add_cell_from_nodes_and_type(line[0], &cell, cell_type, degree);
615            }
616
617            line_n += num_elements_in_block;
618        }
619    }
620
621    fn import_from_gmsh_binary_v4(
622        &mut self,
623        mut reader: BufReader<File>,
624        data_size: usize,
625        is_le: bool,
626    ) {
627        let mut line = String::new();
628        let mut buf = Vec::new();
629
630        macro_rules! read_exact {
631            ($size: expr, $msg: expr) => {{
632                buf.resize($size, 0);
633                reader.read_exact(&mut buf).expect($msg);
634            }};
635        }
636
637        const GMSH_INT_SIZE: usize = 4;
638
639        loop {
640            let Ok(num_bytes) = reader.read_line(&mut line) else {
641                continue;
642            };
643
644            // EOF reached.
645            if num_bytes == 0 {
646                break;
647            }
648
649            match line.trim() {
650                // Load all nodes.
651                "$Nodes" => {
652                    read_exact!(4 * data_size, "Unable to read node section info");
653                    let [num_entity_blocks, _num_nodes, _min_node_tag, _max_node_tag] = buf
654                        .chunks(data_size)
655                        .map(|i| parse::<usize>(i, data_size, is_le))
656                        .collect::<Vec<_>>()[..]
657                    else {
658                        panic!("Could not parse node section info")
659                    };
660
661                    for _ in 0..num_entity_blocks {
662                        read_exact!(3 * GMSH_INT_SIZE, "Unable to read node entity block info");
663                        let [_entity_dim, _entity_tag, parametric] = buf
664                            .chunks(GMSH_INT_SIZE)
665                            .map(|i| parse::<usize>(i, GMSH_INT_SIZE, is_le))
666                            .collect::<Vec<_>>()[..]
667                        else {
668                            panic!("Could not parse node entity block info")
669                        };
670
671                        if parametric == 1 {
672                            unimplemented!("Parametric nodes currently not supported")
673                        }
674
675                        read_exact!(data_size, "Unable to read num nodes in block");
676                        let num_nodes_in_block = parse::<usize>(&buf, data_size, is_le);
677
678                        read_exact!(num_nodes_in_block * data_size, "Unable to read node tags");
679                        let tags = buf
680                            .chunks(data_size)
681                            .map(|i| parse::<usize>(i, data_size, is_le))
682                            .collect::<Vec<_>>();
683
684                        read_exact!(
685                            3 * num_nodes_in_block * data_size,
686                            "Unable to read node coords"
687                        );
688                        let coords = buf
689                            .chunks(3 * data_size)
690                            .map(|i| {
691                                i.chunks(data_size)
692                                    .map(|j| parse::<T>(j, data_size, is_le))
693                                    .collect::<Vec<_>>()
694                            })
695                            .collect::<Vec<_>>();
696
697                        for (t, c) in izip!(tags, coords) {
698                            self.add_point(t, &c[..self.gdim()]);
699                        }
700                    }
701
702                    line.clear();
703                }
704                // Load all elements.
705                "$Elements" => {
706                    read_exact!(4 * data_size, "Unable to read element section info");
707                    let [
708                        num_entity_blocks,
709                        _num_elements,
710                        _min_element_tag,
711                        _max_element_tag,
712                    ] = buf
713                        .chunks(data_size)
714                        .map(|i| parse::<usize>(i, data_size, is_le))
715                        .collect::<Vec<_>>()[..]
716                    else {
717                        panic!("Could not parse element section info")
718                    };
719
720                    for _ in 0..num_entity_blocks {
721                        read_exact!(
722                            3 * GMSH_INT_SIZE,
723                            "Unable to read element entity block info"
724                        );
725                        let [_entity_dim, _entity_tag, element_type] = buf
726                            .chunks(GMSH_INT_SIZE)
727                            .map(|i| parse::<usize>(i, GMSH_INT_SIZE, is_le))
728                            .collect::<Vec<_>>()[..]
729                        else {
730                            panic!("Could not parse element entity block info")
731                        };
732
733                        read_exact!(data_size, "Unable to read num elements in block");
734                        let num_elements_in_block = parse::<usize>(&buf, data_size, is_le);
735
736                        let (cell_type, degree) = interpret_gmsh_cell(element_type);
737                        let gmsh_perm = get_permutation_to_gmsh(cell_type, degree);
738
739                        for _ in 0..num_elements_in_block {
740                            read_exact!(data_size, "Unable to read element tag");
741                            let tag = parse::<usize>(&buf, data_size, is_le);
742
743                            read_exact!(data_size * gmsh_perm.len(), "Unable to read node tags");
744                            let node_tags = buf
745                                .chunks(data_size)
746                                .map(|i| parse::<usize>(i, data_size, is_le))
747                                .collect::<Vec<_>>();
748
749                            let mut cell = vec![0; node_tags.len()];
750                            for (i, j) in gmsh_perm.iter().enumerate() {
751                                cell[*j] = node_tags[i];
752                            }
753
754                            self.add_cell_from_nodes_and_type(tag, &cell, cell_type, degree);
755                        }
756                    }
757                    line.clear();
758                }
759                _ => {
760                    line.clear();
761                }
762            }
763        }
764    }
765}
766
767#[cfg(test)]
768mod test {
769    use super::*;
770    use crate::{
771        MixedGridBuilder, SingleElementGridBuilder,
772        shapes::regular_sphere,
773        traits::{Builder, Topology},
774    };
775    use approx::*;
776
777    #[test]
778    fn test_regular_sphere_gmsh_io() {
779        let g = regular_sphere::<f64>(2);
780        g.export_as_gmsh("_test_io_sphere.msh");
781    }
782
783    #[test]
784    fn test_export_quads() {
785        let mut b = SingleElementGridBuilder::<f64>::new(3, (ReferenceCellType::Quadrilateral, 1));
786        b.add_point(0, &[0.0, 0.0, 0.0]);
787        b.add_point(1, &[1.0, 0.0, 0.0]);
788        b.add_point(2, &[0.0, 1.0, 0.0]);
789        b.add_point(3, &[1.0, 1.0, 0.0]);
790        b.add_point(4, &[0.0, 0.0, 1.0]);
791        b.add_point(5, &[1.0, 0.0, 1.0]);
792        b.add_point(6, &[0.0, 1.0, 1.0]);
793        b.add_point(7, &[1.0, 1.0, 1.0]);
794        b.add_cell(0, &[0, 2, 1, 3]);
795        b.add_cell(1, &[0, 1, 4, 5]);
796        b.add_cell(2, &[0, 4, 2, 6]);
797        b.add_cell(3, &[1, 3, 5, 7]);
798        b.add_cell(4, &[2, 6, 3, 7]);
799        b.add_cell(5, &[4, 5, 6, 7]);
800        let g = b.create_grid();
801        g.export_as_gmsh("_test_io_cube.msh");
802    }
803
804    #[test]
805    fn test_export_tetrahedra() {
806        let mut b = SingleElementGridBuilder::<f64>::new(3, (ReferenceCellType::Tetrahedron, 1));
807        b.add_point(0, &[0.0, 0.0, 0.0]);
808        b.add_point(1, &[1.0, 0.0, 0.0]);
809        b.add_point(2, &[0.0, 1.0, 0.0]);
810        b.add_point(3, &[1.0, 1.0, 0.0]);
811        b.add_point(4, &[0.0, 0.0, 1.0]);
812        b.add_point(5, &[1.0, 0.0, 1.0]);
813        b.add_point(6, &[0.0, 1.0, 1.0]);
814        b.add_point(7, &[1.0, 1.0, 1.0]);
815        b.add_cell(0, &[0, 1, 5, 7]);
816        b.add_cell(1, &[0, 2, 6, 7]);
817        b.add_cell(2, &[0, 4, 5, 7]);
818        b.add_cell(3, &[0, 1, 3, 7]);
819        b.add_cell(4, &[0, 2, 3, 7]);
820        b.add_cell(5, &[0, 4, 6, 7]);
821        let g = b.create_grid();
822        g.export_as_gmsh("_test_io_tetrahedra.msh");
823    }
824
825    #[test]
826    fn test_hexahedra() {
827        let mut b = SingleElementGridBuilder::<f64>::new(3, (ReferenceCellType::Hexahedron, 1));
828        b.add_point(0, &[0.0, 0.0, 0.0]);
829        b.add_point(1, &[1.0, 0.0, 0.0]);
830        b.add_point(2, &[0.0, 1.0, 0.0]);
831        b.add_point(3, &[1.0, 1.0, 0.0]);
832        b.add_point(4, &[0.0, 0.0, 1.0]);
833        b.add_point(5, &[1.0, 0.0, 1.0]);
834        b.add_point(6, &[0.0, 1.0, 1.0]);
835        b.add_point(7, &[1.0, 1.0, 1.0]);
836        b.add_point(8, &[0.0, 0.0, 2.0]);
837        b.add_point(9, &[1.0, 0.0, 2.0]);
838        b.add_point(10, &[0.0, 1.0, 2.0]);
839        b.add_point(11, &[1.0, 1.0, 1.5]);
840        b.add_cell(1, &[0, 1, 2, 3, 4, 5, 6, 7]);
841        b.add_cell(2, &[4, 5, 6, 7, 8, 9, 10, 11]);
842        let g = b.create_grid();
843        g.export_as_gmsh("_test_io_hexahedra.msh");
844
845        let mut b = SingleElementGridBuilder::<f64>::new(3, (ReferenceCellType::Hexahedron, 1));
846        b.import_from_gmsh("_test_io_hexahedra.msh");
847        let g2 = b.create_grid();
848
849        let mut p1 = [0.0; 3];
850        let mut p2 = [0.0; 3];
851        for (v1, v2) in izip!(
852            g.entity_iter(ReferenceCellType::Point),
853            g2.entity_iter(ReferenceCellType::Point)
854        ) {
855            for (pt1, pt2) in izip!(v1.geometry().points(), v2.geometry().points()) {
856                pt1.coords(&mut p1);
857                pt2.coords(&mut p2);
858                for (c1, c2) in izip!(&p1, &p2) {
859                    assert_relative_eq!(c1, c2, epsilon = 1e-10);
860                }
861            }
862        }
863        for (h1, h2) in izip!(
864            g.entity_iter(ReferenceCellType::Hexahedron),
865            g2.entity_iter(ReferenceCellType::Hexahedron)
866        ) {
867            for (v1, v2) in izip!(
868                h1.topology().sub_entity_iter(ReferenceCellType::Point),
869                h2.topology().sub_entity_iter(ReferenceCellType::Point)
870            ) {
871                assert_eq!(v1, v2);
872            }
873        }
874    }
875
876    #[test]
877    fn test_high_order_triangles() {
878        let mut b = SingleElementGridBuilder::<f64>::new(3, (ReferenceCellType::Triangle, 5));
879        b.add_point(0, &[0.0, 0.0, 0.0]);
880        b.add_point(1, &[5.0, 0.0, 0.2]);
881        b.add_point(2, &[0.0, 5.0, 0.4]);
882        b.add_point(3, &[4.0, 1.0, -0.1]);
883        b.add_point(4, &[3.0, 2.0, 0.2]);
884        b.add_point(5, &[2.0, 3.0, -0.3]);
885        b.add_point(6, &[1.0, 4.0, -0.5]);
886        b.add_point(7, &[0.0, 1.0, 0.6]);
887        b.add_point(8, &[0.0, 2.0, 0.2]);
888        b.add_point(9, &[0.0, 3.0, 0.1]);
889        b.add_point(10, &[0.0, 4.0, -0.2]);
890        b.add_point(11, &[1.0, 0.0, -0.3]);
891        b.add_point(12, &[2.0, 0.0, -0.4]);
892        b.add_point(13, &[3.0, 0.0, -0.5]);
893        b.add_point(14, &[4.0, 0.0, -0.2]);
894        b.add_point(15, &[1.0, 1.0, 0.1]);
895        b.add_point(16, &[2.0, 1.0, 0.1]);
896        b.add_point(17, &[3.0, 1.0, 0.1]);
897        b.add_point(18, &[2.0, 1.0, 0.2]);
898        b.add_point(19, &[2.0, 2.0, 0.1]);
899        b.add_point(20, &[3.0, 1.0, 0.1]);
900        b.add_cell(
901            1,
902            &[
903                0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
904            ],
905        );
906        let g = b.create_grid();
907        g.export_as_gmsh("_test_io_high_order_triangle.msh");
908
909        let mut b = SingleElementGridBuilder::<f64>::new(3, (ReferenceCellType::Triangle, 5));
910        b.import_from_gmsh("_test_io_high_order_triangle.msh");
911        let g2 = b.create_grid();
912
913        let mut p1 = [0.0; 3];
914        let mut p2 = [0.0; 3];
915        for (v1, v2) in izip!(
916            g.entity_iter(ReferenceCellType::Point),
917            g2.entity_iter(ReferenceCellType::Point)
918        ) {
919            v1.geometry().points().next().unwrap().coords(&mut p1);
920            v2.geometry().points().next().unwrap().coords(&mut p2);
921        }
922        for (v1, v2) in izip!(
923            g.entity_iter(ReferenceCellType::Point),
924            g2.entity_iter(ReferenceCellType::Point)
925        ) {
926            v1.geometry().points().next().unwrap().coords(&mut p1);
927            v2.geometry().points().next().unwrap().coords(&mut p2);
928            for (c1, c2) in izip!(&p1, &p2) {
929                assert_relative_eq!(c1, c2, epsilon = 1e-10);
930            }
931        }
932        for (h1, h2) in izip!(
933            g.entity_iter(ReferenceCellType::Triangle),
934            g2.entity_iter(ReferenceCellType::Triangle)
935        ) {
936            for (v1, v2) in izip!(
937                h1.topology().sub_entity_iter(ReferenceCellType::Point),
938                h2.topology().sub_entity_iter(ReferenceCellType::Point)
939            ) {
940                assert_eq!(v1, v2);
941            }
942        }
943    }
944
945    #[test]
946    fn test_export_mixed() {
947        let mut b = MixedGridBuilder::<f64>::new(2);
948        b.add_point(0, &[0.0, 0.0]);
949        b.add_point(1, &[1.0, 0.0]);
950        b.add_point(2, &[0.0, 1.0]);
951        b.add_point(3, &[1.0, 1.0]);
952        b.add_point(4, &[2.0, 0.5]);
953        b.add_point(5, &[1.6, 0.9]);
954        b.add_point(6, &[1.0, 0.5]);
955        b.add_point(7, &[1.6, 0.1]);
956
957        b.add_cell(0, (ReferenceCellType::Quadrilateral, 1, &[0, 1, 2, 3]));
958        b.add_cell(1, (ReferenceCellType::Triangle, 2, &[1, 4, 3, 5, 6, 7]));
959        let g = b.create_grid();
960        g.export_as_gmsh("_test_io_mixed.msh");
961
962        let mut b = MixedGridBuilder::<f64>::new(2);
963        b.import_from_gmsh("_test_io_mixed.msh");
964        let _g = b.create_grid();
965    }
966
967    #[test]
968    fn test_import_mixed() {
969        let mut b = MixedGridBuilder::<f64>::new(2);
970        b.import_from_gmsh("meshes/mixed.msh");
971        let _g = b.create_grid();
972    }
973    #[test]
974    fn test_import_mixed_triangle() {
975        let mut b = MixedGridBuilder::<f64>::new(3);
976        b.import_from_gmsh("meshes/sphere_triangle.msh");
977        let _g = b.create_grid();
978    }
979
980    #[test]
981    fn test_import_triangle() {
982        let mut b = SingleElementGridBuilder::<f64>::new(3, (ReferenceCellType::Triangle, 1));
983        b.import_from_gmsh("meshes/sphere_triangle.msh");
984        let _g = b.create_grid();
985    }
986
987    #[test]
988    fn test_import_quadrilateral() {
989        let mut b = SingleElementGridBuilder::<f64>::new(3, (ReferenceCellType::Quadrilateral, 1));
990        b.import_from_gmsh("meshes/cube_quadrilateral.msh");
991        let _g = b.create_grid();
992    }
993
994    #[test]
995    fn test_import_tetrahedron() {
996        let mut b = SingleElementGridBuilder::<f64>::new(3, (ReferenceCellType::Tetrahedron, 1));
997        b.import_from_gmsh("meshes/cube_tetrahedron.msh");
998        let _g = b.create_grid();
999    }
1000
1001    #[test]
1002    #[should_panic]
1003    fn test_import_wrong_cell() {
1004        let mut b = SingleElementGridBuilder::<f64>::new(3, (ReferenceCellType::Quadrilateral, 1));
1005        b.import_from_gmsh("meshes/cube_tetrahedron.msh");
1006        let _g = b.create_grid();
1007    }
1008
1009    #[test]
1010    #[should_panic]
1011    fn test_import_wrong_degree() {
1012        let mut b = SingleElementGridBuilder::<f64>::new(3, (ReferenceCellType::Triangle, 2));
1013        b.import_from_gmsh("meshes/sphere_triangle.msh");
1014        let _g = b.create_grid();
1015    }
1016
1017    #[test]
1018    fn test_import_triangle_bin() {
1019        let mut b = SingleElementGridBuilder::<f64>::new(3, (ReferenceCellType::Triangle, 1));
1020        b.import_from_gmsh("meshes/sphere_triangle_bin.msh");
1021        let _g = b.create_grid();
1022    }
1023
1024    #[test]
1025    fn test_import_quadrilateral_bin() {
1026        let mut b = SingleElementGridBuilder::<f64>::new(3, (ReferenceCellType::Quadrilateral, 1));
1027        b.import_from_gmsh("meshes/cube_quadrilateral_bin.msh");
1028        let _g = b.create_grid();
1029    }
1030
1031    #[test]
1032    fn test_import_tetrahedron_bin() {
1033        let mut b = SingleElementGridBuilder::<f64>::new(3, (ReferenceCellType::Tetrahedron, 1));
1034        b.import_from_gmsh("meshes/cube_tetrahedron_bin.msh");
1035        let _g = b.create_grid();
1036    }
1037
1038    #[test]
1039    #[should_panic]
1040    fn test_import_wrong_cell_bin() {
1041        let mut b = SingleElementGridBuilder::<f64>::new(3, (ReferenceCellType::Quadrilateral, 1));
1042        b.import_from_gmsh("meshes/cube_tetrahedron_bin.msh");
1043        let _g = b.create_grid();
1044    }
1045
1046    #[test]
1047    #[should_panic]
1048    fn test_import_wrong_degree_bin() {
1049        let mut b = SingleElementGridBuilder::<f64>::new(3, (ReferenceCellType::Triangle, 2));
1050        b.import_from_gmsh("meshes/sphere_triangle_bin.msh");
1051        let _g = b.create_grid();
1052    }
1053}
ndgrid/io/gmsh.rs

ndgrid/io/
gmsh.rs
