Add source plane interface (#55)

* Allow usage of custom source interface

* Use initial ray weight from source

* Add resetSource and test interface

* Make mean free path particle specific

include/viennaray/rayParticle.hpp

@@ -54,6 +54,10 @@ template <typename NumericType> class rayAbstractParticle {
   /// Set the power of the cosine source distribution for this particle.
   virtual NumericType getSourceDistributionPower() const = 0;
 
+  // Set the mean free path of the particle. If the mean free path is negative,
+  // the mean free path is infinite.
+  virtual NumericType getMeanFreePath() const = 0;
+
   /// Set the number of required data vectors for this particle to
   /// collect data. If an empty vector is returned, no local data will be
   /// provided
@@ -92,6 +96,7 @@ class rayParticle : public rayAbstractParticle<NumericType> {
                    rayRNG &Rng) override { // collect data for this hit
   }
   virtual NumericType getSourceDistributionPower() const override { return 1.; }
+  virtual NumericType getMeanFreePath() const override { return -1.; }
   virtual std::vector<std::string> getLocalDataLabels() const override {
     return {};
   }

include/viennaray/raySource.hpp

@@ -3,17 +3,12 @@
 #include <rayRNG.hpp>
 #include <rayUtil.hpp>
 
-template <typename Derived> class raySource {
-protected:
-  raySource() = default;
-  ~raySource() = default;
-
+template <typename NumericType, int D> class raySource {
 public:
-  Derived &derived() { return static_cast<Derived &>(*this); }
-  const Derived &derived() const { return static_cast<const Derived &>(*this); }
+  virtual ~raySource() = default;
 
-  void fillRay(RTCRay &ray, const size_t idx, rayRNG &RngState) const {
-    derived().fillRay(ray, idx, RngState);
-  }
-  size_t getNumPoints() const { return derived().getNumPoints(); }
+  virtual rayPair<rayTriple<NumericType>>
+  getOriginAndDirection(const size_t idx, rayRNG &RngState) const = 0;
+  virtual size_t getNumPoints() const = 0;
+  virtual NumericType getInitialRayWeight(const size_t idx) const { return 1.; }
 };

include/viennaray/raySourceGrid.hpp

@@ -3,7 +3,7 @@
 #include <raySource.hpp>
 
 template <typename NumericType, int D>
-class raySourceGrid : public raySource<raySourceGrid<NumericType, D>> {
+class raySourceGrid : public raySource<NumericType, D> {
 public:
   raySourceGrid(std::vector<rayTriple<NumericType>> &sourceGrid,
                 NumericType cosinePower,
@@ -13,30 +13,15 @@ class raySourceGrid : public raySource<raySourceGrid<NumericType, D>> {
         secondDir_(traceSettings[2]), minMax_(traceSettings[3]),
         posNeg_(traceSettings[4]), ee_(((NumericType)2) / (cosinePower + 1)) {}
 
-  void fillRay(RTCRay &ray, const size_t idx, rayRNG &RngState) const {
+  rayPair<rayTriple<NumericType>>
+  getOriginAndDirection(const size_t idx, rayRNG &RngState) const override {
     auto origin = sourceGrid_[idx % numPoints_];
     auto direction = getDirection(RngState);
 
-#ifdef ARCH_X86
-    reinterpret_cast<__m128 &>(ray) =
-        _mm_set_ps(1e-4f, (float)origin[2], (float)origin[1], (float)origin[0]);
-
-    reinterpret_cast<__m128 &>(ray.dir_x) = _mm_set_ps(
-        0.0f, (float)direction[2], (float)direction[1], (float)direction[0]);
-#else
-    ray.org_x = (float)origin[0];
-    ray.org_y = (float)origin[1];
-    ray.org_z = (float)origin[2];
-    ray.tnear = 1e-4f;
-
-    ray.dir_x = (float)direction[0];
-    ray.dir_y = (float)direction[1];
-    ray.dir_z = (float)direction[2];
-    ray.tnear = 0.0f;
-#endif
+    return {origin, direction};
   }
 
-  size_t getNumPoints() const { return numPoints_; }
+  size_t getNumPoints() const override { return numPoints_; }
 
 private:
   rayTriple<NumericType> getDirection(rayRNG &RngState) const {

include/viennaray/raySourceRandom.hpp

@@ -3,7 +3,7 @@
 #include <raySource.hpp>
 
 template <typename NumericType, int D>
-class raySourceRandom : public raySource<raySourceRandom<NumericType, D>> {
+class raySourceRandom : public raySource<NumericType, D> {
   using boundingBoxType = rayPair<rayTriple<NumericType>>;
 
 public:
@@ -19,7 +19,8 @@ class raySourceRandom : public raySource<raySourceRandom<NumericType, D>> {
         customDirection_(customDirection), orthonormalBasis_(orthonormalBasis) {
   }
 
-  void fillRay(RTCRay &ray, const size_t idx, rayRNG &RngState) const {
+  rayPair<rayTriple<NumericType>>
+  getOriginAndDirection(const size_t idx, rayRNG &RngState) const override {
     auto origin = getOrigin(RngState);
     rayTriple<NumericType> direction;
     if (customDirection_) {
@@ -28,26 +29,10 @@ class raySourceRandom : public raySource<raySourceRandom<NumericType, D>> {
       direction = getDirection(RngState);
     }
 
-#ifdef ARCH_X86
-    reinterpret_cast<__m128 &>(ray) =
-        _mm_set_ps(1e-4f, (float)origin[2], (float)origin[1], (float)origin[0]);
-
-    reinterpret_cast<__m128 &>(ray.dir_x) = _mm_set_ps(
-        0.0f, (float)direction[2], (float)direction[1], (float)direction[0]);
-#else
-    ray.org_x = (float)origin[0];
-    ray.org_y = (float)origin[1];
-    ray.org_z = (float)origin[2];
-    ray.tnear = 1e-4f;
-
-    ray.dir_x = (float)direction[0];
-    ray.dir_y = (float)direction[1];
-    ray.dir_z = (float)direction[2];
-    ray.time = 0.0f;
-#endif
+    return {origin, direction};
   }
 
-  size_t getNumPoints() const { return numPoints_; }
+  size_t getNumPoints() const override { return numPoints_; }
 
 private:
   rayTriple<NumericType> getOrigin(rayRNG &RngState) const {

include/viennaray/rayTrace.hpp

@@ -38,9 +38,10 @@ template <class NumericType, int D> class rayTrace {
     std::array<rayTriple<NumericType>, 3> orthonormalBasis;
     if (usePrimaryDirection_)
       orthonormalBasis = rayInternal::getOrthonormalBasis(primaryDirection_);
-    auto raySource = raySourceRandom<NumericType, D>(
-        boundingBox, pParticle_->getSourceDistributionPower(), traceSettings,
-        geometry_.getNumPoints(), usePrimaryDirection_, orthonormalBasis);
+    if (!pSource_)
+      pSource_ = std::make_unique<raySourceRandom<NumericType, D>>(
+          boundingBox, pParticle_->getSourceDistributionPower(), traceSettings,
+          geometry_.getNumPoints(), usePrimaryDirection_, orthonormalBasis);
 
     auto localDataLabels = pParticle_->getLocalDataLabels();
     if (!localDataLabels.empty()) {
@@ -52,10 +53,10 @@ template <class NumericType, int D> class rayTrace {
       }
     }
 
-    rayTraceKernel tracer(device_, geometry_, boundary, raySource, pParticle_,
-                          dataLog_, numberOfRaysPerPoint_, numberOfRaysFixed_,
-                          useRandomSeeds_, calcFlux_, lambda_, runNumber_++,
-                          hitCounter_, RTInfo_);
+    rayTraceKernel tracer(device_, geometry_, boundary, std::move(pSource_),
+                          pParticle_, dataLog_, numberOfRaysPerPoint_,
+                          numberOfRaysFixed_, useRandomSeeds_, calcFlux_,
+                          runNumber_++, hitCounter_, RTInfo_);
     tracer.setTracingData(&localData_, pGlobalData_);
     tracer.apply();
 
@@ -121,6 +122,16 @@ template <class NumericType, int D> class rayTrace {
     }
   }
 
+  /// Set a custom source for the ray tracing. Per default a random source is
+  /// set up. The source has to be a user defined object that has to interface
+  /// the raySource class.
+  void setSource(std::unique_ptr<raySource<NumericType, D>> source) {
+    pSource_ = std::move(source);
+  }
+
+  /// Reset the source to the default random source.
+  void resetSource() { pSource_.reset(); }
+
   /// Set the number of rays per geometry point.
   /// The total number of rays, that are traced, is the set number set here
   /// times the number of points in the geometry.
@@ -151,8 +162,6 @@ template <class NumericType, int D> class rayTrace {
     usePrimaryDirection_ = true;
   }
 
-  void setMeanFreePath(const NumericType lambda) { lambda_ = lambda; }
-
   /// Set whether random seeds for the internal random number generators
   /// should be used.
   void setUseRandomSeeds(const bool useRand) { useRandomSeeds_ = useRand; }
@@ -371,21 +380,26 @@ template <class NumericType, int D> class rayTrace {
 
 private:
   RTCDevice device_;
+
   rayGeometry<NumericType, D> geometry_;
   std::unique_ptr<rayAbstractParticle<NumericType>> pParticle_ = nullptr;
+  std::unique_ptr<raySource<NumericType, D>> pSource_ = nullptr;
+
   size_t numberOfRaysPerPoint_ = 1000;
   size_t numberOfRaysFixed_ = 0;
   NumericType diskRadius_ = 0;
   NumericType gridDelta_ = 0;
+
   rayBoundaryCondition boundaryConditions_[D] = {};
   rayTraceDirection sourceDirection_ = rayTraceDirection::POS_Z;
   rayTriple<NumericType> primaryDirection_ = {0.};
+
   bool usePrimaryDirection_ = false;
   bool useRandomSeeds_ = false;
   size_t runNumber_ = 0;
   bool calcFlux_ = true;
   bool checkError_ = true;
-  NumericType lambda_ = -1.;
+
   rayHitCounter<NumericType> hitCounter_;
   rayTracingData<NumericType> localData_;
   rayTracingData<NumericType> *pGlobalData_ = nullptr;