Merge pull request #10 from sjkelly/sjk/dihedral

Sjk/dihedral

Author: sjkelly

src/quality.jl

@@ -8,9 +8,9 @@ function triqual(p1, p2, p3)
     d13 = p3 - p1
     d23 = p3 - p2
     n   = cross(d12,d13)
-    vol = sqrt(sum(n.^2))/2
+    vol = sqrt(sum(n.^2))
     den = dot(d12,d12) + dot(d13,d13) + dot(d23,d23)
-    return sqrt(3)*4*vol/den
+    return sqrt(3)*2*vol/den
 end
 
 function triangle_qualities(p,tets)
@@ -32,4 +32,56 @@ end
 
 function triangle_qualities!(tris::Vector,qualities::Vector,p,tets)
     triangle_qualities!(tris,Set{eltype(tris)}(),qualities,p,tets)
-end
+end
+
+const ⋅ = dot
+const × = cross
+
+function dihedral(p0,p1,p2,p3)\
+    b1 = p1 - p0
+    b2 = p2 - p1
+    b3 = p3 - p2
+
+    abs(atan(((b1×b2)×(b2×b3))⋅normalize(b2), (b1×b2)⋅(b2×b3)))
+end
+
+"""
+    Compute dihedral angles within a tetrahedra
+    radians
+"""
+function dihedral_angles(p,t)
+    AT = eltype(eltype(p))
+    nangs = length(t)*6
+    a = fill(zero(AT), nangs)
+    for i in 1:length(t)
+        t1, t2, t3, t4 = t[i]
+        a[i*6-5] = dihedral(p[t1],p[t2],p[t3],p[t4])
+        a[i*6-4] = dihedral(p[t1],p[t2],p[t4],p[t3])
+        a[i*6-3] = dihedral(p[t2],p[t1],p[t4],p[t3])
+        a[i*6-2] = dihedral(p[t2],p[t3],p[t4],p[t1])
+        a[i*6-1] = dihedral(p[t2],p[t1],p[t3],p[t4])
+        a[i*6]   = dihedral(p[t3],p[t1],p[t2],p[t4])
+    end
+    a
+end
+
+"""
+    Compute the minimum dihedral angle within a tetrahedra
+    radians
+"""
+function min_dihedral_angles(p,t)
+    AT = eltype(eltype(p))
+    nangs = length(t)
+    a = fill(zero(AT), nangs)
+    for i in 1:length(t)
+        t1, t2, t3, t4 = t[i]
+        d = (dihedral(p[t1],p[t2],p[t3],p[t4]),
+             dihedral(p[t1],p[t2],p[t4],p[t3]),
+             dihedral(p[t2],p[t1],p[t4],p[t3]),
+             dihedral(p[t2],p[t3],p[t4],p[t1]),
+             dihedral(p[t2],p[t1],p[t3],p[t4]),
+             dihedral(p[t3],p[t1],p[t2],p[t4]))
+        a[i] = minimum(d)
+    end
+    a
+end

test/runtests.jl

@@ -85,3 +85,15 @@ end
         @test isapprox(getproperty(s,fn), getproperty(stat_04,fn))
     end
 end
+
+@testset "dihedral metrics" begin
+    d(p) = sqrt(sum(p.^2))-1
+    p,t,_ = distmesh(d,HUniform(),0.2)
+    all_angs =  DistMesh.dihedral_angles(p,t)
+    min_angs =  DistMesh.min_dihedral_angles(p,t)
+    ax = extrema(all_angs)
+    mx = extrema(min_angs)
+    @test ax[1] == mx[1]
+    @test ax == (0.023502688273828173, 3.104396619996953)
+    @test mx == (0.023502688273828173, 1.2044902180168893)
+end