Merge pull request #16 from pshriwise/particle-tool

Particle Tool

CMakeLists.txt

@@ -5,6 +5,7 @@ option(XDG_ENABLE_MOAB    "Enable support for the MOAB mesh library"    ON)
 option(XDG_ENABLE_MFEM    "Enable support for the MFEM mesh library"    OFF)
 option(XDG_ENABLE_LIBMESH "Enable support for the libMesh mesh library" OFF)
 option(XDG_BUILD_TESTS    "Enable C++ unit testing"                     ON)
+option(XDG_BUILD_TOOLS    "Enable tools and miniapps"                   ON)
 
 # Set version numbers
 set(XDG_VERSION_MAJOR 0)
@@ -25,7 +26,7 @@ find_package(MOAB REQUIRED HINTS ${MOAB_DIR})
 endif()
 
 # use Embree for CPU ray tracing
-find_package(embree 3.6.1 REQUIRED)
+find_package(embree 3.0.0...4.0.0 REQUIRED)
 if (NOT ${EMBREE_VERSION} VERSION_GREATER 3.6.0)
   message(FATAL_ERROR "XDG requires Embree v3.6.1 or higher.")
 endif()
@@ -76,10 +77,18 @@ else()
   target_compile_definitions(xdg PUBLIC XDG_EMBREE3)
 endif()
 
+if (${CMAKE_BUILD_TYPE} MATCHES "Debug")
+  target_compile_definitions(xdg PUBLIC XDG_DEBUG)
+endif()
+
 if (XDG_BUILD_TESTS)
   add_subdirectory("tests")
 endif()
 
+if (XDG_BUILD_TOOLS)
+  add_subdirectory("tools")
+endif()
+
 target_include_directories(xdg
   PUBLIC
   $<INSTALL_INTERFACE:include>

include/xdg/constants.h

@@ -9,6 +9,16 @@ namespace xdg {
 
 constexpr double INFTY {std::numeric_limits<double>::max()};
 
+#ifdef XDG_DEBUG
+  // from Embree, if the floating point value of Tfar is larger than this value,
+  // it is considered overflow by the internal hit verification function. This
+  // is only enabled when Embree is compiled in debug mode. I think a
+  // corresponding behavior makes sense to keep here for now
+  constexpr double INFTYF {1.844E18f};
+#else
+  constexpr double INFTYF {std::numeric_limits<float>::max()};
+#endif
+
 // Whether information pertains to a surface or volume
 enum class GeometryType {
  SURFACE = 2,

include/xdg/mesh_manager_interface.h

@@ -62,7 +62,7 @@ class MeshManager {
 
   virtual void add_surface_to_volume(MeshID volume, MeshID surface, Sense sense, bool overwrite=false) = 0;
 
-  MeshID next_volume(MeshID surface, MeshID current_volume) const;
+  MeshID next_volume(MeshID current_volume, MeshID surface) const;
 
   // Methods
   MeshID next_volume_id() const;
@@ -74,8 +74,11 @@ class MeshManager {
   // Metadata methods
   virtual void parse_metadata() = 0;
 
-  virtual Property get_volume_property(MeshID volume, PropertyType type) const = 0;
-  virtual Property get_surface_property(MeshID surface, PropertyType type) const = 0;
+  bool volume_has_property(MeshID volume, PropertyType type) const;
+  bool surface_has_property(MeshID surface, PropertyType type) const;
+
+  Property get_volume_property(MeshID volume, PropertyType type) const;
+  Property get_surface_property(MeshID surface, PropertyType type) const;
 
   // Accessors
   const std::vector<MeshID>& volumes() const { return volumes_; }

include/xdg/moab/mesh_manager.h

@@ -81,11 +81,6 @@ class MOABMeshManager : public MeshManager {
   // Metadata
   void parse_metadata() override;
 
-  // TODO: move to mesh_manger_interface???
-  Property get_volume_property(MeshID volume, PropertyType type) const override;
-
-  Property get_surface_property(MeshID surface, PropertyType type) const override;
-
   // Other
   MeshLibrary mesh_library() const override {return MeshLibrary::MOAB; }
 

include/xdg/ray.h

@@ -3,6 +3,7 @@
 #define _XDG_RAY_H
 
 #include <set>
+#include <vector>
 
 #include "xdg/vec3da.h"
 
@@ -13,6 +14,9 @@
 
 namespace xdg {
 
+// forward declaration
+class TriangleRef;
+
 // TO-DO: there should be a few more double elements here (barycentric coords)
 
 /*! Stucture that is an extension of Embree's RTCRay with
@@ -23,7 +27,8 @@ struct RTCDRay: RTCRay {
 
   RTCDRay() {
     this->tnear = 0.0;
-    this->tfar = INFTY;
+    this->tfar = INFTYF;
+    this->mask = -1;
   }
 
   //! \brief Set both the single and double precision versions of the ray origin
@@ -64,7 +69,7 @@ struct RTCDRay: RTCRay {
 
   //! \brief Set both the single and double precision versions of the ray max distance
   void set_tfar(double d) {
-    tfar = d;
+    tfar = std::min(d, INFTYF);
     dtfar = d;
   }
 
@@ -86,10 +91,14 @@ struct RTCDHit : RTCHit {
   RTCDHit() {
     this->geomID = RTC_INVALID_GEOMETRY_ID;
     this->primID = RTC_INVALID_GEOMETRY_ID;
+    this->Ng_x = 0.0;
+    this->Ng_y = 0.0;
+    this->Ng_z = 0.0;
   }
 
   // data members
   const PrimitiveRef* primitive_ref {nullptr}; //!< Pointer to the primitive reference for this hit
+  MeshID surface {ID_NONE}; //!< ID of the surface this hit belongs to
   Vec3da dNg; //!< Double precision version of the primitive normal
 };
 
@@ -115,7 +124,7 @@ struct RTCDPointQuery : RTCPointQuery {
 
   //! \brief Set both the single and double precision versions of the query radius
   void set_radius(double rad) {
-    radius = rad;
+    radius = std::min(rad, INFTYF);
     dradius = rad;
   }
 

include/xdg/ray_tracing_interface.h

@@ -11,6 +11,8 @@
 #include "xdg/primitive_ref.h"
 #include "xdg/geometry_data.h"
 
+
+
 namespace xdg
 {
 
@@ -30,17 +32,15 @@ class RayTracer {
   TreeID register_volume(const std::shared_ptr<MeshManager> mesh_manager, MeshID volume);
 
   // Query Methods
-
   bool point_in_volume(TreeID scene,
                        const Position& point,
                        const Direction* direction = nullptr,
-                       const std::vector<MeshID>* exclude_primitives = nullptr);
+                       const std::vector<MeshID>* exclude_primitives = nullptr) const;
 
-  void ray_fire(TreeID scene,
-                const Position& origin,
-                const Direction& direction,
-                double& distance,
-                const std::vector<MeshID>* exclude_primitives = nullptr);
+  std::pair<double, MeshID> ray_fire(TreeID scene,
+                                     const Position& origin,
+                                     const Direction& direction,
+                                     std::vector<MeshID>* const exclude_primitives = nullptr);
 
   void closest(TreeID scene,
                const Position& origin,
@@ -59,7 +59,7 @@ class RayTracer {
 // Accessors
   int num_registered_scenes() const { return scenes_.size(); }
 
-  const GeometryUserData& geometry_data(MeshID surface) const { return user_data_map_.at(surface_to_geometry_map_.at(surface)); }
+  const std::shared_ptr<GeometryUserData>& geometry_data(MeshID surface) const { return user_data_map_.at(surface_to_geometry_map_.at(surface)); }
 
 // Data members
 private:
@@ -74,7 +74,7 @@ class RayTracer {
   RTCScene gloabal_scene_;
 
   // Internal Embree Mappings
-  std::unordered_map<RTCGeometry, GeometryUserData> user_data_map_;
+  std::unordered_map<RTCGeometry, std::shared_ptr<GeometryUserData>> user_data_map_;
 
   // Internal parameters
   double numerical_precision_ {1e-3};

include/xdg/vec3da.h

@@ -150,7 +150,7 @@ __forceinline std::ostream& operator <<(std::ostream &os, Vec3da  const& v) {
 inline bool lower(const Vec3da& a, const Vec3da& b)
 {
   for (int i = 0; i < 3; i++)
-    if (a[i] != b[i])
+  if (a[i] != b[i])
       return a[i] < b[i];
   return false;
 }
@@ -160,6 +160,15 @@ using Vertex = Vec3da;
 using Position = Vec3da;
 using Direction = Vec3da;
 
+inline Direction rand_dir() {
+  double theta = drand48() * 2.0 * M_PI;
+  double u = (1.0 - drand48()) - 1.0;
+  double phi = acos(u);
+  return Direction(sin(phi) * cos(theta), sin(phi) * sin(theta), cos(phi)).normalize();
+
+}
+
+
 } // end namespace xdg
 
 #endif

include/xdg/xdg.h

@@ -29,15 +29,14 @@ MeshID find_volume(const Position& point,
                    const Direction& direction) const;
 
 bool point_in_volume(MeshID volume,
-                     const Position& point,
-                     const Direction* direction = nullptr,
-                     const std::vector<MeshID>* exclude_primitives = nullptr) const;
+                          const Position point,
+                          const Direction* direction = nullptr,
+                          const std::vector<MeshID>* exclude_primitives = nullptr) const;
 
-void ray_fire(MeshID volume,
-              const Position& origin,
-              const Direction& direction,
-              double& distance,
-              const std::vector<MeshID>* exclude_primitives = nullptr) const;
+std::pair<double, MeshID> ray_fire(MeshID volume,
+                                   const Position& origin,
+                                   const Direction& direction,
+                                   std::vector<MeshID>* const exclude_primitives = nullptr) const;
 
 void closest(MeshID volume,
               const Position& origin,
@@ -57,6 +56,7 @@ Direction surface_normal(MeshID surface,
                          Position point,
                          const std::vector<MeshID>* exclude_primitives = nullptr) const;
 
+
   // Geometric Measurements
   double measure_volume(MeshID volume) const;
   double measure_surface_area(MeshID surface) const;
@@ -75,21 +75,19 @@ Direction surface_normal(MeshID surface,
   const std::shared_ptr<MeshManager>& mesh_manager() const {
     return mesh_manager_;
   }
-
 // Private methods
 private:
   double _triangle_volume_contribution(const PrimitiveRef& triangle) const;
   double _triangle_area_contribution(const PrimitiveRef& triangle) const;
 
-private:
+// Data members
   const std::shared_ptr<RayTracer> ray_tracing_interface_ {std::make_shared<RayTracer>()};
   std::shared_ptr<MeshManager> mesh_manager_ {nullptr};
 
   std::map<MeshID, TreeID> volume_to_scene_map_;  //<! Map from mesh volume to embree scene
   std::map<MeshID, TreeID> surface_to_scene_map_; //<! Map from mesh surface to embree scnee
   std::map<MeshID, RTCGeometry> surface_to_geometry_map_; //<! Map from mesh surface to embree geometry
   TreeID gloabal_scene_;
-
 };
 
 }

src/mesh_manager_interface.cpp

@@ -40,7 +40,33 @@ MeshID MeshManager::next_surface_id() const
     return *std::max_element(surfaces().begin(), surfaces().end()) + 1;
 }
 
-MeshID MeshManager::next_volume(MeshID surface, MeshID current_volume) const
+bool
+MeshManager::volume_has_property(MeshID volume, PropertyType type) const
+{
+  return volume_metadata_.count({volume, type}) > 0;
+}
+
+bool
+MeshManager::surface_has_property(MeshID surface, PropertyType type) const
+{
+  return surface_metadata_.count({surface, type}) > 0;
+}
+
+Property
+MeshManager::get_volume_property(MeshID volume, PropertyType type) const
+{
+  return volume_metadata_.at({volume, type});
+}
+
+Property
+MeshManager::get_surface_property(MeshID surface, PropertyType type) const
+{
+  if (surface_metadata_.count({surface, type}) == 0)
+    return {PropertyType::BOUNDARY_CONDITION, "transmission"};
+  return surface_metadata_.at({surface, type});
+}
+
+MeshID MeshManager::next_volume(MeshID current_volume, MeshID surface) const
 {
   auto parent_vols = this->get_parent_volumes(surface);
 

src/moab/mesh_manager.cpp

@@ -299,18 +299,6 @@ MOABMeshManager::get_volume_surfaces(MeshID volume) const
   return surface_ids;
 }
 
-Property
-MOABMeshManager::get_volume_property(MeshID volume, PropertyType type) const
-{
-  return volume_metadata_.at({volume, type});
-}
-
-Property
-MOABMeshManager::get_surface_property(MeshID surface, PropertyType type) const
-{
-  return surface_metadata_.at({surface, type});
-}
-
 void
 MOABMeshManager::parse_metadata()
 {