Merge pull request #15 from jlchan/jc/abctorst

add `abctorst` and `rsttoabc` for Tri, Tet, and Pyr elements

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)
 

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

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

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