Export to VTK support

Project.toml

@@ -5,6 +5,7 @@ version = "0.1.0"
 [deps]
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192"
 
 [weakdeps]
 Gmsh = "705231aa-382f-11e9-3f0c-b7cb4346fdeb"

src/Inti.jl

@@ -4,6 +4,7 @@ const PROJECT_ROOT = pkgdir(Inti)
 
 using StaticArrays
 using Printf
+using WriteVTK
 
 # helper functions
 include("utils.jl")
@@ -19,6 +20,9 @@ include("entities.jl")
 include("domain.jl")
 include("mesh.jl")
 
+# IO
+include("vtk.jl")
+
 # # integral operators
 # include("integral_potentials.jl")
 # include("integral_operators.jl")

src/mesh.jl

@@ -76,6 +76,10 @@ Return the element types present in the `msh`.
 """
 element_types(msh::LagrangeMesh) = keys(msh.etype2mat)
 
+nodes(msh::LagrangeMesh) = msh.nodes
+elements(msh::LagrangeMesh) = msh.etype2mat
+ent2tags(msh::LagrangeMesh) = msh.ent2tags
+
 """
     dom2elt(m::LagrangeMesh,Ω,E)
 

src/reference_interpolation.jl

@@ -30,11 +30,15 @@ function jacobian(f, x)
     return interface_method(f)
 end
 
+# TODO Should we use SType here?
 domain(::ReferenceInterpolant{D,T}) where {D,T} = D()
 domain(::Type{<:ReferenceInterpolant{D,T}}) where {D,T} = D()
 return_type(::ReferenceInterpolant{D,T}) where {D,T} = T
-domain_dimension(t::ReferenceInterpolant) = domain(t) |> center |> length
+return_type(::Type{<:ReferenceInterpolant{D,T}}) where {D,T} = T
+domain_dimension(t::ReferenceInterpolant{D,T}) where {D,T} = domain(t) |> center |> length
+domain_dimension(t::Type{<:ReferenceInterpolant{D,T}}) where {D,T} = domain(t) |> center |> length
 range_dimension(el::ReferenceInterpolant{R,T}) where {R,T} = domain(el) |> center |> el |> length
+range_dimension(el::Type{<:ReferenceInterpolant{R,T}}) where {R,T} = domain(el) |> center |> el |> length
 
 """
     struct LagrangeElement{D,Np,T} <: ReferenceInterpolant{D,T}

src/vtk.jl

@@ -0,0 +1,144 @@
+@info "including vtk.jl"
+
+using WriteVTK
+
+"""
+    vtk_mesh_file(mesh::LagrangeMesh[, Ω::Domain], name::String)
+
+Creates a `VTK` file (.vtu) with name `name` containing the mesh information.
+It is possible to export only a `Domain` (i.e. only a part of the mesh).
+
+Note that all the heavy lifting is done by the package `WriteVTK`.
+We refer to its documentation for more information.
+
+Warning: the output is not an actual file (on disk).
+To save it, simply write:
+
+    vtk_mesh_file(mesh, "my_mesh") |> vtk_save
+
+(Noting that you will need to be `using WriteVTK` to do so)
+
+To add some data (scalar or vector values at the mesh nodes or mesh cells)
+for visualization purposes, you can do for instance:
+
+    vtkfile = vtk_mesh_file(mesh, name)
+    vtkfile["my_point_data", VTKPointData()] = pdata
+    vtkfile["my_cell_data", VTKCellData()] = cdata
+    vtk_save(vtkfile)
+
+It is possible also to export a partition `Ωs::Vector{Domain}` using
+`Multiblock` files (.vtm), for instance like so
+
+    vtmfile = vtk_multiblock(name)
+    for (Ω, pdata) in zip(Ωs, pdatas)
+        vtkfile = vtk_mesh_file(mesh, Ω, name)
+        vtkfile["my_point_data", VTKPointData()] = pdata
+    end
+    vtk_save(vtmfile)
+
+To save a sequence of solutions (time steps, iterations), simply append
+the number of the element to the file name.
+Paraview will recognize the sequence automatically.
+"""
+function vtk_mesh_file(mesh::LagrangeMesh, name::String)
+    points = _vtk_points(mesh)
+    cells = _vtk_cells(mesh)
+    return vtk_grid(name * ".vtu", points, cells)
+end
+function vtk_mesh_file(mesh::LagrangeMesh, Ω::Domain, name::String)
+    points = _vtk_points(mesh)
+    cells = _vtk_cells(mesh, Ω)
+    return vtk_grid(name * ".vtu", points, cells)
+end
+
+"""
+    _vtk_points(mesh::LagrangeMesh)
+
+Creates the matrix of nodes in the format required by `WriteVTK`.
+"""
+function _vtk_points(mesh::LagrangeMesh)
+    vtx = zeros(Float64, 3, length(nodes(mesh)))
+    for (i, nd) in enumerate(nodes(mesh))
+        vtx[1:ambient_dimension(mesh), i] = nd
+    end
+    return vtx
+end
+
+"""
+    _vtk_cells(mesh::LagrangeMesh, E::DataType)
+    _vtk_cells(mesh::LagrangeMesh, Ω::Domain)
+
+Creates the vector of all elements contained in the mesh in the format
+required by `WriteVTK` for a particular element type `E<:AbstractElement`
+or a `Domain` instance.
+"""
+function _vtk_cells end
+
+function _vtk_cells(tags, E::DataType)
+    vtk_cell_type, ind = etype_to_vtk_cell_type[E]
+    return [MeshCell(vtk_cell_type, tags[ind, i]) for i in 1:size(tags, 2)]
+end
+function _vtk_cells(mesh::LagrangeMesh, Ω::Domain)
+    cells = MeshCell[]
+    # Loop on `ElementaryEntity`
+    for ω in Ω
+        # Loop on `AbstractElement`
+        for (E, ind) in ent2tags(mesh)[ω]
+            # Subset corresponding to the `ElementaryEntity`
+            tags = elements(mesh)[E][:, ind]
+            append!(cells, _vtk_cells(tags, E))
+        end
+    end
+    return cells
+end
+function _vtk_cells(mesh::LagrangeMesh)
+    cells = MeshCell[]
+    # Loop on `AbstractElement`
+    for (E, tags) in elements(mesh)
+        # Export only the cells of the largest geometrical dimension
+        if domain_dimension(E) == ambient_dimension(mesh)
+            append!(cells, _vtk_cells(tags, E))
+        end
+    end
+    return cells
+end
+
+"""
+    const etype_to_vtk_cell_type
+
+Dictionary mapping internal element types to a tuple containing:
+    - the corresponding `WriteVTK` cell types (following the convention
+    chosen by `VTK`, see below);
+    - the indices in the `elements` column that defines the element.
+    This is because we want to support the export of more than just the flat
+    elements available in the `VTK` specification, hence which may require
+    a conversion of some sort.
+
+See VTK specification [Fig. 2] on
+[http://www.vtk.org/VTK/img/file-formats.pdf](http://www.vtk.org/VTK/img/file-formats.pdf)
+
+- VTK_VERTEX (=1)
+- VTK_POLY_VERTEX (=2)
+- VTK_LINE (=3)
+- VTK_POLY_LINE (=4)
+- VTK_TRIANGLE (=5)
+- VTK_TRIANGLE_STRIP (=6)
+- VTK_POLYGON (=7)
+- VTK_PIXEL (=8)
+- VTK_QUAD (=9)
+- VTK_TETRA (=10)
+- VTK_VOXEL (=11)
+- VTK_HEXAHEDRON (=12)
+- VTK_WEDGE (=13)
+- VTK_PYRAMID (=14)
+"""
+const etype_to_vtk_cell_type = Dict(SVector{3,Float64} => (VTKCellTypes.VTK_VERTEX,
+                                                           collect(1:1)),
+                                    LagrangeLine{2,SVector{3,Float64}} => (VTKCellTypes.VTK_LINE,
+                                                                           collect(1:2)),
+                                    LagrangeTriangle{3,SVector{2,Float64}} => (VTKCellTypes.VTK_TRIANGLE,
+                                                                               collect(1:3)),
+                                    LagrangeTriangle{3,SVector{3,Float64}} => (VTKCellTypes.VTK_TRIANGLE,
+                                                                               collect(1:3)),
+                                    LagrangeTetrahedron{4,SVector{3,Float64}} => (VTKCellTypes.VTK_TETRA,
+                                                                                  collect(1:4)))

test/gmsh_test.jl

@@ -13,7 +13,7 @@ using Test
     Ω   = Inti.gmsh_import_domain(;dim=2)
     msh = Inti.gmsh_import_mesh(Ω;dim=2)
     gmsh.finalize()
-    return true
+    true == true
 end
 
 @test begin
@@ -27,5 +27,5 @@ end
     Ω   = Inti.gmsh_import_domain(;dim=3)
     msh = Inti.gmsh_import_mesh(Ω;dim=3)
     gmsh.finalize()
-    return true
+    true == true
 end

test/vtkIO_test.jl

@@ -0,0 +1,28 @@
+using Test
+using Gmsh
+using Inti
+using WriteVTK
+
+@test begin
+# This test should simply not throw an error
+# TODO Replace the gmsh code with reading from a vtk file; currently tests writing
+    gmsh.initialize()
+    gmsh.option.setNumber("General.Verbosity", 2)
+    gmsh.model.add("Sphere")
+    gmsh.model.occ.addSphere(0,0,0,1)
+    gmsh.model.occ.synchronize()
+    gmsh.model.mesh.generate(3)
+    ents = gmsh.model.getEntities()
+    Ω = Inti.gmsh_import_domain(;dim=3)
+    M = Inti.gmsh_import_mesh(Ω;dim=3)
+    vtk_save(Inti.vtk_mesh_file(M, joinpath(Inti.PROJECT_ROOT, "test", "ball")))
+    rm(joinpath(Inti.PROJECT_ROOT, "test", "ball.vtu"))
+    vtk_save(Inti.vtk_mesh_file(M, Ω, joinpath(Inti.PROJECT_ROOT, "test", "ball")))
+    rm(joinpath(Inti.PROJECT_ROOT, "test", "ball.vtu"))
+    vtk_save(Inti.vtk_mesh_file(M,
+                               Inti.external_boundary(Ω),
+                               joinpath(Inti.PROJECT_ROOT, "test", "sphere")))
+    rm(joinpath(@__DIR__, "sphere.vtu"))
+    gmsh.finalize()
+    true == true
+end