decoupling normal triangle

src/charts.jl

@@ -12,10 +12,40 @@ struct Simplex{U,D,C,N,T}
     normals::SVector{C,SVector{U,T}}
     volume::T
 end
-
+normal(t::Simplex{3,2,1,3,<:Number}) = t.normals[1]
 dimtype(splx::Simplex{U,D}) where {U,D} = Val{D}
+"""
+    permute_simplex(simplex,permutation)
+
+Permutation is a Vector v which sets the v[i]-th vertex at the i-th place.
+
+Return Simplex with permuted vertices list, tangents are recalculated, normal is kept the same
+"""
+function permute_vertices(s::Simplex{3,2,1,3,T},permutation::Union{Vector{P},SVector{P}}) where {T,P}
+    vert = SVector{3,SVector{3,T}}(s.vertices[permutation])
+    simp = simplex(vert)
+    Simplex(simp.vertices,simp.tangents,s.normals,simp.volume)
+end
+
+"""
+    flip_normal(simplex)
+
+Flips the normal of the simplex. Only on triangles embedded in 3D space
+"""
+flip_normal(t::Simplex{3,2,1,3,<:Number}) = Simplex(t.vertices,t.tangents,-t.normals,t.volume)
+
+"""
+    flip_normal(simplex, sign)
 
+Flips the normal of the simplex if sign is -1. Only on triangles embedded in 3D space
+"""
+function flip_normal(t::Simplex{3,2,1,3,<:Number},sign::Int)
+    sign == 1 && return t
+    return flip_normal(t)
+    end
+export flip_normal
 
+tangents(s::Simplex{3,2,1,3,<:Number},i) = s.tangents[i]
 """
     coordtype(simplex)
 
@@ -126,7 +156,7 @@ simplex(vertices...) = simplex(SVector((vertices...,)))
     :(simplex(SVector{$D1,$P}($xp)))
 end
 
-normal(s::Simplex{3,2,1,3,T}) where {T} = normalize(cross(s[1]-s[3],s[2]-s[3]))
+# normal(s::Simplex{3,2,1,3,T}) where {T} = normalize(cross(s[1]-s[3],s[2]-s[3]))
 
 
 function _normals(tangents, ::Type{Val{1}})
@@ -341,4 +371,12 @@ Returns a matrix whose columns are the tangents of the simplex `splx`.
 """
 tangents(splx::Simplex) = hcat((splx.tangents)...)
 
-vertices(splx::Simplex) = hcat((splx.vertices)...)
+vertices(splx::Simplex) = hcat((splx.vertices)...)
+
+"""
+    verticeslist(simplex)
+
+Returns the vertices as a list.
+"""
+verticeslist(splx::Simplex) = splx.vertices
+export verticeslist

src/intersect.jl

@@ -36,16 +36,44 @@ The output inherits the orientation from the first operand
 """
 function intersection(p1::Simplex{U,2,C,3}, p2::Simplex{U,2,C,3}) where {U,C}
   pq = sutherlandhodgman(p1.vertices, p2.vertices)
-  [ simplex(pq[1], pq[i], pq[i+1]) for i in 2:length(pq)-1 ]
+  return [ simplex(pq[1], pq[i], pq[i+1]) for i in 2:length(pq)-1 ]
 end
 
+function intersection(p1::Simplex{3,2,1,3}, p2::Simplex{3,2,1,3};tol=eps()) 
+    pq = sutherlandhodgman(p1.vertices, p2.vertices)
+    nonoriented_simplexes = [ simplex(pq[1], pq[i], pq[i+1]) for i in 2:length(pq)-1 ]
+    nonoriented_simplexes = nonoriented_simplexes[volume.(nonoriented_simplexes).>tol]
+    signs = Int.(sign.(dot.(normal.(nonoriented_simplexes),Ref(normal(p1)))))
+    flip_normal.(nonoriented_simplexes,signs)
+  end
+
 
 function intersection(p1::Simplex{U,3,C,4}, p2::Simplex{U,3,C,4}) where {U,C}
   @assert overlap(p1,p2)
   volume(p1) <= volume(p2) ? [p1] : [p2]
 end
 
+"""
+    intersection2(triangleA, triangleB)
+
+returns intersection in which the first operand inhirits orientation from first argument
+and second argument inhirits orientation of second argument.
+"""
+function intersection2(p1::Simplex, p2::Simplex)
+    a = intersection(p1,p2)
+    return [(i,i) for i in a]
+end
 
+function intersection2(p1::Simplex{3,2,1,3}, p2::Simplex{3,2,1,3}; tol=eps())
+    pq = sutherlandhodgman(p1.vertices, p2.vertices)
+    nonoriented_simplexes = [ simplex(pq[1], pq[i], pq[i+1]) for i in 2:length(pq)-1 ]
+    nonoriented_simplexes = nonoriented_simplexes[volume.(nonoriented_simplexes).>tol]
+    signs1 = Int.(sign.(dot.(normal.(nonoriented_simplexes),Ref(normal(p1)))))
+
+    signs2 = Int.(sign.(dot.(normal.(nonoriented_simplexes),Ref(normal(p2)))))
+    return [(flip_normal(s,signs1[i]),flip_normal(s,signs2[i])) for (i,s) in enumerate(nonoriented_simplexes)]
+end
+export intersection2
 
 """
     intersectline(a,b,p,q)

src/mesh.jl

@@ -139,7 +139,7 @@ or vertices not appearing in any cell. In other words the following is not neces
     numvertices(mesh) == numcells(skeleton(mesh,0))
 ```
 """
-numvertices(m::Mesh) = length(m.vertices)
+numvertices(m::Mesh) = length(vertices(m))
 
 
 
@@ -148,7 +148,7 @@ numvertices(m::Mesh) = length(m.vertices)
 
 Returns the number of cells in the mesh.
 """
-numcells(m::AbstractMesh) = length(m.faces)
+numcells(m::AbstractMesh) = length(cells(m))
 
 
 

src/neighborhood.jl

@@ -26,7 +26,7 @@ function parametric(mp::MeshPointNM)
     mp.bary
 end
 
-chart(nbd::MeshPointNM) = mp.patch
+chart(nbd::MeshPointNM) = nbd.patch
 
 "Return the barycentric coordinates of `mp`"
 barycentric(mp::MeshPointNM) = SVector(mp.bary[1], mp.bary[2], 1-mp.bary[1]-mp.bary[2])
@@ -80,3 +80,4 @@ to parameter `(1/(D+1), 1/(D+1), ...)` where `D` is the simplex dimension.
     uv = ones(T,D)/(D+1)
     :(neighborhood(p, $uv))
 end
+

src/quadpoints.jl

@@ -80,7 +80,6 @@ this method used `quadpoints(chart, rule)` to retrieve the points and weights fo
 a certain quadrature rule over `chart`.
 """
 function quadpoints(f, charts, rules)
-
     pw = quadpoints(charts[1], rules[1])
     P = typeof(pw[1][1])
     W = typeof(pw[1][2])

src/subd_chart.jl

@@ -30,3 +30,4 @@ function getelementVertices(chart::subd_chart{T}) where {T}
     verticecoords3[3] = verticescoords[chart.N-6]
     return verticecoords3
 end
+verticeslist(chart::subd_chart) = chart.vertices

test/runtests.jl

@@ -36,6 +36,7 @@ include("test_sh_intersection.jl")
 include("test_isinside.jl")
 include("test_jctweld.jl")
 include("test_isinclosure.jl")
+include("test_permute_vertices.jl")
 
 include("test_trgauss.jl")
 include("test_chartquad.jl")

test/test_permute_vertices.jl

@@ -0,0 +1,34 @@
+using CompScienceMeshes
+using Test
+
+splx = simplex(
+    point(0,0,0),
+    point(1,0,0),
+    point(0,1,0)
+);
+
+@test normal(splx) ≈ point(0,0,1)
+
+splx1 = CompScienceMeshes.permute_vertices(splx, [2,1,3])
+@test normal(splx1) ≈ point(0,0,1)
+
+splx2 = CompScienceMeshes.flip_normal(splx)
+@test normal(splx2) ≈ point(0,0,-1)
+
+for i in 1:2
+    @test CompScienceMeshes.tangents(splx2,i) ≈ CompScienceMeshes.tangents(splx,i)
+end
+
+splx3 = simplex(
+    point(1,0.5,0),
+    point(-1,0.5,0),
+    point(-1,2.5,0)
+)
+
+isct = CompScienceMeshes.intersection2(splx, splx3)
+@test length(isct) == 1
+
+splx4, splx5 = isct[1]
+@test volume(splx4) ≈ volume(splx5)
+@test normal(splx4) ≈ -normal(splx5)
+