A number of fixes

Manifest.toml

@@ -882,6 +882,12 @@ git-tree-sha1 = "6d105d40e30b635cfed9d52ec29cf456e27d38f8"
 uuid = "f57f5aa1-a3ce-4bc8-8ab9-96f992907883"
 version = "0.3.12"
 
+[[deps.PackageCompiler]]
+deps = ["Artifacts", "LazyArtifacts", "Libdl", "Pkg", "Printf", "RelocatableFolders", "TOML", "UUIDs"]
+git-tree-sha1 = "c497e2bb9c2127a411b74dbff56b11f258d67d12"
+uuid = "9b87118b-4619-50d2-8e1e-99f35a4d4d9d"
+version = "2.0.9"
+
 [[deps.Packing]]
 deps = ["GeometryBasics"]
 git-tree-sha1 = "1155f6f937fa2b94104162f01fa400e192e4272f"

Project.toml

@@ -21,6 +21,7 @@ KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
 LightXML = "9c8b4983-aa76-5018-a973-4c85ecc9e179"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
+PackageCompiler = "9b87118b-4619-50d2-8e1e-99f35a4d4d9d"
 PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"

README.md

@@ -1,15 +1,17 @@
 # Geom4hep
 An exercise for evaluation of the Julia language for a Geometry modeling package used in HEP (VecGeom). The current functionality:
 
-- Basic shapes: `Box`, `Tube`, `Trd`, `Cone`
+- Basic shapes: `Box`, `Tube`, `Trd`, `Cone`, `CutTube`, `Polycone`, `Boolean` 
 - Materials: `Isotope`, `Element`, `Material`
 - Essential functions to navigate on them (`distanceToIn`, `distanceToOut`, ...)
-- Hierarchical geometry models with `Volume` and `PlacedVolumes` 
-- Very basic navigator (linear search among all daughters)
+- Construction of hierarchical geometry models with `Volume` and `PlacedVolumes` 
+- Very basic navigator `TrivialNavigator` (linear search among all daughters)
+- Bounding Volume Hierarchy accelerated navigator `BVHNavigator`
 - 3D graphics support for displaying the shapes and geometry models using the [Makie](https://makie.juliaplots.org/stable/) package
 - GDML reader to input geometries
 - Benchmark functions for shapes
 - A number of examples to demonstrate the available functionality 
+- Tests with using GPUs (CUDA.jl) 
 
 ## Running the unit tests
 A number of unit tests is provided under /test for each of the geometrical shapes. To run all of them do:

examples/DrawDetector.jl

@@ -29,7 +29,7 @@ world = processGDML("examples/boxes.gdml")
 using GLMakie
 fig = Figure()
 scene = LScene(fig[1, 1])
-draw(scene, world)
+draw!(scene, world)
 display(fig)
 
 

examples/XRay.jl

@@ -54,28 +54,34 @@ end
 
 using GLMakie
 
-if abspath(PROGRAM_FILE) == @__FILE__
+#if abspath(PROGRAM_FILE) == @__FILE__
     #-----build detector---------------------------------------
-    full = processGDML("examples/trackML.gdml", Float64)
-    volume = full.daughters[2].volume.daughters[1].volume
+    #full = processGDML("examples/trackML.gdml", Float64)
+    #volume = full[2,1]
+
+    full = processGDML("examples/cms2018.gdml", Float64)
+    volume = full[1,7,1]
+
     world = getWorld(volume)
     nav = BVHNavigator(world)
     #nav = TrivialNavigator(world)
+
     #-----Generate and Plot results----------------------------
     @printf "generating x-projection\n"
     rx = generateXRay(nav, world, 1e6, 1)
     @printf "generating y-projection\n"
     ry = generateXRay(nav, world, 1e6, 2)
     @printf "generating z-projection\n"
     rz = generateXRay(nav, world, 1e6, 3)
-    limits = (0, max(maximum(rx[3]), maximum(ry[3]), maximum(rz[3])))
+    limits = (min(minimum(rx[3]), minimum(ry[3]), minimum(rz[3])), max(maximum(rx[3]), maximum(ry[3]), maximum(rz[3])))
+
     #----Plot the results--------------------------------------
-    fig = Figure(resolution = (1000, 1000))
+    fig = Figure(resolution = (2000, 2000))
     @printf "ploting the results\n"
     heatmap!(Axis(fig[1, 1], title = "X direction"), rx..., colormap=:grayC, colorrange=limits)
     heatmap!(Axis(fig[2, 1], title = "Y direction"), ry..., colormap=:grayC, colorrange=limits)
-    heatmap!(Axis(fig[1, 2], title = "Z direction"), rz..., colormap=:grayC, colorrange=(0,maximum(rz[3])))
-    draw!(LScene(fig[2, 2]), volume; wireframe=true, maxlevel=3)
+    heatmap!(Axis(fig[1, 2], title = "Z direction"), rz..., colormap=:grayC, colorrange=limits)
+    draw!(LScene(fig[2, 2]), volume; wireframe=true, maxlevel=2)
     #display(fig)
-    save("trackML.png", fig)
-end
+    save("cms2018.png", fig)
+#end

src/BVH.jl

@@ -1,6 +1,9 @@
 #---Boxed volumes used to build the Bounding Volume Hierarchy---------------------------------------
 #   BoxedVolumes wraps around a set of PlacedVolumes indices, and the AABBs (axis aligned bounding 
 #   boxes) of the corresponding daughters, identified by min/max coordinates
+#
+#   Algorithm adapted from https://github.com/jonalm/BoundingVolumeHierarchies.jl
+#---------------------------------------------------------------------------------------------------
 
 struct BoxedVolumes{T<:AbstractFloat}
     min::Matrix{T}
@@ -95,6 +98,7 @@ BVH(b::AABB{T}, inds::Vector{Int64}) where T = BVH{T}(b, (), inds)
 function Base.show(io::IO, bvh::BVH{T}) where T
     print(io, "BVH{$T}",(aabb=bvh.aabb, leaves=_leaves(bvh)))
 end
+
 #--- _Branch and _Leaf are only used to construct the BVH recursively
 struct _Branch{T} data::T; end
 struct _Leaf{T} data::T; end
@@ -122,7 +126,6 @@ buildBVH(pvols::Vector{PlacedVolume{T}}) where T = buildBVH(pvols, BVHParam{T}()
 resetBVH(bvh::BVH, vertices, faces) = nothing
 
 #---Returns an iterator over all subsets of PlacedVolume indices such that `f(x::AABB) == true` for all parent axis aligned bounding boxes `x`.
-
 function _pvolindices!(ind::Vector{Int64}, f::Function, b::BVH{T}) where T
     isempty(b.children) && return append!(ind, b.indices)
     for c in b.children

src/Boolean.jl

@@ -10,9 +10,9 @@ function Boolean(op::Symbol, left::AbstractShape{T}, right::AbstractShape{T}, pl
     @assert op in (:Union, :Subtraction, :Intersection) "Boolean suppoted operations are: :Union, :Subtraction, :Intersection"
     Boolean{T,typeof(left),typeof(right),op}(left, right, place)
 end
-#function BooleanUnion(left::AbstractShape{T}, right::AbstractShape{T}, place::Transformation3D{T}=one(Transformation3D{T})) where T<:AbstractFloat
-#    Boolean{T,typeof(left),typeof(right),:Union}(left, right, place)
-#end
+function BooleanUnion(left::AbstractShape{T}, right::AbstractShape{T}, place::Transformation3D{T}=one(Transformation3D{T})) where T<:AbstractFloat
+    Boolean{T,typeof(left),typeof(right),:Union}(left, right, place)
+end
 function Subtraction(left::AbstractShape{T}, right::AbstractShape{T}, place::Transformation3D{T}=one(Transformation3D{T})) where T<:AbstractFloat
     Boolean{T,typeof(left),typeof(right),:Subtraction}(left, right, place)
 end
@@ -29,17 +29,10 @@ function Base.show(io::IO, shape::Boolean{T, SL, SR, OP}) where {T,SL,SR,OP}
     print(io, "Boolean{$OP}",(left=left, right=right, placement=placement))
 end
 
-function GeometryBasics.coordinates(shape::Boolean{T, SL, SR, OP}, facets=36) where {T,SL,SR,OP}
+function GeometryBasics.mesh(shape::Boolean{T, SL, SR, OP}) where {T,SL,SR,OP}
     (; left, right, transformation) = shape
-    coors = (coordinates(left, facets), coordinates(right, facets))
-    transform = (one(Transformation3D{T}), transformation)
-    return (c * transform[i] for i in 1:2 for c in coors[i])
-end
-
-function GeometryBasics.faces(shape::Boolean{T, SL, SR, OP}, facets=36) where {T,SL,SR,OP}
-    (; left, right) = shape
-    offset = length(coordinates(left, facets))
-    append!(faces(left, facets), [t .+ offset for t in faces(right, facets)])
+    rmesh = mesh(right)
+    merge([mesh(left), Mesh(map(c -> Point3{T}(c * transformation), coordinates(rmesh)), faces(rmesh))])
 end
 
 #---Basic functions---------------------------------------------------------------------------------
@@ -133,13 +126,15 @@ function distanceToOut(shape::Boolean{T, SL, SR, :Union}, point::Point3{T}, dir:
     if positionA != kOutside  # point inside A
         while(true)
             distA = distanceToOut(left, point, dir)
-            dist += (distA > 0 && distA < Inf) ? distA : 0 + kPushTolerance(T)
+            dist += (distA > 0 && distA < Inf) ? distA : 0
+            dist += kPushTolerance(T)
             npoint = point + dist * dir
             lpoint = transformation * npoint
             if inside(right, lpoint) != kOutside # B could be overlapping with A -- and/or connecting A to another part of A
                 ldir = transformation * dir
                 distB = distanceToOut(right, lpoint, ldir)
-                dist += (distB > 0 && distB < Inf) ? distB : 0 + kPushTolerance(T)
+                dist += (distB > 0 && distB < Inf) ? distB : 0
+                dist += kPushTolerance(T)
                 npoint = point + dist * dir
                 if inside(left, npoint) == kOutside
                     break
@@ -156,12 +151,14 @@ function distanceToOut(shape::Boolean{T, SL, SR, :Union}, point::Point3{T}, dir:
             ldir = transformation * dir
             while(true)
                 distB = distanceToOut(right, lpoint, ldir)
-                dist += (distB > 0 && distB < Inf) ? distB : 0 + kPushTolerance(T)
+                dist += (distB > 0 && distB < Inf) ? distB : 0
+                dist += kPushTolerance(T) # Give a push
                 npoint = point + dist * dir
                 if inside(left, npoint) != kOutside # A could be overlapping with B -- and/or connecting B to another part of B
                     ldir = transformation * dir
                     distA = distanceToOut(left, npoint, dir)
-                    dist += (distA > 0 && distA < Inf) ? distA : 0 + kPushTolerance(T)
+                    dist += (distA > 0 && distA < Inf) ? distA : 0
+                    dist += kPushTolerance(T)
                     npoint = point + dist * dir
                     lpoint = transformation * npoint
                     if inside(right, lpoint) == kOutside
@@ -267,7 +264,7 @@ function distanceToIn(shape::Boolean{T, SL, SR, :Subtraction}, point::Point3{T},
             # propagate to outside of '- / RightShape'
             d1 = distanceToOut(right, lpoint, ldir)
             dist += (d1 >= 0 && d1 < Inf) ? d1 + kPushTolerance(T) : 0
-            npoint = point + dist * dir
+            npoint = point + (dist + kPushTolerance(T)) * dir
             lpoint = transformation * npoint
             # now master outside 'B'; check if inside 'A'
             inside(left, npoint) == kInside && distanceToOut(left, npoint, dir) > kPushTolerance(T) && return dist
@@ -285,8 +282,8 @@ function distanceToIn(shape::Boolean{T, SL, SR, :Subtraction}, point::Point3{T},
         end
 
         #   propagate to '-'
-        dist += (d1 >= 0 && d1 < Inf) ? d1 + kPushTolerance(T) : 0
-        point + dist * dir
+        dist += (d1 >= 0 && d1 < Inf) ? d1 : 0
+        npoint = point + (dist + kPushTolerance(T)) * dir
         lpoint = transformation * npoint
         inRight = true
     end

src/Box.jl

@@ -109,8 +109,8 @@ function GeometryBasics.coordinates(box::Box{T}, facets=6) where {T<:AbstractFlo
 end
 
 function GeometryBasics.faces(box::Box{T}, facets=6) where {T<:AbstractFloat}
-    return TriangleFace{Int64}[(1,2,4), (4,3,1), (1,3,7), (7,5,1), (1,5,6), (6,2,1), 
-                               (8,6,5), (5,7,8), (8,7,3), (3,4,8), (8,4,2), (2,6,8)] 
+    iface = ((1,5,6,2),(3,4,8,7),(1,3,7,5),(2,6,8,4),(1,2,4,3),(5,6,8,7))
+    (QuadFace{Int64}(f...) for f in iface)
 end
 
 function GeometryBasics.normals(box::Box{T}, facets=24) where {T<:AbstractFloat}

src/Cone.jl

@@ -501,41 +501,34 @@ function GeometryBasics.coordinates(cone::Cone{T}, facets=36) where {T<:Abstract
     nbf = Int(facets / 2)            # Number of faces
     nbv = issector ? nbf + 1 : nbf   # Number of vertices
     nbc = ishollow ? nbv : 1         # Number of centers
-
-    indexes = Vector{TriangleFace{Int}}()
+    indexes = Vector{QuadFace{Int64}}()
     for j in 1:nbf
         a,b = 2j-1, 2j
         c,d = !issector && j == nbf ? (1, 2) : (2j+1, 2j+2) 
-        push!(indexes, (a,b,d))
-        push!(indexes, (d,c,a))
+        push!(indexes, (a,b,d,c))
         if ishollow
             a′,b′ = 2j-1+2nbv, 2j+2nbv
             c′,d′ = !issector && j == nbf ? (2nbv+1, 2nbv+2) : (2j+1+2nbv, 2j+2+2nbv)
             # inner wall
-            push!(indexes, (a′,d′,b′))
-            push!(indexes, (d′,a′,c′))
+            push!(indexes, (a′,b′,d′,c′))
             # top
-            push!(indexes, (a, c ,a′))
-            push!(indexes, (c, c′,a′))
+            push!(indexes, (c, c′, a′, a))
             # bottom
-            push!(indexes, (b, b′, d))
-            push!(indexes, (b′,d′, d))
+            push!(indexes, (b, b′, d′, d))
         else
             a′,b′ = 2nbv+1, 2nbv+2
             # top
-            push!(indexes, (a′,a, c))
+            push!(indexes, (a′,a, c, c))
             # bottom
-            push!(indexes, (b′,d, b))
+            push!(indexes, (b′,d, b, b))
         end
     end
     if issector
         # wedge walls
         a, b, c, d  = ( 1, 2, 2nbv-1, 2nbv)
         a′,b′,c′,d′ = ishollow ? (2nbv+1, 2nbv+2, 4nbv-1, 4nbv ) : (2nbv+1, 2nbv+2, 2nbv+1, 2nbv+2)
-        push!(indexes, (a,  b, a′))
-        push!(indexes, (b,  b′,a′))
-        push!(indexes, (c′, d′,c ))
-        push!(indexes, (d′, d, c ))
+        push!(indexes, (a,  b, b′, a′))
+        push!(indexes, (c′, d′, d, c ))
     end
     return indexes
 end

src/Drawing.jl

@@ -4,21 +4,23 @@ using Printf
 
 import GeometryBasics: coordinates, faces, mesh, Mesh
 
-colors = colormap("Grays", 8)
+colors = colormap("Grays", 16)
 
 global drawn = 0
 
+GeometryBasics.mesh(shape::AbstractShape{T}) where T = GeometryBasics.mesh(Tesselation(shape, 64), facetype=typeof(first(GeometryBasics.faces(shape))))
+
 #---Draw a Volume---------------------------------------------------------------
 function draw!(s::LScene, vol::Volume{T}, t::Transformation3D{T}, level::Int64, wireframe::Bool, bvh::Bool, maxlevel::Int64) where T
-    m = GeometryBasics.mesh(Tesselation(vol.shape, 64), facetype=typeof(first(GeometryBasics.faces(vol.shape))))
+    m = GeometryBasics.mesh(vol.shape)
     if ! isone(t)
         points = GeometryBasics.coordinates(m)
         faces  = GeometryBasics.faces(m)
         map!(c -> c * t, points, points)
         m = GeometryBasics.Mesh(points, faces)
     end
     if wireframe
-        wireframe!(s, m, color=colors[9-level], visible = level == 1 ? false : true)
+        wireframe!(s, m, color=colors[16-level], visible = level == 1 ? false : true)
     else
         mesh!(s, m, color=colors[level], transparency=true, ambient=0.7, visible = level == 1 ? false : true)
     end
@@ -52,7 +54,7 @@ end
 
 #---Draw a Shape---------------------------------------------------------------
 function draw!(s::LScene, shape::AbstractShape; wireframe::Bool=false)
-    m = GeometryBasics.mesh(Tesselation(shape, 64), facetype=typeof(first(faces(shape))))
+    m = GeometryBasics.mesh(shape)
     if wireframe
         wireframe!(s, m)
     else

src/GDML.jl

@@ -262,7 +262,7 @@ end
 function getTransformation(dicts::GDMLDicts{T}, element::XMLElement) where T
     (; positions, rotations) = dicts
     position = Vector3{T}(0,0,0)
-    rotation = RotXYZ{T}(0,0,0)
+    rotation = RotZYX{T}(0,0,0)
     for c in child_nodes(element)
         if is_elementnode(c)
             aa = attributes_dict(XMLElement(c))
@@ -273,9 +273,9 @@ function getTransformation(dicts::GDMLDicts{T}, element::XMLElement) where T
                                       parse(T, aa["z"]) * unit)
             elseif name(c) == "rotation"
                 unit = eval(Meta.parse(aa["unit"]))
-                rotation = RotXYZ{T}(parse(T, aa["x"]) * unit, 
-                                     parse(T, aa["y"]) * unit,
-                                     parse(T, aa["z"]) * unit)
+                rotation = RotZYX{T}(parse(T, aa["z"]) * unit, 
+                                     parse(T, aa["y"]) * unit, 
+                                     parse(T, aa["x"]) * unit)
             elseif name(c) == "positionref"
                 position = positions[aa["ref"]]
             elseif name(c) == "rotationref"

src/Geom4hep.jl

@@ -3,7 +3,7 @@ module Geom4hep
 export Point3, Point2, Vector3, Vector2, nonzero
 export AbstractShape, AbstractMaterial
 export Shape, NoShape, Box, Trd, TBox, TTrd, Tube, Wedge, isInside, isOutside, Cone, Polycone, getNz, getNSections, getSectionIndex
-export CutTube, Plane, Boolean, Trap
+export CutTube, Plane, Boolean, Trap, BooleanUnion
 export getindex, capacity, surface, extent, normal, distanceToOut,  distanceToIn, inside, safetyToOut, safetyToIn
 export Material, Isotope, Element
 export Transformation3D, RotMatrix3, RotXYZ, one, isone, transform, hasrotation, hastranslation, inv, lmul!
@@ -13,7 +13,7 @@ export AbstractNavigator, TrivialNavigator, BVHNavigator, NavigatorState, comput
 export kTolerance
 export processGDML
 export Triangle, Intersection, intersect, distanceToPlane
-export Tesselation, coordinates, faces, normals
+export Tesselation, coordinates, faces, normals, mesh
 
 using StaticArrays, GeometryBasics, LinearAlgebra, Rotations, AbstractTrees
 include("BasicTypes.jl")
@@ -35,6 +35,6 @@ include("Navigators.jl")
 include("Drawing.jl")
 include("GDML.jl")
 include("Benchmark.jl")
-include("CuGeom.jl")
+#include("CuGeom.jl")  # PackageCompiler fails?
 
 end # module

src/Navigators.jl

@@ -177,12 +177,16 @@ function computeStep!(state::NavigatorState{T,NAV}, gpoint::Point3{T}, gdir::Vec
     #---If didn't hit any daughter return distance to out
     if idx == 0
         step = distanceToOut(volume.shape, lpoint, ldir)
-        if step >= -kTolerance(T)
+        if step > -kTolerance(T)
             popOut!(state)
+        elseif step < 0
+            shape = volume.shape
+            volstack = state.volstack
+            error("Negative distanceToOut. \nstep = $step \nshape = $shape \nvolstack = $volstack \npoint = $lpoint \ndir =   $ldir" )
         end
     else
     #---We hit a daughter, push it into the stack
-        step += kTolerance(T) # to ensure that do not stay in the surface of the daughter
+        step += kPushTolerance(T)  # to ensure that do not stay in the surface of the daughter
         pvol = volume.daughters[idx]
         pushIn!(state, pvol)
     end