environment tree

Documentation/Examples/geometry_example.cc

@@ -39,10 +39,10 @@ int main() {
   std::cout << "p2-p1 norm^2: " << norm << std::endl;
 
   Sphere s(p1, 10_m); // define a sphere around a point with a radius
-  std::cout << "p1 inside s:  " << s.isInside(p2) << std::endl;
+  std::cout << "p1 inside s:  " << s.Contains(p2) << std::endl;
 
   Sphere s2(p1, 3_um); // another sphere
-  std::cout << "p1 inside s2: " << s2.isInside(p2) << std::endl;
+  std::cout << "p1 inside s2: " << s2.Contains(p2) << std::endl;
 
   // let's try parallel projections:
   auto const v1 = Vector<length_d>(root, {1_m, 1_m, 0_m});

Environment/BaseEnvironment.h

-#ifndef _include_EnvironmentBase_h
-#define _include_EnvironmentBase_h
-
-#include <vector>
-#include <memory>
-#include <algorithm>
-#include <corsika/geometry/Sphere.h>
-#include <corsika/geometry/Point.h>
-
-namespace corsika::environment {
-
-class BaseEnvironment {
-    std::vector<std::unique_ptr<BaseEnvironment>> childNodes, excludedOverlapNodes;
-    corsika::geometry::Sphere shape;
-    
-public:
-    bool ShapeContains(corsika::geometry::Point const& p) const {
-        return shape.Contains(p);
-    }
-
-    BaseEnvironment* GetContainingNode(corsika::geometry::Point const& p) const {
-        if (!shape.Contains(p))
-        {
-            return nullptr;
-        }
-            
-        auto const predicate = [&] (auto const& s) {return bool(s->GetContainingNode(p));};
-        
-        auto const childFountIt = std::find_if(childNodes.cbegin(), childNodes.cend(), predicate);
-        
-        if (childFoundIt == childNodes.cend())
-        {
-            if (auto const excludedIt = std::find_if(excludedOverlapNodes.cbegin(), excludedOverlapNodes.cend(), predicate); excluded == excludedOverlapNodes.cend())
-            {
-                return *this;
-            }
-            else
-            {
-                return excludedIt;
-            }
-        }
-        
-        return (*childFoundIt)->GetContainingNode(p);
-    }
-
-    /**
-     * Get mass density at given point \a p. It is the caller's responsibility
-     * to call this method from the appropriate node in whose domain \a p is located.
-     */
-    virtual MassDensityType GetDensity(corsika::geometry::Point const& p) = 0;
-    
-    virtual std::pair<corsika::particles::Code const*, corsika::particles::Code const*> GetTargetComposition(Point const& p)
-};
-
-}
-
-#endif

Environment/VolumeTreeNode.h

+#ifndef _include_VolumeTreeNode_H
+#define _include_VolumeTreeNode_H
+
+namespace corsika::environment {
+
+class Empty {}; //<! intended for usage as default template arguments
+
+template <typename IModelProperties=Empty>
+class VolumeTreeNode {
+public:
+    using VTNUPtr = std::unique_ptr<VolumeTreeNode<IModelProperties>>;
+    using IMPSharedPtr = std::shared_ptr<IModelProperties>;
+    
+    VolumeTreeNode(VolUPtr pVolume = nullptr) : fGeoVolume(std::move(pVolume)) { }
+
+    bool Contains(Point const& p) const { return fGeoVolume->Contains(p); }
+    
+    VolumeTreeNode<IModelProperties> const* Excludes(Point const& p) const {
+        auto exclContainsIter = std::find_if(fExcludedNodes.cbegin(), fExcludedNodes.cend(),
+            [&] (auto const& s) {return bool(s->Contains(p));});
+        
+        return exclContainsIter != fExcludedNodes.cend() ? *exclContainsIter : nullptr;
+    }
+
+    /** returns a pointer to the sub-VolumeTreeNode which is "responsible" for the given
+     * \class Point \arg p, or nullptr iff \arg p is not contained in this volume.
+     */
+    VolumeTreeNode<IModelProperties> const* GetContainingNode(Point const& p) const {
+        if (!Contains(p))
+        {
+            return nullptr;
+        }
+        
+        if (auto const childContainsIter = std::find_if(fChildNodes.cbegin(), fChildNodes.cend(),
+            [&] (auto const& s) {return bool(s->Contains(p));});
+            childContainsIter == fChildNodes.cend()) // not contained in any of the children
+        {
+            if (auto const exclContainsIter = Excludes(p)) // contained in any excluded nodes
+            {
+                return exclContainsIter->GetContainingNode(p);
+            }
+            else
+            {
+                return this;
+            }
+        }
+        else 
+        {
+            return (*childContainsIter)->GetContainingNode(p);
+        }
+    }
+    
+    void addChild(VTNUPtr pChild) {
+        pChild->fParentNode = this;
+        fChildNodes.push_back(std::move(pChild));
+        // It is a bad idea to return an iterator to the inserted element
+        // because it might get invalidated when the vector needs to grow
+        // later and the caller won't notice.
+    }
+    
+    void excludeOverlapWith(VTNUPtr const& pNode) {        
+        fExcludedNodes.push_back(pNode.get());
+    }
+    
+    auto GetParent() const {
+        return fParentNode;
+    };
+    
+    auto const& GetVolume() const {
+        return *fGeoVolume;
+    }
+    
+    auto const& GetModelProperties() const {
+        return *fModelProperties;
+    }
+    
+    template <typename ModelProperties, typename... Args>
+    auto SetModelProperties(Args&&... args) {
+        static_assert(std::is_base_of_v<IModelProperties, ModelProperties>, "unusable type provided");
+        
+        fModelProperties = std::make_unique<ModelProperties>(std::forward<Args>(args)...);
+    }
+    
+    auto SetModelProperties(IMPSharedPtr ptr) {        
+        fModelProperties = IMPSharedPtr;
+    }
+    
+    template <class MediumType, typename... Args>
+    static auto CreateMedium(Args&&... args) {
+        static_assert(std::is_base_of_v<IMediumModel, MediumType>, "unusable type provided, needs to be derived from \"IMediumModel\"");
+        
+        return std::make_shared<MediumType>(std::forward<Args>(args)...);
+    }
+        
+    // factory methods for creation of nodes
+    template <class VolumeType, typename... Args>
+    static auto CreateNode(Args&&...args) {
+        static_assert(std::is_base_of_v<Volume, VolumeType>, "unusable type provided, needs to be derived from \"Volume\"");
+        
+        return std::make_unique<VolumeTreeNode<IModelProperties>>(std::make_unique<VolumeType>(std::forward<Args>(args)...));
+    }
+    
+private:
+    std::vector<VTNUPtr> fChildNodes;
+    std::vector<VolumeTreeNode<IModelProperties> const*> fExcludedNodes;
+    VolumeTreeNode<IModelProperties> const* fParentNode = nullptr;
+    VolUPtr fGeoVolume;
+    std::shared_ptr<IModelProperties> fModelProperties;    
+};
+
+}
+
+#endif

Framework/Geometry/CMakeLists.txt

@@ -9,6 +9,7 @@ set (
   Vector.h
   Point.h
   Sphere.h
+  Volume.h
   CoordinateSystem.h
   Helix.h
   BaseVector.h

Framework/Geometry/Sphere.h

 #ifndef _include_SPHERE_H_
 #define _include_SPHERE_H_
 
+#include <corsika/geometry/Volume.h>
 #include <corsika/geometry/Point.h>
 #include <corsika/units/PhysicalUnits.h>
 
 namespace corsika::geometry {
 
-  class Sphere {
-    Point center;
-    LengthType const radius;
+  class Sphere : public Volume {
+    Point const fCenter;
+    LengthType const fRadius;
 
   public:
     Sphere(Point const& pCenter, LengthType const pRadius)
-        : center(pCenter)
-        , radius(pRadius) {}
+        : fCenter(pCenter)
+        , fRadius(pRadius) {}
 
     //! returns true if the Point p is within the sphere
-    auto isInside(Point const& p) const {
-      return radius * radius > (center - p).squaredNorm();
+    bool Contains(Point const& p) const override {
+      return fRadius * fRadius > (fCenter - p).squaredNorm();
     }
   };
 

Framework/Geometry/Volume.h

+#ifndef _include_VOLUME_H_
+#define _include_VOLUME_H_
+
+#include <corsika/geometry/Point.h>
+
+namespace corsika::geometry {
+
+  class Volume {
+  
+  public:
+    //! returns true if the Point p is within the volume
+    virtual bool Contains(Point const& p) const = 0;
+    
+    virtual ~Volume() = default;
+  };
+
+} // namespace corsika::geometry
+
+#endif

Framework/Geometry/testGeometry.cc

@@ -116,8 +116,8 @@ TEST_CASE("Sphere") {
   Sphere sphere(center, 5_m);
 
   SECTION("isInside") {
-    REQUIRE_FALSE(sphere.isInside(Point(rootCS, {100_m, 0_m, 0_m})));
-    REQUIRE(sphere.isInside(Point(rootCS, {2_m, 3_m, 4_m})));
+    REQUIRE_FALSE(sphere.Contains(Point(rootCS, {100_m, 0_m, 0_m})));
+    REQUIRE(sphere.Contains(Point(rootCS, {2_m, 3_m, 4_m})));
   }
 }