Added fin tapering option (#63)

* Fin tapter angle in 2D

* Fix surface mesh printing in hole etching example

* Enable primary direction in cell set tracing

* Read average filling fraction

* csTracign point source

* Tapered fin in 3D

* Added fin unit test

examples/holeEtching/holeEtching.py

@@ -53,10 +53,10 @@
 process.setProcessDuration(params["processTime"])  # seconds
 
 # print initial surface
-geometry.saveSurface(filename="initial.vtp", addMaterialIds=True)
+geometry.saveSurfaceMesh(filename="initial.vtp", addMaterialIds=True)
 
 # run the process
 process.apply()
 
 # print final surface
-geometry.saveSurface(filename="final.vtp", addMaterialIds=True)
+geometry.saveSurfaceMesh(filename="final.vtp", addMaterialIds=True)

examples/volumeModel/volumeModel.cpp

@@ -34,12 +34,8 @@ int main(int argc, char *argv[]) {
       params.meanFreePath /* damage ion mean free path */,
       psMaterial::None /*mask material (no mask)*/);
 
-  psProcess<NumericType, D> process;
-  process.setDomain(geometry);
-  process.setProcessModel(model);
-  process.setProcessDuration(0.); // apply only damage model
-
-  process.apply();
+  psProcess<NumericType, D>(geometry, model, 0.)
+      .apply(); // apply only damage model
 
   geometry->getCellSet()->writeVTU("DamageModel.vtu");
 }

include/cellSet/csDenseCellSet.hpp

@@ -195,6 +195,23 @@ template <class T, int D> class csDenseCellSet {
     return getFillingFractions()->at(idx);
   }
 
+  T getAverageFillingFraction(const std::array<T, 3> &point,
+                              const T radius) const {
+    T sum = 0.;
+    int count = 0;
+    for (int i = 0; i < numberOfCells; i++) {
+      auto &cell = cellGrid->template getElements<(1 << D)>()[i];
+      auto node = cellGrid->getNodes()[cell[0]];
+      for (int j = 0; j < D; j++)
+        node[j] += gridDelta / 2.;
+      if (csUtil::distance(node, point) < radius) {
+        sum += fillingFractions->at(i);
+        count++;
+      }
+    }
+    return sum / count;
+  }
+
   int getIndex(const std::array<T, 3> &point) { return findIndex(point); }
 
   std::vector<T> *getScalarData(std::string name) {

include/cellSet/csPointSource.hpp

@@ -0,0 +1,40 @@
+#pragma once
+
+#include <csUtil.hpp>
+#include <raySource.hpp>
+
+template <typename NumericType, int D>
+class csPointSource : public raySource<NumericType, D> {
+  const unsigned mNumPoints;
+  const csTriple<NumericType> origin;
+  const csTriple<NumericType> direction;
+
+public:
+  csPointSource(csTriple<NumericType> passedOrigin,
+                csTriple<NumericType> passedDirection,
+                std::array<int, 5> &pTraceSettings, const size_t pNumPoints)
+      : origin(passedOrigin), direction(passedDirection),
+        mNumPoints(pNumPoints) {}
+
+  void fillRay(RTCRay &ray, const size_t idx, rayRNG &RngState) override final {
+#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
+  }
+
+  size_t getNumPoints() const override final { return mNumPoints; }
+};

include/cellSet/csTracing.hpp

@@ -3,6 +3,7 @@
 #include <embree3/rtcore.h>
 
 #include <csDenseCellSet.hpp>
+#include <csPointSource.hpp>
 #include <csTracingKernel.hpp>
 #include <csTracingParticle.hpp>
 
@@ -25,10 +26,16 @@ template <class T, int D> class csTracing {
   size_t mNumberOfRaysFixed = 1000;
   T mGridDelta = 0;
   rayBoundaryCondition mBoundaryConditions[D] = {};
-  rayTraceDirection mSourceDirection = rayTraceDirection::POS_Z;
+  const rayTraceDirection mSourceDirection =
+      D == 2 ? rayTraceDirection::POS_Y : rayTraceDirection::POS_Z;
   bool mUseRandomSeeds = true;
+  bool usePrimaryDirection = false;
+  rayTriple<T> primaryDirection = {0.};
   size_t mRunNumber = 0;
   int excludeMaterialId = -1;
+  bool usePointSource = false;
+  csTriple<T> pointSourceOrigin = {0.};
+  csTriple<T> pointSourceDirection = {0.};
 
 public:
   csTracing() : mDevice(rtcNewDevice("hugepages=1")) {
@@ -55,15 +62,37 @@ template <class T, int D> class csTracing {
                                       traceSettings);
 
     std::array<rayTriple<T>, 3> orthoBasis;
-    auto raySource = raySourceRandom<T, D>(
-        boundingBox, mParticle->getSourceDistributionPower(), traceSettings,
-        mGeometry.getNumPoints(), false, orthoBasis);
+    if (usePrimaryDirection) {
+      orthoBasis = rayInternal::getOrthonormalBasis(primaryDirection);
+      psLogger::getInstance()
+          .addInfo("Using primary direction: " +
+                   std::to_string(primaryDirection[0]) + " " +
+                   std::to_string(primaryDirection[1]) + " " +
+                   std::to_string(primaryDirection[2]))
+          .print();
+    }
 
-    auto tracer = csTracingKernel<T, D>(
-        mDevice, mGeometry, boundary, raySource, mParticle,
-        mNumberOfRaysPerPoint, mNumberOfRaysFixed, mUseRandomSeeds,
-        mRunNumber++, cellSet, excludeMaterialId - 1);
-    tracer.apply();
+    if (usePointSource) {
+      auto raySource =
+          csPointSource<T, D>(pointSourceOrigin, pointSourceDirection,
+                              traceSettings, mGeometry.getNumPoints());
+
+      csTracingKernel<T, D>(mDevice, mGeometry, boundary, raySource, mParticle,
+                            mNumberOfRaysPerPoint, mNumberOfRaysFixed,
+                            mUseRandomSeeds, mRunNumber++, cellSet,
+                            excludeMaterialId - 1)
+          .apply();
+    } else {
+      auto raySource = raySourceRandom<T, D>(
+          boundingBox, mParticle->getSourceDistributionPower(), traceSettings,
+          mGeometry.getNumPoints(), usePrimaryDirection, orthoBasis);
+
+      csTracingKernel<T, D>(mDevice, mGeometry, boundary, raySource, mParticle,
+                            mNumberOfRaysPerPoint, mNumberOfRaysFixed,
+                            mUseRandomSeeds, mRunNumber++, cellSet,
+                            excludeMaterialId - 1)
+          .apply();
+    }
 
     averageNeighborhood();
     boundary.releaseGeometry();
@@ -73,6 +102,13 @@ template <class T, int D> class csTracing {
     cellSet = passedCellSet;
   }
 
+  void setPointSource(const csTriple<T> &passedOrigin,
+                      const csTriple<T> &passedDirection) {
+    usePointSource = true;
+    pointSourceOrigin = passedOrigin;
+    pointSourceDirection = passedDirection;
+  }
+
   template <typename ParticleType>
   void setParticle(std::unique_ptr<ParticleType> &p) {
     static_assert(std::is_base_of<csAbstractParticle<T>, ParticleType>::value &&
@@ -90,6 +126,11 @@ template <class T, int D> class csTracing {
     mNumberOfRaysFixed = 0;
   }
 
+  void setPrimaryDirection(const rayTriple<T> pPrimaryDirection) {
+    primaryDirection = pPrimaryDirection;
+    usePrimaryDirection = true;
+  }
+
   void setExcludeMaterialId(int passedId) { excludeMaterialId = passedId; }
 
   lsSmartPointer<csDenseCellSet<T, D>> getCellSet() const { return cellSet; }

include/cellSet/csTracingKernel.hpp

@@ -216,18 +216,12 @@ template <typename T, int D> class csTracingKernel {
           rayHit.ray.time = 0.0f;
 #endif
         } while (reflect);
-
-        if (psLogger::getLogLevel() >= 3)
-          psUtils::printProgress(idx, mNumRays);
       } // end ray tracing for loop
 
 #pragma omp critical
       myCellSet->mergePath(path, mNumRays);
     } // end parallel section
 
-    if (psLogger::getLogLevel() >= 3)
-      std::cout << std::endl;
-
     rtcReleaseGeometry(rtcGeometry);
     rtcReleaseGeometry(rtcBoundary);
   }