From 43235ecbd396d087cbc5038f47f18f673008c075 Mon Sep 17 00:00:00 2001 From: Jesse Chan <1156048+jlchan@users.noreply.github.com> Date: Sat, 9 Dec 2023 21:07:16 -0600 Subject: [PATCH] Make triangle and tet node orderings consistent (#153) * simplifying test * fixing vertex orderings * splitting wedge-pyr MeshData tests out * fix triangulate tests * fix wedge mesh ordering * updating hardcoded values in VTK tests * update tests * bump compat for NodesAndModes to 1 --- Project.toml | 2 +- src/mesh/gmsh_utilities.jl | 4 +- src/mesh/simple_meshes.jl | 17 ++-- src/mesh/triangulate_utils.jl | 5 +- test/MeshData_tests.jl | 90 +---------------- test/MeshData_wedge_pyr_tests.jl | 164 +++++++++++++++++++++++++++++++ test/reference_elem_tests.jl | 4 +- test/runtests.jl | 5 +- test/write_vtk_tests.jl | 4 +- 9 files changed, 186 insertions(+), 109 deletions(-) create mode 100644 test/MeshData_wedge_pyr_tests.jl diff --git a/Project.toml b/Project.toml index 89ae6440..8030717b 100644 --- a/Project.toml +++ b/Project.toml @@ -34,7 +34,7 @@ ConstructionBase = "1" FillArrays = "0.13, 1" HDF5 = "0.16, 0.17" Kronecker = "0.5" -NodesAndModes = "0.9" +NodesAndModes = "1" PathIntersections = "0.1" RecipesBase = "1" RecursiveArrayTools = "2" diff --git a/src/mesh/gmsh_utilities.jl b/src/mesh/gmsh_utilities.jl index db958d7e..51a471ea 100644 --- a/src/mesh/gmsh_utilities.jl +++ b/src/mesh/gmsh_utilities.jl @@ -205,7 +205,6 @@ function read_Gmsh_2D_v4(filename::String, options::MeshImportOptions) block_line_start = block_line_start + num_elem_in_block + 1 end - EToV = EToV[:, vec([1 3 2])] # permute for Gmsh ordering EToV = correct_negative_Jacobians!((VX, VY), EToV) if grouping @@ -320,7 +319,6 @@ function read_Gmsh_2D_v2(filename::String) end end - EToV = EToV[:, vec([1 3 2])] # permute for Gmsh ordering EToV = correct_negative_Jacobians!((VX, VY), EToV) return (VX, VY), EToV @@ -331,7 +329,7 @@ end # Computes the area of a triangle given `tri`, which is a tuple of three points (vectors), # using the [Shoelace_formula](https://en.wikipedia.org/wiki/Shoelace_formula). function compute_triangle_area(tri) - B, A, C = tri + A, B, C = tri return 0.5 * (A[1] * (B[2] - C[2]) + B[1] * (C[2] - A[2]) + C[1] * (A[2] - B[2])) end diff --git a/src/mesh/simple_meshes.jl b/src/mesh/simple_meshes.jl index ffa92075..36aef802 100644 --- a/src/mesh/simple_meshes.jl +++ b/src/mesh/simple_meshes.jl @@ -29,8 +29,8 @@ function uniform_mesh(elem::Tri, Kx, Ky) id2 = id(ex + 1, ey) id3 = id(ex + 1, ey + 1) id4 = id(ex, ey + 1) - EToV[2*sk-1, :] = [id1 id3 id2] - EToV[2*sk, :] = [id3 id1 id4] + EToV[2*sk-1, :] = [id1 id2 id3] + EToV[2*sk, :] = [id3 id4 id1] sk += 1 end end @@ -123,15 +123,15 @@ function uniform_mesh(elem::Hex, Nx, Ny, Nz) end -function uniform_mesh(elem::Wedge, Kx, Ky, Kz) +function uniform_mesh(::Wedge, Kx, Ky, Kz) (VY, VX, VZ) = meshgrid(LinRange(-1, 1, Ky + 1), LinRange(-1, 1, Kx + 1), LinRange(-1, 1, Kz + 1)) sk = 1 - shift = (Kx+1)*(Ky+1) + shift = (Kx + 1) * (Ky + 1) + id(ex, ey, ez) = ex + (ey - 1) * (Kx + 1) + (ez - 1) * (shift) # index function EToV = zeros(Int, 2 * Kx * Ky * Kz, 6) for ey = 1:Ky for ex = 1:Kx for ez = 1:Kz - id(ex, ey, ez) = ex + (ey - 1) * (Kx + 1) + (ez - 1) * (shift)# index function id1 = id(ex, ey, ez) id2 = id(ex + 1, ey, ez) id3 = id(ex, ey + 1, ez) @@ -140,14 +140,13 @@ function uniform_mesh(elem::Wedge, Kx, Ky, Kz) id6 = id(ex + 1, ey, ez + 1) id7 = id(ex, ey + 1, ez + 1) id8 = id(ex + 1, ey + 1, ez + 1) - EToV[2*sk-1, :] = [id1 id5 id4 id8 id2 id6] - EToV[2*sk, :] = [id1 id5 id3 id7 id4 id8] + EToV[2 * sk - 1, :] = [id1 id2 id4 id5 id6 id8] + EToV[2 * sk, :] = [id1 id4 id3 id5 id8 id7] sk += 1 end end end - EToV .= EToV[:, SVector(1, 3, 5, 2, 4, 6)] - return (VX[:], VY[:], VZ[:]), EToV + return vec.((VX, VY, VZ)), EToV end # split each cube into 6 pyramids. Pyramid 1: diff --git a/src/mesh/triangulate_utils.jl b/src/mesh/triangulate_utils.jl index f189af7f..918c16e9 100644 --- a/src/mesh/triangulate_utils.jl +++ b/src/mesh/triangulate_utils.jl @@ -19,7 +19,6 @@ Computes `VX`,`VY`,`EToV` from a `TriangulateIO` object. function triangulateIO_to_VXYEToV(triout::TriangulateIO) VX,VY = (triout.pointlist[i,:] for i = 1:size(triout.pointlist,1)) EToV = permutedims(triout.trianglelist) - Base.swapcols!(EToV,2,3) # to match MeshData ordering return (VX, VY), Matrix{Int}(EToV) end @@ -37,8 +36,8 @@ function get_boundary_face_labels(triout::TriangulateIO, rd::RefElemData{2, Tri} for (f,boundary_face) in enumerate(boundary_faces) element = (boundary_face - 1) ÷ rd.Nfaces + 1 face = (boundary_face - 1) % rd.Nfaces + 1 - vertices_on_face = sort(md.EToV[element,rd.fv[face]]) - tag_id = findfirst(c->view(segmentlist,:,c)==vertices_on_face,axes(segmentlist,2)) + vertices_on_face = sort(md.EToV[element, rd.fv[face]]) + tag_id = findfirst(c -> view(segmentlist,:,c) == vertices_on_face,axes(segmentlist, 2)) boundary_face_tags[f] = triout.segmentmarkerlist[tag_id] end return boundary_face_tags, boundary_faces diff --git a/test/MeshData_tests.jl b/test/MeshData_tests.jl index 39f86075..d939f90d 100644 --- a/test/MeshData_tests.jl +++ b/test/MeshData_tests.jl @@ -208,11 +208,11 @@ end end +# Note: we test wedge and pyr elements in a separate file approx_elem_types_to_test = [(Polynomial(), Hex()), (SBP(), Hex()), (Polynomial(), Tet()), - (Polynomial(), Wedge()), - (Polynomial(), Pyr())] + ] @testset "3D MeshData tests" begin @testset "$approximation_type $element_type MeshData initialization" for (approximation_type, element_type) in approx_elem_types_to_test tol = 5e2*eps() @@ -288,90 +288,4 @@ approx_elem_types_to_test = [(Polynomial(), Hex()), @test isempty(md.mapB) end - - @testset "TensorProductWedge MeshData" begin - element_type = Wedge() - tol = 5e2*eps() - @testset "Degree $tri_grad triangle" for tri_grad = [2, 3] - @testset "Degree $line_grad line" for line_grad = [2, 3] - line = RefElemData(Line(), line_grad) - tri = RefElemData(Tri(), tri_grad) - tensor = TensorProductWedge(tri, line) - - rd = RefElemData(element_type, tensor) - K1D = 2 - md = MeshData(uniform_mesh(element_type, K1D)..., rd) - (; wq, Dr, Ds, Dt, Vq, Vf, wf ) = rd - Nfaces = length(rd.fv) - (; x, y, z, xq, yq, zq, wJq, xf, yf, zf, K ) = md - (; rxJ, sxJ, txJ, ryJ, syJ, tyJ, rzJ, szJ, tzJ, J ) = md - (; nxJ, nyJ, nzJ, sJ ) = md - (; FToF, mapM, mapP, mapB ) = md - - @test StartUpDG._short_typeof(rd.approximation_type) == "TensorProductWedge{Polynomial, Polynomial}" - @test typeof(md.mesh_type) <: StartUpDG.VertexMappedMesh{<:typeof(rd.element_type)} - @test md.x == md.xyz[1] - @test md.y == md.xyz[2] - @test md.z == md.xyz[3] - - # check positivity of Jacobian - # @show J[1,:] - @test all(J .> 0) - h = estimate_h(rd, md) - @test h <= 2 / K1D + tol - - # check differentiation - u = @. x^2 + 2 * x * y - y^2 + x * y * z - dudx_exact = @. (2*x + 2*y + y*z) - dudy_exact = @. 2*x - 2*y + x*z - dudz_exact = @. x*y - dudr,duds,dudt = (D->D*u).((Dr, Ds, Dt)) - dudx = @. (rxJ * dudr + sxJ * duds + txJ * dudt) / J - dudy = @. (ryJ * dudr + syJ * duds + tyJ * dudt) / J - dudz = @. (rzJ * dudr + szJ * duds + tzJ * dudt) / J - @test dudx ≈ dudx_exact - @test dudy ≈ dudy_exact - @test dudz ≈ dudz_exact - - # check volume integration - @test Vq * x ≈ xq - @test Vq * y ≈ yq - @test Vq * z ≈ zq - @test diagm(wq) * (Vq * J) ≈ wJq - @test abs(sum(xq .* wJq)) < tol - @test abs(sum(yq .* wJq)) < tol - @test abs(sum(zq .* wJq)) < tol - - # check surface integration - @test Vf * x ≈ xf - @test Vf * y ≈ yf - @test Vf * z ≈ zf - @test abs(sum(diagm(wf) * nxJ)) < tol - @test abs(sum(diagm(wf) * nyJ)) < tol - @test abs(sum(diagm(wf) * nzJ)) < tol - @test md.nx .* md.Jf ≈ md.nxJ - @test md.ny .* md.Jf ≈ md.nyJ - @test md.nz .* md.Jf ≈ md.nzJ - - # check connectivity and boundary maps - u = @. (1-x) * (1+x) * (1-y) * (1+y) * (1-z) * (1+z) - uf = Vf * u - @test uf ≈ uf[mapP] - @test norm(uf[mapB]) < tol - - # check periodic node connectivity maps - md_periodic = make_periodic(md, (true, true, true)) - @test md_periodic.mapP != md.mapP # check that the node mapping actually changed - - u = @. sin(pi * (.5 + x)) * sin(pi * (.5 + y)) * sin(pi * (.5 + z)) - (; mapP ) = md_periodic - uf = Vf * u - @test uf ≈ uf[mapP] - - md = MeshData(uniform_mesh(rd.element_type, K1D)..., rd; is_periodic=true) - @test isempty(md.mapB) - - end - end - end end diff --git a/test/MeshData_wedge_pyr_tests.jl b/test/MeshData_wedge_pyr_tests.jl new file mode 100644 index 00000000..e3a139c6 --- /dev/null +++ b/test/MeshData_wedge_pyr_tests.jl @@ -0,0 +1,164 @@ +approx_elem_types_to_test = [(Polynomial(), Wedge()), + (Polynomial(), Pyr())] + +@testset "3D MeshData tests for wedges and pyramids" begin + @testset "$approximation_type $element_type MeshData initialization" for (approximation_type, element_type) in approx_elem_types_to_test + tol = 5e2*eps() + + N = 3 + K1D = 2 + rd = RefElemData(element_type, approximation_type, N) + md = MeshData(uniform_mesh(element_type, K1D)..., rd) + (; wq, Dr, Ds, Dt, Vq, Vf, wf ) = rd + Nfaces = length(rd.fv) + (; x, y, z, xq, yq, zq, wJq, xf, yf, zf, K ) = md + (; rxJ, sxJ, txJ, ryJ, syJ, tyJ, rzJ, szJ, tzJ, J ) = md + (; nxJ, nyJ, nzJ, sJ ) = md + (; FToF, mapM, mapP, mapB ) = md + + @test typeof(md.mesh_type) <: StartUpDG.VertexMappedMesh{<:typeof(rd.element_type)} + @test md.x == md.xyz[1] + + # check positivity of Jacobian + @test all(J .> 0) + h = estimate_h(rd, md) + @test h <= 2 / K1D + tol + + # check differentiation + u = @. x^2 + 2 * x * y - y^2 + x * y * z + dudx_exact = @. 2*x + 2*y + y*z + dudy_exact = @. 2*x - 2*y + x*z + dudz_exact = @. x*y + dudr,duds,dudt = (D->D*u).((Dr, Ds, Dt)) + dudx = @. (rxJ * dudr + sxJ * duds + txJ * dudt) / J + dudy = @. (ryJ * dudr + syJ * duds + tyJ * dudt) / J + dudz = @. (rzJ * dudr + szJ * duds + tzJ * dudt) / J + @test dudx ≈ dudx_exact + @test dudy ≈ dudy_exact + @test dudz ≈ dudz_exact + + # check volume integration + @test Vq * x ≈ xq + @test Vq * y ≈ yq + @test Vq * z ≈ zq + @test diagm(wq) * (Vq * J) ≈ wJq + @test abs(sum(xq .* wJq)) < tol + @test abs(sum(yq .* wJq)) < tol + @test abs(sum(zq .* wJq)) < tol + + # check surface integration + @test Vf * x ≈ xf + @test Vf * y ≈ yf + @test Vf * z ≈ zf + @test abs(sum(diagm(wf) * nxJ)) < tol + @test abs(sum(diagm(wf) * nyJ)) < tol + @test abs(sum(diagm(wf) * nzJ)) < tol + @test md.nx .* md.Jf ≈ md.nxJ + @test md.ny .* md.Jf ≈ md.nyJ + @test md.nz .* md.Jf ≈ md.nzJ + + # check connectivity and boundary maps + u = @. (1-x) * (1+x) * (1-y) * (1+y) * (1-z) * (1+z) + uf = Vf * u + @test uf ≈ uf[mapP] + @test norm(uf[mapB]) < tol + + # check periodic node connectivity maps + md_periodic = make_periodic(md, (true, true, true)) + @test md_periodic.mapP != md.mapP # check that the node mapping actually changed + + u = @. sin(pi * (.5 + x)) * sin(pi * (.5 + y)) * sin(pi * (.5 + z)) + (; mapP ) = md_periodic + uf = Vf * u + @test uf ≈ uf[mapP] + + md = MeshData(uniform_mesh(rd.element_type, K1D)..., rd; is_periodic=true) + @test isempty(md.mapB) + + end + + @testset "TensorProductWedge MeshData" begin + element_type = Wedge() + tol = 5e2*eps() + @testset "Degree $tri_grad triangle" for tri_grad = [2, 3] + @testset "Degree $line_grad line" for line_grad = [2, 3] + line = RefElemData(Line(), line_grad) + tri = RefElemData(Tri(), tri_grad) + tensor = TensorProductWedge(tri, line) + + rd = RefElemData(element_type, tensor) + K1D = 2 + md = MeshData(uniform_mesh(element_type, K1D)..., rd) + (; wq, Dr, Ds, Dt, Vq, Vf, wf ) = rd + Nfaces = length(rd.fv) + (; x, y, z, xq, yq, zq, wJq, xf, yf, zf, K ) = md + (; rxJ, sxJ, txJ, ryJ, syJ, tyJ, rzJ, szJ, tzJ, J ) = md + (; nxJ, nyJ, nzJ, sJ ) = md + (; FToF, mapM, mapP, mapB ) = md + + @test StartUpDG._short_typeof(rd.approximation_type) == "TensorProductWedge{Polynomial, Polynomial}" + @test typeof(md.mesh_type) <: StartUpDG.VertexMappedMesh{<:typeof(rd.element_type)} + @test md.x == md.xyz[1] + @test md.y == md.xyz[2] + @test md.z == md.xyz[3] + + # check positivity of Jacobian + @test all(J .> 0) + h = estimate_h(rd, md) + @test h <= 2 / K1D + tol + + # check differentiation + u = @. x^2 + 2 * x * y - y^2 + x * y * z + dudx_exact = @. (2*x + 2*y + y*z) + dudy_exact = @. 2*x - 2*y + x*z + dudz_exact = @. x*y + dudr,duds,dudt = (D->D*u).((Dr, Ds, Dt)) + dudx = @. (rxJ * dudr + sxJ * duds + txJ * dudt) / J + dudy = @. (ryJ * dudr + syJ * duds + tyJ * dudt) / J + dudz = @. (rzJ * dudr + szJ * duds + tzJ * dudt) / J + @test dudx ≈ dudx_exact + @test dudy ≈ dudy_exact + @test dudz ≈ dudz_exact + + # check volume integration + @test Vq * x ≈ xq + @test Vq * y ≈ yq + @test Vq * z ≈ zq + @test diagm(wq) * (Vq * J) ≈ wJq + @test abs(sum(xq .* wJq)) < tol + @test abs(sum(yq .* wJq)) < tol + @test abs(sum(zq .* wJq)) < tol + + # check surface integration + @test Vf * x ≈ xf + @test Vf * y ≈ yf + @test Vf * z ≈ zf + @test abs(sum(diagm(wf) * nxJ)) < tol + @test abs(sum(diagm(wf) * nyJ)) < tol + @test abs(sum(diagm(wf) * nzJ)) < tol + @test md.nx .* md.Jf ≈ md.nxJ + @test md.ny .* md.Jf ≈ md.nyJ + @test md.nz .* md.Jf ≈ md.nzJ + + # check connectivity and boundary maps + u = @. (1-x) * (1+x) * (1-y) * (1+y) * (1-z) * (1+z) + uf = Vf * u + @test uf ≈ uf[mapP] + @test norm(uf[mapB]) < tol + + # check periodic node connectivity maps + md_periodic = make_periodic(md, (true, true, true)) + @test md_periodic.mapP != md.mapP # check that the node mapping actually changed + + u = @. sin(pi * (.5 + x)) * sin(pi * (.5 + y)) * sin(pi * (.5 + z)) + (; mapP ) = md_periodic + uf = Vf * u + @test uf ≈ uf[mapP] + + md = MeshData(uniform_mesh(rd.element_type, K1D)..., rd; is_periodic=true) + @test isempty(md.mapB) + + end + end + end +end diff --git a/test/reference_elem_tests.jl b/test/reference_elem_tests.jl index 7d487087..9349dd90 100644 --- a/test/reference_elem_tests.jl +++ b/test/reference_elem_tests.jl @@ -24,8 +24,8 @@ @test abs(sum(rd.wf)) ≈ 6 @test abs(sum(rd.wf .* rd.nrJ)) + abs(sum(rd.wf .* rd.nsJ)) < tol @test rd.Pq * rd.Vq ≈ I - Vfp = vandermonde(Line(), N, quad_nodes(Line(), N)[1]) / vandermonde(Line(), N, nodes(Line(), N)) - rstf = (x->Vfp * x[reshape(rd.Fmask, rd.Nfq ÷ rd.Nfaces, rd.Nfaces)]).(rd.rst) + Vf = vandermonde(Line(), N, quad_nodes(Line(), N)[1]) / vandermonde(Line(), N, nodes(Line(), N)) + rstf = (x->Vf * x[reshape(rd.Fmask, :, rd.Nfaces)]).(rd.rst) @test all(vec.(rstf) .≈ rd.rstf) @test StartUpDG.eigenvalue_inverse_trace_constant(rd) ≈ inverse_trace_constant(rd) @test propertynames(rd)[1] == :element_type diff --git a/test/runtests.jl b/test/runtests.jl index af4b0d34..60303f2c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -8,13 +8,14 @@ using SummationByPartsOperators using StartUpDG -include("write_vtk_tests.jl") +include("write_vtk_tests.jl") include("named_array_partition_tests.jl") -include("triangulate_tests.jl") +include("triangulate_tests.jl") include("reference_elem_tests.jl") include("multidim_sbp_tests.jl") include("SummationByPartsOperatorsExt_tests.jl") include("MeshData_tests.jl") +include("MeshData_wedge_pyr_tests.jl") include("boundary_util_tests.jl") include("hybrid_mesh_tests.jl") include("noncon_mesh_tests.jl") diff --git a/test/write_vtk_tests.jl b/test/write_vtk_tests.jl index e2bcdde8..d52eb3ed 100644 --- a/test/write_vtk_tests.jl +++ b/test/write_vtk_tests.jl @@ -11,7 +11,9 @@ deg_one_order(::Hex) = [-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 -1.0 -1.0 1.0 1.0 1.0 1.0] - 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_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)