diff --git a/src/pyr_element.jl b/src/pyr_element.jl index b21a363..4caa3bf 100644 --- a/src/pyr_element.jl +++ b/src/pyr_element.jl @@ -63,16 +63,30 @@ end """ abctorst(elem::Pyr,a,b,c) -Converts from Stroud coordinates (a,b,c) on [-1,1]^3 to reference element -coordinates (r,s,t). +Converts from Stroud coordinates (a,b,c) on [-1,1]^3 to reference +element coordinates (r,s,t). """ -function abctorst(elem::Pyr,a,b,c) - r = @. .5 * (1+a) * (1-c) - 1 - s = @. .5 * (1+b) * (1-c) - 1 +function abctorst(elem::Pyr, a, b, c) + r = @. .5 * (1 + a) * (1 - c) - 1 + s = @. .5 * (1 + b) * (1 - c) - 1 t = @. c return r, s, t end +""" + rsttoabc(elem::Pyr,a,b,c) + +Converts from reference element coordinates (r,s,t) to Stroud +coordinates (a,b,c) on [-1,1]^3. +""" +function rsttoabc(::Pyr, r, s, t) + a = @. 2 * (r + 1) / (1 - t) - 1 + b = @. 2 * (s + 1) / (1 - t) - 1 + c = @. t + return a, b, c +end + + """ nodes(elem::Pyr,N) diff --git a/src/tet_element.jl b/src/tet_element.jl index 49ed5fb..fca6067 100644 --- a/src/tet_element.jl +++ b/src/tet_element.jl @@ -125,7 +125,20 @@ function jaskowiec_sukumar_quad_nodes(elem::Tet, N) return rstw[:,1], rstw[:,2], rstw[:,3], rstw[:,4] end -function stroud_quad_nodes(elem::Tet,N) +function abctorst(::Tet, a, b, c) + r = @. .5 * (1+a) * .5 * (1-b) * (1-c) - 1 + s = @. .5 * (1+b) * (1-c) - 1 + t = copy(c) + return r, s, t +end + +function rsttoabc(::Tet, r, s, t) + a = @. 2 * (1+r) / (-s-t) - 1 + b = @. 2 * (1+s) / (1-t) - 1 + c = copy(t) + return a, b, c +end + function stroud_quad_nodes(elem::Tet, N) cubA,cubWA = gauss_quad(0,0,N) cubB,cubWB = gauss_quad(1,0,N) @@ -134,12 +147,10 @@ function stroud_quad_nodes(elem::Tet, N) a,b,c = vec.(meshgrid(cubA,cubB,cubC)) wa,wb,wc = vec.(meshgrid(cubWA,cubWB,cubWC)) - r = @. .5*(1+a)*.5*(1-b)*(1-c)-1 - s = @. .5*(1+b)*(1-c)-1 - t = c - w = @. wa*wb*wc - w = (4/3)*w./sum(w) - return r,s,t,w + r, s, t = abctorst(elem, a, b, c) + w = @. wa * wb * wc + w = (4/3) * w ./ sum(w) + return r, s, t, w end quad_nodes(elem::Tet, N) = quad_nodes_tet(2 * N) @@ -149,9 +160,7 @@ function basis(elem::Tet, N, r, s, t) V, Vr, Vs, Vt = ntuple(x -> zeros(length(r), Np), 4) - a = @. 2*(1+r)/(-s-t)-1 - b = @. 2*(1+s)/(1-t)-1 - c = t + a, b, c = rsttoabc(elem, r, s, t) tol = 1e-14 ida = @. abs(s+t) < tol diff --git a/src/tri_element.jl b/src/tri_element.jl index 807ea2d..fa0866e 100644 --- a/src/tri_element.jl +++ b/src/tri_element.jl @@ -63,7 +63,22 @@ function rstoab(r, s, tol = 1e-12) a[n] = -1 end end - return a, s + b = copy(s) + return a, b +end + +rstoab(::Tri, r, s, kwargs...) = rstoab(r, s, kwargs...) + +""" + abtors(Tri(), r, s, tol = 1e-12) + +Converts from polynomial basis evaluation coordinates (a,b) on the +domain [-1,1]^2 to reference bi-unit right triangle coordinate (r,s). +""" +function abtors(::Tri, a, b) + r = @. 0.5 * (a + 1) * (1 - b) - 1 + s = copy(b) + return r, s end """ @@ -92,8 +107,7 @@ function stroud_quad_nodes(::Tri, N) cubA, cubB = vec.(meshgrid(cubA, cubB)) - r = @. 0.5 * (1+cubA) * (1-cubB) - 1 - s = cubB + r, s = abtors(Tri(), cubA, cubB) w = 0.5 * cubWB * (cubWA') return vec.((r, s, w)) end diff --git a/test/runtests.jl b/test/runtests.jl index 2bf1593..9064434 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -72,6 +72,14 @@ area(elem::Tri) = 2.0 num_faces = (elem == Tri()) ? 3 : 4 @test length(NodesAndModes.face_vertices(elem)) == num_faces + if elem == Tri() + a, b = NodesAndModes.abtors(Tri(), r, s) + r2, s2 = NodesAndModes.rstoab(Tri(), a, b) + r3, s3 = NodesAndModes.rstoab(a, b) + @test r ≈ r2 ≈ r3 + @test s ≈ s2 ≈ s3 + end + # check for Kronecker structure if elem == Quad() r, s = nodes(elem, N) @@ -157,6 +165,12 @@ num_faces(::Hex) = 6 @test Dr * r ≈ one.(r) @test Ds * s ≈ one.(s) @test Dt * t ≈ one.(t) + + a, b, c = NodesAndModes.rsttoabc(Pyr(), rq, sq, tq) + r2, s2, t2 = NodesAndModes.abctorst(Pyr(), a, b, c) + @test r2 ≈ rq + @test s2 ≈ sq + @test t2 ≈ tq end end