Adding surface and volume measurement methods (+basic test)

CMakeLists.txt

@@ -41,13 +41,15 @@ if (NOT fmt_FOUND)
 endif()
 
 list(APPEND xdg_sources
-src/closest.cpp
+src/geometry/measure.cpp
+src/geometry/plucker.cpp
+src/geometry/closest.cpp
 src/error.cpp
 src/mesh_manager_interface.cpp
-src/plucker.cpp
 src/ray_tracing_interface.cpp
 src/triangle_ref.cpp
 src/util/str_utils.cpp
+src/xdg.cpp
 )
 
 if (XDG_ENABLE_MOAB)

include/xdg/geometry/measure.h

@@ -0,0 +1,16 @@
+#ifndef _XDG_AREA_H
+#define _XDG_AREA_H
+
+#include "xdg/vec3da.h"
+
+namespace xdg {
+
+double triangle_volume_contribution(const std::array<Vertex, 3>& vertices);
+double triangle_volume_contribution(const Vertex& v0, const Vertex& v1, const Vertex& v2);
+
+double triangle_area(const std::array<Vertex, 3>& vertices);
+double triangle_area(const Vertex& v0, const Vertex& v1, const Vertex& v2);
+
+}
+
+#endif // include guard

include/xdg/mesh_manager_interface.h

@@ -56,6 +56,8 @@ class MeshManager {
 
   virtual std::vector<MeshID> get_volume_surfaces(MeshID volume) const = 0;
 
+  virtual Sense surface_sense(MeshID surface, MeshID volume) const = 0;
+
   virtual MeshID create_volume() = 0;
 
   virtual void add_surface_to_volume(MeshID volume, MeshID surface, Sense sense, bool overwrite=false) = 0;

include/xdg/moab/mesh_manager.h

@@ -76,6 +76,8 @@ class MOABMeshManager : public MeshManager {
 
   std::vector<MeshID> get_volume_surfaces(MeshID volume) const override;
 
+  Sense surface_sense(MeshID surface, MeshID volume) const override;
+
   // Metadata
   void parse_metadata() override;
 

include/xdg/xdg.h

@@ -11,29 +11,45 @@ namespace xdg {
 class XDG {
 
 public:
+// Constructors
+  XDG() = default;
+
+  XDG(std::shared_ptr<MeshManager> mesh_manager) :
+    mesh_manager_(mesh_manager) {}
+
 
 // Methods
   void prepare_raytracer() {
-    ray_tracing_interface_->register_all_volumes(mesh_manager_interface_);
+    ray_tracing_interface_->register_all_volumes(mesh_manager_);
   }
 
+  double measure_volume(MeshID volume) const;
+
+  double measure_surface_area(MeshID surface) const;
+  double measure_volume_area(MeshID surface) const;
+
 // Mutators
-  void set_mesh_manager_interface(std::shared_ptr<MeshManager> mesh_manager_interface) {
-    mesh_manager_interface_ = mesh_manager_interface;
+  void set_mesh_manager_interface(std::shared_ptr<MeshManager> mesh_manager) {
+    mesh_manager_ = mesh_manager;
   }
 
 // Accessors
   const RayTracer* ray_tracing_interface() const {
     return ray_tracing_interface_.get();
   }
 
-  const MeshManager* mesh_manager_interface() const {
-    return mesh_manager_interface_.get();
+  const MeshManager* mesh_manager() const {
+    return mesh_manager_.get();
   }
 
+// Private methods
+private:
+  double _triangle_volume_contribution(const TriangleRef& triangle) const;
+  double _triangle_area_contribution(const TriangleRef& triangle) const;
+
 private:
   const std::shared_ptr<RayTracer> ray_tracing_interface_ {std::make_shared<RayTracer>()};
-  std::shared_ptr<MeshManager> mesh_manager_interface_ {nullptr};
+  std::shared_ptr<MeshManager> mesh_manager_ {nullptr};
 };
 
 }

src/geometry/measure.cpp

@@ -0,0 +1,23 @@
+#include "xdg/geometry/measure.h"
+
+namespace xdg {
+  double triangle_volume_contribution(const std::array<Vertex, 3>& vertices)
+  {
+    return triangle_volume_contribution(vertices[0], vertices[1], vertices[2]);
+  }
+
+  double triangle_volume_contribution(const Vertex& v0, const Vertex& v1, const Vertex& v2)
+  {
+    return v0.dot((v1-v0).cross(v2-v0));
+  }
+
+  double triangle_area(const std::array<Vertex, 3>& vertices)
+  {
+    return triangle_area(vertices[0], vertices[1], vertices[2]);
+  }
+
+  double triangle_area(const Vertex& v0, const Vertex& v1, const Vertex& v2)
+  {
+    return 0.5 * (v1-v0).cross(v2-v0).length();
+  }
+}

src/moab/mesh_manager.cpp

@@ -258,6 +258,16 @@ MOABMeshManager::get_parent_volumes(MeshID surface) const
   return this->surface_senses(surface);
 }
 
+Sense
+MOABMeshManager::surface_sense(MeshID surface, MeshID volume) const
+{
+  auto sense_data = this->get_parent_volumes(surface);
+  if (sense_data.first == volume) return Sense::FORWARD;
+  else if (sense_data.second == volume) return Sense::REVERSE;
+  else fatal_error("Volume {} is not a parent of surface {}", volume, surface);
+  return Sense::UNSET;
+}
+
 std::vector<moab::EntityHandle>
 MOABMeshManager::_ents_of_dim(int dim) const {
   std::array<moab::Tag, 1> tags = {geometry_dimension_tag_};

src/ray_tracing_interface.cpp

@@ -258,12 +258,10 @@ Direction RayTracer::get_normal(MeshID surface,
     if (triangle_ref.surface_id != surface) {
       fatal_error("Point {} was closest to surface {}, not surface {}, in volume {}.", point, triangle_ref.surface_id, surface, surface_vols.first);
     }
-
     element = triangle_ref.triangle_id;
-
   }
 
-  // set the normal based on the triangle
+  // return the normal of the selected triangle
   return mesh_manager->triangle_normal(element);
 }
 

src/xdg.cpp

@@ -0,0 +1,52 @@
+#include <vector>
+
+#include "xdg/xdg.h"
+
+#include "xdg/constants.h"
+#include "xdg/geometry/measure.h"
+namespace xdg {
+
+double XDG::measure_volume(MeshID volume) const
+{
+  double volume_total {0.0};
+
+  auto surfaces = mesh_manager()->get_volume_surfaces(volume);
+
+  std::vector<Sense> surface_senses;
+  for (auto surface : surfaces) {
+    surface_senses.push_back(mesh_manager()->surface_sense(surface, volume));
+  }
+
+  for (int i = 0; i < surfaces.size(); ++i) {
+    MeshID& surface = surfaces[i];
+    double surface_contribution {0.0};
+    auto triangles = mesh_manager()->get_surface_elements(surface);
+    for (auto triangle : triangles) {
+      surface_contribution += triangle_volume_contribution(mesh_manager()->triangle_vertices(triangle));
+    }
+    if (surface_senses[i] == Sense::REVERSE) surface_contribution *= -1.0;
+    volume_total += surface_contribution;
+  }
+
+  return volume_total / 6.0;
+}
+
+double XDG::measure_surface_area(MeshID surface) const
+{
+  double area {0.0};
+  for (auto triangle : mesh_manager()->get_surface_elements(surface)) {
+    area += triangle_area(mesh_manager()->triangle_vertices(triangle));
+  }
+  return area;
+}
+
+double XDG::measure_volume_area(MeshID volume) const
+{
+  double area {0.0};
+  for (auto surface : mesh_manager()->get_volume_surfaces(volume)) {
+    area += measure_surface_area(surface);
+  }
+  return area;
+}
+
+} // namespace xdg

tests/CMakeLists.txt

@@ -10,6 +10,7 @@ test_occluded
 test_ray_fire
 test_point_in_volume
 test_normal
+test_measure
 test_xdg_interface
 )
 

tests/mesh_mock.h

@@ -103,6 +103,10 @@ class MeshMock : public MeshManager {
     return {0, 1, 2, 3, 4, 5};
   }
 
+  Sense surface_sense(MeshID surface, MeshID volume) const override {
+    return Sense::FORWARD;
+  }
+
   virtual MeshID create_volume() override {
     fatal_error("MockMesh does not support create_volume()");
     return ID_NONE;

tests/test_bvh.cpp

@@ -1,5 +1,3 @@
-
-
 // for testing
 #include <catch2/catch_test_macros.hpp>
 

tests/test_measure.cpp

@@ -0,0 +1,36 @@
+#include <memory>
+#include <vector>
+
+// for testing
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_floating_point.hpp>
+
+#include "xdg/xdg.h"
+
+// xdg test includes
+#include "mesh_mock.h"
+#include "util.h"
+
+using namespace xdg;
+
+TEST_CASE("Test Mesh Mock")
+{
+  std::shared_ptr<MeshManager> mm = std::make_shared<MeshMock>();
+  mm->init(); // this should do nothing, but its good practice to call it
+
+  XDG xdg{mm};
+
+  double volume = xdg.measure_volume(mm->volumes()[0]);
+  REQUIRE_THAT(volume, Catch::Matchers::WithinAbs(693., 1e-6));
+
+  double area = xdg.measure_volume_area(mm->volumes()[0]);
+  REQUIRE_THAT(area, Catch::Matchers::WithinAbs(478., 1e-6));
+
+  std::vector<double> surface_areas = {63., 63., 99., 99., 77., 77.};
+
+  for (int i = 0; i < mm->surfaces().size(); ++i) {
+    double area = xdg.measure_surface_area(mm->surfaces()[i]);
+    REQUIRE_THAT(area, Catch::Matchers::WithinAbs(surface_areas[i], 1e-6));
+  }
+
+}

tests/test_xdg_interface.cpp

@@ -12,5 +12,5 @@ TEST_CASE("XDG Interface")
 {
   std::shared_ptr<XDG> xdg = std::make_shared<XDG>();
   REQUIRE(xdg->ray_tracing_interface() != nullptr);
-  REQUIRE(xdg->mesh_manager_interface() == nullptr);
+  REQUIRE(xdg->mesh_manager() == nullptr);
 }