Finishing concrete implementation.

include/openmc/mesh.h

@@ -18,6 +18,7 @@
 #ifdef DAGMC
 #include "moab/Core.hpp"
 #include "moab/AdaptiveKDTree.hpp"
+#include "moab/Matrix3.hpp"
 #endif
 
 namespace openmc {
@@ -26,11 +27,13 @@ namespace openmc {
 // Global variables
 //==============================================================================
 
+
 class Mesh;
 
 namespace model {
 
 extern std::vector<std::unique_ptr<Mesh>> meshes;
+
 extern std::unordered_map<int32_t, int32_t> mesh_map;
 
 } // namespace model
@@ -50,14 +53,6 @@ public:
   virtual void bins_crossed(const Particle* p, std::vector<int>& bins,
                             std::vector<double>& lengths) const = 0;
 
-  //! Check where a line segment intersects the mesh and if it intersects at all
-  //
-  //! \param[in,out] r0 In: starting position, out: intersection point
-  //! \param[in] r1 Ending position
-  //! \param[out] ijk Indices of the mesh bin containing the intersection point
-  //! \return Whether the line segment connecting r0 and r1 intersects mesh
-  virtual bool intersects(Position& r0, Position r1, int* ijk) const = 0;
-
   //! Write mesh data to an HDF5 group
   //
   //! \param[in] group HDF5 group
@@ -69,14 +64,6 @@ public:
   //! \return Mesh bin
   virtual int get_bin(Position r) const = 0;
 
-  //! Count number of bank sites in each mesh bin / energy bin
-  //
-  //! \param[in] bank Array of bank sites
-  //! \param[out] Whether any bank sites are outside the mesh
-  //! \return Array indicating number of sites in each mesh/energy bin
-  virtual xt::xarray<double> count_sites(const std::vector<Particle::Bank>& bank,
-                                         bool* outside) const = 0;
-
   int id_ {-1};  //!< User-specified ID
   int n_dimension_; //!< Number of dimensions
 };
@@ -281,6 +268,7 @@ private:
 #ifdef DAGMC
 
 class UnstructuredMesh : public Mesh {
+public:
   UnstructuredMesh() { };
   UnstructuredMesh(pugi::xml_node);
 
@@ -294,24 +282,20 @@ class UnstructuredMesh : public Mesh {
 
   bool intersects(Position& r0, Position r1, int* ijk);
 
+  bool point_in_tet(const Position& r, moab::EntityHandle tet) const;
+
   //! Write mesh data to an HDF5 group
   //
   //! \param[in] group HDF5 group
-  void to_hdf5(hid_t group);
+  void to_hdf5(hid_t group) const;
 
   //! Get bin at a given position in space
   //
   //! \param[in] r Position to get bin for
   //! \return Mesh bin
-  int get_bin(Position r);
+  int get_bin(Position r) const;
 
-  //! Count number of bank sites in each mesh bin / energy bin
-  //
-  //! \param[in] bank Array of bank sites
-  //! \param[out] Whether any bank sites are outside the mesh
-  //! \return Array indicating number of sites in each mesh/energy bin
-  xt::xarray<double> count_sites(const std::vector<Particle::Bank>& bank,
-                                         bool* outside);
+  void compute_barycentric_data(const moab::Range& all_tets);
 
   int get_bin_from_ent_handle(moab::EntityHandle eh) const;
 
@@ -325,6 +309,7 @@ private:
   moab::EntityHandle meshset_;
   std::shared_ptr<moab::Interface> mbi_;
   std::unique_ptr<moab::AdaptiveKDTree> kdtree_;
+  std::vector<moab::Matrix3> baryc_data_;
 };
 
 #endif

src/mesh.cpp

@@ -32,6 +32,7 @@ namespace openmc {
 
 namespace model {
 
+
 std::vector<std::unique_ptr<Mesh>> meshes;
 std::unordered_map<int32_t, int32_t> mesh_map;
 
@@ -1605,6 +1606,89 @@ UnstructuredMesh::bins_crossed(const Particle* p, std::vector<int>& bins,
 
 };
 
+
+
+int
+UnstructuredMesh::get_bin(Position r) const {
+  moab::CartVect pnt(r.x, r.y, r.z);
+  moab::AdaptiveKDTreeIter kdtree_iter;
+  moab::ErrorCode rval = kdtree_->point_search(pnt.array(), kdtree_iter);
+  if (rval != moab::MB_SUCCESS) { return -1; }
+
+  moab::EntityHandle leaf = kdtree_iter.handle();
+  moab::Range tets;
+  rval = mbi_->get_entities_by_dimension(leaf, 3, tets);
+  for (const auto& tet : tets) {
+    if (point_in_tet(r, tet)) {
+      return get_bin_from_ent_handle(tet);
+    }
+  }
+
+  return -1;
+}
+
+void
+UnstructuredMesh::compute_barycentric_data(const moab::Range& all_tets) {
+  moab::ErrorCode rval;
+  for (auto& tet : all_tets) {
+    moab::Range verts;
+    rval = mbi_->get_connectivity(&tet, 1, verts);
+    if (rval != moab::MB_SUCCESS) {
+      fatal_error("Failed to get connectivity of tet on umesh: " + filename_);
+    }
+
+    moab::CartVect p[4];
+    rval = mbi_->get_coords(verts, p[0].array());
+    if (rval != moab::MB_SUCCESS) {
+      fatal_error("Failed to get coordinates of a tet in umesh: " + filename_);
+    }
+
+    moab::Matrix3 a(p[1] - p[0], p[2] - p[0], p[3] - p[0], true);
+
+    baryc_data_.push_back(a.transpose().inverse());
+  }
+}
+
+// TODO: write this function
+void
+UnstructuredMesh::to_hdf5(hid_t group) const { }
+
+bool
+UnstructuredMesh::point_in_tet(const Position& r, moab::EntityHandle tet) const {
+
+  moab::ErrorCode rval;
+
+  // get tet vertices
+  moab::Range verts;
+  rval = mbi_->get_connectivity(&tet, 1, verts);
+  if (rval != moab::MB_SUCCESS) {
+    warning("Failed to get vertices of tet in umesh: " + filename_);
+    return false;
+  }
+
+  moab::EntityHandle v_zero = verts[0];
+  moab::CartVect p_zero;
+  rval = mbi_->get_coords(&v_zero, 1, p_zero.array());
+  if (rval != moab::MB_SUCCESS) {
+    warning("Failed to get coordinates of a vertex in umesh: " + filename_);
+    return false;
+  }
+
+  moab::CartVect pos(r.x, r.y, r.z);
+
+  // look up barycentric data
+  int idx = get_bin_from_ent_handle(tet);
+  const moab::Matrix3& a_inv = baryc_data_[idx];
+
+  moab::CartVect bary_coords = a_inv * (pos - p_zero);
+
+  bool in_tet = (bary_coords[0] >= 0 && bary_coords[1] >= 0 && bary_coords[2] >= 0 &&
+                 bary_coords[0] + bary_coords[1] + bary_coords[2] <= 1.);
+
+  return in_tet;
+
+}
+
 int
 UnstructuredMesh::get_bin_from_ent_handle(moab::EntityHandle eh) const {
   auto pos = ehs_.find(eh);
@@ -1638,6 +1722,9 @@ void read_meshes(pugi::xml_node root)
       model::meshes.push_back(std::make_unique<RegularMesh>(node));
     } else if (mesh_type == "rectilinear") {
       model::meshes.push_back(std::make_unique<RectilinearMesh>(node));
+    }
+    else if (mesh_type == "unstructured") {
+      model::meshes.push_back(std::make_unique<UnstructuredMesh>(node));
     } else {
       fatal_error("Invalid mesh type: " + mesh_type);
     }