Implement basic parser for triangular meshes from .MSH version 4 (GMSH)

src/MeshIO.jl

@@ -6,7 +6,8 @@ using ColorTypes
 using Printf
 import FileIO
 
-import FileIO: DataFormat, @format_str, Stream, File, filename, stream, skipmagic
+import FileIO: DataFormat, @format_str, Stream, File, filename, stream,
+    skipmagic, add_format
 import Base.show
 
 
@@ -16,6 +17,7 @@ include("io/ply.jl")
 include("io/stl.jl")
 include("io/obj.jl")
 include("io/2dm.jl")
+include("io/msh.jl")
 
 load(fn::File{format}, MeshType=GLNormalMesh) where {format} = open(fn) do s
     skipmagic(s)

src/io/msh.jl

@@ -0,0 +1,128 @@
+
+function load(fs::Stream{format"MSH"}, MeshType=GLNormalMesh)
+    #GMSH MSH format (version 4)
+    #http://gmsh.info/doc/texinfo/gmsh.html#MSH-file-format
+    io = stream(fs)
+
+    FaceType = facetype(MeshType)
+    VertexType = vertextype(MeshType)
+    faces = FaceType[]
+    nodes = VertexType[]
+    node_tags = Int[]
+
+    while !eof(io)
+        BlockType = _parse_blocktype!(io)
+        if BlockType == MSHNodesBlock()
+            _parse_nodes!(io, nodes, node_tags)
+        elseif BlockType == MSHElementsBlock()
+            _parse_elements!(io, faces)
+        else
+            _skip_block!(io)
+        end
+    end
+
+    _remap_faces!(faces, node_tags)
+
+    return MeshType(nodes, faces)
+end
+
+struct MSHFormatBlock end
+struct MSHNodesBlock end
+struct MSHElementsBlock end
+struct MSHUnknownBlock end
+
+function _parse_blocktype!(io)
+    header = readline(io)
+    if header == "\$MeshFormat"
+        return MSHFormatBlock()
+    elseif header == "\$Nodes"
+        return MSHNodesBlock()
+    elseif header == "\$Elements"
+        return MSHElementsBlock()
+    else
+        return MSHUnknownBlock()
+    end
+end
+
+function _parse_format!(io)
+    version, binary, size = map(parse, (Float64, Int, Int), split(readline(io)))
+    if version < 4
+        error("version $(version[1]) < 4.0 not supported.")
+    elseif binary == 1
+        error("binary format not supported.")
+    end
+    endblock = readline(io)
+    if endblock != "\$EndMeshFormat"
+        error("expected block end tag, got $endblock.")
+    end
+    return version
+end
+
+function _skip_block!(io)
+    while true
+        line = readline(io)
+        if line[1:4] == "\$End"
+            break
+        end
+    end
+    return nothing
+end
+
+function _parse_nodes!(io, nodes, node_tags)
+    numEntityBlocks, numNodes, minNodeTag, maxNodeTag =
+        map(parse, (Int, Int, Int, Int), split(readline(io)))
+    for index_entity in 1:numEntityBlocks
+        entityDim, entityTag, parametric, numNodesInBlock =
+            map(parse, (Int, Int, Int, Int), split(readline(io)))
+        for i in 1:numNodesInBlock
+            push!(node_tags, parse(eltype(node_tags), readline(io)))
+        end
+        for i in 1:numNodesInBlock
+            x, y, z = map(parse, (Float64, Float64, Float64), split(readline(io)))
+            push!(nodes, eltype(nodes)(x, y, z))
+        end
+    end
+    endblock = readline(io)
+    if endblock != "\$EndNodes"
+        error("expected end block tag, got $endblock")
+    end
+    return nodes, node_tags
+end
+
+function _parse_elements!(io, faces::Vector{T}) where T <: Triangle
+    numEntityBlocks, numElements, minElementTag, maxElementTag =
+        map(parse, (Int, Int, Int, Int), split(readline(io)))
+    for index_entity in 1:numEntityBlocks
+        entityDim, entityTag, elementType, numElementsInBlock =
+            map(parse, (Int, Int, Int, Int), split(readline(io)))
+        if elementType == 2 # Triangles
+            for i in 1:numElementsInBlock
+                tag, n1, n2, n3 =
+                    map(parse, (Int, Int, Int, Int), split(readline(io)))
+                push!(faces, eltype(faces)(n1, n2, n3))
+            end
+        else
+            # for now we ignore all other elements (points, lines, hedrons, etc)
+            for i in 1:numElementsInBlock
+                readline(io)
+            end
+        end
+    end
+    endblock = readline(io)
+    if endblock != "\$EndElements"
+        error("expected end block tag, got $endblock")
+    end
+    return faces
+end
+
+function _remap_faces!(faces, node_tags)
+    node_map = indexin(1:maximum(node_tags), node_tags)
+    for (i, face) in enumerate(faces)
+        faces[i] = eltype(faces)(
+            node_map[face[1]],
+            node_map[face[2]],
+            node_map[face[3]]
+        )
+    end
+    return faces
+end

test/runtests.jl

@@ -158,6 +158,13 @@ end
             #@test length(vertices(msh)) == 2248
             #@test length(normals(msh)) == 2248
         end
+        @testset "GMSH" begin
+            msh = load(joinpath(tf, "cube.msh"))
+            @test typeof(msh) == GLNormalMesh
+            @test length(faces(msh)) == 24
+            @test length(vertices(msh)) == 14
+            test_face_indices(msh)
+        end
     end
 end
 

test/testfiles/cube.msh

@@ -0,0 +1,164 @@
+$MeshFormat
+4.1 0 8
+$EndMeshFormat
+$Entities
+8 12 6 1
+1 0 0 1 0 
+2 0 0 0 0 
+3 0 1 1 0 
+4 0 1 0 0 
+5 1 0 1 0 
+6 1 0 0 0 
+7 1 1 1 0 
+8 1 1 0 0 
+1 -1e-07 -1e-07 -9.999999994736442e-08 1e-07 1e-07 1.0000001 0 2 2 -1 
+2 -1e-07 -9.999999994736442e-08 0.9999999000000001 1e-07 1.0000001 1.0000001 0 2 1 -3 
+3 -1e-07 0.9999999000000001 -9.999999994736442e-08 1e-07 1.0000001 1.0000001 0 2 4 -3 
+4 -1e-07 -9.999999994736442e-08 -1e-07 1e-07 1.0000001 1e-07 0 2 2 -4 
+5 0.9999999000000001 -1e-07 -9.999999994736442e-08 1.0000001 1e-07 1.0000001 0 2 6 -5 
+6 0.9999999000000001 -9.999999994736442e-08 0.9999999000000001 1.0000001 1.0000001 1.0000001 0 2 5 -7 
+7 0.9999999000000001 0.9999999000000001 -9.999999994736442e-08 1.0000001 1.0000001 1.0000001 0 2 8 -7 
+8 0.9999999000000001 -9.999999994736442e-08 -1e-07 1.0000001 1.0000001 1e-07 0 2 6 -8 
+9 -9.999999994736442e-08 -1e-07 -1e-07 1.0000001 1e-07 1e-07 0 2 2 -6 
+10 -9.999999994736442e-08 -1e-07 0.9999999000000001 1.0000001 1e-07 1.0000001 0 2 1 -5 
+11 -9.999999994736442e-08 0.9999999000000001 -1e-07 1.0000001 1.0000001 1e-07 0 2 4 -8 
+12 -9.999999994736442e-08 0.9999999000000001 0.9999999000000001 1.0000001 1.0000001 1.0000001 0 2 3 -7 
+1 -1e-07 -9.999999994736442e-08 -9.999999994736442e-08 1e-07 1.0000001 1.0000001 0 4 1 2 -3 -4 
+2 0.9999999000000001 -9.999999994736442e-08 -9.999999994736442e-08 1.0000001 1.0000001 1.0000001 0 4 5 6 -7 -8 
+3 -9.999999994736442e-08 -1e-07 -9.999999994736442e-08 1.0000001 1e-07 1.0000001 0 4 9 5 -10 -1 
+4 -9.999999994736442e-08 0.9999999000000001 -9.999999994736442e-08 1.0000001 1.0000001 1.0000001 0 4 11 7 -12 -3 
+5 -9.999999994736442e-08 -9.999999994736442e-08 -1e-07 1.0000001 1.0000001 1e-07 0 4 4 11 -8 -9 
+6 -9.999999994736442e-08 -9.999999994736442e-08 0.9999999000000001 1.0000001 1.0000001 1.0000001 0 4 2 12 -6 -10 
+1 -9.999999994736442e-08 -9.999999994736442e-08 -9.999999994736442e-08 1.0000001 1.0000001 1.0000001 0 6 1 2 3 4 5 6 
+$EndEntities
+$Nodes
+27 14 1 14
+0 1 0 1
+1
+0 0 1
+0 2 0 1
+2
+0 0 0
+0 3 0 1
+3
+0 1 1
+0 4 0 1
+4
+0 1 0
+0 5 0 1
+5
+1 0 1
+0 6 0 1
+6
+1 0 0
+0 7 0 1
+7
+1 1 1
+0 8 0 1
+8
+1 1 0
+1 1 0 0
+1 2 0 0
+1 3 0 0
+1 4 0 0
+1 5 0 0
+1 6 0 0
+1 7 0 0
+1 8 0 0
+1 9 0 0
+1 10 0 0
+1 11 0 0
+1 12 0 0
+2 1 0 1
+9
+0 0.5 0.5
+2 2 0 1
+10
+1 0.5 0.5
+2 3 0 1
+11
+0.5 0 0.5
+2 4 0 1
+12
+0.5 1 0.5
+2 5 0 1
+13
+0.5 0.5 0
+2 6 0 1
+14
+0.5 0.5 1
+3 1 0 0
+$EndNodes
+$Elements
+26 44 1 44
+0 1 15 1
+1 1 
+0 2 15 1
+2 2 
+0 3 15 1
+3 3 
+0 4 15 1
+4 4 
+0 5 15 1
+5 5 
+0 6 15 1
+6 6 
+0 7 15 1
+7 7 
+0 8 15 1
+8 8 
+1 1 1 1
+9 2 1 
+1 2 1 1
+10 1 3 
+1 3 1 1
+11 4 3 
+1 4 1 1
+12 2 4 
+1 5 1 1
+13 6 5 
+1 6 1 1
+14 5 7 
+1 7 1 1
+15 8 7 
+1 8 1 1
+16 6 8 
+1 9 1 1
+17 2 6 
+1 10 1 1
+18 1 5 
+1 11 1 1
+19 4 8 
+1 12 1 1
+20 3 7 
+2 1 2 4
+21 2 1 9 
+22 1 3 9 
+23 4 2 9 
+24 3 4 9 
+2 2 2 4
+25 6 10 5 
+26 5 10 7 
+27 8 10 6 
+28 7 10 8 
+2 3 2 4
+29 1 2 11 
+30 5 1 11 
+31 2 6 11 
+32 6 5 11 
+2 4 2 4
+33 3 12 4 
+34 7 12 3 
+35 4 12 8 
+36 8 12 7 
+2 5 2 4
+37 2 4 13 
+38 6 2 13 
+39 4 8 13 
+40 8 6 13 
+2 6 2 4
+41 1 14 3 
+42 5 14 1 
+43 3 14 7 
+44 7 14 5 
+$EndElements