Skip to content

Commit

Permalink
formatting/updating tests
Browse files Browse the repository at this point in the history
  • Loading branch information
@jlchan
jlchan committed Jun 24, 2021
1 parent 763a4d6 commit c6155d4
Show file tree
Hide file tree
Showing 5 changed files with 181 additions and 121 deletions.
2 changes: 1 addition & 1 deletion test/MeshData_tests.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -262,7 +262,7 @@ end
@test md.x == md.xyz[1]

# check positivity of Jacobian
# @show J[1,:]

@test all(J .> 0)
# h = estimate_h(rd,md)
# @test h ≈ 2/K1D
Expand Down
4 changes: 2 additions & 2 deletions test/geometric_term_tests.jl
View file
Original file line number Diff line number Diff line change
@@ -1,9 +1,9 @@
@testset "Geometric terms for $elem elements" for elem in [Tri() Quad() Hex()]
@testset "Geometric terms for $elem elements" for elem in [Tri() Quad() Hex() Tet()]
tol = 5e2*eps()
N = 3
rd = RefElemData(elem,N)
geofacs = geometric_factors(rd.rst...,rd.Drst...)
if elem != Hex()
if ndims(elem)==2
rxJ,sxJ,ryJ,syJ,J = geofacs
@test all(rxJ .≈ 1)
@test norm(sxJ) < tol
Expand Down
176 changes: 112 additions & 64 deletions test/reference_elem_tests.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,72 +3,120 @@
tol = 5e2*eps()

N = 2
@testset "Interval" begin
rd = RefElemData(Line(),N)
@test rd.r == rd.rst[1]
@test rd.Np == length(rd.r)
@test abs(sum(rd.rq.*rd.wq)) < tol
@test rd.nrJ ≈ [-1,1]
@test rd.Pq*rd.Vq ≈ I
@test rd.r[rd.Fmask[:]] ≈ rd.rf
@test invoke(inverse_trace_constant,Tuple{RefElemData},rd) ≈ inverse_trace_constant(rd)
@test propertynames(rd)[1] == :elementType
end

#####
##### interval
#####
rd = RefElemData(Line(),N)
@test rd.r == rd.rst[1]
@test rd.Np == length(rd.r)
@test abs(sum(rd.rq.*rd.wq)) < tol
@test rd.nrJ ≈ [-1,1]
@test rd.Pq*rd.Vq ≈ I
@test rd.r[rd.Fmask[:]] ≈ rd.rf
@test invoke(inverse_trace_constant,Tuple{RefElemData},rd) ≈ inverse_trace_constant(rd)
@test propertynames(rd)[1] == :elementType
@testset "Triangle" begin
rd = RefElemData(Tri(),N)
@test rd.r == rd.rst[1]
@test rd.Np == length(rd.r)
@test rd.Nq == length(rd.rq)
@test abs(sum(rd.wq)) ≈ 2
@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)
@test all(vec.(rstf) .≈ rd.rstf)
@test invoke(inverse_trace_constant,Tuple{RefElemData},rd) ≈ inverse_trace_constant(rd)
@test propertynames(rd)[1] == :elementType
end

#####
##### triangles
#####
rd = RefElemData(Tri(),N)
@test rd.r == rd.rst[1]
@test rd.Np == length(rd.r)
@test rd.Nq == length(rd.rq)
@test abs(sum(rd.wq)) ≈ 2
@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)
@test all(vec.(rstf) .≈ rd.rstf)
@test invoke(inverse_trace_constant,Tuple{RefElemData},rd) ≈ inverse_trace_constant(rd)
@test propertynames(rd)[1] == :elementType
@testset "Quad" begin
rd = RefElemData(Quad(),N)
@test rd.r == rd.rst[1]
@test rd.Np == length(rd.r)
@test rd.Nq == length(rd.rq)
@test abs(sum(rd.wq)) ≈ 4
@test abs(sum(rd.wf)) ≈ 8
@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)
@test all(vec.(rstf) .≈ rd.rstf)
@test invoke(inverse_trace_constant,Tuple{RefElemData},rd) ≈ inverse_trace_constant(rd)
end

#####
##### quads
#####
rd = RefElemData(Quad(),N)
@test rd.r == rd.rst[1]
@test rd.Np == length(rd.r)
@test rd.Nq == length(rd.rq)
@test abs(sum(rd.wq)) ≈ 4
@test abs(sum(rd.wf)) ≈ 8
@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)
@test all(vec.(rstf) .≈ rd.rstf)
@test invoke(inverse_trace_constant,Tuple{RefElemData},rd) ≈ inverse_trace_constant(rd)
@testset "DGSEM quad" begin
rd = RefElemData(Quad(),SBP(),N)
@test rd.r == rd.rst[1]
@test rd.Np == length(rd.r)
@test rd.Nq == length(rd.rq)
@test abs(sum(rd.wq)) ≈ 4
@test abs(sum(rd.wf)) ≈ 8
@test abs(sum(rd.wf .* rd.nrJ)) + abs(sum(rd.wf .* rd.nsJ)) < tol
@test rd.Pq*rd.Vq ≈ I
@test all(vec.(rd.rstf) .≈ (x->getindex(x,rd.Fmask)).(rd.rst))
@test invoke(inverse_trace_constant,Tuple{RefElemData},rd) ≈ inverse_trace_constant(rd)
end

#####
##### hexes
#####
rd = RefElemData(Hex(),N)
@test propertynames(rd)[1] == :elementType
@test rd.t == rd.rst[3]
@test rd.tf == rd.rstf[3]
@test rd.tq == rd.rstq[3]
@test rd.rp == rd.rstp[1]
@test rd.sp == rd.rstp[2]
@test rd.tp == rd.rstp[3]
@test rd.Np == length(rd.r)
@test rd.Nq == length(rd.rq)
@test abs(sum(rd.wq)) ≈ 8
@test abs(sum(rd.wf)) ≈ 6*4
@test abs(sum(rd.wf .* rd.nrJ)) < tol
@test abs(sum(rd.wf .* rd.nsJ)) < tol
@test abs(sum(rd.wf .* rd.ntJ)) < tol
@test rd.Pq*rd.Vq ≈ I
@test invoke(inverse_trace_constant,Tuple{RefElemData},rd) ≈ inverse_trace_constant(rd)
# TODO: test interpolation of Fmask matches rd.rstf.
@testset "Hex" begin
rd = RefElemData(Hex(),N)
@test propertynames(rd)[1] == :elementType
@test rd.t == rd.rst[3]
@test rd.tf == rd.rstf[3]
@test rd.tq == rd.rstq[3]
@test rd.rp == rd.rstp[1]
@test rd.sp == rd.rstp[2]
@test rd.tp == rd.rstp[3]
@test rd.Np == length(rd.r)
@test rd.Nq == length(rd.rq)
@test abs(sum(rd.wq)) ≈ 8
@test abs(sum(rd.wf)) ≈ 6*4
@test abs(sum(rd.wf .* rd.nrJ)) < tol
@test abs(sum(rd.wf .* rd.nsJ)) < tol
@test abs(sum(rd.wf .* rd.ntJ)) < tol
@test rd.Pq*rd.Vq ≈ I
@test invoke(inverse_trace_constant,Tuple{RefElemData},rd) ≈ inverse_trace_constant(rd)
# TODO: test interpolation of Fmask matches rd.rstf.
end

@testset "DGSEM Hex" begin
rd = RefElemData(Hex(),SBP(),N)
@test propertynames(rd)[1] == :elementType
@test rd.t == rd.rst[3]
@test rd.tf == rd.rstf[3]
@test rd.tq == rd.rstq[3]
@test rd.rp == rd.rstp[1]
@test rd.sp == rd.rstp[2]
@test rd.tp == rd.rstp[3]
@test rd.Np == length(rd.r)
@test rd.Nq == length(rd.rq)
@test abs(sum(rd.wq)) ≈ 8
@test abs(sum(rd.wf)) ≈ 6*4
@test abs(sum(rd.wf .* rd.nrJ)) < tol
@test abs(sum(rd.wf .* rd.nsJ)) < tol
@test abs(sum(rd.wf .* rd.ntJ)) < tol
@test rd.Pq*rd.Vq ≈ I
@test invoke(inverse_trace_constant,Tuple{RefElemData},rd) ≈ inverse_trace_constant(rd)
end

@testset "Tet" begin
rd = RefElemData(Tet(),N)
@test propertynames(rd)[1] == :elementType
@test rd.t == rd.rst[3]
@test rd.tf == rd.rstf[3]
@test rd.tq == rd.rstq[3]
@test rd.rp == rd.rstp[1]
@test rd.sp == rd.rstp[2]
@test rd.tp == rd.rstp[3]
@test rd.Np == length(rd.r)
@test rd.Nq == length(rd.rq)
@test abs(sum(rd.wq)) ≈ 4/3
@test abs(sum(rd.wf)) ≈ 8
@test abs(sum(rd.wf .* rd.nrJ)) < tol
@test abs(sum(rd.wf .* rd.nsJ)) < tol
@test abs(sum(rd.wf .* rd.ntJ)) < tol
@test rd.Pq*rd.Vq ≈ I
@test invoke(inverse_trace_constant,Tuple{RefElemData},rd) ≈ inverse_trace_constant(rd)
end
end
43 changes: 25 additions & 18 deletions test/sbp_tests.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,24 +3,31 @@
tol = 200*eps()

N = 3
f(N,r,s) = r^N + s^N # function for testing

@testset "Hicken's 2N triangular SBP operators" begin
rd = RefElemData(Tri(),SBP(), N, quadrature_strength=2*N)
@unpack rq,sq,wq = rd
rq2,sq2,wq2 = quad_nodes(Tri(),2*N)
@test sum(wq.*f.(2*N,rq,sq)) ≈ sum(wq2.*f.(2*N,rq2,sq2))
@test rd.Dr*rd.r.^N ≈ N*rd.r.^(N-1)
@test rd.Ds*rd.s.^N ≈ N*rd.s.^(N-1)
@test norm(rd.Dr*rd.s + rd.Ds*rd.r) < tol
end

rd = RefElemData(Tri(),SBP(),N,quadrature_strength=2*N)
@unpack rq,sq,wq = rd
rq2,sq2,wq2 = quad_nodes(Tri(),2*N)
f(N,r,s) = r^N + s^N
@test sum(wq.*f.(2*N,rq,sq)) ≈ sum(wq2.*f.(2*N,rq2,sq2))
@test rd.Dr*rd.r.^N ≈ N*rd.r.^(N-1)
@test rd.Ds*rd.s.^N ≈ N*rd.s.^(N-1)
@test norm(rd.Dr*rd.s + rd.Ds*rd.r) < tol
@testset "Kubatko's Legendre face node triangular SBP operators" begin
rd = RefElemData(Tri(), SBP(), N, quadrature_strength=2*N-1, quad_rule_face=:Legendre)
@unpack rq,sq,wq = rd
rq2,sq2,wq2 = quad_nodes(Tri(),2*N-1)
@test sum(wq.*f.(2*N-1,rq,sq)) ≈ sum(wq2.*f.(2*N-1,rq2,sq2))
@test rd.Dr*rd.r.^N ≈ N*rd.r.^(N-1)
@test rd.Ds*rd.s.^N ≈ N*rd.s.^(N-1)
@test norm(rd.Dr*rd.s + rd.Ds*rd.r) < tol
end

@testset "Warning for N=6 Kubatko Lobatto SBP nodes" begin
@test_logs (:warn,"N=6 SBP operators with quadrature strength 2N-1 and Lobatto face nodes may require very small timesteps.") RefElemData(Tri(),SBP(),6,quadrature_strength=11)
@test_throws SystemError RefElemData(Tri(), SBP(), 7,quadrature_strength=14)
end

@test_logs (:warn,"N=6 SBP operators with quadrature strength 2N-1 and Lobatto face nodes may require very small timesteps.") RefElemData(Tri(),SBP(),6,quadrature_strength=11)
@test_throws SystemError RefElemData(Tri(),SBP(),7,quadrature_strength=14)

rd = RefElemData(Tri(), SBP(), N, quadrature_strength=2*N-1, quad_rule_face=:Legendre)
@unpack rq,sq,wq = rd
rq2,sq2,wq2 = quad_nodes(Tri(),2*N-1)
@test sum(wq.*f.(2*N-1,rq,sq)) ≈ sum(wq2.*f.(2*N-1,rq2,sq2))
@test rd.Dr*rd.r.^N ≈ N*rd.r.^(N-1)
@test rd.Ds*rd.s.^N ≈ N*rd.s.^(N-1)
@test norm(rd.Dr*rd.s + rd.Ds*rd.r) < tol
end
77 changes: 41 additions & 36 deletions test/triangulate_tests.jl
View file
Original file line number Diff line number Diff line change
@@ -1,44 +1,49 @@
@testset "Mesh, timestep, and Triangulate utils" begin
tol = 5e2*eps()
@testset "Timestep utils" begin
rk4a,rk4b,rk4c = ck45()
@test rk4c[1] ≈ 0.0 && rk4c[end] ≈ 1.0
@test rk4a[1] ≈ 0.0
@test all(rk4b .> 0)
end

rk4a,rk4b,rk4c = ck45()
@test rk4c[1] ≈ 0.0 && rk4c[end] ≈ 1.0
@test rk4a[1] ≈ 0.0
@test all(rk4b .> 0)
@testset "Gmsh reading" begin
VX,VY,EToV = readGmsh2D("squareCylinder2D.msh")
@test size(EToV)==(3031,3)
end

VX,VY,EToV = readGmsh2D("squareCylinder2D.msh")
@test size(EToV)==(3031,3)
@testset "Triangulate utils/example meshes" begin
# test triangulate
meshIO = square_domain()
VX,VY,EToV = triangulateIO_to_VXYEToV(meshIO)
rd = RefElemData(Tri(),2)
md = MeshData(VX,VY,EToV,rd)
@test size(EToV,1)==md.num_elements==620
@test length(VX)==length(VY)==338
@test sort(unique(get_node_boundary_tags(meshIO,rd,md)))==[0,1,2,3,4]

# test triangulate
meshIO = square_domain()
VX,VY,EToV = triangulateIO_to_VXYEToV(meshIO)
rd = RefElemData(Tri(),2)
md = MeshData(VX,VY,EToV,rd)
@test size(EToV,1)==md.num_elements==620
@test length(VX)==length(VY)==338
@test sort(unique(get_node_boundary_tags(meshIO,rd,md)))==[0,1,2,3,4]
h = .1
meshIO = square_hole_domain(h)
VX,VY,EToV = triangulateIO_to_VXYEToV(meshIO)
rd = RefElemData(Tri(),2)
md = MeshData(VX,VY,EToV,rd)
@test size(EToV,1)==md.num_elements==598
@test length(VX)==length(VY)==327
@test sort(unique(get_node_boundary_tags(meshIO,rd,md)))==[0,1,2]
meshIO2 = refine(meshIO,h)
@test sort(unique(get_node_boundary_tags(meshIO2,rd,MeshData(triangulateIO_to_VXYEToV(meshIO2)...,rd))))==[0,1,2]

h = .1
meshIO = square_hole_domain(h)
VX,VY,EToV = triangulateIO_to_VXYEToV(meshIO)
rd = RefElemData(Tri(),2)
md = MeshData(VX,VY,EToV,rd)
@test size(EToV,1)==md.num_elements==598
@test length(VX)==length(VY)==327
@test sort(unique(get_node_boundary_tags(meshIO,rd,md)))==[0,1,2]
meshIO2 = refine(meshIO,h)
@test sort(unique(get_node_boundary_tags(meshIO2,rd,MeshData(triangulateIO_to_VXYEToV(meshIO2)...,rd))))==[0,1,2]
meshIO = scramjet()
VX,VY,EToV = triangulateIO_to_VXYEToV(meshIO)
rd = RefElemData(Tri(),2)
md = MeshData(VX,VY,EToV,rd)
@test size(EToV,1)==md.num_elements==1550
@test length(VX)==length(VY)==871
@test sort(unique(get_node_boundary_tags(meshIO,rd,md)))==[0,1,2,3]

meshIO = scramjet()
VX,VY,EToV = triangulateIO_to_VXYEToV(meshIO)
rd = RefElemData(Tri(),2)
md = MeshData(VX,VY,EToV,rd)
@test size(EToV,1)==md.num_elements==1550
@test length(VX)==length(VY)==871
@test sort(unique(get_node_boundary_tags(meshIO,rd,md)))==[0,1,2,3]

boundary_faces = tag_boundary_faces(meshIO, rd, md, Dict(:wall=>1, :inflow=>2, :outflow=>3))
@test length(boundary_faces[:wall])==173
@test length(boundary_faces[:inflow])==14
@test length(boundary_faces[:outflow])==5
boundary_faces = tag_boundary_faces(meshIO, rd, md, Dict(:wall=>1, :inflow=>2, :outflow=>3))
@test length(boundary_faces[:wall])==173
@test length(boundary_faces[:inflow])==14
@test length(boundary_faces[:outflow])==5
end
end

0 comments on commit c6155d4

Please sign in to comment.