HOHQMesh curved tet mesh reader (#96)

* add fixed mesh * add test for tets * fix tet reader number of faces * start work on HOHQMesh tet reader * add curved tet reader * remove TODO
jlchan · Jul 27, 2023 · 8696c8c · 8696c8c
1 parent 07bc9cc
commit 8696c8c
Show file tree

Hide file tree

Showing 3 changed files with 39,761 additions and 39,682 deletions.
diff --git a/src/mesh/hohqmesh_utilities.jl b/src/mesh/hohqmesh_utilities.jl
@@ -20,6 +20,7 @@ vertex_reordering(::Tri) = SVector(1, 3, 2)
 # for more details on HOHQMesh's face ordering. 
 get_HOHQMesh_to_StartUpDG_face_ordering(::Quad) = SVector(3, 2, 4, 1) 
 get_HOHQMesh_to_StartUpDG_face_ordering(::Tri) = SVector(1, 2, 3) 
+get_HOHQMesh_to_StartUpDG_face_ordering(::Tet) = 1:4
 
 function get_HOHQMesh_ids(::Quad, curved_faces, f_HOHQMesh, f, polydeg) 
     active_face_offset = max.(0, cumsum(curved_faces) .- 1) * (polydeg + 1)
@@ -167,6 +168,83 @@ function MeshData(hmd::HOHQMeshData{3}, rd::RefElemData{3, <:Hex})
     return @set md_curved.mesh_type = HOHQMeshType(hmd, boundary_faces)
 end
 
+vertex_reordering(::Tet) = SVector(1, 2, 3, 4)
+function get_HOHQMesh_ids(::Tet, curved_faces, f_HOHQMesh, f, polydeg) 
+    num_face_nodes = (polydeg + 1)^2
+    return (1:num_face_nodes) .+ (f_HOHQMesh - 1) * num_face_nodes #active_face_offset[f_HOHQMesh]
+end
+
+function abtors(a, b)
+    s = b
+    r = @. (1 - s) * 0.5 * (a + 1) - 1
+    return r, s
+end
+
+# for tetrahedral meshes
+function MeshData(hmd::HOHQMeshData{3}, rd::RefElemData{3, <:Tet})
+    (; VXYZ, EToV) = hmd    
+    md = MeshData(VXYZ, EToV[:, vertex_reordering(rd.element_type)], rd)
+
+    # interpolation matrix from chebyshev_to_lobatto nodes
+    a_chebyshev = vec([-cos(j * pi / hmd.polydeg) for j in 0:hmd.polydeg, i in 0:hmd.polydeg])
+    b_chebyshev = vec([-cos(j * pi / hmd.polydeg) for i in 0:hmd.polydeg, j in 0:hmd.polydeg])
+    r_chebyshev, s_chebyshev = abtors(a_chebyshev, b_chebyshev)
+
+    chebyshev_to_lobatto = 
+        vandermonde(face_type(rd.element_type), hmd.polydeg, nodes(face_type(rd.element_type), rd.N)...) / vandermonde(face_type(rd.element_type), hmd.polydeg, r_chebyshev, s_chebyshev)
+
+    # permute rows of chebyshev_to_lobatto according to the permutation between StartUpDG rd.Fmask 
+    # ordering  and nodes(face_type(Tet())) ordering
+    Fmask = reshape(rd.Fmask, :, 4) # 4 = number of tet faces
+    rf, sf = nodes(face_type(rd.element_type), rd.N)
+    rFmask, sFmask = rd.r[Fmask[:, 4]], rd.s[Fmask[:, 4]] # use face 4 nodes for which t = -1
+    permuted_face_ids = match_coordinate_vectors((rf, sf), (rFmask, sFmask))    
+    chebyshev_to_lobatto = chebyshev_to_lobatto[permuted_face_ids, :]                
+
+    warp_face_nodes_to_volume_nodes = 
+        face_basis(rd.element_type, rd.N, rd.rst...) / face_basis(rd.element_type, rd.N, rd.r[rd.Fmask], rd.s[rd.Fmask], rd.t[rd.Fmask])            
+
+    HOHQMesh_to_StartUpDG_face_ordering = get_HOHQMesh_to_StartUpDG_face_ordering(rd.element_type)
+    curved_face_coordinates = ntuple(_ -> similar(reshape(md.x[rd.Fmask, 1], :, rd.num_faces)), 3)
+    x, y, z = copy.(md.xyz)
+    for curved_elem in hmd.curved_elements
+        (; element, curved_faces, curved_edge_coordinates) = curved_elem
+
+        # initialize face coordinates as linear coordinates
+        curved_face_coordinates[1] .= reshape(md.x[rd.Fmask, element], :, rd.num_faces)
+        curved_face_coordinates[2] .= reshape(md.y[rd.Fmask, element], :, rd.num_faces)
+        curved_face_coordinates[3] .= reshape(md.z[rd.Fmask, element], :, rd.num_faces)
+
+        for (f_HOHQMesh, f) in enumerate(HOHQMesh_to_StartUpDG_face_ordering)            
+            if curved_faces[f_HOHQMesh] == 1
+                # incorporate any permutations or face reorderings into indices
+                ids_HOHQMesh = get_HOHQMesh_ids(rd.element_type, curved_faces, f_HOHQMesh, f, hmd.polydeg)
+
+                curved_lobatto_coordinates = chebyshev_to_lobatto * curved_edge_coordinates[ids_HOHQMesh, :]
+                curved_face_coordinates[1][:, f] .= curved_lobatto_coordinates[:, 1]
+                curved_face_coordinates[2][:, f] .= curved_lobatto_coordinates[:, 2]
+                curved_face_coordinates[3][:, f] .= curved_lobatto_coordinates[:, 3]
+            end
+        end        
+
+        x[:, element] .= warp_face_nodes_to_volume_nodes * vec(curved_face_coordinates[1])
+        y[:, element] .= warp_face_nodes_to_volume_nodes * vec(curved_face_coordinates[2])
+        z[:, element] .= warp_face_nodes_to_volume_nodes * vec(curved_face_coordinates[3])
+    end
+
+    md_curved = MeshData(rd, md, x, y, z) 
+
+    # find boundary tags
+    boundary_tags = hmd.boundary_tags[:, 1:num_faces(rd.element_type)] # ignore `#undef` entries in last columns 5 and 6
+    boundary_strings = unique(boundary_tags)
+    deleteat!(boundary_strings, findfirst(boundary_strings .=== "---")) # remove non-boundary faces
+    face_indices = [map(x -> global_face_index(x, rd.num_faces, HOHQMesh_to_StartUpDG_face_ordering), 
+                        findall(boundary_tags .== name)) for name in boundary_strings]
+    boundary_faces = NamedTuple(Pair.(Symbol.(boundary_strings), face_indices))
+
+    return @set md_curved.mesh_type = HOHQMeshType(hmd, boundary_faces)
+end
+
 function read_HOHQMesh(filename::String)
     f = open(filename)
     lines = readlines(f)
@@ -185,7 +263,7 @@ end
 struct CurvedHOHQMeshElement
     element::Int
     curved_faces::Vector{Int}
-    curved_edge_coordinates::Matrix{Float64}
+    curved_edge_coordinates::Matrix{Float64} # TODO: change this to curved_face_coordinates for 2D/3D consistency 
 end
 
 # for hexes and quads

diff --git a/test/hohqmesh_tests.jl b/test/hohqmesh_tests.jl
@@ -40,8 +40,9 @@
     @test sum(md.wJq) ≈ 1.0
 
     # Tet meshes
-    filename = "testset_HOHQMesh_meshes/MSMappedTet4P4.mesh"
+    filename = "testset_HOHQMesh_meshes/TetMesh44.mesh"
     @test_nowarn hmd = read_HOHQMesh(filename, Tet())
-    # rd = RefElemData(Tet(), 4)
-    # md = MeshData(hmd, rd)
+    rd = RefElemData(Tet(), 4)
+    md = MeshData(hmd, rd)
+    @test all(md.J .> 0)
 end