Skip to content

Commit

Permalink
Merge branch 'dev' into jc/fix_RefElemData_dispatch
Browse files Browse the repository at this point in the history
  • Loading branch information
jlchan authored Dec 18, 2023
2 parents 61ff4d4 + 43235ec commit effa35f
Show file tree
Hide file tree
Showing 9 changed files with 186 additions and 109 deletions.
2 changes: 1 addition & 1 deletion Project.toml
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand Down
4 changes: 1 addition & 3 deletions src/mesh/gmsh_utilities.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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

Expand Down
17 changes: 8 additions & 9 deletions src/mesh/simple_meshes.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand All @@ -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:
Expand Down
5 changes: 2 additions & 3 deletions src/mesh/triangulate_utils.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand All @@ -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
Expand Down
90 changes: 2 additions & 88 deletions test/MeshData_tests.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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()
Expand Down Expand Up @@ -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
164 changes: 164 additions & 0 deletions test/MeshData_wedge_pyr_tests.jl
Original file line number Diff line number Diff line change
@@ -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
4 changes: 2 additions & 2 deletions test/reference_elem_tests.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
5 changes: 3 additions & 2 deletions test/runtests.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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")
Expand Down
Loading

0 comments on commit effa35f

Please sign in to comment.