add tet() functions for VTK

src/mesh/mesh_visualization.jl

@@ -73,12 +73,12 @@ RecipesBase.@recipe function f(m::VertexMeshPlotter{2})
 end
 
 """
-    MeshData_to_vtk(md, rd, dim, data, dataname, datatype, filename, write_data = false, equi_dist_nodes = true)
+    MeshData_to_vtk(md, rd, data, dataname, datatype, filename, write_data = false, equi_dist_nodes = true)
 
 Translate the given mesh into a vtk-file.
 `md` holds a `MeshData` object
 `rd` holds a reference element data/`RefElemData` object. 
-`data` holds an array of arrays (of size `num_nodes` by `num_elements`) with plotting data
+`data` holds an array of matrices (of size `num_nodes` by `num_elements`) with plotting data
 `dataname` is an array of strings with name of the associated data
 `write_data`, flag if data should be written or not (e.g., if data is not written, only the mesh will be saved as output)
 `equi_dist_nodes` flag if points should be interpolated to equidstant nodes
@@ -114,12 +114,12 @@ function MeshData_to_vtk(md::MeshData, rd::RefElemData{DIM}, data, dataname, fil
 end
 
 """
-MeshData_to_vtk(md, rd, dim, data, dataname, datatype, filename, write_data = false, equi_dist_nodes = true)
+MeshData_to_vtk(md, rd, data, dataname, datatype, filename, write_data = false, equi_dist_nodes = true)
 
 Translate the given mesh into a vtk-file.
 `md` holds a `MeshData` object
 `rd` holds a reference element data/`RefElemData` of a TensorProductWedge
-`data` holds an array of arrays (of size `num_nodes` by `num_elements`) with plotting data
+`data` holds an array of matrices (of size `num_nodes` by `num_elements`) with plotting data
 `dataname` is an array of strings with name of the associated data
 `write_data`, flag if data should be written or not (e.g., if data is not written, only the mesh will be saved as output)
 `equi_dist_nodes` flag if points should be interpolated to equidstant nodes

src/mesh/vtk_helper.jl

@@ -223,6 +223,94 @@
     return coords
 end
 
+"""
+    tet_vtk_order(corner_verts, order, dim, skip = false)
+
+Compute the coordinates of a VTK_LAGRANGE_TETRAHEDRON of a quad of order `order`
+defined by the coordinates of the vertices given in `corner_verts`. `dim` is the
+dimension of the coordinates given.
+
+Inspired by: https://github.com/ju-kreber/paraview-scripts/blob/master/node_ordering.py
+"""
+function tet_vtk_order(corner_verts, order, dim, skip = false)
+
+    @assert order >= 0 "`order` must be non-negative."
+
+    coords = copy(corner_verts) 
+
+    if order == 0
+        return coords
+    end
+
+    # Edges.
+    num_nodes_on_edge = order - 1
+    edges = SVector(
+      # VTK ordering. https://www.kitware.com/modeling-arbitrary-order-lagrange-finite-elements-in-the-visualization-toolkit/
+      (1,2), (2, 3), (3, 1), (1, 4), (2, 4), (3, 4)
+    )
+    for (frm, to) in edges
+        tmp = n_verts_between(num_nodes_on_edge, corner_verts[:, frm], corner_verts[:, to])
+        tmp_vec = Vector{Float64}(undef, dim)
+        for i in 1:num_nodes_on_edge
+            for j in 1:dim
+                tmp_vec[j] = tmp[j][i]
+            end
+            coords = hcat(coords, tmp_vec)
+        end
+    end
+
+    # Faces.
+    tri_faces = SVector(
+        # VTK ordering. 
+        (1, 2, 4), # left
+        (3, 4, 2), # right
+        (1, 4, 3), # front
+        (1, 3, 2), # back
+        )
+    for indices in tri_faces
+        tmp_vec = Vector{Float64}(undef, dim)
+        tri_nodes = Matrix{Float64}(undef, dim, 0)
+        for j in indices
+            tri_nodes = hcat(tri_nodes, corner_verts[:, j])
+        end
+        face_coords = triangle_vtk_order(tri_nodes, order, 3, true)
+        if length(face_coords) > 0
+            for i in 1:size(face_coords)[2]
+                for j in 1:dim
+                    tmp_vec[j] = face_coords[j,i]
+                end
+                coords = hcat(coords, tmp_vec)
+            end
+        end
+    end
+
+    # Volume. Recursively find nodes in volume.
+    e_12 = (corner_verts[:,2] - corner_verts[:, 1]) ./ order
+    e_13 = (corner_verts[:,3] - corner_verts[:, 1]) ./ order
+    e_14 = (corner_verts[:,4] - corner_verts[:, 1]) ./ order
+    e_23 = (corner_verts[:,3] - corner_verts[:, 2]) ./ order
+    e_24 = (corner_verts[:,4] - corner_verts[:, 2]) ./ order
+    e_34 = (corner_verts[:,4] - corner_verts[:, 3]) ./ order
+
+    # Find inner corners of reduced order tetrahedron
+    a_1 = corner_verts[:, 1] + e_12 + e_13 + e_14
+    a_2 = corner_verts[:, 2] - e_12 + e_23 + e_24
+    a_3 = corner_verts[:, 3] - e_13 - e_23 + e_34
+    a_4 = corner_verts[:, 4] - e_34 - e_24 - e_14
+
+    interior_verts = hcat(a_1, a_2, a_3, a_4)
+    if order > 4
+        tmp_coords = tet_vtk_order(interior_verts, order - 4, dim)
+        coords = hcat(coords, tmp_coords)
+    # Order 4 tetrahedra will have one node in the center, we do this to account for not having an implementation for N = 0.
+    elseif order == 4
+        center = (corner_verts[:, 1] + corner_verts[:, 2] + corner_verts[:, 3] + corner_verts[:, 4]) / 4
+        coords = hcat(coords, center)
+    end
+
+    return coords
+end
+
 """
     wedge_vtk_order(corner_verts, order, dim)
 
@@ -362,6 +450,20 @@
     return wedge_vtk_order(wedge_vtk_vertices, order, 3)
 end
 
+"""
+    vtk_order(elem::Tet, order)
+
+Construct all node-points of a VTK_LAGRANGE_TETRAHEDRON of order `order`. The corner-nodes are
+given by the reference-tet used by StartUpDG
+"""
+function vtk_order(::Tet, order)
+    tet_sud_vertices = permutedims(hcat(nodes(Tet(), 1)...))
+    perm = SVector(1, 2, 3, 4)
+    tet_vtk_vertices = tet_sud_vertices[:, perm]
+    return tet_vtk_order(tet_vtk_vertices, order, 3)
+end
+
+
 """
     SUD_to_vtk_order(rd::RefElemData, dim)
 
@@ -413,5 +515,10 @@
     return VTKCellTypes.VTK_LAGRANGE_WEDGE
 end
 
-
-
+"""
+    type_to_vtk(elem::Tet)
+    return the VTK-type
+"""
+function type_to_vtk(elem::Tet)
+    return VTKCellTypes.VTK_LAGRANGE_TETRAHEDRON
+end

test/write_vtk_tests.jl

@@ -14,10 +14,13 @@
     deg_one_order(::Wedge) = [-1.0  -1.0   1.0  -1.0  -1.0   1.0
                               -1.0   1.0  -1.0  -1.0   1.0  -1.0
                               -1.0  -1.0  -1.0   1.0   1.0   1.0]
-    deg_zero_order(elem::Union{Quad, Hex, Wedge}) = deg_one_order(elem)
+    deg_one_order(::Tet) = [-1.0  1.0 -1.0 -1.0;
+                            -1.0 -1.0  1.0 -1.0;
+                            -1.0 -1.0 -1.0  1.0]
+    deg_zero_order(elem::Union{Quad, Hex, Wedge, Tet}) = deg_one_order(elem)
 
 
-    @testset "VTKWriter test for $elem" for elem in [Tri(), Quad(), Hex(), Wedge()]
+    @testset "VTKWriter test for $elem" for elem in [Tri(), Quad(), Hex(), Wedge(), Tet()]
         N = 2 # test only N=2 for CI time
         @testset "Write Mesh" begin
             rd = RefElemData(elem, N)