From a2aef1c8635f626d31931a8545d2617fa4d73950 Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Wed, 9 Aug 2017 22:27:23 +0300 Subject: [PATCH] Add functions to calculate polygon clips in 3D --- docs/src/api.md | 1 + src/Mortar3D.jl | 1 + src/calculate_polygon_clip.jl | 182 ++++++++++++++++++++++++++++ test/runtests.jl | 1 + test/test_calculate_polygon_clip.jl | 53 ++++++++ 5 files changed, 238 insertions(+) create mode 100644 src/calculate_polygon_clip.jl create mode 100644 test/test_calculate_polygon_clip.jl diff --git a/docs/src/api.md b/docs/src/api.md index f64477a..2e21f26 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -10,6 +10,7 @@ DocTestSetup = quote using Mortar3D using Mortar3D: calculate_normals using Mortar3D: project_from_plane_to_surface, project_from_surface_to_plane + using Mortar3D: approx_in, is_vertex_inside_polygon, check_orientation!, calculate_polygon_clip end ``` diff --git a/src/Mortar3D.jl b/src/Mortar3D.jl index 930cc13..cc36bc7 100644 --- a/src/Mortar3D.jl +++ b/src/Mortar3D.jl @@ -4,4 +4,5 @@ module Mortar3D include("calculate_normals.jl") include("calculate_projections.jl") +include("calculate_polygon_clip.jl") end diff --git a/src/calculate_polygon_clip.jl b/src/calculate_polygon_clip.jl new file mode 100644 index 0000000..ea3a0f0 --- /dev/null +++ b/src/calculate_polygon_clip.jl @@ -0,0 +1,182 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/Mortar3D.jl/blob/master/LICENSE + +""" + approx_in(q, P) + +Test does vector `P` contain approximately `q`. This function uses isapprox() +internally to make boolean test. + +# Example + +```jldoctest +P = Vector[[1.0, 1.0], [2.0, 2.0]] +q = [1.0, 1.0] + eps(Float64) +in(q, P), approx_in(q, P) + +# output + +(false, true) +``` +""" +function approx_in(q, P; rtol=1.0e-9, atol=1.0e-9) + for p in P + if isapprox(q, p; rtol=rtol, atol=atol) + return true + end + end + return false +end + +""" + is_vertex_inside_polygon(q, P; atol=1.0e-3) + +Test is vertex `q` inside polygon `P` with given tolerance. + +# Example + +```jldoctest +P = Vector[[0.0,0.0,0.0], [1.0,0.0,0.0], [0.0,1.0,0.0]] +q1 = [ 0.1, 0.0, 0.0] +q2 = [-0.1, 0.1, 0.0] +is_vertex_inside_polygon(q1, P), is_vertex_inside_polygon(q2, P) + +# output + +(true, false) + +``` +""" +function is_vertex_inside_polygon(q, P; atol=1.0e-3) + N = length(P) + a = 0.0 + for i=1:N + A = P[i] - q + B = P[mod(i,N)+1] - q + c = norm(A)*norm(B) + isapprox(c, 0.0; atol=atol) && return true + cosa = dot(A,B)/c + isapprox(cosa, 1.0; atol=atol) && return false + isapprox(cosa, -1.0; atol=atol) && return true + a += acos(cosa) + end + return isapprox(a, 2*pi; atol=atol) +end + +""" + check_orientation!(P, n) + +Given polygon `P` and normal direction `n`, ensure that polygon vertices are +ordered in counter clock wise direction with respect to surface normal. It is +assumed that polygon is convex. + +# Example + +Unit triangle, normal in z-direction: + +```jldoctest +P = Vector[[0.0, 0.0, 0.0], [0.0, 1.0, 0.0], [1.0, 0.0, 0.0]] +n = [0.0, 0.0, 1.0] +check_orientation!(P, n) + +# output + +3-element Array{Array{T,1} where T,1}: + [1.0, 0.0, 0.0] + [0.0, 0.0, 0.0] + [0.0, 1.0, 0.0] +``` +""" +function check_orientation!(P, n) + C = mean(P) + np = length(P) + s = [dot(n, cross(P[i]-C, P[mod(i+1,np)+1]-C)) for i=1:np] + all(s .< 0) && return + # debug("polygon not in ccw order, fixing") + # project points to new orthogonal basis Q and sort there + t1 = (P[1]-C)/norm(P[1]-C) + t2 = cross(n, t1) + Q = [n t1 t2] + sort!(P, lt=(A, B) -> begin + A_proj = Q'*(A-C) + B_proj = Q'*(B-C) + a = atan2(A_proj[3], A_proj[2]) + b = atan2(B_proj[3], B_proj[2]) + return a > b + end) + return P +end + +""" + calculate_polygon_clip(xs, xm, n) + +Given two polygons `xs` and `xm`, calculate polygon clipping `P`. It is +assumed that clipping is done in direction `n`. + +# Example + +```jldoctest +xs = Vector[[0, 0, 0], [1, 0, 0], [0, 1, 0]] +xm = Vector[[-1/4, 1/2, 0], [1/2, -1/4, 0], [3/4, 3/4, 0]] +n = [0, 0, 1] +P = calculate_polygon_clip(xs, xm, n) + +# output + +6-element Array{Array{T,1} where T,1}: + [0.65, 0.35, 0.0] + [0.5625, 0.0, 0.0] + [0.25, 0.0, 0.0] + [0.0, 0.25, 0.0] + [0.0, 0.5625, 0.0] + [0.35, 0.65, 0.0] + +``` +""" +function calculate_polygon_clip{T}(xs::Vector{T}, xm::Vector{T}, n::T) + # objective: search does line xm1 - xm2 clip xs + nm = length(xm) + ns = length(xs) + P = T[] + + # 1. test is master point inside slave, if yes, add to clip + for i=1:nm + if is_vertex_inside_polygon(xm[i], xs) + push!(P, xm[i]) + end + end + + # 2. test is slave point inside master, if yes, add to clip + for i=1:ns + # FIXME: point of this? + # approx_in(xs[i], P) && continue + if is_vertex_inside_polygon(xs[i], xm) + push!(P, xs[i]) + end + end + + for i=1:nm + # 2. find possible intersection + xm1 = xm[i] + xm2 = xm[mod(i,nm)+1] + for j=1:ns + xs1 = xs[j] + xs2 = xs[mod(j,ns)+1] + tnom = dot(cross(xm1-xs1, xm2-xm1), n) + tdenom = dot(cross(xs2-xs1, xm2-xm1), n) + isapprox(tdenom, 0) && continue + t = tnom/tdenom + (0 <= t <= 1) || continue + q = xs1 + t*(xs2 - xs1) + # FIXME: point of this? + # if approx_in(q, P) + # continue + # end + if is_vertex_inside_polygon(q, xm) + push!(P, q) + end + end + end + check_orientation!(P, n) + return P +end diff --git a/test/runtests.jl b/test/runtests.jl index e05e613..11641c5 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -8,6 +8,7 @@ const to = TimerOutput() test_files = String[] push!(test_files, "test_calculate_normals.jl") push!(test_files, "test_calculate_projections.jl") +push!(test_files, "test_calculate_polygon_clip.jl") @testset "Mortar3D.jl" begin for fn in test_files diff --git a/test/test_calculate_polygon_clip.jl b/test/test_calculate_polygon_clip.jl new file mode 100644 index 0000000..a01bdee --- /dev/null +++ b/test/test_calculate_polygon_clip.jl @@ -0,0 +1,53 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/Mortar3D.jl/blob/master/LICENSE + +using Base.Test + +using Mortar3D: approx_in, + is_vertex_inside_polygon, + calculate_polygon_clip + +@testset "approx in" begin + P = Vector[[1.0, 1.0], [2.0, 2.0]] + q1 = [1.0, 1.0] + eps(Float64) + q2 = [1.0, 1.0] + 0.1 + @test approx_in(q1, P) + @test !approx_in(q2, P) +end + +@testset "vertex inside polygon" begin + P = Vector[[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]] + q1 = [ 0.1, 0.1, 0.0] + q2 = [-0.1, 0.1, 0.0] + @test is_vertex_inside_polygon(q1, P) == true + @test is_vertex_inside_polygon(q2, P) == false +end + +@testset "calculate polygon clip general case" begin + xs = Vector[[0, 0, 0], [1, 0, 0], [0, 1, 0]] + xm = Vector[[-1/4, 1/2, 0], [1/2, -1/4, 0], [3/4, 3/4, 0]] + n = [0, 0, 1] + P = calculate_polygon_clip(xs, xm, n) + @test isapprox(P[6], [7/20, 13/20, 0]) + @test isapprox(P[5], [0, 9/16, 0]) + @test isapprox(P[4], [0, 1/4, 0]) + @test isapprox(P[3], [1/4, 0, 0]) + @test isapprox(P[2], [9/16, 0, 0]) + @test isapprox(P[1], [13/20, 7/20, 0]) +end + +@testset "calculate polygon clip, master node inside slave" begin + xs = Vector[[0.0,0.0,0.0], [1.0,0.0,0.0], [0.0,1.0,0.0]] + xm = Vector[[0.0,-0.9,0.0], [0.5,0.1,0.0], [1.0,-0.9,0.0]] + n = [0.0, 0.0, 1.0] + P = calculate_polygon_clip(xs, xm, n) + @test length(P) == 3 +end + +@testset "calculate polygon clip, slave node inside master" begin + xm = Vector[[0.0,0.0,0.0], [1.0,0.0,0.0], [0.0,1.0,0.0]] + xs = Vector[[0.0,-0.9,0.0], [0.5,0.1,0.0], [1.0,-0.9,0.0]] + n = [0.0, 0.0, 1.0] + P = calculate_polygon_clip(xs, xm, n) + @test length(P) == 3 +end