Skip to content

Commit

Permalink
Merge pull request #451 from JuliaPhysics/polycone
Browse files Browse the repository at this point in the history
Add `Polycone`
  • Loading branch information
fhagemann authored Jan 24, 2025
2 parents 286ba04 + b92f3e2 commit 22bcccb
Show file tree
Hide file tree
Showing 13 changed files with 273 additions and 29 deletions.
2 changes: 1 addition & 1 deletion Project.toml
Original file line number Diff line number Diff line change
Expand Up @@ -63,7 +63,7 @@ LinearAlgebra = "<0.0.1, 1"
OrderedCollections = "1.1"
ParallelProcessingTools = "0.4"
Parameters = "0.12"
PolygonOps = "0.1"
PolygonOps = "0.1.1"
Polynomials = "2, 3, 4"
ProgressMeter = "1.5"
RadiationDetectorSignals = "0.3.5"
Expand Down
16 changes: 16 additions & 0 deletions docs/src/man/csg.md
Original file line number Diff line number Diff line change
Expand Up @@ -73,6 +73,22 @@ cone = CSG.Geometry(T, cfn)
plot(cone)
````

#### Polycone


````@example primitives
cfn = joinpath(path_to_example_primitives_config_files, "Polycone.yaml")
print(open(f -> read(f, String), cfn))
````

Load the primitive from the configuration file via `CSG.Geometry`

````@example primitives
polycone = CSG.Geometry(T, cfn)
plot(polycone)
````


### Ellipsoid
#### Sphere

Expand Down
Original file line number Diff line number Diff line change
@@ -1,26 +1,7 @@
geometry:
difference:
- tube:
r: 35
h: 80
origin:
z: 40
- cone:
r:
bottom:
from: 35
to: 36
top:
from: 23.71
to: 36
h: 64
origin:
z: 52
- tube:
r: 5
h: 80
origin:
z: 65
polycone:
r: [0, 35, 35, 24.42, 5, 5, 0]
z: [0, 0, 20, 80, 80, 25, 25]
contacts:
- material: HPGe
id: 1
Expand Down
3 changes: 3 additions & 0 deletions examples/example_primitive_files/Polycone.yaml
Original file line number Diff line number Diff line change
@@ -0,0 +1,3 @@
polycone:
r: [0, 35, 35, 25, 5, 5, 0, 0]
z: [0, 0, 20, 80, 80, 25, 25, 0]
32 changes: 31 additions & 1 deletion ext/Geant4/io_gdml.jl
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
using SolidStateDetectors
using SolidStateDetectors: Cylindrical, Cartesian, World, Simulation, SSDFloat
using SolidStateDetectors.ConstructiveSolidGeometry: CSGUnion, CSGDifference, CSGIntersection,
Box, Cone, Ellipsoid, Torus, RegularPrism, show_CSG_tree, AbstractVolumePrimitive, AbstractGeometry,
Box, Cone, Ellipsoid, Torus, RegularPrism, Polycone, AbstractVolumePrimitive, AbstractGeometry,
AbstractConstructiveGeometry, CylindricalPoint, origin, rotation
using IntervalSets
using OrderedCollections: OrderedDict
Expand Down Expand Up @@ -171,6 +171,14 @@ function has_volume(e::RegularPrism, v::Bool = false)
return true
end

function has_volume(p::Polycone, v::Bool = false)
if isapprox(PolygonOps.area(tuple.(p.r, p.z)), 0)
v && @warn "Polycone: The points passed result in a zero-volume Polycone"
return false
end
return true
end


# Check whether Boolean solid has any volume at all
function has_volume(e::AbstractConstructiveGeometry, v::Bool = false)
Expand Down Expand Up @@ -374,6 +382,28 @@ function parse_geometry(e::RegularPrism{T, <:Any, N}, x_solids::XMLElement, x_de
end
end

function parse_geometry(p::Polycone, x_solids::XMLElement, x_define::XMLElement, id::Integer, pf::AbstractString, v::Bool)::Nothing
if has_volume(p, v)
y = new_child(x_solids, "genericPolycone")

for i in eachindex(p.r)
rzpoint = new_child(y, "rzpoint")
set_attributes(rzpoint, OrderedDict(
"r" => p.r[i],
"z" => p.z[i]
))
end

set_attributes(y, OrderedDict(
"name" => pf * string(id),
"startphi" => 0,
"deltaphi" => 360,
"lunit" => SolidStateDetectors.internal_length_unit,
"aunit" => "deg"
))
end
end


function parse_geometry(e::World{T, 3, Cylindrical}, x_solids::XMLElement, x_define::XMLElement, id::Integer, pf::AbstractString, v::Bool)::Nothing where {T}
w = new_child(x_solids, "tube")
Expand Down
1 change: 1 addition & 0 deletions ext/SolidStateDetectorsGeant4Ext.jl
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,7 @@ using Geant4.SystemOfUnits
using LightXML
using LinearAlgebra
using Parameters
using PolygonOps
using ProgressMeter
using RadiationDetectorSignals
using StaticArrays
Expand Down
6 changes: 3 additions & 3 deletions src/ConstructiveSolidGeometry/ConstructiveSolidGeometry.jl
Original file line number Diff line number Diff line change
Expand Up @@ -50,9 +50,9 @@ module ConstructiveSolidGeometry

abstract type AbstractConstructiveGeometry{T} <: AbstractGeometry{T} end

_csg_convert_args(eltype::Type{T}, r::Real) where T = convert(T, r)
_csg_convert_args(eltype::Type{T}, r::Tuple) where T = broadcast(x -> _csg_convert_args(T, x), r)
_csg_convert_args(eltype::Type{T}, r::Nothing) where T = nothing
_csg_convert_args(::Type{T}, r::Real) where T = convert(T, r)
_csg_convert_args(::Type{T}, r::Tuple) where T = broadcast(x -> _csg_convert_args(T, x), r)
_csg_convert_args(::Type{T}, r::Nothing) where T = nothing

_csg_get_promoted_eltype(::Type{T}) where {T <: AbstractArray} = eltype(T)
_csg_get_promoted_eltype(::Type{T}) where {T <: Real} = T
Expand Down
3 changes: 2 additions & 1 deletion src/ConstructiveSolidGeometry/IO.jl
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,7 @@ const CSG_dict = Dict{String, Any}(
"sphere" => Ellipsoid,
"box" => Box,
"torus" => Torus,
"polycone" => Polycone,
"TriangularPrism" => TriangularPrism,
"QuadranglePrism" => QuadranglePrism,
"PentagonalPrism" => PentagonalPrism,
Expand Down Expand Up @@ -59,7 +60,7 @@ end
T(to_internal_units(float(x)))
end
@inline _parse_value(::Type{T}, s::String, unit::Unitful.Units) where {T} = _parse_value(T, uparse(s), unit)
@inline _parse_value(::Type{T}, a::AbstractVector, unit::Unitful.Units) where {T} = _parse_value.(T, a, unit)
@inline _parse_value(::Type{T}, a::Union{<:AbstractVector, <:Tuple}, unit::Unitful.Units) where {T} = _parse_value.(T, a, unit)

# parses dictionary entries of type {"from": ..., "to": ... } to a Tuple of the interval boundaries
function _parse_interval_from_to(::Type{T}, dict::AbstractDict, unit::Unitful.Units)::Tuple{T,T} where {T}
Expand Down
184 changes: 184 additions & 0 deletions src/ConstructiveSolidGeometry/VolumePrimitives/Polycone.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,184 @@
"""
struct Polyone{T,CO,N,TP} <: AbstractVolumePrimitive{T, CO}
Volume primitive describing a [Polycone](@ref), similar to the G4Polycone defined in Geant4.
## Parametric types
* `T`: Precision type.
* `CO`: Describes whether the surface belongs to the primitive.
It can be `ClosedPrimitive`, i.e. the surface points belong to the primitive,
or `OpenPrimitive`, i.e. the surface points do not belong to the primitive.
* `N`: Integer describing the number of corners.
* `TP`: Type of the angular range `φ`.
* `TP == Nothing`: Full 2π Cone.
## Fields
* `r::NTuple{N,T}`: `r`-coordinates of the corners of the polycone.
* `z::NTuple{N,T}`: `z`-coordinates of the corners of the polycone.
* `φ::TP`: Range in polar angle `φ` over which the `Cone` extends (in radians).
* `origin::CartesianPoint{T}`: The position of the center of the `Cone` at the middle height.
* `rotation::SMatrix{3,3,T,9}`: Matrix that describes a rotation of the `Cone` around its `origin`.
## Definition in Configuration File
A `Polycone` is defined in the configuration file as part of the `geometry` field
of an object through the field `polycone`.
Example definitions of a polycone looks like this:
```yaml
polycone:
r: [0, 35, 35, 24.42, 5, 5, 0, 0]
z: [0, 0, 20, 80, 80, 25, 25, 0]
origin:
z: 1.0 # => origin = [0.0, 0.0, 1.0]
```
This is a polycone describing the semiconductor of the example inverted coaxial detector.
See also [Constructive Solid Geometry (CSG)](@ref).
"""
struct Polycone{T,CO,N,TP<:Nothing} <: AbstractVolumePrimitive{T,CO}
r::NTuple{N,T}
z::NTuple{N,T}
φ::TP

origin::CartesianPoint{T}
rotation::SMatrix{3,3,T,9}

function Polycone{T,CO}(r::Union{<:AbstractArray, <:Tuple}, z::Union{<:AbstractArray, <:Tuple}, φ, origin, rotation) where {T, CO}
nr = length(r)
nz = length(z)
nr != nz && throw(ConfigFileError("In PolyCone: r and z must have the same length."))
# convert to Tuple and determine the length N
_r, _z = if first(r) != last(r) || first(z) != last(z)
nr += 1
nz += 1
(r..., first(r)), (z..., first(z))
else
(r...,), (z...,)
end
# sort the points counter-clockwise in the r-z-plane
if PolygonOps.area(tuple.(_r,_z)) < 0
_r, _z = reverse(_r), reverse!(_z)
end
new{T,CO,nr,typeof(φ)}(_r, _z, φ, origin, rotation)
end
end

#Type conversion happens here
function Polycone{T}(CO, r, z, φ, origin, rotation) where {T}
_r = _csg_convert_args.(T, r)
_z = _csg_convert_args.(T, z)
= _csg_convert_args(T, φ)
Polycone{T,CO}(_r, _z, _φ, origin, rotation)
end

#Type promotion happens here
function Polycone(CO, r::TR, z::TZ, φ::TP, origin::PT, rotation::ROT) where {TR, TZ, TP, PT, ROT}
eltypes = _csg_get_promoted_eltype.((TR, TZ, TP, PT, ROT))
T = float(promote_type(eltypes...))
Polycone{T}(CO, r, z, φ, origin, rotation)
end

function Polycone(::Type{CO}=ClosedPrimitive;
r = (0,1,1,0),
z = (0,0,1,1),
φ = nothing,
origin = zero(CartesianPoint{Int}),
rotation = one(SMatrix{3, 3, Int, 9})
) where {CO}
Polycone(CO, r, z, φ, origin, rotation)
end

function Polycone{T}(::Type{CO}=ClosedPrimitive;
r = (0.,1.,1.,0.),
z = (0.,0.,1.,1.),
φ = nothing,
origin = zero(CartesianPoint{Float64}),
rotation = one(SMatrix{3, 3, Float64, 9})
) where {T, CO}
Polycone{T}(CO, r, z, φ, origin, rotation)
end

Polycone{T,CO,N,TP}( c::Polycone{T,CO,N,TP}; COT = CO,
origin::CartesianPoint{T} = c.origin,
rotation::SMatrix{3,3,T,9} = c.rotation) where {T, CO<:Union{ClosedPrimitive, OpenPrimitive}, N, TP} =
Polycone{T,COT}(c.r, c.z, c.φ, origin, rotation)


####################################################################
####################################################################


function _inpolygon(pt::Tuple{T,T}, polygon::NTuple{N,Tuple{T,T}}; csgtol::T = csg_default_tol(T)) where {N,T <: Real}
# use the normal PolygonOps.inpolygon method to check if the point is inside
PolygonOps.inpolygon(pt, polygon, in = true, on = true, out = false) && return true
# if it is flagged as outside, check if the distance to any of the edges is smaller than csgtol
rp,zp = pt
@inbounds for i in Base.OneTo(N-1)
r1,z1 = polygon[i]
r2,z2 = polygon[i+1]
r1 == r2 && z1 == z2 && continue
s = clamp(((zp - z1) * (z2 - z1) + (rp - r1) * (r2 - r1)) / ((z2 - z1)^2 + (r2 - r1)^2), zero(T), one(T))
if hypot(rp - r1 - s * (r2 - r1), zp - z1 - s * (z2 - z1)) <= csgtol
return true
end
end
# if not, then it is really outside
return false
end

function _in(pt::CartesianPoint{T}, c::Polycone{<:Any, ClosedPrimitive}; csgtol::T = csg_default_tol(T)) where {T}
return (isnothing(c.φ) || _in_angular_interval_closed(atan(pt.y, pt.x), c.φ, csgtol = csgtol / r)) &&
_inpolygon((hypot(pt.x, pt.y), pt.z), tuple.(c.r, c.z); csgtol)
end

function _in(pt::CartesianPoint{T}, c::Polycone{<:Any, OpenPrimitive}; csgtol::T = csg_default_tol(T)) where {T}
return (isnothing(c.φ) || _in_angular_interval_open(atan(pt.y, pt.x), c.φ, csgtol = csgtol / r)) &&
PolygonOps.inpolygon((hypot(pt.x, pt.y), pt.z), tuple.(c.r, c.z), in = true, on = iszero(pt.x) && iszero(pt.y), out = false)
end

function surfaces(c::Polycone{T,<:Any,N,Nothing}) where {T,N}
s = []
@inbounds for i in Base.OneTo(N-1)
r1::T = c.r[i]
r2::T = c.r[i+1]
## skip the surfaces on the z-axis
iszero(r1) && iszero(r2) && continue
z1::T = c.z[i]
z2::T = c.z[i+1]
origin = _transform_into_global_coordinate_system(CartesianPoint{T}(zero(T), zero(T), (z1+z2)/2), c)
vol = if z1 == z2
EllipticalSurface{T}(r = r1 <= r2 ? (r1,r2) : (r2,r1), origin = origin, rotation = r1 <= r2 ? c.rotation : -c.rotation * RotZ{T}(π))
else
FullConeMantle{T, z1 < z2 ? (:inwards) : (:outwards)}(z1 < z2 ? (r1,r2) : (r2,r1), c.φ, abs(z1-z2)/2, origin, c.rotation)
end
push!(s, vol)
end
(s...,)
end


function Geometry(::Type{T}, ::Type{Polycone}, dict::AbstractDict, input_units::NamedTuple, transformations::Transformations{T}) where {T}
length_unit = input_units.length
origin = get_origin(T, dict, length_unit)
rotation = get_rotation(T, dict, input_units.angle)
haskey(dict, "r") || throw(ConfigFileError("Polycone needs entry `r`."))
haskey(dict, "z") || throw(ConfigFileError("Polycone needs entry `z`."))
r = _parse_value(T, dict["r"], length_unit)
z = _parse_value(T, dict["z"], length_unit)
polycone = Polycone{T, ClosedPrimitive}(r, z, nothing, origin, rotation)
transform(polycone, transformations)
end


function Dictionary(c::Polycone{T})::OrderedDict{String, Any} where {T}
dict = OrderedDict{String, Any}()
@assert isnothing(c.φ) "Polycone needs `φ` field to be nothing."
dict["r"] = c.r
dict["z"] = c.z
if c.origin != zero(CartesianVector{T}) dict["origin"] = c.origin end
if c.rotation != one(SMatrix{3,3,T,9}) dict["rotation"] = Dictionary(c.rotation) end
OrderedDict{String, Any}("polycone" => dict)
end

extremum(c::Polycone{T}) where {T} = hypot(max(abs.(extrema(c.r))...), max(abs.(extrema(c.z))...))
Original file line number Diff line number Diff line change
@@ -1,3 +1,9 @@
# All volume primitives need to have the following methods defined:
# * _in(::CartesianPoint, ::VolumePrimitive) (in the local coordinate system)
# * surfaces(::VolumePrimitive) --> all surfaces of the volume primitive as Tuple
# * Geometry(...) --> a method to read the geometry from a config file
# * Dictionary(...) --> a method to convert the primitive to a config

ClosedPrimitive(p::VP) where {VP <:AbstractVolumePrimitive} = VP(p, COT = ClosedPrimitive)
OpenPrimitive(p::VP) where {VP <: AbstractVolumePrimitive} = VP(p, COT = OpenPrimitive)
switchClosedOpen(p::AbstractVolumePrimitive{<:Any,ClosedPrimitive}) = OpenPrimitive(p)
Expand All @@ -19,4 +25,5 @@ include("Box.jl")
include("Cone.jl")
include("Ellipsoid.jl")
include("RegularPrism.jl")
include("Polycone.jl")
include("Torus.jl")
5 changes: 5 additions & 0 deletions test/ConstructiveSolidGeometry/CSG_IO.jl
Original file line number Diff line number Diff line change
Expand Up @@ -30,6 +30,11 @@ example_primitive_dir = joinpath(@__DIR__, "../../examples/example_primitive_fil
@test ellipsoid_full_sphere isa CSG.FullSphere{T}
end

@testset "Polycone" begin
polycone = Geometry(T, joinpath(example_primitive_dir, "Polycone.yaml"))
@test polycone isa CSG.Polycone{T}
end

@testset "RegularPrism" begin
hexagon = Geometry(T, joinpath(example_primitive_dir, "RegularPrism_hexagon.yaml"))
@test hexagon isa CSG.HexagonalPrism{T}
Expand Down
16 changes: 16 additions & 0 deletions test/ConstructiveSolidGeometry/CSG_primitives.jl
Original file line number Diff line number Diff line change
Expand Up @@ -164,6 +164,22 @@ no_translations = (rotation = one(SMatrix{3, 3, T, 9}), translation = zero(Carte
rot_torus = @inferred CSG.Torus{Float64}=π/2=nothing,rotation=SMatrix{3}(0,0,-1,0,1,0,1,0,0) * SMatrix{3}(cos/4),sin/4),0,-sin/4),cos/4),0,0,0,1))
@test tuple_torus ==rot_torus
end
@testset "Polycone" begin
polycone1 = CSG.Polycone(CSG.ClosedPrimitive,r=Float32[0,2,4,0,0],z=Float32[0,1,2,3,0],origin=zero(CartesianPoint{Float16}),rotation=one(SMatrix{3, 3, Float16, 9}))
polycone2 = CSG.Polycone{Float32}(r = [0,2,4,0], z = [0,1,2,3]) # omit last entry
@test polycone1 == polycone2

dict = Dict("polycone" => Dict(
"r" => (0,2,4,0,0),
"z" => (0,1,2,3,0),
))
polycone3 = Geometry(T,dict,default_units,no_translations)
@test dict == Dictionary(polycone3)

polycone_open = CSG.Polycone(CSG.OpenPrimitive,r=Float32[0,2,4,0,0],z=Float32[0,1,2,3,0],origin=zero(CartesianPoint{Float16}),rotation=one(SMatrix{3, 3, Float16, 9}))
@test in(CartesianPoint{Float32}(4,0,2), polycone1)
@test !in(CartesianPoint{Float32}(4,0,2), polycone_open)
end
@testset "Ellipsoid" begin
ellip1 = @inferred CSG.Ellipsoid(CSG.ClosedPrimitive,r=1f0,origin = zero(CartesianPoint{Float16}),rotation = one(SMatrix{3, 3, Float16, 9}))
ellip2 = @inferred CSG.Ellipsoid{Float32}(r=1.0)
Expand Down
Loading

0 comments on commit 22bcccb

Please sign in to comment.