Merge pull request #53 from vpuri3/vp/vectorize_freps

WIP: vectorize freps

.gitignore

@@ -6,3 +6,6 @@
 *.obj
 .cache/*
 Manifest*
+
+*.swp
+*.DS_Store

Project.toml

@@ -21,10 +21,10 @@ julia = "1.6"
 
 [extras]
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
+FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
+MeshIO = "7269a6da-0436-5bbc-96c2-40638cbb6118"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
-MeshIO = "7269a6da-0436-5bbc-96c2-40638cbb6118"
-FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
 
 [targets]
 test = ["Test", "Combinatorics", "StaticArrays", "MeshIO", "FileIO"]

src/frep.jl

@@ -1,41 +1,93 @@
 # http://en.wikipedia.org/wiki/Function_representation
+#----------------------------------
 
-function _radius(a,b,r)
-    if abs(a-b) >= r
-        return min(a,b)
-    else
-        return b+r*sin(pi/4+asin((a-b)/(r*sqrt(2))))-r
-    end
+"""
+ Overloading CoordinateTransformations.Translation for vectorized implementation.
+"""
+function (trans::Translation{V})(x::AbstractMatrix) where {V}
+    x .+ trans.translation
 end
 
-function FRep(p::MapContainer{3,T,P}, v) where {T,P}
-    FRep(p.primitive, p.inv(v))
+#function FRep(p::MapContainer{N,T,P}, v) where {N,T,P}
+#    FRep(p.primitive, p.inv(v))
+#end
+
+function FRep(p::MapContainer{N,T,P}, v...) where {N,T,P}
+    FRep(p.primitive, p.inv(vcat(v...)))
 end
+#----------------------------------
 
+function FRep(u::CSGUnion, v...)
+    min.(FRep(u.left, v...),FRep(u.right, v...))
+end
 
-function FRep(p::MapContainer{2,T,P}, v) where {T,P}
-    FRep(p.primitive, p.inv(SVector(v[1],v[2])))
+function FRep(u::CSGDiff, v...)
+    max.(FRep(u.left, v...), -FRep(u.right, v...))
 end
 
-function FRep(p::Sphere, v)
-    x = v[1]
-    y = v[2]
-    z = v[3]
-    sqrt(x*x + y*y + z*z) - p.radius
+function FRep(u::CSGIntersect, v...)
+    max.(FRep(u.left, v...), FRep(u.right, v...))
+end
+#----------------------------------
+
+function FRep(p::Sphere,v::AbstractVector)
+    norm(v) - p.radius
+end
+
+function FRep(p::Sphere,x::AbstractMatrix,y::AbstractMatrix,z::AbstractMatrix)
+    @assert size(x) == size(y) == size(z)
+    @. sqrt(x*x + y*y + z*z) - p.radius
 end
 
-function FRep(p::Cylinder, v)
+function FRep(p::Sphere,xyz::AbstractMatrix)
+    @assert size(xyz)[1] == 3
+    r2 = sum(xyz .* xyz,dims=1)
+    @. sqrt(r2) - p.radius
+end
+#----------------------------------
+
+function FRep(p::Cylinder,v::AbstractVector)
     x = v[1]
     y = v[2]
     z = v[3]
     max(-z+p.bottom, z-p.height-p.bottom, sqrt(x*x + y*y) - p.radius)
 end
 
-function FRep(p::Circle, v)
+function FRep(p::Cylinder,x::AbstractMatrix,y::AbstractMatrix,z::AbstractMatrix)
+    @assert size(x) == size(y) == size(z)
+    r = @. sqrt(x*x + y*y) - p.radius
+    a = @. -z+p.bottom
+    b = @.  z-(p.height+p.bottom)
+    max.(a, b, r)
+end
+
+function FRep(p::Cylinder,xyz::AbstractMatrix)
+    @assert size(xyz)[1] == 3
+
+    x = xyz[1,:]'
+    y = xyz[2,:]'
+    z = xyz[3,:]'
+    fr = FRep(p,x,y,z)
+end
+#----------------------------------
+
+function FRep(p::Circle, v::AbstractVector)
     norm(v) - p.radius
 end
 
-function FRep(p::Cuboid, v)
+function FRep(p::Circle,x::AbstractMatrix,y::AbstractMatrix)
+    @assert size(x) == size(y)
+    @. sqrt(x*x + y*y) - p.radius
+end
+
+function FRep(p::Circle,xy::AbstractMatrix)
+    @assert size(xy)[1] == 2
+    r2 = sum(xy .* xy,dims=1)
+    @. sqrt(r2) - p.radius
+end
+#----------------------------------
+
+function FRep(p::Cuboid, v::AbstractVector)
     x = v[1]
     y = v[2]
     z = v[3]
@@ -46,7 +98,30 @@ function FRep(p::Cuboid, v)
         -z+lbz, z-dz-lbz)
 end
 
-function FRep(p::Square, v)
+function FRep(p::Cuboid,x::AbstractMatrix,y::AbstractMatrix,z::AbstractMatrix)
+    @assert size(x) == size(y) == size(z)
+    dx, dy, dz = p.dimensions
+    lbx, lby,lbz = p.lowercorner
+    a = @. -x + lbx
+    b = @.  x - (dx+lbx)
+    c = @. -y + lby
+    d = @.  y - (dy+lby)
+    e = @. -z + lbz
+    f = @.  z - (dz+lbz)
+    max.(a,b,c,d,e,f)
+end
+
+function FRep(p::Cuboid,xyz::AbstractMatrix)
+    @assert size(xyz)[1] == 3
+
+    x = xyz[1,:]'
+    y = xyz[2,:]'
+    z = xyz[3,:]'
+    fr = FRep(p,x,y,z)
+end
+#----------------------------------
+
+function FRep(p::Square, v::AbstractVector)
     x = v[1]
     y = v[2]
     dx, dy = p.dimensions
@@ -55,21 +130,47 @@ function FRep(p::Square, v)
         -y+lby, y-dy-lby)
 end
 
-function FRep(u::CSGUnion, v)
-    min(FRep(u.left, v),FRep(u.right, v))
-end
+function FRep(p::Square,x::AbstractMatrix,y::AbstractMatrix)
+    @assert size(x) == size(y)
+    dx, dy = p.dimensions
+    lbx, lby = p.lowercorner
+    a = @. -x + lbx
+    b = @.  x - (dx+lbx)
+    c = @. -y + lby
+    d = @.  y - (dy+lby)
 
-function FRep(u::CSGDiff, v)
-    max(FRep(u.left, v), -FRep(u.right, v))
+    max.(a,b,c,d)
 end
 
-function FRep(u::CSGIntersect, v)
-    max(FRep(u.left, v), FRep(u.right, v))
+function FRep(p::Square,xy::AbstractMatrix)
+    @assert size(xy)[1] == 2
+
+    x = xy[1,:]'
+    y = xy[2,:]'
+    fr = FRep(p,x,y)
 end
+#----------------------------------
 
 function FRep(s::Shell, v)
     r = FRep(s.primitive, v)
-    max(r, -r-s.distance)
+    max.(r, -r .- s.distance)
+end
+
+#----------------------------------
+function _radius(a::Real,b::Real,r::Real)
+    if abs(a-b) >= r
+        return min(a,b)
+    else
+        return b+r*sin(pi/4+asin((a-b)/(r*sqrt(2))))-r
+    end
+end
+
+function _radius(a::AbstractVector,b::AbstractVector,c::AbstractVector)
+    cond = @. (abs(a-b) >= r) # bit-vector
+    ans1 = @. min(a,b)
+    ans2 = @. b+r*sin(pi/4+asin((a-b)/(r*sqrt(2))))-r
+
+    @. cond*ans1 + (1-cond)*ans2
 end
 
 function FRep(u::RadiusedCSGUnion, v)
@@ -78,8 +179,9 @@ function FRep(u::RadiusedCSGUnion, v)
     r = u.radius
     _radius(a,b,r)
 end
+#----------------------------------
 
-function FRep(p::Piping{T}, v) where {T}
+function FRep(p::Piping{T}, v::AbstractVector) where {T}
     num_pts = length(p.points)
 
     val = typemax(T)
@@ -101,10 +203,42 @@ function FRep(p::Piping{T}, v) where {T}
     val - p.radius
 end
 
-function FRep(p::LinearExtrude, v)
+function FRep(p::Piping{T}, x::AbstractMatrix, y::AbstractMatrix, z::AbstractMatrix) where {T}
+
+    xyz = vcat(x,y,z)
+    collect(FRep(p,xyz[:,i]) for i in 1:size(xyz)[2])'
+end
+
+function FRep(p::Piping{T}, xyz::AbstractMatrix) where {T}
+    @warn "F-Reps for Piping is not vectorized. switching to non-vectorized implementation" 
+
+    collect(FRep(p,xyz[:,i]) for i in 1:size(xyz)[2])'
+end
+#----------------------------------
+
+function FRep(p::LinearExtrude, v::AbstractVector)
     x = v[1]
     y = v[2]
     z = v[3]
-    r = FRep(p.primitive, v)
-    max(max(-z,z-p.distance), r)
+    r = FRep(p.primitive, [x,y])
+    max.(max.(-z,z-p.distance), r)
+end
+
+function FRep(p::LinearExtrude,x::AbstractMatrix,y::AbstractMatrix,z::AbstractMatrix)
+    @assert size(x) == size(y) == size(z)
+
+    xy = vcat(x,y)
+    r = FRep(p.primitive,xy)
+    @. max(max(-z,z-p.distance), r)
+end
+
+function FRep(p::LinearExtrude,xyz::AbstractMatrix)
+    @assert size(xyz)[1] == 3
+
+    x = xyz[1,:]'
+    y = xyz[2,:]'
+    z = xyz[3,:]'
+    fr = FRep(p,x,y,z)
+
 end
+#----------------------------------

test/frep.jl

@@ -0,0 +1,67 @@
+#
+@testset "FReps" begin
+    N = 10
+    xy  = rand(2,N)
+    xyz = rand(3,N)
+
+    x = rand(1,N)
+    y = rand(1,N)
+    z = rand(1,N)
+    #----------------------------------
+
+    p  = Circle(1.0)
+
+    @test FRep(p,xy[:,1]) |> size == ()
+    @test FRep(p,xy)      |> size == (1,N)
+    @test FRep(p,x,y)     |> size == (1,N)
+    #----------------------------------
+
+    p = Square([1.0,1.0],center=true)
+
+    @test FRep(p,xy[:,1]) |> size == ()
+    @test FRep(p,xy)      |> size == (1,N)
+    @test FRep(p,x,y)     |> size == (1,N)
+    #----------------------------------
+
+    p = Sphere(1.0)
+
+    @test FRep(p,xyz[:,1]) |> size == ()
+    @test FRep(p,xyz)      |> size == (1,N)
+    @test FRep(p,x,y,z)    |> size == (1,N)
+    #----------------------------------
+
+    p = Cuboid([1.0,1.0,1.0],center=true)
+
+    @test FRep(p,xyz[:,1]) |> size == ()
+    @test FRep(p,xyz)      |> size == (1,N)
+    @test FRep(p,x,y,z)    |> size == (1,N)
+    #----------------------------------
+
+    p = Cylinder(1.0,1.0,1.0)
+
+    @test FRep(p,xyz[:,1]) |> size == ()
+    @test FRep(p,xyz)      |> size == (1,N)
+    @test FRep(p,x,y,z)    |> size == (1,N)
+    #----------------------------------
+
+    p = LinearExtrude(Circle(1.0),2.0)
+
+    @test FRep(p,xyz[:,1]) |> size == ()
+    @test FRep(p,xyz)      |> size == (1,N)
+    @test FRep(p,x,y,z)    |> size == (1,N)
+    #----------------------------------
+    c = translate([1.0,0.0])Circle(1.0)
+    s = translate([0.0,1.0])Square([3.0,3.0];center=true)
+    p = diff(s,c)
+
+    @test FRep(p,xy[:,1]) |> size == ()
+    @test FRep(p,xy)      |> size == (1,N)
+    @test FRep(p,x,y)     |> size == (1,N)
+    #----------------------------------
+    p = Piping(1.0, [[0,0,0],[10,0,0],[10,10,0],[10,10,10],[5,5,5]])
+
+    @test FRep(p,xyz[:,1]) |> size == ()
+    @test FRep(p,xyz)      |> size == (1,N)
+    @test FRep(p,x,y,z)    |> size == (1,N)
+    #----------------------------------
+end

test/runtests.jl

@@ -17,3 +17,4 @@ include("2d.jl")
 include("hyperrectangles.jl")
 include("meshes.jl")
 include("examples.jl")
+include("frep.jl")