diff --git a/src/charts.jl b/src/charts.jl index 03d8ccd..861f5ee 100644 --- a/src/charts.jl +++ b/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)...) \ No newline at end of file +vertices(splx::Simplex) = hcat((splx.vertices)...) + +""" + verticeslist(simplex) + +Returns the vertices as a list. +""" +verticeslist(splx::Simplex) = splx.vertices +export verticeslist \ No newline at end of file diff --git a/src/intersect.jl b/src/intersect.jl index 92cd380..c224330 100644 --- a/src/intersect.jl +++ b/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) diff --git a/src/mesh.jl b/src/mesh.jl index eea46ff..5a43d59 100644 --- a/src/mesh.jl +++ b/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)) diff --git a/src/neighborhood.jl b/src/neighborhood.jl index e107e2c..c56c28a 100644 --- a/src/neighborhood.jl +++ b/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 + diff --git a/src/quadpoints.jl b/src/quadpoints.jl index 1b54e07..3d789d1 100644 --- a/src/quadpoints.jl +++ b/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]) diff --git a/src/subd_chart.jl b/src/subd_chart.jl index 351058c..bddf5af 100644 --- a/src/subd_chart.jl +++ b/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 \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 69c255e..8a2372f 100644 --- a/test/runtests.jl +++ b/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") diff --git a/test/test_permute_vertices.jl b/test/test_permute_vertices.jl new file mode 100644 index 0000000..a095eff --- /dev/null +++ b/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) +