Condensing filter translation into single commit.

include/openmc/mesh.h

@@ -71,18 +71,28 @@ public:
 
   //! Determine which bins were crossed by a particle
   //
-  //! \param[in] p Particle to check
+  //! \param[in] last_r Previous position of the particle
+  //! \param[in] r Current position of the particle
+  //! \param[in] u Particle direction
   //! \param[out] bins Bins that were crossed
   //! \param[out] lengths Fraction of tracklength in each bin
-  virtual void bins_crossed(const Particle& p, std::vector<int>& bins,
+  virtual void bins_crossed(Position last_r,
+                            Position r,
+                            const Direction& u,
+                            std::vector<int>& bins,
                             std::vector<double>& lengths) const = 0;
 
   //! Determine which surface bins were crossed by a particle
   //
-  //! \param[in] p Particle to check
+  //! \param[in] r0 Previous position of the particle
+  //! \param[in] r1 Current position of the particle
+  //! \param[in] u Particle direction
   //! \param[out] bins Surface bins that were crossed
   virtual void
-  surface_bins_crossed(const Particle& p, std::vector<int>& bins) const = 0;
+  surface_bins_crossed(Position r0,
+                       Position r1,
+                       const Direction& u,
+                       std::vector<int>& bins) const = 0;
 
   //! Get bin at a given position in space
   //
@@ -137,7 +147,10 @@ public:
 
   int n_surface_bins() const override;
 
-  void bins_crossed(const Particle& p, std::vector<int>& bins,
+  void bins_crossed(Position last_r,
+                    Position r,
+                    const Direction& u,
+                    std::vector<int>& bins,
                     std::vector<double>& lengths) const override;
 
   //! Count number of bank sites in each mesh bin / energy bin
@@ -219,9 +232,11 @@ public:
   RegularMesh(pugi::xml_node node);
 
   // Overriden methods
-
-  void surface_bins_crossed(const Particle& p, std::vector<int>& bins)
-  const override;
+  void
+  surface_bins_crossed(Position r0,
+                       Position r1,
+                       const Direction& u,
+                       std::vector<int>& bins) const override;
 
   int get_index_in_direction(double r, int i) const override;
 
@@ -259,9 +274,11 @@ public:
   RectilinearMesh(pugi::xml_node node);
 
   // Overriden methods
-
-  void surface_bins_crossed(const Particle& p, std::vector<int>& bins)
-  const override;
+  void
+  surface_bins_crossed(Position r0,
+                       Position r1,
+                       const Direction& u,
+                       std::vector<int>& bins) const override;
 
   int get_index_in_direction(double r, int i) const override;
 
@@ -325,7 +342,7 @@ public:
   std::string bin_label(int bin) const override;
 
   void surface_bins_crossed(const Particle& p,
-                            std::vector<int>& bins) const override;
+                            std::vector<int>& bins) const;
 
   void to_hdf5(hid_t group) const override;
 
@@ -348,10 +365,18 @@ public:
   MOABMesh(pugi::xml_node);
   MOABMesh(const std::string& filename);
 
-  void bins_crossed(const Particle& p,
+  void bins_crossed(Position last_r,
+                    Position r,
+                    const Direction& u,
                     std::vector<int>& bins,
                     std::vector<double>& lengths) const override;
 
+  void
+  surface_bins_crossed(Position r0,
+                       Position r1,
+                       const Direction& u,
+                       std::vector<int>& bins) const override;
+
   int get_bin(Position r) const;
 
   int n_bins() const override;
@@ -496,10 +521,18 @@ public:
   LibMesh(const std::string& filename);
 
   // Methods
-  void bins_crossed(const Particle& p,
+  void bins_crossed(Position last_r,
+                    Position r,
+                    const Direction& u,
                     std::vector<int>& bins,
                     std::vector<double>& lengths) const override;
 
+  void
+  surface_bins_crossed(Position r0,
+                       Position r1,
+                       const Direction& u,
+                       std::vector<int>& bins) const override;
+
   int get_bin(Position r) const override;
 
   int n_bins() const override;

include/openmc/tallies/filter_mesh.h

@@ -4,6 +4,7 @@
 #include <cstdint>
 
 #include "openmc/tallies/filter.h"
+#include "openmc/position.h"
 
 namespace openmc {
 
@@ -42,11 +43,22 @@ public:
 
   virtual void set_mesh(int32_t mesh);
 
+  virtual void set_translation(const Position& translation);
+
+  virtual void set_translation(const double translation[3]);
+
+  virtual const Position& translation() const {return translation_;}
+
+  virtual bool translated() const {return translated_;}
+
+
 protected:
   //----------------------------------------------------------------------------
   // Data members
 
   int32_t mesh_;
+  bool translated_ {false};
+  Position translation_ {0.0, 0.0, 0.0};
 };
 
 } // namespace openmc

src/mesh.cpp

@@ -581,17 +581,15 @@ bool StructuredMesh::intersects_3d(Position& r0, Position r1, int* ijk) const
   return min_dist < INFTY;
 }
 
-void StructuredMesh::bins_crossed(const Particle& p, std::vector<int>& bins,
-  std::vector<double>& lengths) const
+void StructuredMesh::bins_crossed(Position last_r,
+                                  Position r,
+                                  const Direction& u,
+                                  std::vector<int>& bins,
+                                  std::vector<double>& lengths) const
 {
   // ========================================================================
   // Determine where the track intersects the mesh and if it intersects at all.
 
-  // Copy the starting and ending coordinates of the particle.
-  Position last_r {p.r_last_};
-  Position r {p.r()};
-  Direction u {p.u()};
-
   // Compute the length of the entire track.
   double total_distance = (r - last_r).norm();
   if (total_distance == 0.0) return;
@@ -797,17 +795,12 @@ double RegularMesh::negative_grid_boundary(int* ijk, int i) const
   return lower_left_[i] + (ijk[i] - 1) * width_[i];
 }
 
-void RegularMesh::surface_bins_crossed(const Particle& p,
-                                       std::vector<int>& bins) const
+void
+RegularMesh::surface_bins_crossed(Position r0,
+                                  Position r1,
+                                  const Direction& u,
+                                  std::vector<int>& bins) const
 {
-  // ========================================================================
-  // Determine if the track intersects the tally mesh.
-
-  // Copy the starting and ending coordinates of the particle.
-  Position r0 {p.r_last_current_};
-  Position r1 {p.r()};
-  Direction u {p.u()};
-
   // Determine indices for starting and ending location.
   int n = n_dimension_;
   std::vector<int> ijk0(n), ijk1(n);
@@ -1107,17 +1100,11 @@ int RectilinearMesh::set_grid()
   return 0;
 }
 
-void RectilinearMesh::surface_bins_crossed(const Particle& p,
+void RectilinearMesh::surface_bins_crossed(Position r0,
+                                           Position r1,
+                                           const Direction& u,
                                            std::vector<int>& bins) const
 {
-  // ========================================================================
-  // Determine if the track intersects the tally mesh.
-
-  // Copy the starting and ending coordinates of the particle.
-  Position r0 {p.r_last_current_};
-  Position r1 {p.r()};
-  Direction u {p.u()};
-
   // Determine indices for starting and ending location.
   int ijk0[3], ijk1[3];
   bool start_in_mesh;
@@ -1701,14 +1688,13 @@ MOABMesh::intersect_track(const moab::CartVect& start,
 }
 
 void
-MOABMesh::bins_crossed(const Particle& p,
-                       std::vector<int>& bins,
-                       std::vector<double>& lengths) const
+MOABUnstructuredMesh::bins_crossed(Position last_r,
+                                  Position r,
+                                  const Direction& u,
+                                  std::vector<int>& bins,
+                                  std::vector<double>& lengths) const
 {
-  Position last_r{p.r_last_};
-  Position r{p.r()};
-  Direction u{p.u()};
-  u /= u.norm();
+  // u /= u.norm();
   moab::CartVect r0(last_r.x, last_r.y, last_r.z);
   moab::CartVect r1(r.x, r.y, r.z);
   moab::CartVect dir(u.x, u.y, u.z);
@@ -1833,6 +1819,16 @@ double MOABMesh::tet_volume(moab::EntityHandle tet) const
   return 1.0 / 6.0 * (((p[1] - p[0]) * (p[2] - p[0])) % (p[3] - p[0]));
 }
 
+void MOABUnstructuredMesh::surface_bins_crossed(Position r0,
+                                                Position r1,
+                                                const Direction& u,
+                                                std::vector<int>& bins) const
+{
+
+  // TODO: Implement triangle crossings here
+  throw std::runtime_error{"Unstructured mesh surface tallies are not implemented."};
+}
+
 int MOABMesh::get_bin(Position r) const
 {
   moab::EntityHandle tet = get_tet(r);
@@ -2212,8 +2208,18 @@ int LibMesh::n_bins() const
   return m_->n_elem();
 }
 
-int LibMesh::n_surface_bins() const
+void
+LibMesh::surface_bins_crossed(Position r0,
+                              Position r1,
+                              const Direction& u,
+                              std::vector<int>& bins) const
 {
+  // TODO: Implement triangle crossings here
+  throw std::runtime_error{"Unstructured mesh surface tallies are not implemented."};
+}
+
+int
+LibMesh::n_surface_bins() const {
   int n_bins = 0;
   for (int i = 0; i < this->n_bins(); i++) {
     const libMesh::Elem& e = get_element_from_bin(i);
@@ -2300,7 +2306,9 @@ void LibMesh::write(const std::string& filename) const
 }
 
 void
-LibMesh::bins_crossed(const Particle& p,
+LibMesh::bins_crossed(Position last_r,
+                      Position r,
+                      const Direction& u,
                       std::vector<int>& bins,
                       std::vector<double>& lengths) const
 {

src/tallies/filter_mesh.cpp

@@ -34,14 +34,25 @@ void
 MeshFilter::get_all_bins(const Particle& p, TallyEstimator estimator, FilterMatch& match)
 const
 {
+
+  Position last_r = p.r_last_;
+  Position r = p.r();
+  Position u = p.u();
+
+  // apply translation if present
+  if (translated_) {
+    last_r -= translation_;
+    r -= translation_;
+  }
+
   if (estimator != TallyEstimator::TRACKLENGTH) {
-    auto bin = model::meshes[mesh_]->get_bin(p.r());
+    auto bin = model::meshes[mesh_]->get_bin(r);
     if (bin >= 0) {
       match.bins_.push_back(bin);
       match.weights_.push_back(1.0);
     }
   } else {
-    model::meshes[mesh_]->bins_crossed(p, match.bins_, match.weights_);
+    model::meshes[mesh_]->bins_crossed(last_r, r, u, match.bins_, match.weights_);
   }
 }
 
@@ -66,6 +77,18 @@ MeshFilter::set_mesh(int32_t mesh)
   n_bins_ = model::meshes[mesh_]->n_bins();
 }
 
+void MeshFilter::set_translation(const Position& translation)
+{
+  translated_ = true;
+  translation_ = translation;
+}
+
+void MeshFilter::set_translation(const double translation[3])
+{
+  translated_ = true;
+  translation_ = translation;
+}
+
 //==============================================================================
 // C-API functions
 //==============================================================================
@@ -123,4 +146,27 @@ openmc_mesh_filter_set_mesh(int32_t index, int32_t index_mesh)
   return 0;
 }
 
+extern "C" int
+openmc_mesh_filter_set_translation(int32_t index, double translation[3])
+{
+  // Make sure this is a valid index to an allocated filter
+  if (int err = verify_filter(index)) return err;
+
+  // Get a pointer to the filter and downcast.
+  const auto& filt_base = model::tally_filters[index].get();
+  auto* filt = dynamic_cast<MeshFilter*>(filt_base);
+
+  // Check the filter type
+  if (!filt) {
+    set_errmsg("Tried to set mesh on a non-mesh filter.");
+    return OPENMC_E_INVALID_TYPE;
+  }
+
+  // Set the translation
+  filt->set_translation(translation);
+
+  return 0;
+}
+
+
 } // namespace openmc

src/tallies/filter_meshsurface.cpp

@@ -11,7 +11,10 @@ void
 MeshSurfaceFilter::get_all_bins(const Particle& p, TallyEstimator estimator,
                                 FilterMatch& match) const
 {
-  model::meshes[mesh_]->surface_bins_crossed(p, match.bins_);
+  Position r0 = p.r_last_current_;
+  Position r1 = p.r();
+  Direction u = p.u();
+  model::meshes[mesh_]->surface_bins_crossed(r0, r1, u, match.bins_);
   for (auto b : match.bins_) match.weights_.push_back(1.0);
 }