Rectangle Surface Primitive and Integration with Box

src/ConstructiveSolidGeometry/PointsAndVectors.jl

@@ -36,11 +36,15 @@ end
 
 @inline _in_planar_α(p::PlanarPoint{T}, α::AbstractInterval) where {T} = _in_angular_interval_closed(mod(atan(p.v, p.u), T(2π)), α)
 
+@inline _isapprox_x(p::CartesianPoint{T}, x::Real) where {T} = isapprox(p.x, x, atol = geom_atol_zero(T))
+
 @inline _in_x(p::CartesianPoint, x::Real) = abs(p.x) <= x
 @inline _in_x(p::CartesianPoint, x::AbstractInterval) = p.x in x
 
 @inline _in_planar_u(p::PlanarPoint, u::AbstractInterval) = p.u in u
 
+@inline _isapprox_y(p::CartesianPoint{T}, y::Real) where {T} = isapprox(p.y, y, atol = geom_atol_zero(T))
+
 @inline _in_y(p::CartesianPoint, y::Real) = abs(p.y) <= y
 @inline _in_y(p::CartesianPoint, y::AbstractInterval) = p.y in y
 
@@ -102,9 +106,13 @@ end
 
 @inline _in_φ(p::CylindricalPoint, φ::AbstractInterval) = _in_angular_interval_closed(p.φ, φ)
 
+@inline _isapprox_x(p::CylindricalPoint{T}, x::Real) where {T} = isapprox(p.r * cos(p.φ), x, atol = geom_atol_zero(T))
+
 @inline _in_x(p::CylindricalPoint, x::Real) = abs(p.r * cos(p.φ)) <= x
 @inline _in_x(p::CylindricalPoint, x::AbstractInterval) = p.r * cos(p.φ) in x
 
+@inline _isapprox_y(p::CylindricalPoint{T}, y::Real) where {T} = isapprox(p.r * sin(p.φ), y, atol = geom_atol_zero(T))
+
 @inline _in_y(p::CylindricalPoint, y::Real) = abs(p.r * sin(p.φ)) <= y
 @inline _in_y(p::CylindricalPoint, y::AbstractInterval) = p.r * sin(p.φ) in y
 

src/ConstructiveSolidGeometry/SurfacePrimitives/Rectangle.jl

@@ -0,0 +1,190 @@
+struct Rectangle{T,TN,TL,TW} <: AbstractSurfacePrimitive{T}
+    normal::TN
+    l::TL
+    w::TW
+    loc::T
+    function Rectangle( ::Type{T},
+                        normal::Union{Val{:x}, Val{:y}, Val{:z}},
+                        l::Union{T, <:AbstractInterval{T}},
+                        w::Union{T, <:AbstractInterval{T}},
+                        loc::T) where {T}
+        new{T,typeof(normal),typeof(l),typeof(w)}(normal,l,w,loc)
+    end
+end
+
+#Constructors
+Rectangle(b::Box{T}, normal::Val{:x}, x::Real) where {T} = Rectangle(T, normal, b.y, b.z, T(x))
+Rectangle(b::Box{T}, normal::Val{:y}, y::Real) where {T} = Rectangle(T, normal, b.x, b.z, T(y))
+Rectangle(b::Box{T}, normal::Val{:z}, z::Real) where {T} = Rectangle(T, normal, b.x, b.y, T(z))
+
+function Rectangle(;normal = Val(:z), lMin = -1, lMax = 1, wMin = -1, wMax = 1, loc = 0)
+    T = float(promote_type(typeof.((lMin, lMax, wMin, wMax, loc))...))
+    l = lMax == -lMin ? T(lMax) : T(lMin)..T(lMax)
+    w = wMax == -wMin ? T(wMax) : T(wMin)..T(wMax)
+    Rectangle(T, normal, l, w, T(loc))
+end
+
+Rectangle(normal, lMin, lMax, wMin, wMax, loc) = Rectangle(; normal = normal, lMin = lMin, lMax = lMax, wMin = wMin, wMax = wMax, loc = loc)
+
+function Rectangle(normal::Union{Val{:x}, Val{:y}, Val{:z}}, l::L, w::W, loc::C) where {L<:Real, W<:Real, C<:Real}
+    T = float(promote_type(L,W,C))
+    Rectangle(T, normal, T(l/2), T(w/2), T(loc))
+end
+
+@inline in(p::AbstractCoordinatePoint, r::Rectangle{<:Any, Val{:x}}) = begin
+    _in_y(p, r.l) && _in_z(p, r.w) && _isapprox_x(p, r.loc)
+end
+
+@inline in(p::AbstractCoordinatePoint, r::Rectangle{<:Any, Val{:y}}) = begin
+    _in_x(p, r.l) && _in_z(p, r.w) && _isapprox_y(p, r.loc)
+end
+
+@inline in(p::AbstractCoordinatePoint, r::Rectangle{<:Any, Val{:z}}) = begin
+    _in_x(p, r.l) && _in_y(p, r.w) && _eq_z(p, r.loc)
+end    
+
+get_l_limits(r::Rectangle) = (_left_linear_interval(r.l), _right_linear_interval(r.l))
+get_w_limits(r::Rectangle) = (_left_linear_interval(r.w), _right_linear_interval(r.w))
+
+function sample(r::Rectangle{T, Val{:x}}, Nsamps::NTuple{3,Int} = (2,2,2)) where {T}
+    lMin::T, lMax::T = get_l_limits(r)
+    wMin::T, wMax::T = get_w_limits(r)
+    samples = [
+        CartesianPoint{T}(r.loc,y,z)
+        for y in (Nsamps[2] ≤ 1 ? lMin : range(lMin, lMax, length = Nsamps[2]))
+        for z in (Nsamps[3] ≤ 1 ? wMin : range(wMin, wMax, length = Nsamps[3]))
+    ]
+end
+
+function sample(r::Rectangle{T, Val{:y}}, Nsamps::NTuple{3,Int} = (2,2,2)) where {T}
+    lMin::T, lMax::T = get_l_limits(r)
+    wMin::T, wMax::T = get_w_limits(r)
+    samples = [
+        CartesianPoint{T}(x,r.loc,z)
+        for x in (Nsamps[1] ≤ 1 ? lMin : range(lMin, lMax, length = Nsamps[1]))
+        for z in (Nsamps[3] ≤ 1 ? wMin : range(wMin, wMax, length = Nsamps[3]))
+    ]
+end
+
+function sample(r::Rectangle{T, Val{:z}}, Nsamps::NTuple{3,Int} = (2,2,2)) where {T}
+    lMin::T, lMax::T = get_l_limits(r)
+    wMin::T, wMax::T = get_w_limits(r)
+    samples = [
+        CartesianPoint{T}(x,y,r.loc)
+        for x in (Nsamps[1] ≤ 1 ? lMin : range(lMin, lMax, length = Nsamps[1]))
+        for y in (Nsamps[2] ≤ 1 ? wMin : range(wMin, wMax, length = Nsamps[2]))
+    ]
+end
+
+
+function get_r(r::Rectangle{T, Val{:x}}, φ::T) where {T}
+    atol = geom_atol_zero(T)
+    x::T = r.loc
+    yMin::T, yMax::T = get_l_limits(r)
+    sφ::T, cφ::T = sincos(φ)
+    if x * cφ < 0 return () end
+    tanφ::T = tan(φ)
+    yMin - atol ≤ x * tanφ ≤ yMax + atol ? (x / cφ,) : ()
+end
+
+function get_r(r::Rectangle{T, Val{:y}}, φ::T) where {T}
+    atol = geom_atol_zero(T)
+    y::T = r.loc
+    xMin::T, xMax::T = get_l_limits(r)
+    sφ::T, cφ::T = sincos(φ)
+    if y * sφ < 0 return () end
+    cotφ::T = cot(φ)
+    xMin - atol ≤ y * cotφ ≤ xMax + atol ? (y / sφ,) : ()
+end
+
+function sample(r::Union{Rectangle{T, Val{:x}}, Rectangle{T, Val{:y}}}, g::CylindricalTicksTuple{T}) where {T}
+    samples = [
+        CylindricalPoint{T}(R, φ, z)
+        for φ in g.φ
+        for R in get_r(r,φ)
+        for z in _get_ticks(g.z, get_w_limits(r)...)
+    ]
+end
+
+function get_r_ticks(r::Rectangle{T, Val{:z}}, g::CylindricalTicksTuple{T}, φ::T)::Vector{T} where {T}
+    xMin::T, xMax::T = get_l_limits(r)
+    yMin::T, yMax::T = get_w_limits(r)
+    sφ::T, cφ::T = sincos(φ)
+    tanφ::T = tan(φ)
+    cotφ::T = inv(tanφ)
+    R = sort!(filter!(R -> R > 0,
+        [
+            (xMin ≤ yMin * cotφ ≤ xMax ? yMin / sφ : T(NaN)),
+            (xMin ≤ yMax * cotφ ≤ xMax ? yMax / sφ : T(NaN)),
+            (yMin ≤ xMin * tanφ ≤ yMax ? xMin / cφ : T(NaN)),
+            (yMin ≤ xMax * tanφ ≤ yMax ? xMax / cφ : T(NaN))
+        ]
+    ))
+    ticks = if length(R) == 2
+        _get_ticks(g.r, R...)
+    elseif length(R) == 1
+        _get_ticks(g.r, T(0), R[1])
+    else
+        []
+    end
+end
+
+function sample(r::Rectangle{T, Val{:z}}, g::CylindricalTicksTuple{T}) where {T}
+    samples = [
+        CylindricalPoint{T}(R, φ, z)
+        for z in _get_ticks(g.z, r.loc, r.loc)
+        for φ in g.φ
+        for R in get_r_ticks(r, g, φ)
+    ]
+end
+
+function sample(r::Rectangle{T, Val{:x}}, g::CartesianTicksTuple{T}) where {T}
+    samples = [
+        CartesianPoint{T}(r.loc,y,z)
+        for y in _get_ticks(g.y, get_l_limits(r)...)
+        for z in _get_ticks(g.z, get_w_limits(r)...)
+    ]
+end
+
+function sample(r::Rectangle{T, Val{:y}}, g::CartesianTicksTuple{T}) where {T}
+    samples = [
+        CartesianPoint{T}(x,r.loc,z)
+        for x in _get_ticks(g.x, get_l_limits(r)...)
+        for z in _get_ticks(g.z, get_w_limits(r)...)
+    ]
+end
+
+function sample(r::Rectangle{T, Val{:z}}, g::CartesianTicksTuple{T}) where {T}
+    samples = [
+        CartesianPoint{T}(x,y,r.loc)
+        for x in _get_ticks(g.x, get_l_limits(r)...)
+        for y in _get_ticks(g.y, get_w_limits(r)...)
+    ]
+end
+
+function get_vertices(r::Rectangle{T, Val{:x}}) where {T}
+    lMin::T, lMax::T = get_l_limits(r)
+    wMin::T, wMax::T = get_w_limits(r)
+    (CartesianPoint{T}(r.loc, lMin, wMin),
+    CartesianPoint{T}(r.loc, lMax, wMin),
+    CartesianPoint{T}(r.loc, lMin, wMax),
+    CartesianPoint{T}(r.loc, lMax, wMax))
+end
+
+function get_vertices(r::Rectangle{T, Val{:y}}) where {T}
+    lMin::T, lMax::T = get_l_limits(r)
+    wMin::T, wMax::T = get_w_limits(r)
+    (CartesianPoint{T}(lMin, r.loc, wMin),
+    CartesianPoint{T}(lMax, r.loc, wMin),
+    CartesianPoint{T}(lMin, r.loc, wMax),
+    CartesianPoint{T}(lMax, r.loc, wMax))
+end
+
+function get_vertices(r::Rectangle{T, Val{:z}}) where {T}
+    lMin::T, lMax::T = get_l_limits(r)
+    wMin::T, wMax::T = get_w_limits(r)
+    (CartesianPoint{T}(lMin, wMin, r.loc),
+    CartesianPoint{T}(lMax, wMin, r.loc),
+    CartesianPoint{T}(lMin, wMax, r.loc),
+    CartesianPoint{T}(lMax, wMax, r.loc))
+end

src/ConstructiveSolidGeometry/SurfacePrimitives/SurfacePrimitives.jl

@@ -1,6 +1,7 @@
 include("ConalPlane.jl")
 include("ConeMantle.jl")
 include("CylindricalAnnulus.jl")
+include("Rectangle.jl")
 include("RegularPolygon.jl")
 include("RegularPrismMantle.jl")
 include("SphereMantle.jl")

src/ConstructiveSolidGeometry/VolumePrimitives/Box.jl

@@ -26,7 +26,7 @@ function Box(x::X, y::Y, z::Z) where {X<:Real, Y<:Real, Z<:Real}
     Box( T, T(x/2), T(y/2), T(z/2))
 end
 
-in(p::AbstractCoordinatePoint, b::Box{<:Any, <:Any, <:Any, <:Any}) =
+in(p::AbstractCoordinatePoint, b::Box) =
     _in_x(p, b.x) && _in_y(p, b.y) && _in_z(p, b.z)
 
 
@@ -43,6 +43,26 @@ get_x_limits(b::Box{T}) where {T} = (_left_linear_interval(b.x), _right_linear_i
 get_y_limits(b::Box{T}) where {T} = (_left_linear_interval(b.y), _right_linear_interval(b.y))
 get_z_limits(b::Box{T}) where {T} = (_left_linear_interval(b.z), _right_linear_interval(b.z))
 
+function get_decomposed_surfaces(b::Box{T}) where {T}
+    xMin::T, xMax::T = get_x_limits(b)
+    yMin::T, yMax::T = get_y_limits(b)
+    zMin::T, zMax::T = get_z_limits(b)
+    tol = geom_atol_zero(T)
+    if isapprox(xMin, xMax, atol = tol)
+        return AbstractSurfacePrimitive[Rectangle(b, Val(:x), xMin)]
+    elseif isapprox(yMin, yMax, atol = tol)
+        return AbstractSurfacePrimitive[Rectangle(b, Val(:y), yMin)]
+    elseif isapprox(zMin, zMax, atol = tol)
+        return AbstractSurfacePrimitive[Rectangle(b, Val(:z), zMin)]
+    else
+        return AbstractSurfacePrimitive[
+                                            Rectangle(b, Val(:x), xMin), Rectangle(b, Val(:x), xMax), 
+                                            Rectangle(b, Val(:y), yMin), Rectangle(b, Val(:y), yMax), 
+                                            Rectangle(b, Val(:z), zMin), Rectangle(b, Val(:z), zMax)
+                                        ]
+    end
+end
+
 function sample(b::Box{T}, Nsamps::NTuple{3,Int} = (2,2,2)) where {T}
     xMin::T, xMax::T = get_x_limits(b)
     yMin::T, yMax::T = get_y_limits(b)

src/ConstructiveSolidGeometry/plotting/SurfacePrimitives/ConeMantle.jl

@@ -24,17 +24,18 @@ end
 function mesh(c::ConeMantle{T}; n = 30) where {T <: AbstractFloat}
 
     rbot::T, rtop::T = get_r_limits(c)
-    φMin::T, φMax::T, φ_is_full_2π::Bool = get_φ_limits(c)
+    φMin::T, φMax::T, _ = get_φ_limits(c)
     zMin::T, zMax::T = get_z_limits(c)
     f = (φMax - φMin)/(2π)
     n = Int(ceil(n*f))
     m = (rtop-rbot)/(zMax-zMin)
-    φ = range(φMin, φMax, length = n+1)
+    φrange = range(φMin, φMax, length = n+1)
+    scφrange = sincos.(φrange)
     z = range(zMin, zMax, length = 2)
 
-    X::Array{T,2} = [(m*(z[j]-zMin)+rbot)*cos(φ_i) for φ_i in φ, j in 1:length(z)]
-    Y::Array{T,2} = [(m*(z[j]-zMin)+rbot)*sin(φ_i) for φ_i in φ, j in 1:length(z)]
-    Z::Array{T,2} = [j for i in 1:length(φ), j in z]
+    X::Array{T,2} = [(m*(z[j]-zMin)+rbot)*cφ for (_,cφ) in scφrange, j in 1:length(z)]
+    Y::Array{T,2} = [(m*(z[j]-zMin)+rbot)*sφ for (sφ,_) in scφrange, j in 1:length(z)]
+    Z::Array{T,2} = [j for i in φrange, j in z]
 
     Mesh(X, Y, Z)
 end

src/ConstructiveSolidGeometry/plotting/SurfacePrimitives/Rectangle.jl

@@ -0,0 +1,9 @@
+function get_plot_points(r::Rectangle{T}; n = 30) where {T <: AbstractFloat}
+    v = get_vertices(r)
+    [Vector{CartesianPoint{T}}([v[1], v[2]]), Vector{CartesianPoint{T}}([v[2], v[4]]), Vector{CartesianPoint{T}}([v[3], v[4]]), Vector{CartesianPoint{T}}([v[3], v[1]])]
+end
+
+function mesh(r::Rectangle{T}; n = 30) where {T <: AbstractFloat}
+    vertices = [get_vertices(r)...]
+    Mesh(reshape(map(p -> p.x, vertices), (2,2)),reshape(map(p -> p.y, vertices), (2,2)),reshape(map(p -> p.z, vertices), (2,2)))
+ end

src/ConstructiveSolidGeometry/plotting/SurfacePrimitives/SurfacePrimitives.jl

@@ -1,6 +1,7 @@
 include("ConalPlane.jl")
 include("ConeMantle.jl")
 include("CylindricalAnnulus.jl")
+include("Rectangle.jl")
 include("RegularPolygon.jl")
 include("RegularPrismMantle.jl")
 include("SphereMantle.jl")

test/ConstructiveSolidGeometry/CSG_test.jl

@@ -8,7 +8,7 @@ using StaticArrays
 using SolidStateDetectors.ConstructiveSolidGeometry:
     CartesianPoint, CartesianVector, CylindricalPoint, scale,
     _in_angular_interval_closed, _in_angular_interval_open,
-    ConalPlane, ConeMantle, CylindricalAnnulus, RegularHexagon, HexagonalPrismMantle, RegularPolygon, RegularPrismMantle,
+    ConalPlane, ConeMantle, CylindricalAnnulus, Rectangle, RegularHexagon, HexagonalPrismMantle, RegularPolygon, RegularPrismMantle,
     SphereMantle, ToroidalAnnulus, TorusMantle,
     Tube, Cone, Torus, Box, Sphere, HexagonalPrism
 
@@ -177,6 +177,17 @@ using SolidStateDetectors.ConstructiveSolidGeometry:
                     @test !(CartesianPoint{T}(1, 0, -0.1) in prism)
                     @test !(CartesianPoint{T}(1, 0, 1.1) in prism)
                 end
+                @testset "Rectangle" begin
+                    rectangle = Rectangle(Val(:z), T(1.0), T(2.0), T(3.0)) # x from -0.5..0.5, y from -1.0..1.0, z = 3
+                    @test !(CartesianPoint{T}(0, 0, 0) in rectangle)
+                    @test CartesianPoint{T}(rectangle.l, rectangle.w, rectangle.loc) in rectangle
+                    @test CartesianPoint{T}(-rectangle.l, -rectangle.w, rectangle.loc) in rectangle
+                    @test CartesianPoint{T}(rectangle.l, -rectangle.w, rectangle.loc) in rectangle
+                    @test CartesianPoint{T}(-rectangle.l, rectangle.w, rectangle.loc) in rectangle
+                    @test CartesianPoint{T}(0, rectangle.w, rectangle.loc) in rectangle
+                    @test CartesianPoint{T}(-rectangle.l, 0, rectangle.loc) in rectangle
+                    @test !(CartesianPoint{T}(rectangle.l, -rectangle.w, 1.1*rectangle.loc) in rectangle)
+                end
             end
 
             @testset "Volume Primitives" begin

test/comparison_to_analytic_solutions.jl

@@ -1,7 +1,7 @@
 ϵ0 = SolidStateDetectors.ϵ0 * u"F / m"
 e = SolidStateDetectors.elementary_charge * u"C"
 
-#=
+
 @testset "Infinite Parallel Plate Capacitor" begin
     sim = Simulation{T}(SSD_examples[:InfiniteParallelPlateCapacitor])
     calculate_electric_potential!(sim, 
@@ -10,7 +10,7 @@ e = SolidStateDetectors.elementary_charge * u"C"
     )
     calculate_electric_field!(sim)
     BV_true = SolidStateDetectors._get_abs_bias_voltage(sim.detector) 
-    Δd = (sim.detector.contacts[2].geometry.x[2] - sim.detector.contacts[1].geometry.x[2]) * u"m"
+    Δd = (sim.detector.contacts[2].decomposed_surfaces[1].loc - sim.detector.contacts[1].decomposed_surfaces[1].loc) * u"m"
     Δx = (sim.detector.world.intervals[2].right - sim.detector.world.intervals[2].left) * u"m"
     A = Δx * Δx 
     V = A * Δd
@@ -24,7 +24,7 @@ e = SolidStateDetectors.elementary_charge * u"C"
         @test isapprox(C_ssd, C_true, rtol = 0.001) 
     end
 end
-=#
+
 
 struct DummyImpurityDensity{T} <: SolidStateDetectors.AbstractImpurityDensity{T} end
 

test/runtests.jl

@@ -69,19 +69,19 @@ end
         end
         @test isapprox( signalsum, T(2), atol = 5e-4 )
     end
-    #=
     @testset "Simulate example detector: CGD" begin
         sim = Simulation(SSD_examples[:CGD])
         SolidStateDetectors.apply_initial_state!(sim, ElectricPotential)
-        # simulate!(sim, max_refinements = 0, verbose = true)
-        # evt = Event([CartesianPoint{T}(5e-3, 5e-3, 5e-3)])
-        # simulate!(evt, sim)
-        # signalsum = T(0)
-        # for i in 1:length(evt.waveforms)
-        #     signalsum += abs(evt.waveforms[i].value[end])
-        # end
-        # @test isapprox( signalsum, T(2), atol = 1e-2 )
+        simulate!(sim, max_refinements = 0, verbose = true)
+        evt = Event([CartesianPoint{T}(5e-3, 5e-3, 5e-3)])
+        simulate!(evt, sim)
+        signalsum = T(0)
+        for i in 1:length(evt.waveforms)
+            signalsum += abs(evt.waveforms[i].value[end])
+        end
+        @test isapprox( signalsum, T(2), atol = 1e-2 )
     end
+    #=
     @testset "Simulate example detector: Spherical" begin
         sim = Simulation(SSD_examples[:Spherical])
         SolidStateDetectors.apply_initial_state!(sim, ElectricPotential)