diff --git a/src/CompScienceMeshes.jl b/src/CompScienceMeshes.jl index 817372a..deea6f7 100644 --- a/src/CompScienceMeshes.jl +++ b/src/CompScienceMeshes.jl @@ -38,6 +38,7 @@ export connectivity, cellpairs # marked for deprecation export translate, translate!, rotate, rotate! export fliporientation!, fliporientation export weld, union +export permutate # mesh refinement export barycentric_refinement, bisecting_refinement @@ -142,6 +143,7 @@ include("primitives/primitives.jl") include("baryref.jl") include("subdivision.jl") include("weld.jl") +include("permutatie_mesh.jl") include("mapper.jl") include("restrict.jl") diff --git a/src/charts.jl b/src/charts.jl index b74ead2..beb761d 100644 --- a/src/charts.jl +++ b/src/charts.jl @@ -217,7 +217,7 @@ function _normals(tangents::SVector{2,SVector{2,T}}, ::Type{Val{0}}) where {T} t = tangents[1] s = tangents[2] - v = (t[1]*s[2] - t[2]*s[1])/2 + v = abs(t[1]*s[2] - t[2]*s[1])/2 # n[3] = tangents[1] × tangents[2] # l = norm(n) diff --git a/src/mesh.jl b/src/mesh.jl index f801417..b1f1b2c 100644 --- a/src/mesh.jl +++ b/src/mesh.jl @@ -331,7 +331,7 @@ Base.getindex(m::AbstractMesh, I::Vector{Int}) = Mesh(vertices(m), cells(m)[I]) Returns the boundary of `mesh` as a mesh of lower dimension. """ -function boundary(mesh) +function boundary(mesh, include_interior_points=true) D = dimension(mesh) @@ -366,7 +366,22 @@ function boundary(mesh) end resize!(bnd_edges, i-1) - bnd = Mesh(vertices(mesh), bnd_edges) + if include_interior_points + return Mesh(vertices(mesh), bnd_edges) + else + originaltonew = Dict{Int,Int}() + sizehint!(originaltonew, i-1) + bnd_edges_new = Vector{I}(undef, length(bnd_edges)) + for (e,edge) in enumerate(bnd_edges) + get!(originaltonew, edge.indices[1], e) + end + for (e,edge) in enumerate(bnd_edges) + bnd_edges_new[e] =SimplexGraph( + get(originaltonew, edge.indices[1], -1), + get(originaltonew, edge.indices[2], -1)) + end + return Mesh(vertices(mesh)[collect(keys(originaltonew))], bnd_edges_new) + end end diff --git a/src/overlap.jl b/src/overlap.jl index bef6432..193e165 100644 --- a/src/overlap.jl +++ b/src/overlap.jl @@ -106,3 +106,34 @@ function overlap(p::Simplex{3,3,0,4,T}, q::Simplex{3,3,0,4,T}) where T all(0+tol .<= u .<= 1-tol) && return true return false end + + +function overlap(p::Simplex{2,2,0,3,T}, q::Simplex{2,2,0,3,T}) where T + + # tol = sqrt(eps(T)) + tol = 1e3 * eps(T) + + # Are the patches in the same plane? + u1 = q.tangents[1] + u2 = q.tangents[2] + v = p.vertices[1] - q.vertices[2] + + for i in 1:3 + a = p.vertices[mod1(i+1,3)] + b = p.vertices[mod1(i+2,3)] + c = p.vertices[i] + t = b - a + m = StaticArrays.SVector{2,T}(t[2],-t[1]) + + sp = zeros(T,3); sp[i] = dot(c-a, m) + sq = T[dot(q.vertices[j]-a, m) for j in 1:3] + + minp, maxp = extrema(sp) + minq, maxq = extrema(sq) + + maxq <= minp + tol && return false + maxp <= minq + tol && return false + end + + return true +end \ No newline at end of file diff --git a/src/permutatie_mesh.jl b/src/permutatie_mesh.jl new file mode 100644 index 0000000..b74dfca --- /dev/null +++ b/src/permutatie_mesh.jl @@ -0,0 +1,113 @@ +""" +permutate_mesh permutate the vertices of a mesh, while keeping the same cells. + +Permutation is represented as a Vector v which sets the v[i]-th vertex at the i-th place. + +mesh can be permutated from: +- Vector or SVector +- Permutations.Permutation +- another mesh with same universedimension +- another mesh with different universedimension, only works for combination of udims 2,3. + +In the last two, control that the other mesh is a submesh of the starting mesh. In either direct + + +""" + +function permutate(mesh::Mesh, σ::Permutations.Permutation) + mesh.vertices = mesh.vertices[σ.data] + for (i, a) in enumerate(mesh.faces) + ind = σ'.data[a.indices] + mesh.faces[i] = SimplexGraph(ind...) + end +end + +function permutate(mesh::Mesh, σ::Union{Vector{<:Integer}, StaticArrays.SVector{Int32, <:Integer}}) + @assert numvertices(mesh) == length(σ) + permutate(mesh, Permutations.Permutation(σ)) +end + +function permutate(X::Mesh{U}, Y::Mesh{U}) where {U} + tol = sqrt(eps(coordtype(X))) + permut = Vector{Int32}() + temp = collect(1:numvertices(X)) + for p in Y.vertices + index = findfirst(isapprox(p;atol = tol), X.vertices) + @assert !isnothing(index) + + temp[index] = 0 + append!(permut, index) + end + + for i in temp + if i !=0 + push!(permut, i) + end end + permutate(X, Permutations.Permutation(permut)) + X +end + +function permutate(X::Mesh{2}, Y::Mesh{3}) + tol = sqrt(eps(coordtype(X))) + + # assert that the z-coordinate for the vertices in X is constant. + + permut = Vector{Int32}() + temp = collect(1:numvertices(X)) + + for p in Y.vertices + q = zeros(eltype(X.vertices[begin]), 2) + q[begin] = p[begin] + q[2] = p[2] + + + index = findfirst(isapprox(q;atol = tol), X.vertices) + @assert !isnothing(index) # vertex from Y exist in X + @assert temp[index] != 0 # vertex from Y unique in X + + temp[index] = 0 + append!(permut, index) + end + + for i in temp + if i !=0 + push!(permut, i) + end end + + permutate(X, Permutations.Permutation(permut)) + X +end + +function permutate(X::Mesh{3}, Y::Mesh{2}) + tol = sqrt(eps(coordtype(X))) + + # assert that the z-coordinate for the vertices in X is constant. + for v in X.vertices + @assert v[3] == X.vertices[1][3] + end + + permut = Vector{Int32}() + temp = collect(1:numvertices(X)) + + for p in Y.vertices + q = zeros(eltype(X.vertices[begin]), 3) + q[begin] = p[begin] + q[2] = p[2] + q[3] = X.vertices[1][3] + + index = findfirst(isapprox(q;atol = tol), X.vertices) + @assert !isnothing(index) # vertex from Y exist in X + @assert temp[index] != 0 # vertex from Y unique in X + + temp[index] = 0 + append!(permut, index) + end + + for i in temp + if i !=0 + push!(permut, i) + end end + + permutate(X, Permutations.Permutation(permut)) + X +end diff --git a/src/plotlyjs_glue.jl b/src/plotlyjs_glue.jl index 9fc87d2..d57f0e3 100644 --- a/src/plotlyjs_glue.jl +++ b/src/plotlyjs_glue.jl @@ -44,6 +44,44 @@ function __init__() return s end + @eval function patch(Γ::CompScienceMeshes.AbstractMesh{2}, fcr=nothing; + caxis=nothing, showscale=true, color="red", kwargs...) + + v = vertexarray(Γ) + c = cellarray(Γ) + + x = v[:,1]; y = v[:,2]; z = zeros(length(y)) + i = c[:,1].-1; j = c[:,2].-1; k = c[:,3].-1 + + + if fcr != nothing + m, M = extrema(fcr) + if caxis != nothing + m, M = caxis + end + + s = PlotlyBase.mesh3d(; + x=x, y=y, z=z, + i=i, j=j, k=k, + intensitymode="cell", + intensity=fcr, + colorscale="Viridis", + showscale=showscale, + cmin=m, + cmax=M, + kwargs... + ) + else + s = PlotlyBase.mesh3d(; + x=x, y=y, z=z, + i=i, j=j, k=k, + color=color, + kwargs... + ) + end + return s + end + @eval function patch(a::Vector{<:Simplex}; kwargs...) vertices = reduce(vcat, [v.vertices for v in a]) faces = collect(SVector(3*(i-1)+1, 3*(i-1)+2, 3*(i-1)+3) for i in 1:length(a)) diff --git a/test/runtests.jl b/test/runtests.jl index f279d15..ff0c0c9 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -32,12 +32,14 @@ include("test_geometry.jl") include("test_patches.jl") include("test_submesh.jl") include("test_overlap.jl") +include("test_overlap2D.jl") include("test_intersect.jl") include("test_sh_intersection.jl") include("test_isinside.jl") include("test_jctweld.jl") include("test_isinclosure.jl") include("test_permute_vertices.jl") +include("test_permutate_mesh.jl") include("test_trgauss.jl") include("test_trdunavant.jl") diff --git a/test/test_mesh.jl b/test/test_mesh.jl index 749375d..a4f5e18 100644 --- a/test/test_mesh.jl +++ b/test/test_mesh.jl @@ -82,4 +82,23 @@ T = Float64 Σᵀ = connectivity(edges, faces, identity) Σ = connectivity(faces, edges, identity) @test norm(Σᵀ - Σ', Inf) == 0 + + rectangle = meshrectangle(T(1.0), T(1.0), T(0.5)); + bnd = boundary(rectangle, true) + + @test dimension(bnd) == 1 + @test numvertices(bnd) == 9 + @test numcells(bnd) == 8 + + bnd = boundary(rectangle, false) + + @test dimension(bnd) == 1 + @test numvertices(bnd) == 8 + @test numcells(bnd) == 8 + + + m = meshrectangle(T(1.0),T(1.0),T(0.5)) + CompScienceMeshes.rotate!(m, [1,0,0]*T(π/2)) + + @test CompScienceMeshes.isoriented(boundary(m, false)) # end \ No newline at end of file diff --git a/test/test_overlap2D.jl b/test/test_overlap2D.jl new file mode 100644 index 0000000..942402e --- /dev/null +++ b/test/test_overlap2D.jl @@ -0,0 +1,32 @@ +using Test +using CompScienceMeshes +using StaticArrays + +p = simplex( + point(0,0), + point(1,0), + point(0,1)) + +q1 = simplex( + point(0.6, 0.6), + point(1.6, 0.6), + point(0.6, 1.6), +) + +q2 = simplex( + point(0.4, 0.4), + point(1.4, 0.4), + point(0.4, 1.4), +) + +@test overlap(p, q1) == false +@test overlap(p, q2) == true + + +# test a case where the segments are: +# not of unit length +# colinear and opposite +# meet in a common point +ch1 = simplex(point(1/3,0), point(1/3,1/3)) +ch2 = simplex(point(1/3,1/3), point(1/3,2/3)) +@test !overlap(ch1, ch2) \ No newline at end of file diff --git a/test/test_patches.jl b/test/test_patches.jl index 403c156..1c5b63d 100644 --- a/test/test_patches.jl +++ b/test/test_patches.jl @@ -1,6 +1,7 @@ using Test using CompScienceMeshes +# universe dimension 3 for T in [Float32, Float64] local mesh = meshrectangle(T(1.0), T(1.0), T(1.0)) # local faces = skeleton(mesh, 2) @@ -19,3 +20,22 @@ for T in [Float32, Float64] point(T, 0.0, 0.0, -1.0)] @test p.volume == T(0.5) end + +# universe dimension 2 +for T in [Float32, Float64] + local mesh = meshrectangle(T(1.0), T(1.0), T(1.0),2) + # local faces = skeleton(mesh, 2) + # local verts = cellvertices(mesh, 1) + p = chart(mesh, cells(mesh)[1]) + #p = simplex(verts) + + @test p.vertices == [ + point(T, 0.0, 0.0), + point(T, 0.0, 1.0), + point(T, 1.0, 0.0)] + @test p.tangents == [ + point(T, -1.0, 0.0), + point(T, -1.0, 1.0)] + @test p.normals == [] + @test p.volume == T(0.5) +end diff --git a/test/test_permutate_mesh.jl b/test/test_permutate_mesh.jl new file mode 100644 index 0000000..495ff98 --- /dev/null +++ b/test/test_permutate_mesh.jl @@ -0,0 +1,72 @@ +using Test +using CompScienceMeshes +import Permutations + +# check if permutated mesh has the same vertices and that the faces are the same. +function check_mesh(X, X_check) + # test same vertices + @test numvertices(X) == numvertices(X_check) + @test Set(X.vertices) == Set(X_check.vertices) + + # test that each face has the same points in identical order + @test numcells(X) == numcells(X_check) + for i in 1:numcells(X) + @test X.vertices[X.faces[i].indices] == X_check.vertices[X_check.faces[i].indices] + end +end + +# same universe dimension +for u in 2:3 + local X = meshrectangle(1.0, 1.0, 0.5, u) + local X_check = meshrectangle(1.0, 1.0, 0.5, u) + local Y = meshrectangle(0.5, 0.5, 0.5, u) + + local σ = Permutations.Permutation(reverse(collect(1:numvertices(X)))) + + permutate(X,σ.data) + check_mesh(X, X_check) + + permutate(X, σ') + @test X.vertices == X_check.vertices + @test X.faces == X_check.faces + + permutate(X, Y) + check_mesh(X, X_check) + + @test X.vertices[begin: begin+numvertices(Y)-1] == Y.vertices +end + + + +# different universe dimension +let + local X = meshrectangle(1.0, 1.0, 0.5, 2) + local X_check = meshrectangle(1.0, 1.0, 0.5, 2) + local Y = meshrectangle(0.5, 0.5, 0.5, 3) + + permutate(X, Y) + check_mesh(X, X_check) + @test X.vertices[begin: begin+numvertices(Y)-1] == map(x->x[1:2], Y.vertices) + + + permutate(X, X_check) + + translate!(Y, point(0,0,1)) + permutate(X, Y) + check_mesh(X, X_check) + @test X.vertices[begin: begin+numvertices(Y)-1] == map(x->x[1:2], Y.vertices) +end + +let + local X = meshrectangle(1.0, 1.0, 0.5, 3) + local X_check = meshrectangle(1.0, 1.0, 0.5, 3) + local Y = meshrectangle(0.5, 0.5, 0.5, 2) + translate!(X, point(0,0,1)) + translate!(X_check, point(0,0,1)) + + permutate(X, Y) + check_mesh(X, X_check) + + @test map(x->x[1:2], X.vertices[begin: begin+numvertices(Y)-1]) == Y.vertices + +end \ No newline at end of file