Add GPU examples

CMakeLists.txt

@@ -93,6 +93,7 @@ include("cmake/cpm.cmake")
 CPMAddPackage(
   NAME ViennaCore
   VERSION 1.6.1
+  GIT_TAG main
   GIT_REPOSITORY "https://github.com/ViennaTools/ViennaCore"
   OPTIONS "VIENNACORE_USE_GPU ${VIENNARAY_USE_GPU}")
 

gpu/examples/CMakeLists.txt

@@ -1,4 +1,14 @@
 project(ViennaRay-GPU_Examples)
 
-add_gpu_executable(trenchGPU target_name trench.cpp)
+add_custom_target(ViennaRay-GPU_Examples ALL)
+
+add_gpu_executable(GPU_trenchTriangles target_name trenchTriangles.cpp)
 configure_file(Resources/trenchMesh.dat ${CMAKE_CURRENT_BINARY_DIR}/trenchMesh.dat COPYONLY)
+add_dependencies(ViennaRay-GPU_Examples ${target_name})
+
+add_gpu_executable(GPU_trenchDisks target_name trenchDisks.cpp)
+configure_file(${CMAKE_SOURCE_DIR}/examples/trench/trenchGrid3D.dat ${CMAKE_CURRENT_BINARY_DIR}/trenchGrid3D.dat COPYONLY)
+add_dependencies(ViennaRay-GPU_Examples ${target_name})
+
+add_gpu_executable(GPU_trenchLines target_name trenchLines.cpp)
+add_dependencies(ViennaRay-GPU_Examples ${target_name})

gpu/examples/trenchDisks.cpp

@@ -0,0 +1,65 @@
+#include <raygTraceDisk.hpp>
+
+#include <omp.h>
+
+using namespace viennaray;
+
+int main() {
+
+  omp_set_num_threads(16);
+  constexpr int D = 3;
+  using NumericType = float;
+  Logger::setLogLevel(LogLevel::DEBUG);
+
+  auto context = DeviceContext::createContext("../../lib/ptx", 0);
+  // relative to build directory
+
+  // Read stored geometry grid
+  NumericType gridDelta;
+  std::vector<VectorType<NumericType, D>> points;
+  std::vector<VectorType<NumericType, D>> normals;
+  rayInternal::readGridFromFile("trenchGrid3D.dat", gridDelta, points, normals);
+
+  gpu::DiskMesh mesh;
+  mesh.nodes = points;
+  mesh.normals = normals;
+  mesh.gridDelta = static_cast<float>(gridDelta);
+  computeBoundingBox(mesh);
+
+  std::vector<int> materialIds(mesh.nodes.size(), 7);
+  for (int i = mesh.nodes.size() / 2; i < mesh.nodes.size(); ++i) {
+    materialIds[i] = 1;
+  }
+
+  gpu::Particle<NumericType> particle;
+  particle.direction = {0.0f, 0.0f, -1.0f};
+  particle.name = "Particle";
+  particle.sticking = 1.f;
+  particle.dataLabels = {"particleFlux"};
+  particle.materialSticking[7] = 1.f;
+  particle.materialSticking[1] = .1f;
+
+  std::unordered_map<std::string, unsigned int> pMap = {{"Particle", 0}};
+  std::vector<gpu::CallableConfig> cMap = {
+      {0, gpu::CallableSlot::COLLISION, "__direct_callable__particleCollision"},
+      {0, gpu::CallableSlot::REFLECTION,
+       "__direct_callable__particleReflection"}};
+
+  gpu::TraceDisk<NumericType, D> tracer(context);
+  tracer.setGeometry(mesh);
+  tracer.setMaterialIds(materialIds);
+  tracer.setCallables("CallableWrapper", context->modulePath);
+  tracer.setParticleCallableMap({pMap, cMap});
+  tracer.setNumberOfRaysPerPoint(1000);
+  tracer.insertNextParticle(particle);
+  tracer.prepareParticlePrograms();
+
+  tracer.apply();
+
+  std::vector<float> flux(mesh.nodes.size());
+  tracer.getFlux(flux.data(), 0, 0, 1);
+
+  rayInternal::writeVTK<float, D>("trenchResult.vtk", points, flux);
+
+  return 0;
+}

gpu/examples/trenchLines.cpp

@@ -0,0 +1,60 @@
+#include <raygTraceLine.hpp>
+
+#include <omp.h>
+
+using namespace viennaray;
+
+int main() {
+
+  omp_set_num_threads(16);
+  constexpr int D = 3;
+  using NumericType = float;
+  Logger::setLogLevel(LogLevel::DEBUG);
+
+  auto context = DeviceContext::createContext("../../lib/ptx", 0);
+  // relative to build directory
+
+  // Read stored geometry grid
+  NumericType gridDelta;
+  std::vector<VectorType<NumericType, D>> points;
+
+  gpu::LineMesh mesh;
+  mesh.nodes = points;
+  mesh.gridDelta = static_cast<float>(gridDelta);
+  computeBoundingBox(mesh);
+
+  std::vector<int> materialIds(mesh.lines.size(), 7);
+  for (int i = mesh.lines.size() / 2; i < mesh.lines.size(); ++i) {
+    materialIds[i] = 1;
+  }
+
+  gpu::Particle<NumericType> particle;
+  particle.direction = {0.0f, 0.0f, -1.0f};
+  particle.name = "Particle";
+  particle.sticking = 1.f;
+  particle.dataLabels = {"particleFlux"};
+  particle.materialSticking[7] = 1.f;
+  particle.materialSticking[1] = .1f;
+
+  std::unordered_map<std::string, unsigned int> pMap = {{"Particle", 0}};
+  std::vector<gpu::CallableConfig> cMap = {
+      {0, gpu::CallableSlot::COLLISION, "__direct_callable__particleCollision"},
+      {0, gpu::CallableSlot::REFLECTION,
+       "__direct_callable__particleReflection"}};
+
+  gpu::TraceLine<NumericType, D> tracer(context);
+  tracer.setGeometry(mesh);
+  tracer.setMaterialIds(materialIds);
+  tracer.setCallables("CallableWrapper", context->modulePath);
+  tracer.setParticleCallableMap({pMap, cMap});
+  tracer.setNumberOfRaysPerPoint(1000);
+  tracer.insertNextParticle(particle);
+  tracer.prepareParticlePrograms();
+
+  tracer.apply();
+
+  std::vector<float> flux(mesh.lines.size());
+  tracer.getFlux(flux.data(), 0, 0, 1);
+
+  return 0;
+}

gpu/examples/trench.cpp → gpu/examples/trenchTriangles.cpp

@@ -1,12 +1,8 @@
-#include <raygMesh.hpp>
-#include <raygTraceDisk.hpp>
-#include <raygTraceLine.hpp>
 #include <raygTraceTriangle.hpp>
 
-#include <fstream>
 #include <omp.h>
 
-#define COUNT_RAYS
+// #define COUNT_RAYS
 
 using namespace viennaray;
 
@@ -43,7 +39,6 @@ int main(int argc, char **argv) {
   gpu::TraceTriangle<NumericType, D> tracer(context);
   tracer.setGeometry(mesh);
   tracer.setMaterialIds(materialIds);
-  tracer.setPipeline("GeneralPipeline", context->modulePath);
   tracer.setCallables("CallableWrapper", context->modulePath);
   tracer.setParticleCallableMap({pMap, cMap});
   tracer.setNumberOfRaysPerPoint(100);
@@ -63,8 +58,10 @@ int main(int argc, char **argv) {
   std::vector<float> flux(mesh.triangles.size());
   tracer.getFlux(flux.data(), 0, 0);
 
+#ifdef COUNT_RAYS
   rayCountBuffer.download(&rayCount, 1);
   std::cout << "Trace count: " << rayCount << std::endl;
+#endif
 
   tracer.freeBuffers();
 

gpu/include/raygDiskGeometry.hpp

@@ -42,15 +42,15 @@ template <int D> struct DiskGeometry {
       launchParams.source.maxPoint[1] = mesh.maximumExtent[1];
       launchParams.source.planeHeight = mesh.maximumExtent[2] + 2 * mesh.radius;
     }
-    launchParams.numElements = mesh.points.size();
+    launchParams.numElements = mesh.nodes.size();
 
     // 2 inputs: one for the geometry, one for the boundary
     std::array<OptixBuildInput, 2> diskInput{};
     std::array<uint32_t, 2> diskInputFlags{};
 
     // ------------------- geometry input -------------------
     // upload the model to the device: the builder
-    geometryPointBuffer.allocUpload(mesh.points);
+    geometryPointBuffer.allocUpload(mesh.nodes);
     geometryNormalBuffer.allocUpload(mesh.normals);
 
     // create local variables, because we need a *pointer* to the
@@ -59,10 +59,10 @@ template <int D> struct DiskGeometry {
     CUdeviceptr d_geoNormals = geometryNormalBuffer.dPointer();
 
     // AABB build input
-    std::vector<OptixAabb> aabb(mesh.points.size());
+    std::vector<OptixAabb> aabb(mesh.nodes.size());
 
-    for (size_t i = 0; i < mesh.points.size(); ++i) {
-      Vec3Df C = mesh.points[i];
+    for (size_t i = 0; i < mesh.nodes.size(); ++i) {
+      Vec3Df C = mesh.nodes[i];
       Vec3Df N = mesh.normals[i];
       Normalize(N);
 
@@ -89,7 +89,7 @@ template <int D> struct DiskGeometry {
     diskInput[0] = {};
     diskInput[0].type = OPTIX_BUILD_INPUT_TYPE_CUSTOM_PRIMITIVES;
     diskInput[0].customPrimitiveArray.aabbBuffers = &d_aabb;
-    diskInput[0].customPrimitiveArray.numPrimitives = mesh.points.size();
+    diskInput[0].customPrimitiveArray.numPrimitives = mesh.nodes.size();
 
     uint32_t diskInput_flags[1] = {OPTIX_GEOMETRY_FLAG_NONE};
     diskInput[0].customPrimitiveArray.flags = diskInput_flags;
@@ -104,7 +104,7 @@ template <int D> struct DiskGeometry {
     // ------------------------- boundary input -------------------------
     auto boundaryMesh = makeBoundary(mesh);
     // upload the model to the device: the builder
-    boundaryPointBuffer.allocUpload(boundaryMesh.points);
+    boundaryPointBuffer.allocUpload(boundaryMesh.nodes);
     boundaryNormalBuffer.allocUpload(boundaryMesh.normals);
 
     // create local variables, because we need a *pointer* to the
@@ -113,9 +113,9 @@ template <int D> struct DiskGeometry {
     CUdeviceptr d_boundNormals = boundaryNormalBuffer.dPointer();
 
     // AABB build input for boundary disks
-    std::vector<OptixAabb> aabbBoundary(boundaryMesh.points.size());
-    for (size_t i = 0; i < boundaryMesh.points.size(); ++i) {
-      Vec3Df C = boundaryMesh.points[i];
+    std::vector<OptixAabb> aabbBoundary(boundaryMesh.nodes.size());
+    for (size_t i = 0; i < boundaryMesh.nodes.size(); ++i) {
+      Vec3Df C = boundaryMesh.nodes[i];
       Vec3Df N = boundaryMesh.normals[i];
       Normalize(N);
 
@@ -136,8 +136,7 @@ template <int D> struct DiskGeometry {
     diskInput[1] = {};
     diskInput[1].type = OPTIX_BUILD_INPUT_TYPE_CUSTOM_PRIMITIVES;
     diskInput[1].customPrimitiveArray.aabbBuffers = &d_aabbBoundary;
-    diskInput[1].customPrimitiveArray.numPrimitives =
-        boundaryMesh.points.size();
+    diskInput[1].customPrimitiveArray.numPrimitives = boundaryMesh.nodes.size();
 
     diskInput[1].customPrimitiveArray.flags = diskInput_flags;
     diskInput[1].customPrimitiveArray.numSbtRecords = 1;
@@ -235,37 +234,37 @@ template <int D> struct DiskGeometry {
     // xmin - back 0
     Vec3Df xminPoint = {bbMin[0], (bbMin[1] + bbMax[1]) / 2,
                         (bbMin[2] + bbMax[2]) / 2};
-    boundaryMesh.points.push_back(xminPoint);
+    boundaryMesh.nodes.push_back(xminPoint);
     boundaryMesh.normals.push_back({1, 0, 0});
 
     // xmax - front 1
     Vec3Df xmaxPoint = {bbMax[0], (bbMin[1] + bbMax[1]) / 2,
                         (bbMin[2] + bbMax[2]) / 2};
-    boundaryMesh.points.push_back(xmaxPoint);
+    boundaryMesh.nodes.push_back(xmaxPoint);
     boundaryMesh.normals.push_back({-1, 0, 0});
 
     // ymin - left 2
     Vec3Df yminPoint = {(bbMin[0] + bbMax[0]) / 2, bbMin[1],
                         (bbMin[2] + bbMax[2]) / 2};
-    boundaryMesh.points.push_back(yminPoint);
+    boundaryMesh.nodes.push_back(yminPoint);
     boundaryMesh.normals.push_back({0, 1, 0});
 
     // ymax - right 3
     Vec3Df ymaxPoint = {(bbMin[0] + bbMax[0]) / 2, bbMax[1],
                         (bbMin[2] + bbMax[2]) / 2};
-    boundaryMesh.points.push_back(ymaxPoint);
+    boundaryMesh.nodes.push_back(ymaxPoint);
     boundaryMesh.normals.push_back({0, -1, 0});
 
     // zmin - bottom 4
     Vec3Df zminPoint = {(bbMin[0] + bbMax[0]) / 2, (bbMin[1] + bbMax[1]) / 2,
                         bbMin[2]};
-    boundaryMesh.points.push_back(zminPoint);
+    boundaryMesh.nodes.push_back(zminPoint);
     boundaryMesh.normals.push_back({0, 0, 1});
 
     // zmax - top 5
     Vec3Df zmaxPoint = {(bbMin[0] + bbMax[0]) / 2, (bbMin[1] + bbMax[1]) / 2,
                         bbMax[2]};
-    boundaryMesh.points.push_back(zmaxPoint);
+    boundaryMesh.nodes.push_back(zmaxPoint);
     boundaryMesh.normals.push_back({0, 0, -1});
 
     return boundaryMesh;

gpu/include/raygMesh.hpp

@@ -16,45 +16,42 @@ struct LineMesh {
 
   Vec3Df minimumExtent;
   Vec3Df maximumExtent;
-  float gridDelta;
+  float gridDelta = 0.f;
 };
 
 struct TriangleMesh {
-  std::vector<Vec3Df> vertices;
+  std::vector<Vec3Df> nodes;
   std::vector<Vec3D<unsigned>> triangles;
 
   Vec3Df minimumExtent;
   Vec3Df maximumExtent;
-  float gridDelta;
+  float gridDelta = 0.f;
 };
 
 struct DiskMesh {
-  std::vector<Vec3Df> points;
+  std::vector<Vec3Df> nodes;
   std::vector<Vec3Df> normals;
 
   Vec3Df minimumExtent;
   Vec3Df maximumExtent;
-  float radius;
-  float gridDelta;
-};
-
-struct SphereMesh {
-  std::vector<Vec3Df> vertices;
-  std::vector<float> radii;
-
-  Vec3Df minimumExtent;
-  Vec3Df maximumExtent;
-  float gridDelta;
+  float radius = 0.f;
+  float gridDelta = 0.f;
 };
 
-struct OrientedPointCloud {
-  std::vector<Vec3Df> vertices;
-  std::vector<Vec3Df> normals;
-
-  Vec3Df minimumExtent;
-  Vec3Df maximumExtent;
-  float gridDelta;
-};
+template <class MeshType> void computeBoundingBox(MeshType &mesh) {
+  if (mesh.nodes.empty())
+    return;
+  mesh.minimumExtent = mesh.nodes[0];
+  mesh.maximumExtent = mesh.nodes[0];
+  for (const auto &p : mesh.nodes) {
+    for (int d = 0; d < 3; ++d) {
+      if (p[d] < mesh.minimumExtent[d])
+        mesh.minimumExtent[d] = p[d];
+      if (p[d] > mesh.maximumExtent[d])
+        mesh.maximumExtent[d] = p[d];
+    }
+  }
+}
 
 inline TriangleMesh readMeshFromFile(const std::string &fileName) {
   TriangleMesh mesh;
@@ -90,13 +87,12 @@ inline TriangleMesh readMeshFromFile(const std::string &fileName) {
   size_t numTriangles;
   dataFile >> numTriangles;
 
-  mesh.vertices.resize(numPoints);
+  mesh.nodes.resize(numPoints);
   mesh.triangles.resize(numTriangles);
   for (size_t i = 0; i < numPoints; ++i) {
     dataFile >> id;
     assert(id == "n");
-    dataFile >> mesh.vertices[i][0] >> mesh.vertices[i][1] >>
-        mesh.vertices[i][2];
+    dataFile >> mesh.nodes[i][0] >> mesh.nodes[i][1] >> mesh.nodes[i][2];
   }
   for (size_t i = 0; i < numTriangles; ++i) {
     dataFile >> id;