From 15dc389d1460a77e55dedba8730e689942baed30 Mon Sep 17 00:00:00 2001
From: Matthew Scroggs <matthew.w.scroggs@gmail.com>
Date: Mon, 14 Jul 2025 16:57:43 +0100
Subject: [PATCH 1/7] implement import from gmsh

---
 src/grid/local_grid/single_element/builder.rs |  7 ++
 src/io/gmsh.rs                                | 98 ++++++++++++++++++-
 src/traits.rs                                 |  2 +-
 src/traits/builder.rs                         |  3 +
 src/traits/io.rs                              |  2 +-
 src/traits/io/gmsh.rs                         | 16 ++-
 src/traits/io/ron.rs                          |  4 +-
 7 files changed, 125 insertions(+), 7 deletions(-)

diff --git a/src/grid/local_grid/single_element/builder.rs b/src/grid/local_grid/single_element/builder.rs
index ccdc203..074bd6a 100644
--- a/src/grid/local_grid/single_element/builder.rs
+++ b/src/grid/local_grid/single_element/builder.rs
@@ -100,6 +100,13 @@ impl<T: RealScalar> Builder for SingleElementGridBuilder<T> {
         }
     }
 
+    fn add_cell_from_nodes_and_type(&mut self, id: usize, nodes: &[usize], cell_type: ReferenceCellType, cell_degree: usize) {
+        if (cell_type, cell_degree) != self.element_data {
+            panic!("Invalid cell type.");
+        }
+        self.add_cell(id, nodes);
+    }
+
     fn create_grid(&self) -> SingleElementGrid<T, CiarletElement<T, IdentityMap>> {
         let cell_vertices =
             self.extract_vertices(&self.cells, &[self.element_data.0], &[self.element_data.1]);
diff --git a/src/io/gmsh.rs b/src/io/gmsh.rs
index 4fbe29e..a088461 100644
--- a/src/io/gmsh.rs
+++ b/src/io/gmsh.rs
@@ -1,9 +1,10 @@
 //! Gmsh I/O
 
-use crate::traits::{Entity, Geometry, GmshExport, Grid, Point, Topology};
+use crate::traits::{Builder, Entity, Geometry, GmshExport, Grid, Point, Topology, GmshImport};
 use itertools::izip;
 use ndelement::types::ReferenceCellType;
 use num::Zero;
+use std::str::FromStr;
 
 fn get_permutation_to_gmsh(cell_type: ReferenceCellType, degree: usize) -> Vec<usize> {
     match cell_type {
@@ -81,6 +82,23 @@ fn get_gmsh_cell(cell_type: ReferenceCellType, degree: usize) -> usize {
     }
 }
 
+fn interpret_gmsh_cell(gmsh_cell: usize) -> (ReferenceCellType, usize) {
+    match gmsh_cell {
+        2 => (ReferenceCellType::Triangle, 1),
+        9 => (ReferenceCellType::Triangle, 2),
+        21 => (ReferenceCellType::Triangle, 3),
+        23 => (ReferenceCellType::Triangle, 4),
+        25 => (ReferenceCellType::Triangle, 5),
+        3 => (ReferenceCellType::Quadrilateral, 1),
+        10 => (ReferenceCellType::Quadrilateral, 2),
+        4 => (ReferenceCellType::Tetrahedron, 1),
+        5 => (ReferenceCellType::Hexahedron, 1),
+        _ => {
+            panic!("Unsupported cell type.");
+        }
+    }
+}
+
 impl<G: Grid<EntityDescriptor = ReferenceCellType>> GmshExport for G {
     fn to_gmsh_string(&self) -> String {
         let tdim = self.topology_dim();
@@ -156,10 +174,79 @@ impl<G: Grid<EntityDescriptor = ReferenceCellType>> GmshExport for G {
     }
 }
 
+/// Get a section from a gmsh string
+fn gmsh_section(s: &str, section: &str) -> String {
+    let a = s.split(&format!("${section}\n")).collect::<Vec<_>>();
+    if a.len() == 0 {
+        panic!("Section not found: {section}");
+    }
+    String::from(a[1].split(&format!("\n$End{section}")).collect::<Vec<_>>()[0])
+}
+
+impl<T: FromStr, B: Builder<T=T, EntityDescriptor = ReferenceCellType>> GmshImport for B {
+    fn import_from_gmsh_string(&mut self, s: String) {
+        let format = gmsh_section(&s, "MeshFormat");
+        // Check msh file version
+        let [version, ascii_mode, _binary_mode] = format.split(" ").collect::<Vec<_>>()[..] else { panic!("Unrecognised gmsh format"); };
+        if version != "4.1" {
+            unimplemented!("Unsupported gmsh file version");
+        }
+        if ascii_mode != "0" {
+            unimplemented!("Non-ASCII gmsh files currently not supported");
+        }
+        // Load nodes
+        let nodes = gmsh_section(&s, "Nodes");
+        let nodes = nodes.lines().collect::<Vec<_>>();
+
+        let [num_entity_blocks, num_nodes, _min_node_tag, _max_node_tag] = nodes[0].split(" ").map(|i| i.parse::<usize>().unwrap()).collect::<Vec<_>>()[..] else { panic!("Unrecognised gmsh format"); };
+
+        let mut line_n = 1;
+        for _ in 0..num_entity_blocks {
+            let [
+                _entity_dim, _entity_tag, parametric, num_nodes_in_block
+            ] = nodes[line_n].split(" ").map(|i| i.parse::<usize>().unwrap()).collect::<Vec<_>>()[..] else { panic!("Unrecognised gmsh format"); };
+            if parametric == 1 {
+                unimplemented!("Parametric nodes currently not supported")
+            }
+            line_n += 1;
+            let tags = &nodes[line_n..line_n + num_nodes_in_block];
+            let coords = &nodes[line_n + num_nodes_in_block..line_n + 2 * num_nodes_in_block];
+            for (t, c) in izip!(tags, coords) {
+                self.add_point(t.parse::<usize>().unwrap(), &c.split(" ").map(|i| if let Ok(j) = T::from_str(i) {j} else {panic!("Could not parse coordinate");}).collect::<Vec<_>>());
+            }
+            line_n += num_nodes_in_block + 2;
+        }
+
+
+        // Load elements
+        let elements = gmsh_section(&s, "Elements");
+        let elements = elements.lines().collect::<Vec<_>>();
+
+        let [num_entity_blocks, num_elements, _min_element_tag, _max_element_tag] = nodes[0].split(" ").map(|i| i.parse::<usize>().unwrap()).collect::<Vec<_>>()[..] else { panic!("Unrecognised gmsh format"); };
+        
+        let mut line_n = 1;
+        for _ in 0..num_entity_blocks {
+            let [
+                _entity_dim, _entity_tag, element_type, num_elements_in_block
+            ] = elements[line_n].split(" ").map(|i| i.parse::<usize>().unwrap()).collect::<Vec<_>>()[..] else { panic!("Unrecognised gmsh format"); };
+            let (cell_type, degree) = interpret_gmsh_cell(element_type);
+            line_n += 1;
+            for line in &elements[line_n..line_n + num_elements_in_block] {
+                let line = line.split(" ").map(|i| i.parse::<usize>().unwrap()).collect::<Vec<_>>();
+                // TODO: permute
+                self.add_cell_from_nodes_and_type(line[0], &line[1..], cell_type, degree);
+            }
+
+            line_n += num_elements_in_block;
+        }
+    }
+}
+
 #[cfg(test)]
 mod test {
     use super::*;
-    use crate::{shapes::regular_sphere, traits::Builder, SingleElementGridBuilder};
+    use crate::{shapes::regular_sphere, traits::Builder, SingleElementGridBuilder, SingleElementGrid};
+    use ndelement::{ciarlet::CiarletElement, map::IdentityMap};
 
     #[test]
     fn test_regular_sphere_gmsh_io() {
@@ -228,4 +315,11 @@ mod test {
         let g = b.create_grid();
         g.export_as_gmsh("_test_io_hexahedra.msh");
     }
+
+    #[test]
+    fn test_gmsh_import_triangle() {
+        let mut b = SingleElementGridBuilder::<f64>::new(3, (ReferenceCellType::Triangle, 1));
+        b.import_from_gmsh("meshes/sphere_triangle.msh");
+        let g = b.create_grid();
+    }
 }
diff --git a/src/traits.rs b/src/traits.rs
index f0911b1..6bfc49b 100644
--- a/src/traits.rs
+++ b/src/traits.rs
@@ -15,7 +15,7 @@ pub use geometry_map::GeometryMap;
 pub use grid::{Grid, ParallelGrid};
 #[cfg(feature = "serde")]
 pub(crate) use io::ConvertToSerializable;
-pub use io::GmshExport;
+pub use io::{GmshExport, GmshImport};
 #[cfg(feature = "serde")]
 pub use io::{RONExport, RONExportParallel, RONImport, RONImportParallel};
 pub use topology::Topology;
diff --git a/src/traits/builder.rs b/src/traits/builder.rs
index 54bca79..1c996f5 100644
--- a/src/traits/builder.rs
+++ b/src/traits/builder.rs
@@ -25,6 +25,9 @@ pub trait Builder {
     /// Add a cell to the grid
     fn add_cell(&mut self, id: usize, cell_data: Self::CellData<'_>);
 
+    /// Add a cell to the grid
+    fn add_cell_from_nodes_and_type(&mut self, id: usize, nodes: &[usize], cell_type: Self::EntityDescriptor, cell_degree: usize);
+
     /// Create the grid
     fn create_grid(&self) -> Self::Grid;
 
diff --git a/src/traits/io.rs b/src/traits/io.rs
index caabefa..e03e635 100644
--- a/src/traits/io.rs
+++ b/src/traits/io.rs
@@ -2,6 +2,6 @@ mod gmsh;
 #[cfg(feature = "serde")]
 mod ron;
 
-pub use gmsh::GmshExport;
+pub use gmsh::{GmshExport, GmshImport};
 #[cfg(feature = "serde")]
 pub use ron::{ConvertToSerializable, RONExport, RONExportParallel, RONImport, RONImportParallel};
diff --git a/src/traits/io/gmsh.rs b/src/traits/io/gmsh.rs
index c8cd835..d677a0b 100644
--- a/src/traits/io/gmsh.rs
+++ b/src/traits/io/gmsh.rs
@@ -1,7 +1,8 @@
 //! Gmsh I/O
 use std::fs;
+use crate::traits::{Grid, Builder};
 
-pub trait GmshExport {
+pub trait GmshExport: Grid {
     //! Grid export for Gmsh
 
     /// Generate the Gmsh string for a grid
@@ -13,3 +14,16 @@ pub trait GmshExport {
         fs::write(filename, gmsh_s).expect("Unable to write file");
     }
 }
+
+pub trait GmshImport: Builder {
+    //! Grid import for Gmsh
+
+    /// Generate grid from a Gmsh string
+    fn import_from_gmsh_string(&mut self, s: String);
+
+    /// Export as Gmsh
+    fn import_from_gmsh(&mut self, filename: &str) {
+        let content = fs::read_to_string(filename).expect("Unable to read file");
+        self.import_from_gmsh_string(content);
+    }
+}
diff --git a/src/traits/io/ron.rs b/src/traits/io/ron.rs
index fcf7466..f3c4d58 100644
--- a/src/traits/io/ron.rs
+++ b/src/traits/io/ron.rs
@@ -28,10 +28,10 @@ pub trait RONExport: Grid {
 pub trait RONImport: Sized + Grid {
     //! Grid import for RON
 
-    /// Generate the RON string for a grid
+    /// Generate grid from a RON string
     fn from_ron_string(s: String) -> Self;
 
-    /// Export as RON
+    /// Import from RON
     fn import_from_ron(filename: &str) -> Self {
         let content = fs::read_to_string(filename).expect("Unable to read file");
         Self::from_ron_string(content)

From d98db3379679cee3d149e6194a58fd5239e55860 Mon Sep 17 00:00:00 2001
From: Matthew Scroggs <matthew.w.scroggs@gmail.com>
Date: Mon, 14 Jul 2025 17:02:12 +0100
Subject: [PATCH 2/7] more tests

---
 meshes/cube_quadrilateral.msh |  33 +++++
 meshes/cube_tetrahedron.msh   |  33 +++++
 meshes/sphere_triangle.msh    | 271 ++++++++++++++++++++++++++++++++++
 src/io/gmsh.rs                |  47 ++++--
 4 files changed, 375 insertions(+), 9 deletions(-)
 create mode 100644 meshes/cube_quadrilateral.msh
 create mode 100644 meshes/cube_tetrahedron.msh
 create mode 100644 meshes/sphere_triangle.msh

diff --git a/meshes/cube_quadrilateral.msh b/meshes/cube_quadrilateral.msh
new file mode 100644
index 0000000..e2911dd
--- /dev/null
+++ b/meshes/cube_quadrilateral.msh
@@ -0,0 +1,33 @@
+$MeshFormat
+4.1 0 8
+$EndMeshFormat
+$Nodes
+1 8 1 8
+2 1 0 8
+1
+2
+3
+4
+5
+6
+7
+8
+0 0 0
+1 0 0
+0 1 0
+1 1 0
+0 0 1
+1 0 1
+0 1 1
+1 1 1
+$EndNodes
+$Elements
+1 6 1 6
+2 1 3 6
+1 1 3 4 2
+2 1 2 6 5
+3 1 5 7 3
+4 2 4 8 6
+5 3 7 8 4
+6 5 6 8 7
+$EndElements
diff --git a/meshes/cube_tetrahedron.msh b/meshes/cube_tetrahedron.msh
new file mode 100644
index 0000000..372e901
--- /dev/null
+++ b/meshes/cube_tetrahedron.msh
@@ -0,0 +1,33 @@
+$MeshFormat
+4.1 0 8
+$EndMeshFormat
+$Nodes
+1 8 1 8
+3 1 0 8
+1
+2
+3
+4
+5
+6
+7
+8
+0 0 0
+1 0 0
+0 1 0
+1 1 0
+0 0 1
+1 0 1
+0 1 1
+1 1 1
+$EndNodes
+$Elements
+1 6 1 6
+3 1 4 6
+1 1 2 6 8
+2 1 3 7 8
+3 1 5 6 8
+4 1 2 4 8
+5 1 3 4 8
+6 1 5 7 8
+$EndElements
diff --git a/meshes/sphere_triangle.msh b/meshes/sphere_triangle.msh
new file mode 100644
index 0000000..94e943b
--- /dev/null
+++ b/meshes/sphere_triangle.msh
@@ -0,0 +1,271 @@
+$MeshFormat
+4.1 0 8
+$EndMeshFormat
+$Nodes
+1 66 1 66
+2 1 0 66
+1
+2
+3
+4
+5
+6
+7
+8
+9
+10
+11
+12
+13
+14
+15
+16
+17
+18
+19
+20
+21
+22
+23
+24
+25
+26
+27
+28
+29
+30
+31
+32
+33
+34
+35
+36
+37
+38
+39
+40
+41
+42
+43
+44
+45
+46
+47
+48
+49
+50
+51
+52
+53
+54
+55
+56
+57
+58
+59
+60
+61
+62
+63
+64
+65
+66
+0 0 1
+1 0 0
+0 1 0
+-1 0 0
+0 -1 0
+0 0 -1
+0.7071067811865475 0.7071067811865475 0
+0 0.7071067811865475 0.7071067811865475
+0.7071067811865475 0 0.7071067811865475
+-0.7071067811865475 0.7071067811865475 0
+-0.7071067811865475 0 0.7071067811865475
+-0.7071067811865475 -0.7071067811865475 0
+0 -0.7071067811865475 0.7071067811865475
+0.7071067811865475 -0.7071067811865475 0
+0.7071067811865475 0 -0.7071067811865475
+0 0.7071067811865475 -0.7071067811865475
+-0.7071067811865475 0 -0.7071067811865475
+0 -0.7071067811865475 -0.7071067811865475
+0.4082482904638631 0.4082482904638631 0.8164965809277261
+0 0.3826834323650897 0.9238795325112867
+0.3826834323650897 0 0.9238795325112867
+0.816496580927726 0.408248290463863 0.408248290463863
+0.9238795325112867 0 0.3826834323650897
+0.9238795325112867 0.3826834323650897 0
+0.408248290463863 0.816496580927726 0.408248290463863
+0.3826834323650897 0.9238795325112867 0
+0 0.9238795325112867 0.3826834323650897
+-0.4082482904638631 0.4082482904638631 0.8164965809277261
+-0.3826834323650897 0 0.9238795325112867
+-0.408248290463863 0.816496580927726 0.408248290463863
+-0.3826834323650897 0.9238795325112867 0
+-0.816496580927726 0.408248290463863 0.408248290463863
+-0.9238795325112867 0.3826834323650897 0
+-0.9238795325112867 0 0.3826834323650897
+-0.4082482904638631 -0.4082482904638631 0.8164965809277261
+0 -0.3826834323650897 0.9238795325112867
+-0.816496580927726 -0.408248290463863 0.408248290463863
+-0.9238795325112867 -0.3826834323650897 0
+-0.408248290463863 -0.816496580927726 0.408248290463863
+-0.3826834323650897 -0.9238795325112867 0
+0 -0.9238795325112867 0.3826834323650897
+0.4082482904638631 -0.4082482904638631 0.8164965809277261
+0.408248290463863 -0.816496580927726 0.408248290463863
+0.3826834323650897 -0.9238795325112867 0
+0.816496580927726 -0.408248290463863 0.408248290463863
+0.9238795325112867 -0.3826834323650897 0
+0.4082482904638631 0.4082482904638631 -0.8164965809277261
+0.3826834323650897 0 -0.9238795325112867
+0 0.3826834323650897 -0.9238795325112867
+0.408248290463863 0.816496580927726 -0.408248290463863
+0 0.9238795325112867 -0.3826834323650897
+0.816496580927726 0.408248290463863 -0.408248290463863
+0.9238795325112867 0 -0.3826834323650897
+-0.4082482904638631 0.4082482904638631 -0.8164965809277261
+-0.3826834323650897 0 -0.9238795325112867
+-0.816496580927726 0.408248290463863 -0.408248290463863
+-0.9238795325112867 0 -0.3826834323650897
+-0.408248290463863 0.816496580927726 -0.408248290463863
+-0.4082482904638631 -0.4082482904638631 -0.8164965809277261
+0 -0.3826834323650897 -0.9238795325112867
+-0.408248290463863 -0.816496580927726 -0.408248290463863
+0 -0.9238795325112867 -0.3826834323650897
+-0.816496580927726 -0.408248290463863 -0.408248290463863
+0.4082482904638631 -0.4082482904638631 -0.8164965809277261
+0.816496580927726 -0.408248290463863 -0.408248290463863
+0.408248290463863 -0.816496580927726 -0.408248290463863
+$EndNodes
+$Elements
+1 128 1 128
+2 1 2 128
+1 1 21 20
+2 9 19 21
+3 8 20 19
+4 19 20 21
+5 2 24 23
+6 7 22 24
+7 9 23 22
+8 22 23 24
+9 3 27 26
+10 8 25 27
+11 7 26 25
+12 25 26 27
+13 7 25 22
+14 8 19 25
+15 9 22 19
+16 19 22 25
+17 1 20 29
+18 8 28 20
+19 11 29 28
+20 28 29 20
+21 3 31 27
+22 10 30 31
+23 8 27 30
+24 30 27 31
+25 4 34 33
+26 11 32 34
+27 10 33 32
+28 32 33 34
+29 10 32 30
+30 11 28 32
+31 8 30 28
+32 28 30 32
+33 1 29 36
+34 11 35 29
+35 13 36 35
+36 35 36 29
+37 4 38 34
+38 12 37 38
+39 11 34 37
+40 37 34 38
+41 5 41 40
+42 13 39 41
+43 12 40 39
+44 39 40 41
+45 12 39 37
+46 13 35 39
+47 11 37 35
+48 35 37 39
+49 1 36 21
+50 13 42 36
+51 9 21 42
+52 42 21 36
+53 5 44 41
+54 14 43 44
+55 13 41 43
+56 43 41 44
+57 2 23 46
+58 9 45 23
+59 14 46 45
+60 45 46 23
+61 14 45 43
+62 9 42 45
+63 13 43 42
+64 42 43 45
+65 6 49 48
+66 16 47 49
+67 15 48 47
+68 47 48 49
+69 3 26 51
+70 7 50 26
+71 16 51 50
+72 50 51 26
+73 2 53 24
+74 15 52 53
+75 7 24 52
+76 52 24 53
+77 7 52 50
+78 15 47 52
+79 16 50 47
+80 47 50 52
+81 6 55 49
+82 17 54 55
+83 16 49 54
+84 54 49 55
+85 4 33 57
+86 10 56 33
+87 17 57 56
+88 56 57 33
+89 3 51 31
+90 16 58 51
+91 10 31 58
+92 58 31 51
+93 10 58 56
+94 16 54 58
+95 17 56 54
+96 54 56 58
+97 6 60 55
+98 18 59 60
+99 17 55 59
+100 59 55 60
+101 5 40 62
+102 12 61 40
+103 18 62 61
+104 61 62 40
+105 4 57 38
+106 17 63 57
+107 12 38 63
+108 63 38 57
+109 12 63 61
+110 17 59 63
+111 18 61 59
+112 59 61 63
+113 6 48 60
+114 15 64 48
+115 18 60 64
+116 64 60 48
+117 2 46 53
+118 14 65 46
+119 15 53 65
+120 65 53 46
+121 5 62 44
+122 18 66 62
+123 14 44 66
+124 66 44 62
+125 14 66 65
+126 18 64 66
+127 15 65 64
+128 64 65 66
+$EndElements
diff --git a/src/io/gmsh.rs b/src/io/gmsh.rs
index a088461..038141b 100644
--- a/src/io/gmsh.rs
+++ b/src/io/gmsh.rs
@@ -198,7 +198,7 @@ impl<T: FromStr, B: Builder<T=T, EntityDescriptor = ReferenceCellType>> GmshImpo
         let nodes = gmsh_section(&s, "Nodes");
         let nodes = nodes.lines().collect::<Vec<_>>();
 
-        let [num_entity_blocks, num_nodes, _min_node_tag, _max_node_tag] = nodes[0].split(" ").map(|i| i.parse::<usize>().unwrap()).collect::<Vec<_>>()[..] else { panic!("Unrecognised gmsh format"); };
+        let [num_entity_blocks, _num_nodes, _min_node_tag, _max_node_tag] = nodes[0].split(" ").map(|i| i.parse::<usize>().unwrap()).collect::<Vec<_>>()[..] else { panic!("Unrecognised gmsh format"); };
 
         let mut line_n = 1;
         for _ in 0..num_entity_blocks {
@@ -222,7 +222,7 @@ impl<T: FromStr, B: Builder<T=T, EntityDescriptor = ReferenceCellType>> GmshImpo
         let elements = gmsh_section(&s, "Elements");
         let elements = elements.lines().collect::<Vec<_>>();
 
-        let [num_entity_blocks, num_elements, _min_element_tag, _max_element_tag] = nodes[0].split(" ").map(|i| i.parse::<usize>().unwrap()).collect::<Vec<_>>()[..] else { panic!("Unrecognised gmsh format"); };
+        let [num_entity_blocks, _num_elements, _min_element_tag, _max_element_tag] = nodes[0].split(" ").map(|i| i.parse::<usize>().unwrap()).collect::<Vec<_>>()[..] else { panic!("Unrecognised gmsh format"); };
         
         let mut line_n = 1;
         for _ in 0..num_entity_blocks {
@@ -245,8 +245,7 @@ impl<T: FromStr, B: Builder<T=T, EntityDescriptor = ReferenceCellType>> GmshImpo
 #[cfg(test)]
 mod test {
     use super::*;
-    use crate::{shapes::regular_sphere, traits::Builder, SingleElementGridBuilder, SingleElementGrid};
-    use ndelement::{ciarlet::CiarletElement, map::IdentityMap};
+    use crate::{shapes::regular_sphere, traits::Builder, SingleElementGridBuilder};
 
     #[test]
     fn test_regular_sphere_gmsh_io() {
@@ -255,7 +254,7 @@ mod test {
     }
 
     #[test]
-    fn test_gmsh_output_quads() {
+    fn test_export_quads() {
         let mut b = SingleElementGridBuilder::<f64>::new(3, (ReferenceCellType::Quadrilateral, 1));
         b.add_point(0, &[0.0, 0.0, 0.0]);
         b.add_point(1, &[1.0, 0.0, 0.0]);
@@ -276,7 +275,7 @@ mod test {
     }
 
     #[test]
-    fn test_gmsh_output_tetrahedra() {
+    fn test_export_tetrahedra() {
         let mut b = SingleElementGridBuilder::<f64>::new(3, (ReferenceCellType::Tetrahedron, 1));
         b.add_point(0, &[0.0, 0.0, 0.0]);
         b.add_point(1, &[1.0, 0.0, 0.0]);
@@ -296,7 +295,7 @@ mod test {
         g.export_as_gmsh("_test_io_tetrahedra.msh");
     }
     #[test]
-    fn test_gmsh_output_hexahedra() {
+    fn test_export_hexahedra() {
         let mut b = SingleElementGridBuilder::<f64>::new(3, (ReferenceCellType::Hexahedron, 1));
         b.add_point(0, &[0.0, 0.0, 0.0]);
         b.add_point(1, &[1.0, 0.0, 0.0]);
@@ -317,9 +316,39 @@ mod test {
     }
 
     #[test]
-    fn test_gmsh_import_triangle() {
+    fn test_import_triangle() {
         let mut b = SingleElementGridBuilder::<f64>::new(3, (ReferenceCellType::Triangle, 1));
         b.import_from_gmsh("meshes/sphere_triangle.msh");
-        let g = b.create_grid();
+        let _g = b.create_grid();
+    }
+
+    #[test]
+    fn test_import_quadrilateral() {
+        let mut b = SingleElementGridBuilder::<f64>::new(3, (ReferenceCellType::Quadrilateral, 1));
+        b.import_from_gmsh("meshes/cube_quadrilateral.msh");
+        let _g = b.create_grid();
+    }
+
+    #[test]
+    fn test_import_tetrahedron() {
+        let mut b = SingleElementGridBuilder::<f64>::new(3, (ReferenceCellType::Tetrahedron, 1));
+        b.import_from_gmsh("meshes/cube_tetrahedron.msh");
+        let _g = b.create_grid();
+    }
+
+    #[test]
+    #[should_panic]
+    fn test_import_wrong_cell() {
+        let mut b = SingleElementGridBuilder::<f64>::new(3, (ReferenceCellType::Quadrilateral, 1));
+        b.import_from_gmsh("meshes/cube_tetrahedron.msh");
+        let _g = b.create_grid();
+    }
+
+    #[test]
+    #[should_panic]
+    fn test_import_wrong_degree() {
+        let mut b = SingleElementGridBuilder::<f64>::new(3, (ReferenceCellType::Triangle, 2));
+        b.import_from_gmsh("meshes/sphere_triangle.msh");
+        let _g = b.create_grid();
     }
 }

From 985c85b7adead4bccd705a859c0a751ba980d3a4 Mon Sep 17 00:00:00 2001
From: Matthew Scroggs <matthew.w.scroggs@gmail.com>
Date: Mon, 14 Jul 2025 17:52:20 +0100
Subject: [PATCH 3/7] working on gmsh io for higher order meshes

---
 src/geometry/point.rs                         |  18 +-
 .../single_element/entity_geometry.rs         |   4 +-
 src/grid/local_grid/single_element/builder.rs |   8 +-
 src/io/gmsh.rs                                | 181 +++++++++++++++---
 src/traits/builder.rs                         |   8 +-
 src/traits/geometry.rs                        |   3 +
 src/traits/io/gmsh.rs                         |   2 +-
 7 files changed, 187 insertions(+), 37 deletions(-)

diff --git a/src/geometry/point.rs b/src/geometry/point.rs
index eb63eb6..fa550cf 100644
--- a/src/geometry/point.rs
+++ b/src/geometry/point.rs
@@ -4,18 +4,23 @@ use crate::{traits::Point as PointTrait, types::RealScalar};
 /// A points
 #[derive(Debug, Clone, Copy)]
 pub struct Point<'a, T: RealScalar> {
+    index: usize,
     coordinates: &'a [T],
 }
 
 impl<'a, T: RealScalar> Point<'a, T> {
     /// Create new
-    pub fn new(coordinates: &'a [T]) -> Self {
-        Self { coordinates }
+    pub fn new(index: usize, coordinates: &'a [T]) -> Self {
+        Self { index, coordinates }
     }
 }
 impl<T: RealScalar> PointTrait for Point<'_, T> {
     type T = T;
 
+    fn index(&self) -> usize {
+        self.index
+    }
+
     fn dim(&self) -> usize {
         self.coordinates.len()
     }
@@ -28,12 +33,12 @@ impl<T: RealScalar> PointTrait for Point<'_, T> {
 /// Iterator over points
 #[derive(Debug)]
 pub struct PointIter<'a, T: RealScalar> {
-    points: Vec<&'a [T]>,
+    points: Vec<(usize, &'a [T])>,
     index: usize,
 }
 impl<'a, T: RealScalar> PointIter<'a, T> {
     /// Create new
-    pub fn new(points: Vec<&'a [T]>) -> Self {
+    pub fn new(points: Vec<(usize, &'a [T])>) -> Self {
         Self { points, index: 0 }
     }
 }
@@ -43,7 +48,10 @@ impl<'a, T: RealScalar> Iterator for PointIter<'a, T> {
     fn next(&mut self) -> Option<Point<'a, T>> {
         self.index += 1;
         if self.index <= self.points.len() {
-            Some(Point::new(self.points[self.index - 1]))
+            Some(Point::new(
+self.points[self.index - 1].0,
+self.points[self.index - 1].1,
+            ))
         } else {
             None
         }
diff --git a/src/geometry/single_element/entity_geometry.rs b/src/geometry/single_element/entity_geometry.rs
index 605615b..bb0e246 100644
--- a/src/geometry/single_element/entity_geometry.rs
+++ b/src/geometry/single_element/entity_geometry.rs
@@ -55,7 +55,7 @@ impl<T: RealScalar, E: FiniteElement> Geometry for SingleElementEntityGeometry<'
             .unwrap()
         {
             let i = self.geometry.cells()[[*index, self.cell_index]];
-            pts.push(&self.geometry.points().data()[i * gdim..(i + 1) * gdim])
+            pts.push((i, &self.geometry.points().data()[i * gdim..(i + 1) * gdim]))
         }
 
         PointIter::new(pts)
@@ -116,7 +116,7 @@ impl<T: RealScalar, E: FiniteElement> Geometry for SingleElementEntityGeometryBo
             .unwrap()
         {
             let i = self.geometry.cells()[[*index, self.cell_index]];
-            pts.push(&self.geometry.points().data()[i * gdim..(i + 1) * gdim])
+            pts.push((i, &self.geometry.points().data()[i * gdim..(i + 1) * gdim]))
         }
 
         PointIter::new(pts)
diff --git a/src/grid/local_grid/single_element/builder.rs b/src/grid/local_grid/single_element/builder.rs
index 074bd6a..33e088c 100644
--- a/src/grid/local_grid/single_element/builder.rs
+++ b/src/grid/local_grid/single_element/builder.rs
@@ -100,7 +100,13 @@ impl<T: RealScalar> Builder for SingleElementGridBuilder<T> {
         }
     }
 
-    fn add_cell_from_nodes_and_type(&mut self, id: usize, nodes: &[usize], cell_type: ReferenceCellType, cell_degree: usize) {
+    fn add_cell_from_nodes_and_type(
+        &mut self,
+        id: usize,
+        nodes: &[usize],
+        cell_type: ReferenceCellType,
+        cell_degree: usize,
+    ) {
         if (cell_type, cell_degree) != self.element_data {
             panic!("Invalid cell type.");
         }
diff --git a/src/io/gmsh.rs b/src/io/gmsh.rs
index 038141b..aa9bb13 100644
--- a/src/io/gmsh.rs
+++ b/src/io/gmsh.rs
@@ -1,10 +1,11 @@
 //! Gmsh I/O
 
-use crate::traits::{Builder, Entity, Geometry, GmshExport, Grid, Point, Topology, GmshImport};
+use crate::traits::{Builder, Entity, Geometry, GmshExport, GmshImport, Grid, Point, Topology};
 use itertools::izip;
 use ndelement::types::ReferenceCellType;
 use num::Zero;
 use std::str::FromStr;
+use std::collections::HashMap;
 
 fn get_permutation_to_gmsh(cell_type: ReferenceCellType, degree: usize) -> Vec<usize> {
     match cell_type {
@@ -103,7 +104,19 @@ impl<G: Grid<EntityDescriptor = ReferenceCellType>> GmshExport for G {
     fn to_gmsh_string(&self) -> String {
         let tdim = self.topology_dim();
         let gdim = self.geometry_dim();
-        let node_count = self.entity_count(ReferenceCellType::Point);
+
+        let mut points = HashMap::new();
+        for cell in self.cell_iter() {
+            for point in cell.geometry().points() {
+                let mut p = vec![G::T::zero(); gdim];
+                point.coords(&mut p);
+                points.insert(point.index(), p);
+            }
+        }
+        let mut points = points.iter().collect::<Vec<_>>();
+        points.sort_by(|i, j| i.0.cmp(j.0));
+        println!("{points:?}");
+        let node_count = points.len();
 
         let mut gmsh_s = String::from("");
         gmsh_s.push_str("$MeshFormat\n");
@@ -112,12 +125,10 @@ impl<G: Grid<EntityDescriptor = ReferenceCellType>> GmshExport for G {
         gmsh_s.push_str("$Nodes\n");
         gmsh_s.push_str(&format!("1 {node_count} 1 {node_count}\n"));
         gmsh_s.push_str(&format!("{tdim} 1 0 {node_count}\n"));
-        for i in 0..node_count {
-            gmsh_s.push_str(&format!("{}\n", i + 1));
+        for (i, _) in &points {
+            gmsh_s.push_str(&format!("{}\n", *i + 1));
         }
-        let mut coords = vec![G::T::zero(); gdim];
-        for node in self.entity_iter(0) {
-            node.geometry().points().next().unwrap().coords(&mut coords);
+        for (_, coords) in &points {
             for (n, c) in coords.iter().enumerate() {
                 if n != 0 {
                     gmsh_s.push(' ');
@@ -160,9 +171,9 @@ impl<G: Grid<EntityDescriptor = ReferenceCellType>> GmshExport for G {
             for (i, index) in cells.iter().enumerate() {
                 gmsh_s.push_str(&format!("{}", i + 1));
                 let entity = self.entity(tdim, *index).unwrap();
-                let topology = entity.topology();
+                let point_indices = entity.geometry().points().map(|i| i.index()).collect::<Vec<_>>();
                 for j in &gmsh_perm {
-                    gmsh_s.push_str(&format!(" {}", topology.sub_entity(0, *j) + 1));
+                    gmsh_s.push_str(&format!(" {}", point_indices[*j] + 1));
                 }
                 gmsh_s.push('\n');
             }
@@ -177,17 +188,19 @@ impl<G: Grid<EntityDescriptor = ReferenceCellType>> GmshExport for G {
 /// Get a section from a gmsh string
 fn gmsh_section(s: &str, section: &str) -> String {
     let a = s.split(&format!("${section}\n")).collect::<Vec<_>>();
-    if a.len() == 0 {
+    if a.len() <= 1 {
         panic!("Section not found: {section}");
     }
     String::from(a[1].split(&format!("\n$End{section}")).collect::<Vec<_>>()[0])
 }
 
-impl<T: FromStr, B: Builder<T=T, EntityDescriptor = ReferenceCellType>> GmshImport for B {
+impl<T: FromStr, B: Builder<T = T, EntityDescriptor = ReferenceCellType>> GmshImport for B {
     fn import_from_gmsh_string(&mut self, s: String) {
         let format = gmsh_section(&s, "MeshFormat");
         // Check msh file version
-        let [version, ascii_mode, _binary_mode] = format.split(" ").collect::<Vec<_>>()[..] else { panic!("Unrecognised gmsh format"); };
+        let [version, ascii_mode, _binary_mode] = format.split(" ").collect::<Vec<_>>()[..] else {
+            panic!("Unrecognised gmsh format");
+        };
         if version != "4.1" {
             unimplemented!("Unsupported gmsh file version");
         }
@@ -198,13 +211,23 @@ impl<T: FromStr, B: Builder<T=T, EntityDescriptor = ReferenceCellType>> GmshImpo
         let nodes = gmsh_section(&s, "Nodes");
         let nodes = nodes.lines().collect::<Vec<_>>();
 
-        let [num_entity_blocks, _num_nodes, _min_node_tag, _max_node_tag] = nodes[0].split(" ").map(|i| i.parse::<usize>().unwrap()).collect::<Vec<_>>()[..] else { panic!("Unrecognised gmsh format"); };
+        let [num_entity_blocks, _num_nodes, _min_node_tag, _max_node_tag] = nodes[0]
+            .split(" ")
+            .map(|i| i.parse::<usize>().unwrap())
+            .collect::<Vec<_>>()[..]
+        else {
+            panic!("Unrecognised gmsh format");
+        };
 
         let mut line_n = 1;
         for _ in 0..num_entity_blocks {
-            let [
-                _entity_dim, _entity_tag, parametric, num_nodes_in_block
-            ] = nodes[line_n].split(" ").map(|i| i.parse::<usize>().unwrap()).collect::<Vec<_>>()[..] else { panic!("Unrecognised gmsh format"); };
+            let [_entity_dim, _entity_tag, parametric, num_nodes_in_block] = nodes[line_n]
+                .split(" ")
+                .map(|i| i.parse::<usize>().unwrap())
+                .collect::<Vec<_>>()[..]
+            else {
+                panic!("Unrecognised gmsh format");
+            };
             if parametric == 1 {
                 unimplemented!("Parametric nodes currently not supported")
             }
@@ -212,29 +235,57 @@ impl<T: FromStr, B: Builder<T=T, EntityDescriptor = ReferenceCellType>> GmshImpo
             let tags = &nodes[line_n..line_n + num_nodes_in_block];
             let coords = &nodes[line_n + num_nodes_in_block..line_n + 2 * num_nodes_in_block];
             for (t, c) in izip!(tags, coords) {
-                self.add_point(t.parse::<usize>().unwrap(), &c.split(" ").map(|i| if let Ok(j) = T::from_str(i) {j} else {panic!("Could not parse coordinate");}).collect::<Vec<_>>());
+                self.add_point(
+                    t.parse::<usize>().unwrap(),
+                    &c.split(" ")
+                        .map(|i| {
+                            if let Ok(j) = T::from_str(i) {
+                                j
+                            } else {
+                                panic!("Could not parse coordinate");
+                            }
+                        })
+                        .collect::<Vec<_>>(),
+                );
             }
             line_n += num_nodes_in_block + 2;
         }
 
-
         // Load elements
         let elements = gmsh_section(&s, "Elements");
         let elements = elements.lines().collect::<Vec<_>>();
 
-        let [num_entity_blocks, _num_elements, _min_element_tag, _max_element_tag] = nodes[0].split(" ").map(|i| i.parse::<usize>().unwrap()).collect::<Vec<_>>()[..] else { panic!("Unrecognised gmsh format"); };
-        
+        let [num_entity_blocks, _num_elements, _min_element_tag, _max_element_tag] = nodes[0]
+            .split(" ")
+            .map(|i| i.parse::<usize>().unwrap())
+            .collect::<Vec<_>>()[..]
+        else {
+            panic!("Unrecognised gmsh format");
+        };
+
         let mut line_n = 1;
         for _ in 0..num_entity_blocks {
-            let [
-                _entity_dim, _entity_tag, element_type, num_elements_in_block
-            ] = elements[line_n].split(" ").map(|i| i.parse::<usize>().unwrap()).collect::<Vec<_>>()[..] else { panic!("Unrecognised gmsh format"); };
+            let [_entity_dim, _entity_tag, element_type, num_elements_in_block] = elements[line_n]
+                .split(" ")
+                .map(|i| i.parse::<usize>().unwrap())
+                .collect::<Vec<_>>()[..]
+            else {
+                panic!("Unrecognised gmsh format");
+            };
             let (cell_type, degree) = interpret_gmsh_cell(element_type);
+            let gmsh_perm = get_permutation_to_gmsh(cell_type, degree);
+
             line_n += 1;
             for line in &elements[line_n..line_n + num_elements_in_block] {
-                let line = line.split(" ").map(|i| i.parse::<usize>().unwrap()).collect::<Vec<_>>();
-                // TODO: permute
-                self.add_cell_from_nodes_and_type(line[0], &line[1..], cell_type, degree);
+                let line = line
+                    .split(" ")
+                    .map(|i| i.parse::<usize>().unwrap())
+                    .collect::<Vec<_>>();
+                let mut cell = vec![0; line.len() - 1];
+                for (i, j) in gmsh_perm.iter().enumerate() {
+                    cell[*j] = line[i + 1];
+                }
+                self.add_cell_from_nodes_and_type(line[0], &cell, cell_type, degree);
             }
 
             line_n += num_elements_in_block;
@@ -244,6 +295,7 @@ impl<T: FromStr, B: Builder<T=T, EntityDescriptor = ReferenceCellType>> GmshImpo
 
 #[cfg(test)]
 mod test {
+    use approx::*;
     use super::*;
     use crate::{shapes::regular_sphere, traits::Builder, SingleElementGridBuilder};
 
@@ -294,8 +346,9 @@ mod test {
         let g = b.create_grid();
         g.export_as_gmsh("_test_io_tetrahedra.msh");
     }
+
     #[test]
-    fn test_export_hexahedra() {
+    fn test_hexahedra() {
         let mut b = SingleElementGridBuilder::<f64>::new(3, (ReferenceCellType::Hexahedron, 1));
         b.add_point(0, &[0.0, 0.0, 0.0]);
         b.add_point(1, &[1.0, 0.0, 0.0]);
@@ -313,6 +366,80 @@ mod test {
         b.add_cell(2, &[4, 5, 6, 7, 8, 9, 10, 11]);
         let g = b.create_grid();
         g.export_as_gmsh("_test_io_hexahedra.msh");
+
+        let mut b = SingleElementGridBuilder::<f64>::new(3, (ReferenceCellType::Hexahedron, 1));
+        b.import_from_gmsh("_test_io_hexahedra.msh");
+        let g2 = b.create_grid();
+
+        let mut p1 = [0.0; 3];
+        let mut p2 = [0.0; 3];
+        for (v1, v2) in izip!(g.entity_iter(0), g2.entity_iter(0)) {
+            for (pt1, pt2) in izip!(v1.geometry().points(), v2.geometry().points()) {
+                pt1.coords(&mut p1);
+                pt2.coords(&mut p2);
+                for (c1, c2) in izip!(&p1, &p2) {
+                    assert_relative_eq!(c1, c2, epsilon=1e-10);
+                }
+            }
+        }
+        for (h1, h2) in izip!(g.entity_iter(3), g2.entity_iter(3)) {
+            for (v1, v2) in izip!(h1.topology().sub_entity_iter(0), h2.topology().sub_entity_iter(0)) {
+                assert_eq!(v1, v2);
+            }
+        }
+    }
+
+    #[test]
+    fn test_high_order_triangles() {
+        let mut b = SingleElementGridBuilder::<f64>::new(3, (ReferenceCellType::Triangle, 5));
+        b.add_point(0, &[0.0, 0.0, 0.0]);
+        b.add_point(1, &[0.0, 1.0, 0.0]);
+        b.add_point(2, &[0.0, 1.0, 1.0]);
+        b.add_point(3, &[1.0, 0.0, 0.0]);
+        b.add_point(4, &[1.0, 1.0, 0.0]);
+        b.add_point(5, &[1.0, 1.0, 1.0]);
+        b.add_point(6, &[2.0, 0.0, 0.0]);
+        b.add_point(7, &[2.0, 1.0, 0.0]);
+        b.add_point(8, &[2.0, 0.0, 1.0]);
+        b.add_point(9, &[3.0, 0.0, 0.0]);
+        b.add_point(10, &[3.0, 1.0, 0.0]);
+        b.add_point(11, &[3.0, 1.0, 1.0]);
+        b.add_point(12, &[4.0, 0.0, 0.0]);
+        b.add_point(13, &[4.0, 1.0, 0.0]);
+        b.add_point(14, &[4.0, 1.0, 1.0]);
+        b.add_point(15, &[5.0, 0.0, 0.0]);
+        b.add_point(16, &[5.0, 1.0, 0.0]);
+        b.add_point(17, &[5.0, 1.0, 1.0]);
+        b.add_point(18, &[6.0, 0.0, 0.0]);
+        b.add_point(19, &[6.0, 1.0, 0.0]);
+        b.add_point(20, &[6.0, 1.0, 1.0]);
+        b.add_cell(1, &[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20]);
+        let g = b.create_grid();
+        g.export_as_gmsh("_test_io_high_order_triangle.msh");
+
+        let mut b = SingleElementGridBuilder::<f64>::new(3, (ReferenceCellType::Triangle, 5));
+        b.import_from_gmsh("_test_io_high_order_triangle.msh");
+        let g2 = b.create_grid();
+
+        let mut p1 = [0.0; 3];
+        let mut p2 = [0.0; 3];
+        for (v1, v2) in izip!(g.entity_iter(0), g2.entity_iter(0)) {
+            v1.geometry().points().next().unwrap().coords(&mut p1);
+            v2.geometry().points().next().unwrap().coords(&mut p2);
+            println!("{p1:?} {p2:?}");
+        }
+        for (v1, v2) in izip!(g.entity_iter(0), g2.entity_iter(0)) {
+            v1.geometry().points().next().unwrap().coords(&mut p1);
+            v2.geometry().points().next().unwrap().coords(&mut p2);
+            for (c1, c2) in izip!(&p1, &p2) {
+                assert_relative_eq!(c1, c2, epsilon=1e-10);
+            }
+        }
+        for (h1, h2) in izip!(g.entity_iter(3), g2.entity_iter(3)) {
+            for (v1, v2) in izip!(h1.topology().sub_entity_iter(0), h2.topology().sub_entity_iter(0)) {
+                assert_eq!(v1, v2);
+            }
+        }
     }
 
     #[test]
diff --git a/src/traits/builder.rs b/src/traits/builder.rs
index 1c996f5..8a9c449 100644
--- a/src/traits/builder.rs
+++ b/src/traits/builder.rs
@@ -26,7 +26,13 @@ pub trait Builder {
     fn add_cell(&mut self, id: usize, cell_data: Self::CellData<'_>);
 
     /// Add a cell to the grid
-    fn add_cell_from_nodes_and_type(&mut self, id: usize, nodes: &[usize], cell_type: Self::EntityDescriptor, cell_degree: usize);
+    fn add_cell_from_nodes_and_type(
+        &mut self,
+        id: usize,
+        nodes: &[usize],
+        cell_type: Self::EntityDescriptor,
+        cell_degree: usize,
+    );
 
     /// Create the grid
     fn create_grid(&self) -> Self::Grid;
diff --git a/src/traits/geometry.rs b/src/traits/geometry.rs
index 1b45a34..70f224b 100644
--- a/src/traits/geometry.rs
+++ b/src/traits/geometry.rs
@@ -6,6 +6,9 @@ pub trait Point {
     /// Scalar type
     type T: RealScalar;
 
+    /// The point's index
+    fn index(&self) -> usize;
+
     /// Return the dimension of the point.
     fn dim(&self) -> usize;
 
diff --git a/src/traits/io/gmsh.rs b/src/traits/io/gmsh.rs
index d677a0b..a33d011 100644
--- a/src/traits/io/gmsh.rs
+++ b/src/traits/io/gmsh.rs
@@ -1,6 +1,6 @@
 //! Gmsh I/O
+use crate::traits::{Builder, Grid};
 use std::fs;
-use crate::traits::{Grid, Builder};
 
 pub trait GmshExport: Grid {
     //! Grid export for Gmsh

From 3251b16b114da7d81ca8b138762eb4e4aff64537 Mon Sep 17 00:00:00 2001
From: Matthew Scroggs <matthew.w.scroggs@gmail.com>
Date: Tue, 15 Jul 2025 08:13:24 +0100
Subject: [PATCH 4/7] fix test, fmt

---
 src/geometry/point.rs |  4 +--
 src/io/gmsh.rs        | 73 ++++++++++++++++++++++++++-----------------
 2 files changed, 46 insertions(+), 31 deletions(-)

diff --git a/src/geometry/point.rs b/src/geometry/point.rs
index fa550cf..cf10589 100644
--- a/src/geometry/point.rs
+++ b/src/geometry/point.rs
@@ -49,8 +49,8 @@ impl<'a, T: RealScalar> Iterator for PointIter<'a, T> {
         self.index += 1;
         if self.index <= self.points.len() {
             Some(Point::new(
-self.points[self.index - 1].0,
-self.points[self.index - 1].1,
+                self.points[self.index - 1].0,
+                self.points[self.index - 1].1,
             ))
         } else {
             None
diff --git a/src/io/gmsh.rs b/src/io/gmsh.rs
index aa9bb13..414529d 100644
--- a/src/io/gmsh.rs
+++ b/src/io/gmsh.rs
@@ -4,8 +4,8 @@ use crate::traits::{Builder, Entity, Geometry, GmshExport, GmshImport, Grid, Poi
 use itertools::izip;
 use ndelement::types::ReferenceCellType;
 use num::Zero;
-use std::str::FromStr;
 use std::collections::HashMap;
+use std::str::FromStr;
 
 fn get_permutation_to_gmsh(cell_type: ReferenceCellType, degree: usize) -> Vec<usize> {
     match cell_type {
@@ -171,7 +171,11 @@ impl<G: Grid<EntityDescriptor = ReferenceCellType>> GmshExport for G {
             for (i, index) in cells.iter().enumerate() {
                 gmsh_s.push_str(&format!("{}", i + 1));
                 let entity = self.entity(tdim, *index).unwrap();
-                let point_indices = entity.geometry().points().map(|i| i.index()).collect::<Vec<_>>();
+                let point_indices = entity
+                    .geometry()
+                    .points()
+                    .map(|i| i.index())
+                    .collect::<Vec<_>>();
                 for j in &gmsh_perm {
                     gmsh_s.push_str(&format!(" {}", point_indices[*j] + 1));
                 }
@@ -295,9 +299,9 @@ impl<T: FromStr, B: Builder<T = T, EntityDescriptor = ReferenceCellType>> GmshIm
 
 #[cfg(test)]
 mod test {
-    use approx::*;
     use super::*;
     use crate::{shapes::regular_sphere, traits::Builder, SingleElementGridBuilder};
+    use approx::*;
 
     #[test]
     fn test_regular_sphere_gmsh_io() {
@@ -378,12 +382,15 @@ mod test {
                 pt1.coords(&mut p1);
                 pt2.coords(&mut p2);
                 for (c1, c2) in izip!(&p1, &p2) {
-                    assert_relative_eq!(c1, c2, epsilon=1e-10);
+                    assert_relative_eq!(c1, c2, epsilon = 1e-10);
                 }
             }
         }
         for (h1, h2) in izip!(g.entity_iter(3), g2.entity_iter(3)) {
-            for (v1, v2) in izip!(h1.topology().sub_entity_iter(0), h2.topology().sub_entity_iter(0)) {
+            for (v1, v2) in izip!(
+                h1.topology().sub_entity_iter(0),
+                h2.topology().sub_entity_iter(0)
+            ) {
                 assert_eq!(v1, v2);
             }
         }
@@ -393,27 +400,32 @@ mod test {
     fn test_high_order_triangles() {
         let mut b = SingleElementGridBuilder::<f64>::new(3, (ReferenceCellType::Triangle, 5));
         b.add_point(0, &[0.0, 0.0, 0.0]);
-        b.add_point(1, &[0.0, 1.0, 0.0]);
-        b.add_point(2, &[0.0, 1.0, 1.0]);
-        b.add_point(3, &[1.0, 0.0, 0.0]);
-        b.add_point(4, &[1.0, 1.0, 0.0]);
-        b.add_point(5, &[1.0, 1.0, 1.0]);
-        b.add_point(6, &[2.0, 0.0, 0.0]);
-        b.add_point(7, &[2.0, 1.0, 0.0]);
-        b.add_point(8, &[2.0, 0.0, 1.0]);
-        b.add_point(9, &[3.0, 0.0, 0.0]);
-        b.add_point(10, &[3.0, 1.0, 0.0]);
-        b.add_point(11, &[3.0, 1.0, 1.0]);
-        b.add_point(12, &[4.0, 0.0, 0.0]);
-        b.add_point(13, &[4.0, 1.0, 0.0]);
-        b.add_point(14, &[4.0, 1.0, 1.0]);
-        b.add_point(15, &[5.0, 0.0, 0.0]);
-        b.add_point(16, &[5.0, 1.0, 0.0]);
-        b.add_point(17, &[5.0, 1.0, 1.0]);
-        b.add_point(18, &[6.0, 0.0, 0.0]);
-        b.add_point(19, &[6.0, 1.0, 0.0]);
-        b.add_point(20, &[6.0, 1.0, 1.0]);
-        b.add_cell(1, &[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20]);
+        b.add_point(1, &[5.0, 0.0, 0.2]);
+        b.add_point(2, &[0.0, 5.0, 0.4]);
+        b.add_point(3, &[4.0, 1.0, -0.1]);
+        b.add_point(4, &[3.0, 2.0, 0.2]);
+        b.add_point(5, &[2.0, 3.0, -0.3]);
+        b.add_point(6, &[1.0, 4.0, -0.5]);
+        b.add_point(7, &[0.0, 1.0, 0.6]);
+        b.add_point(8, &[0.0, 2.0, 0.2]);
+        b.add_point(9, &[0.0, 3.0, 0.1]);
+        b.add_point(10, &[0.0, 4.0, -0.2]);
+        b.add_point(11, &[1.0, 0.0, -0.3]);
+        b.add_point(12, &[2.0, 0.0, -0.4]);
+        b.add_point(13, &[3.0, 0.0, -0.5]);
+        b.add_point(14, &[4.0, 0.0, -0.2]);
+        b.add_point(15, &[1.0, 1.0, 0.1]);
+        b.add_point(16, &[2.0, 1.0, 0.1]);
+        b.add_point(17, &[3.0, 1.0, 0.1]);
+        b.add_point(18, &[2.0, 1.0, 0.2]);
+        b.add_point(19, &[2.0, 2.0, 0.1]);
+        b.add_point(20, &[3.0, 1.0, 0.1]);
+        b.add_cell(
+            1,
+            &[
+                0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
+            ],
+        );
         let g = b.create_grid();
         g.export_as_gmsh("_test_io_high_order_triangle.msh");
 
@@ -432,11 +444,14 @@ mod test {
             v1.geometry().points().next().unwrap().coords(&mut p1);
             v2.geometry().points().next().unwrap().coords(&mut p2);
             for (c1, c2) in izip!(&p1, &p2) {
-                assert_relative_eq!(c1, c2, epsilon=1e-10);
+                assert_relative_eq!(c1, c2, epsilon = 1e-10);
             }
         }
-        for (h1, h2) in izip!(g.entity_iter(3), g2.entity_iter(3)) {
-            for (v1, v2) in izip!(h1.topology().sub_entity_iter(0), h2.topology().sub_entity_iter(0)) {
+        for (h1, h2) in izip!(g.entity_iter(2), g2.entity_iter(2)) {
+            for (v1, v2) in izip!(
+                h1.topology().sub_entity_iter(0),
+                h2.topology().sub_entity_iter(0)
+            ) {
                 assert_eq!(v1, v2);
             }
         }

From 684b580c358feb2194b0c7b49bd91a20de3f1da5 Mon Sep 17 00:00:00 2001
From: Matthew Scroggs <matthew.w.scroggs@gmail.com>
Date: Tue, 15 Jul 2025 08:13:45 +0100
Subject: [PATCH 5/7] clippy

---
 src/io/gmsh.rs | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/io/gmsh.rs b/src/io/gmsh.rs
index 414529d..2d706b5 100644
--- a/src/io/gmsh.rs
+++ b/src/io/gmsh.rs
@@ -1,6 +1,6 @@
 //! Gmsh I/O
 
-use crate::traits::{Builder, Entity, Geometry, GmshExport, GmshImport, Grid, Point, Topology};
+use crate::traits::{Builder, Entity, Geometry, GmshExport, GmshImport, Grid, Point};
 use itertools::izip;
 use ndelement::types::ReferenceCellType;
 use num::Zero;

From f52d14c2783dfca71d8ceace7e434c9f2e6eec27 Mon Sep 17 00:00:00 2001
From: Matthew Scroggs <matthew.w.scroggs@gmail.com>
Date: Tue, 15 Jul 2025 08:18:09 +0100
Subject: [PATCH 6/7] use

---
 src/io/gmsh.rs | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/io/gmsh.rs b/src/io/gmsh.rs
index 2d706b5..219a030 100644
--- a/src/io/gmsh.rs
+++ b/src/io/gmsh.rs
@@ -300,7 +300,7 @@ impl<T: FromStr, B: Builder<T = T, EntityDescriptor = ReferenceCellType>> GmshIm
 #[cfg(test)]
 mod test {
     use super::*;
-    use crate::{shapes::regular_sphere, traits::Builder, SingleElementGridBuilder};
+    use crate::{shapes::regular_sphere, traits::{Builder, Topology}, SingleElementGridBuilder};
     use approx::*;
 
     #[test]

From 3772867d7ec81dbffe605c7fdf0594b089e0c85c Mon Sep 17 00:00:00 2001
From: Matthew Scroggs <matthew.w.scroggs@gmail.com>
Date: Tue, 15 Jul 2025 08:19:33 +0100
Subject: [PATCH 7/7] fmt

---
 src/io/gmsh.rs | 6 +++++-
 1 file changed, 5 insertions(+), 1 deletion(-)

diff --git a/src/io/gmsh.rs b/src/io/gmsh.rs
index 219a030..0de9b3b 100644
--- a/src/io/gmsh.rs
+++ b/src/io/gmsh.rs
@@ -300,7 +300,11 @@ impl<T: FromStr, B: Builder<T = T, EntityDescriptor = ReferenceCellType>> GmshIm
 #[cfg(test)]
 mod test {
     use super::*;
-    use crate::{shapes::regular_sphere, traits::{Builder, Topology}, SingleElementGridBuilder};
+    use crate::{
+        shapes::regular_sphere,
+        traits::{Builder, Topology},
+        SingleElementGridBuilder,
+    };
     use approx::*;
 
     #[test]