Read meshes in MFEM or VTK format to build Gridap models

[learning-chip]

I am attempting to port MFEM examples to Gridap. They should be useful additions to Gridap's existing tutorial examples.

The first step is to read example mesh files defined in MFEM-supported mesh formats. Most are in MFEM's own format (*.mesh) and some are in VTK format (*.vtk).

Gridap reads Gmsh formats via GridapGmsh.jl, but not other formats. So I need to either:

1. Find a robust tool to convert MFEM *.mesh -> Gmsh *.msh. Can be any standalone tools possibly not in Julia.

• Unfortunately, the versatile meshio does not support MFEM mesh (nschloe/meshio#456).
• MFEM's own IO utilities is able to read MFEM mesh (ReadMFEMMesh()), VTK mesh (ReadVTKMesh()) and GMSH (ReadGmshMesh()), but can only write out VTK format (Mesh::PrintVTK). There is no function to write out GMSH in MFEM; the last-resort is to write WriteGmshMesh() inside MFEM from scratch.

2. Write a parser in Julia from scratch, for MFEM-format mesh files. Then build Gripdap model/grid manually.
3. Use ReadVTK.jl to read VTK-format mesh files, and convert to Gridap's grid representation. The VTK files can be existing files in MFEM's example directory (e.g. beam-hex.vtk), or be converted from MFEM formats (via MFEM's ReadMFEMMesh() + PrintVTK).

Any suggestions on the quickest & most robust way for such mesh conversion? Thanks!

[fverdugo]

Hi @learning-chip, great initiative to port examples from MFEM!

Before deciding which is the best option, I wonder if the .mesh or .vtk files coming from MFEM contain all the info we need. In particular I am thinking about the face labels. E.g., If you don't import face label info you will only be able to apply boundary conditions on the entire boundary.

[learning-chip]

If you don't import face label info you will only be able to apply boundary conditions on the entire boundary

The MFEM mesh file does contain labels, for both boundary and inner elements. This ex2 example has two regions and different boundary conditions on two sides. The mesh file describes the boundary elements by:

boundary
18
3 1 0 1
...
3 1 17 16
1 1 9 0
2 1 8 17

The first column is the "tag" or "label" (1, 2, or 3 here), explained in MFEM's doc:

# Mesh faces/edges on the boundary, e.g. triangles (2)
boundary
<number of boundary elements>
<boundary element attribute> <geometry type> <vertex index 1> ... <vertex index m>
...

[fverdugo]

Then, It seems to contain all info we need. If the .mesh format it is not very complex and it is well documented, perhaps a Julia parser is a good option. If it is short and pure Julia, we can even add it to Gridap directly instead of creating another repo for this.

[learning-chip]

@fverdugo The MFEM *.mesh is a rather structured and simple text file, so I write a parser in Julia, by almost line-by-line translating from MFEM's C++ parser (Mesh::Loader and Mesh::ReadMFEMMesh).

Uncollapse the hidden block for complete parser code: (just to save space)

"""Check if the stream starts with comment_char. If so skip it."""
function skip_comment_lines!(io, comment_char)
    while true
        line = readline(io)
        if length(line) == 0
            break
        end
        if line[1] != comment_char
            break
        end
    end
end


"""Assume the next line is empty, and skip it"""
function skip_empty_line!(io)
    line = readline(io)
    @assert length(line) == 0 "not an empty line"
end


function read_element(io)
    line = readline(io)
    parse.(Int, string.(split(line)))
end


function read_vertex_coordinate(io)
    line = readline(io)
    parse.(Float64, string.(split(line)))
end


"""
Read MFEM mesh format

References:
- https://github.com/mfem/mfem/blob/v4.4/mesh/mesh_readers.cpp#L39-L106
- https://github.com/mfem/mfem/blob/v4.4/mesh/mesh.cpp#L3954
"""
function read_mfem_mesh(filename; has_comments=true)
    mesh = open(filename) do io
        # parse header and comments
        line = readline(io)
        @assert line == "MFEM mesh v1.0" "Only MFEM mesh v1.0 is supported now"

        skip_empty_line!(io)
        if has_comments
            skip_comment_lines!(io, '#')
        end

        # parse "dimension" block
        line = readline(io)
        @assert line == "dimension" "invalid mesh file"

        line = readline(io)
        dim = parse(Int, line)

        skip_empty_line!(io)

        # parse "elements" block
        line = readline(io)
        @assert line == "elements" "invalid mesh file"

        line = readline(io)
        num_of_elements = parse(Int, line)

        elements = []
        for j in 1:num_of_elements
            push!(elements, read_element(io))
        end
        elements = hcat(elements...)

        skip_empty_line!(io)

        # parse "boundary" block
        line = readline(io)
        @assert line == "boundary" "invalid mesh file"

        line = readline(io)
        num_of_bdr_elements = parse(Int, line)

        boundary = []
        for j in 1:num_of_bdr_elements
            push!(boundary, read_element(io))
        end
        boundary = hcat(boundary...)

        skip_empty_line!(io)

        # parse "vertices" block
        line = readline(io)
        @assert line == "vertices" "invalid mesh file"

        line = readline(io)
        num_of_vertices = parse(Int, line)

        line = readline(io)
        space_dim = parse(Int, line)

        vertices = Array{Float64}(undef, space_dim, num_of_vertices)
        for j in 1:num_of_vertices
            vertices[:,j] = read_vertex_coordinate(io)
        end

        # put all information in a dict
        mesh = Dict(
            :dim => dim,
            :num_of_elements => num_of_elements,
            :elements => elements,
            :num_of_bdr_elements => num_of_bdr_elements,
            :boundary => boundary,
            :num_of_vertices => num_of_vertices,
            :space_dim => space_dim,
            :vertices => vertices
        )
        return mesh
    end

    return mesh
end

Use it to read the example star.mesh:

mesh = read_mfem_mesh("./star.mesh")

which gives a dictionary containing all the mesh information:

Dict{Symbol, Any} with 8 entries:
  :num_of_elements     => 20
  :dim                 => 2
  :num_of_vertices     => 31
  :elements            => [1 1 … 1 1; 3 3 … 3 3; … ; 26 27 … 10 25; 14 17 … 25 …
  :boundary            => [1 1 … 1 1; 1 1 … 1 1; 13 12 … 10 8; 2 3 … 25 24]
  :vertices            => [0.0 1.0 … -0.25 0.654509; 0.0 0.0 … -0.76942 -0.4755…
  :num_of_bdr_elements => 20
  :space_dim           => 2

But I am still not sure how to build a valid Gridap DiscreteModel from such vanilla mesh information, and would like to have more guidance :)

My main references are the DiscreteModel class in Gridap core, the GmshDiscreteModels.jl constructor, and the serialized JSON form used in Gridap tutorial.

Some mappings are straightforward:

• Gridap model.grid.node_coordinates corresponds to MFEM mesh[:vertices], for spatial coordinates of vertices
• Gridap model.grid.cell_node_ids corresponds to MFEM mesh[:elements][3:end], to mark the vertice ids for each cell. (mesh[:elements][1] is cell tag to mark possibly different regions -- all 1 for uniform; mesh[:elements][2] is cell shape -- 3 is triangle as documented)

Some are less straighforward:

• face_labeling seems related to mesh[:boundary] for marking boundary cells (to later apply boundary conditions), but the exact mapping of information is not clear to me.
• cell_types probably means the cell shape (mesh[:elements][2])? Or the cell tag [mesh[:elements][1]]?

A Gridap model object actually contains a lot more attributes:

model = DiscreteModelFromFile("./model.json")  # https://github.com/gridap/Tutorials/blob/master/models/model.json

fieldnames(typeof(model))
# (:grid, :grid_topology, :face_labeling)
fieldnames(typeof(model.grid))
# (:node_coordinates, :cell_node_ids, :reffes, :cell_types, :orientation_style, :facet_normal, :cell_map)
fieldnames(typeof(model.grid_topology))
# (:vertex_coordinates, :n_m_to_nface_to_mfaces, :cell_type, :polytopes, :orientation_style)
fieldnames(typeof(model.face_labeling))
# (:d_to_dface_to_entity, :tag_to_entities, :tag_to_name)

Do I need to construct all those attributes in order to build a valid Gridap model/mesh? What is the "minimum required set of attributes"?

Just take this trivial star.mesh as an example. With the Julia parser script to read the mesh file, how would you build a valid Gripdap model mesh? Thanks!

[adigitoleo]

I've just found this package and have previously played around a bit with MFEM. I can't offer much help here yet but am just adding my vote for adopting a formal input mesh spec, and the MFEM spec has some benefits:

• well documented
• simple to write basic meshes by hand (or with simple scripts)
• simple to parse (see the excelent comment above)
• contains all (?) of the required metadata for a complete model
• can be easily visualised using MFEM's standalone glvis tool (using OpenGL so works quite nicely on many platforms)

In fact, the main drawback of MFEM meshes might be the fact that there is no GUI to create them, so with MFEM you have to write them by hand or use the C++ API. Constructing these meshes in the high-level Julia API that you have here would probably be more ergonomic.

@learning-chip I wonder if you had any further success in creating a parser for these meshes?

[adigitoleo]

From poking around a little, I think the type we need to construct would be an Gridap.Geometry.UnstructuredDiscreteModel.

julia> fieldnames(Gridap.Geometry.UnstructuredDiscreteModel)
(:grid, :grid_topology, :face_labeling)

help?> Gridap.Geometry.UnstructuredDiscreteModel
  struct UnstructuredDiscreteModel{Dc,Dp,Tp,B} <: DiscreteModel{Dc,Dp}
    grid::UnstructuredGrid{Dc,Dp,Tp,B}
    grid_topology::UnstructuredGridTopology{Dc,Dp,Tp,B}
    face_labeling::FaceLabeling
  end

The required components are the geometrical grid, the topology information and the face labels.

• Gridap model.grid.node_coordinates corresponds to MFEM mesh[:vertices]

Not sure if it has changed since you wrote this, but I believe this is not correct. MFEM "vertices" probably map to the topological vertex_coordinates. See here.

The rest of that page may help to create the mapping, although it seems there are still some differences. Maybe a collaborative effort with MFEM folks to come up with a complete mesh format would be possible? As far as I can tell, they are also still in active development and may be open to changing their format to suit more general needs.