diff --git a/CMakeLists.txt b/CMakeLists.txt index e9033e6..a5a0e93 100644 --- a/CMakeLists.txt +++ b/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 $ diff --git a/include/xdg/constants.h b/include/xdg/constants.h index 2de9986..0f59d15 100644 --- a/include/xdg/constants.h +++ b/include/xdg/constants.h @@ -9,6 +9,16 @@ namespace xdg { constexpr double INFTY {std::numeric_limits::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::max()}; +#endif + // Whether information pertains to a surface or volume enum class GeometryType { SURFACE = 2, diff --git a/include/xdg/mesh_manager_interface.h b/include/xdg/mesh_manager_interface.h index 1d83cf0..fe14552 100644 --- a/include/xdg/mesh_manager_interface.h +++ b/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& volumes() const { return volumes_; } diff --git a/include/xdg/moab/mesh_manager.h b/include/xdg/moab/mesh_manager.h index 5de4c0d..dba048f 100644 --- a/include/xdg/moab/mesh_manager.h +++ b/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; } diff --git a/include/xdg/ray.h b/include/xdg/ray.h index c2755b2..4d0dee1 100644 --- a/include/xdg/ray.h +++ b/include/xdg/ray.h @@ -3,6 +3,7 @@ #define _XDG_RAY_H #include +#include #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; } diff --git a/include/xdg/ray_tracing_interface.h b/include/xdg/ray_tracing_interface.h index 09ca831..756b555 100644 --- a/include/xdg/ray_tracing_interface.h +++ b/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 mesh_manager, MeshID volume); // Query Methods - bool point_in_volume(TreeID scene, const Position& point, const Direction* direction = nullptr, - const std::vector* exclude_primitives = nullptr); + const std::vector* exclude_primitives = nullptr) const; - void ray_fire(TreeID scene, - const Position& origin, - const Direction& direction, - double& distance, - const std::vector* exclude_primitives = nullptr); + std::pair ray_fire(TreeID scene, + const Position& origin, + const Direction& direction, + std::vector* 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& 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 user_data_map_; + std::unordered_map> user_data_map_; // Internal parameters double numerical_precision_ {1e-3}; diff --git a/include/xdg/vec3da.h b/include/xdg/vec3da.h index 3ea3a07..c082c2e 100644 --- a/include/xdg/vec3da.h +++ b/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 diff --git a/include/xdg/xdg.h b/include/xdg/xdg.h index 0a7da84..197f223 100644 --- a/include/xdg/xdg.h +++ b/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* exclude_primitives = nullptr) const; + const Position point, + const Direction* direction = nullptr, + const std::vector* exclude_primitives = nullptr) const; -void ray_fire(MeshID volume, - const Position& origin, - const Direction& direction, - double& distance, - const std::vector* exclude_primitives = nullptr) const; +std::pair ray_fire(MeshID volume, + const Position& origin, + const Direction& direction, + std::vector* 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* exclude_primitives = nullptr) const; + // Geometric Measurements double measure_volume(MeshID volume) const; double measure_surface_area(MeshID surface) const; @@ -75,13 +75,12 @@ Direction surface_normal(MeshID surface, const std::shared_ptr& 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 ray_tracing_interface_ {std::make_shared()}; std::shared_ptr mesh_manager_ {nullptr}; @@ -89,7 +88,6 @@ Direction surface_normal(MeshID surface, std::map surface_to_scene_map_; // surface_to_geometry_map_; // 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); diff --git a/src/moab/mesh_manager.cpp b/src/moab/mesh_manager.cpp index bb76d49..19be979 100644 --- a/src/moab/mesh_manager.cpp +++ b/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() { diff --git a/src/ray_tracing_interface.cpp b/src/ray_tracing_interface.cpp index 11cc08e..53aa23d 100644 --- a/src/ray_tracing_interface.cpp +++ b/src/ray_tracing_interface.cpp @@ -79,6 +79,8 @@ RayTracer::register_volume(const std::shared_ptr mesh_manager, double max_distance = std::sqrt(dx*dx + dy*dy + dz*dz); double bump = max_distance * std::pow(10, -std::numeric_limits::digits10); + bump = std::max(bump, 1e-03); + bump = 1.0; // create a new geometry for each surface int buffer_start = 0; @@ -89,16 +91,17 @@ RayTracer::register_volume(const std::shared_ptr mesh_manager, unsigned int embree_surface = rtcAttachGeometry(volume_scene, surface_geometry); this->surface_to_geometry_map_[surface] = surface_geometry; - GeometryUserData surface_data; - surface_data.surface_id = surface; - surface_data.mesh_manager = mesh_manager.get(); - surface_data.prim_ref_buffer = tri_ref_ptr + buffer_start; + std::shared_ptr surface_data = std::make_shared(); + surface_data->surface_id = surface; + surface_data->mesh_manager = mesh_manager.get(); + surface_data->prim_ref_buffer = tri_ref_ptr + buffer_start; user_data_map_[surface_geometry] = surface_data; - rtcSetGeometryUserData(surface_geometry, &user_data_map_[surface_geometry]); + // TODO: This could be a problem if user_data_map_ is reallocated? + rtcSetGeometryUserData(surface_geometry, user_data_map_[surface_geometry].get()); for (int i = 0; i < surface_triangles.size(); ++i) { - auto& triangle_ref = surface_data.prim_ref_buffer[i]; + auto& triangle_ref = surface_data->prim_ref_buffer[i]; // triangle_ref.embree_surface = embree_surface; } buffer_start += surface_triangles.size(); @@ -118,7 +121,7 @@ RayTracer::register_volume(const std::shared_ptr mesh_manager, bool RayTracer::point_in_volume(TreeID scene, const Position& point, const Direction* direction, - const std::vector* exclude_primitives) + const std::vector* exclude_primitives) const { RTCDRayHit rayhit; rayhit.ray.set_org(point); @@ -128,6 +131,7 @@ bool RayTracer::point_in_volume(TreeID scene, rayhit.ray.orientation = HitOrientation::ANY; rayhit.ray.set_tfar(INFTY); rayhit.ray.set_tnear(0.0); + if (exclude_primitives != nullptr) rayhit.ray.exclude_primitives = exclude_primitives; { @@ -142,12 +146,11 @@ bool RayTracer::point_in_volume(TreeID scene, return rayhit.ray.ddir.dot(rayhit.hit.dNg) > 0.0; } -void +std::pair RayTracer::ray_fire(TreeID scene, const Position& origin, const Direction& direction, - double& distance, - const std::vector* exclude_primitves) + std::vector* const exclude_primitves) { RTCDRayHit rayhit; // set ray data @@ -170,9 +173,10 @@ RayTracer::ray_fire(TreeID scene, } if (rayhit.hit.geomID == RTC_INVALID_GEOMETRY_ID) - distance = INFTY; + return {INFTY, ID_NONE}; else - distance = rayhit.ray.dtfar; + if (exclude_primitves) exclude_primitves->push_back(rayhit.hit.primitive_ref->primitive_id); + return {rayhit.ray.dtfar, rayhit.hit.surface}; } void RayTracer::closest(TreeID scene, diff --git a/src/triangle_intersect.cpp b/src/triangle_intersect.cpp index 2174cfd..66e4d9d 100644 --- a/src/triangle_intersect.cpp +++ b/src/triangle_intersect.cpp @@ -87,11 +87,14 @@ void TriangleIntersectionFunc(RTCIntersectFunctionNArguments* args) { // zero-out barycentric coords rayhit->hit.u = 0.0; rayhit->hit.v = 0.0; + rayhit->hit.Ng_x = 0.0; + rayhit->hit.Ng_y = 0.0; + rayhit->hit.Ng_z = 0.0; // set the hit information rayhit->hit.geomID = args->geomID; rayhit->hit.primID = args->primID; rayhit->hit.primitive_ref = &primitive_ref; - + rayhit->hit.surface = user_data->surface_id; rayhit->hit.dNg = normal; } diff --git a/src/xdg.cpp b/src/xdg.cpp index 840d3d9..c8c9493 100644 --- a/src/xdg.cpp +++ b/src/xdg.cpp @@ -33,6 +33,15 @@ std::shared_ptr XDG::create(MeshLibrary library) return xdg; } +bool XDG::point_in_volume(MeshID volume, + const Position point, + const Direction* direction, + const std::vector* exclude_primitives) const +{ + TreeID scene = volume_to_scene_map_.at(volume); + return ray_tracing_interface()->point_in_volume(scene, point, direction, exclude_primitives); +} + MeshID XDG::find_volume(const Position& point, const Direction& direction) const { @@ -46,23 +55,14 @@ MeshID XDG::find_volume(const Position& point, return ID_NONE; } -bool XDG::point_in_volume(MeshID volume, - const Position& point, - const Direction* direction, - const std::vector* exclude_primitives) const -{ - TreeID scene = volume_to_scene_map_.at(volume); - return ray_tracing_interface()->point_in_volume(scene, point, direction, exclude_primitives); -} - -void XDG::ray_fire(MeshID volume, +std::pair +XDG::ray_fire(MeshID volume, const Position& origin, const Direction& direction, - double& distance, - const std::vector* exclude_primitives) const + std::vector* const exclude_primitives) const { TreeID scene = volume_to_scene_map_.at(volume); - ray_tracing_interface()->ray_fire(scene, origin, direction, distance, exclude_primitives); + return ray_tracing_interface()->ray_fire(scene, origin, direction, exclude_primitives); } void XDG::closest(MeshID volume, @@ -96,7 +96,7 @@ Direction XDG::surface_normal(MeshID surface, const std::vector* exclude_primitives) const { MeshID element; - if (exclude_primitives != nullptr) { + if (exclude_primitives != nullptr && exclude_primitives->size() > 0) { element = exclude_primitives->back(); } else { auto surface_vols = mesh_manager()->get_parent_volumes(surface); diff --git a/tests/mesh_mock.h b/tests/mesh_mock.h index 927c99e..14fc9b1 100644 --- a/tests/mesh_mock.h +++ b/tests/mesh_mock.h @@ -120,16 +120,6 @@ class MeshMock : public MeshManager { fatal_error("MockMesh does not support parse_metadata()"); } - virtual Property get_volume_property(MeshID volume, PropertyType type) const override { - fatal_error("MockMesh does not support get_volume_property()"); - return Property(); - } - - virtual Property get_surface_property(MeshID surface, PropertyType type) const override { - fatal_error("MockMesh does not support get_surface_property()"); - return Property(); - } - // Other virtual MeshLibrary mesh_library() const override { return MeshLibrary::INTERNAL; } diff --git a/tests/test_moab.cpp b/tests/test_moab.cpp index 6ba9eea..236b41d 100644 --- a/tests/test_moab.cpp +++ b/tests/test_moab.cpp @@ -84,20 +84,24 @@ TEST_CASE("Test Ray Fire MOAB") Position origin {0.0, 0.0, 0.0}; Direction direction {1.0, 0.0, 0.0}; - double intersection_distance {0.0}; + std::pair intersection; - xdg->ray_fire(volume, origin, direction, intersection_distance); + intersection = xdg->ray_fire(volume, origin, direction); // this cube is 10 cm on a side, so the ray should hit the surface at 5 cm - REQUIRE_THAT(intersection_distance, Catch::Matchers::WithinAbs(5.0, 1e-6)); + REQUIRE_THAT(intersection.first, Catch::Matchers::WithinAbs(5.0, 1e-6)); origin = {3.0, 0.0, 0.0}; - xdg->ray_fire(volume, origin, direction, intersection_distance); - REQUIRE_THAT(intersection_distance, Catch::Matchers::WithinAbs(2.0, 1e-6)); + intersection = xdg->ray_fire(volume, origin, direction); + REQUIRE_THAT(intersection.first, Catch::Matchers::WithinAbs(2.0, 1e-6)); origin = {-10.0, 0.0, 0.0}; - xdg->ray_fire(volume, origin, direction, intersection_distance); - REQUIRE_THAT(intersection_distance, Catch::Matchers::WithinAbs(15.0, 1e-6)); + intersection = xdg->ray_fire(volume, origin, direction); + REQUIRE_THAT(intersection.first, Catch::Matchers::WithinAbs(15.0, 1e-6)); + + origin = {0.0, 0.0, 0.0}; + REQUIRE(xdg->point_in_volume(volume, origin)); + } TEST_CASE("TEST XDG Factory Method") diff --git a/tests/test_point_in_volume.cpp b/tests/test_point_in_volume.cpp index 0c00a97..c674b8c 100644 --- a/tests/test_point_in_volume.cpp +++ b/tests/test_point_in_volume.cpp @@ -27,7 +27,7 @@ TEST_CASE("Test Point in Volume") REQUIRE(result == false); // test a point just inside the positive x boundary - point = {5.0 - 1e-06, 0.0, 0.0}; + point = {4.0 - 1e-06, 0.0, 0.0}; result = rti->point_in_volume(volume_tree, point); REQUIRE(result == true); diff --git a/tests/test_ray_fire.cpp b/tests/test_ray_fire.cpp index 1bd6914..1adec46 100644 --- a/tests/test_ray_fire.cpp +++ b/tests/test_ray_fire.cpp @@ -23,41 +23,41 @@ TEST_CASE("Test Ray Fire Mesh Mock") Position origin {0.0, 0.0, 0.0}; Direction direction {1.0, 0.0, 0.0}; - double intersection_distance {0.0}; + std::pair intersection; // fire from the origin toward each face, ensuring that the intersection distances are correct - rti->ray_fire(volume_tree, origin, direction, intersection_distance); - REQUIRE_THAT(intersection_distance, Catch::Matchers::WithinAbs(5.0, 1e-6)); + intersection = rti->ray_fire(volume_tree, origin, direction); + REQUIRE_THAT(intersection.first, Catch::Matchers::WithinAbs(5.0, 1e-6)); direction *= -1; - rti->ray_fire(volume_tree, origin, direction, intersection_distance); - REQUIRE_THAT(intersection_distance, Catch::Matchers::WithinAbs(2.0, 1e-6)); + intersection = rti->ray_fire(volume_tree, origin, direction); + REQUIRE_THAT(intersection.first, Catch::Matchers::WithinAbs(2.0, 1e-6)); direction = {0.0, 1.0, 0.0}; - rti->ray_fire(volume_tree, origin, direction, intersection_distance); - REQUIRE_THAT(intersection_distance, Catch::Matchers::WithinAbs(6.0, 1e-6)); + intersection = rti->ray_fire(volume_tree, origin, direction); + REQUIRE_THAT(intersection.first, Catch::Matchers::WithinAbs(6.0, 1e-6)); direction *= -1; - rti->ray_fire(volume_tree, origin, direction, intersection_distance); - REQUIRE_THAT(intersection_distance, Catch::Matchers::WithinAbs(3.0, 1e-6)); + intersection = rti->ray_fire(volume_tree, origin, direction); + REQUIRE_THAT(intersection.first, Catch::Matchers::WithinAbs(3.0, 1e-6)); direction = {0.0, 0.0, 1.0}; - rti->ray_fire(volume_tree, origin, direction, intersection_distance); - REQUIRE_THAT(intersection_distance, Catch::Matchers::WithinAbs(7.0, 1e-6)); + intersection = rti->ray_fire(volume_tree, origin, direction); + REQUIRE_THAT(intersection.first, Catch::Matchers::WithinAbs(7.0, 1e-6)); direction *= -1; - rti->ray_fire(volume_tree, origin, direction, intersection_distance); - REQUIRE_THAT(intersection_distance, Catch::Matchers::WithinAbs(4.0, 1e-6)); + intersection = rti->ray_fire(volume_tree, origin, direction); + REQUIRE_THAT(intersection.first, Catch::Matchers::WithinAbs(4.0, 1e-6)); // fire from the outside of the cube toward each face, ensuring that the intersection distances are correct // rays should skip entering intersections and intersect with the far side of the cube origin = {-10.0, 0.0, 0.0}; direction = {1.0, 0.0, 0.0}; - rti->ray_fire(volume_tree, origin, direction, intersection_distance); - REQUIRE_THAT(intersection_distance, Catch::Matchers::WithinAbs(15.0, 1e-6)); + intersection = rti->ray_fire(volume_tree, origin, direction); + REQUIRE_THAT(intersection.first, Catch::Matchers::WithinAbs(15.0, 1e-6)); origin = {10.0, 0.0, 0.0}; direction = {-1.0, 0.0, 0.0}; - rti->ray_fire(volume_tree, origin, direction, intersection_distance); - REQUIRE_THAT(intersection_distance, Catch::Matchers::WithinAbs(12.0, 1e-6)); + intersection = rti->ray_fire(volume_tree, origin, direction); + REQUIRE_THAT(intersection.first, Catch::Matchers::WithinAbs(12.0, 1e-6)); } \ No newline at end of file diff --git a/tools/CMakeLists.txt b/tools/CMakeLists.txt new file mode 100644 index 0000000..ecc6051 --- /dev/null +++ b/tools/CMakeLists.txt @@ -0,0 +1,4 @@ + + +add_executable(particle-sim particle_sim.cpp) +target_link_libraries(particle-sim xdg) \ No newline at end of file diff --git a/tools/particle_sim.cpp b/tools/particle_sim.cpp new file mode 100644 index 0000000..57e3985 --- /dev/null +++ b/tools/particle_sim.cpp @@ -0,0 +1,170 @@ +#include +#include +#include +#include + +#include "xdg/error.h" +#include "xdg/mesh_manager_interface.h" +#include "xdg/moab/mesh_manager.h" +#include "xdg/vec3da.h" +#include "xdg/xdg.h" + +using namespace xdg; + +static double mfp {1.0}; + +struct Particle { + +Particle(std::shared_ptr xdg, uint32_t id, bool verbose=true) : verbose_(verbose), xdg_(xdg), id_(id) {} + +template +void log (const std::string& msg, const Params&... fmt_args) { + if (!verbose_) return; + write_message(msg, fmt_args...); +} + +void initialize() { + // TODO: replace with sampling + r_ = {0.0, 0.0, 0.0}; + u_ = {1.0, 0.0, 0.0}; + + volume_ = xdg_->find_volume(r_, u_); +} + +void surf_dist() { + surface_intersection_ = xdg_->ray_fire(volume_, r_, u_, &history_); + if (surface_intersection_.first == 0.0) { + fatal_error("Particle {} stuck at position ({}, {}, {}) on surfacce {}", id_, r_.x, r_.y, r_.z, surface_intersection_.second); + alive_ = false; + return; + } + if (surface_intersection_.second == ID_NONE) { + fatal_error("Particle {} lost in volume {}", id_, volume_); + alive_ = false; + return; + } + log("Intersected surface {} at distance {} ", surface_intersection_.second, surface_intersection_.first); +} + +void sample_collision_distance() { + collision_distance_ = -std::log(1.0 - drand48()) / mfp; +} + +void collide() { + n_events_++; + log("Event {} for particle {}", n_events_, id_); + u_ = rand_dir(); + log("Particle {} collides with material at position ({}, {}, {}), new direction is ({}, {}, {})", id_, r_.x, r_.y, r_.z, u_.z, u_.y, u_.z); + history_.clear(); +} + +void advance() +{ + log("Comparing surface intersection distance {} to collision distance {}", surface_intersection_.first, collision_distance_); + if (collision_distance_ < surface_intersection_.first) { + r_ += collision_distance_ * u_; + log("Particle {} collides with material at position ({}, {}, {}) ", id_, r_.x, r_.y, r_.z); + + } else { + r_ += surface_intersection_.first * u_; + log("Particle {} advances to surface {} at position ({}, {}, {}) ", id_, surface_intersection_.second, r_.x, r_.y, r_.z); + } +} + +void cross_surface() +{ + n_events_++; + log("Event {} for particle {}", n_events_, id_); + // check for the surface boundary condition + if (xdg_->mesh_manager()->get_surface_property(surface_intersection_.second, PropertyType::BOUNDARY_CONDITION).value == "reflecting") { + log("Particle {} reflects off surface {}", id_, surface_intersection_.second); + log("Direction before reflection: ({}, {}, {})", u_.x, u_.y, u_.z); + + Direction normal = xdg_->surface_normal(surface_intersection_.second, r_, &history_); + log("Normal to surface: ({}, {}, {})", normal.x, normal.y, normal.z); + + double proj = dot(normal, u_); + double mag = normal.length(); + normal = normal * (2.0 * proj/mag); + u_ = u_ - normal; + u_ = u_.normalize(); + log("Direction after reflection: ({}, {}, {})", u_.x, u_.y, u_.z); + // reset to last intersection + if (history_.size() > 0) { + log("Resetting particle history to last intersection"); + history_ = {history_.back()}; + } + } else { + volume_ = xdg_->mesh_manager()->next_volume(volume_, surface_intersection_.second); + log("Particle {} enters volume {}", id_, volume_); + if (volume_ == ID_NONE) { + alive_ = false; + return; + } + } +} + +// Data Members +bool verbose_ {true}; +std::shared_ptr xdg_; +uint32_t id_ {0}; +Position r_; +Direction u_; +MeshID volume_ {ID_NONE}; +std::vector history_{}; + +std::pair surface_intersection_ {INFTY, ID_NONE}; +double collision_distance_ {INFTY}; + +int32_t n_events_ {0}; +bool alive_ {true}; +}; + +int main(int argc, char** argv) { + +srand48(42); + +// create a mesh manager +std::shared_ptr xdg = XDG::create(MeshLibrary::MOAB); +const auto& mm = xdg->mesh_manager(); + +std::string filename {argv[1]}; + +mm->load_file(filename); +mm->init(); +mm->parse_metadata(); +xdg->prepare_raytracer(); + +// create a new particle + +const int n_particles {100}; + +const int max_events {1000}; + +bool verbose = true; + + for (int i = 0; i < n_particles; i++) { + write_message("Starting particle {}", i); + Particle p(xdg, i, verbose); + p.initialize(); + while (true) { + p.surf_dist(); + // terminate for leakage + if (!p.alive_) break; + p.sample_collision_distance(); + p.advance(); + if (p.surface_intersection_.first < p.collision_distance_) + p.cross_surface(); + else + p.collide(); + if (!p.alive_) break; + + if (p.n_events_ > max_events) { + write_message("Maximum number of events ({}) reached", max_events); + break; + } + } + } + + return 0; +}