Address #1061 comments

docs/source/io_formats/geometry.rst

@@ -201,7 +201,7 @@ Each ``<cell>`` element can have the following attributes or sub-elements:
 
     .. math::
 
-       \left [ \begin{array}{ccc} \cos\theta \cos\psi & -\cos\theta \sin\psi +
+       \left [ \begin{array}{ccc} \cos\theta \cos\psi & -\cos\phi \sin\psi +
        \sin\phi \sin\theta \cos\psi & \sin\phi \sin\psi + \cos\phi \sin\theta
        \cos\psi \\ \cos\theta \sin\psi & \cos\phi \cos\psi + \sin\phi \sin\theta
        \sin\psi & -\sin\phi \cos\psi + \cos\phi \sin\theta \sin\psi \\

include/openmc/cell.h

@@ -7,10 +7,10 @@
 #include <unordered_map>
 #include <vector>
 
-#include "constants.h"
 #include "hdf5.h"
 #include "pugixml.hpp"
 
+#include "openmc/constants.h"
 #include "openmc/position.h"
 
 

include/openmc/geometry.h

@@ -4,10 +4,15 @@
 #include <cstdint>
 #include <vector>
 
-#include "particle.h"
+#include "openmc/particle.h"
+
 
 namespace openmc {
 
+//==============================================================================
+// Global variables
+//==============================================================================
+
 extern "C" int openmc_root_universe;
 
 extern std::vector<int64_t> overlap_check_count;
@@ -21,6 +26,15 @@ check_cell_overlap(Particle* p);
 
 //==============================================================================
 //! Locate a particle in the geometry tree and set its geometry data fields.
+//!
+//! \param p A particle to be located.  This function will populate the
+//!   geometry-dependent data fields of the particle.
+//! \param search_surf A surface that the particle is expected to be on.  This
+//!   value should be the signed, 1-based index of a surface.  If positive, the
+//!   cells on the positive half-space of the surface will be searched.  If
+//!   negative, the negative half-space will be searched.
+//! \return True if the particle's location could be found and ascribed to a
+//!   valid geometry coordinate stack.
 //==============================================================================
 
 extern "C" bool

include/openmc/position.h

@@ -68,17 +68,9 @@ inline Position operator*(double a, Position b)   { return b *= a; }
 
 inline bool operator==(Position a, Position b)
 {return a.x == b.x && a.y == b.y && a.z == b.z;}
-inline bool operator==(Position a, double b)
-{return a.x == b && a.y == b && a.z == b;}
-inline bool operator==(double a, Position b)
-{return a == b.x && a == b.y && a == b.z;}
 
 inline bool operator!=(Position a, Position b)
 {return a.x != b.x || a.y != b.y || a.z != b.z;}
-inline bool operator!=(Position a, double b)
-{return a.x != b || a.y != b || a.z != b;}
-inline bool operator!=(double a, Position b)
-{return a != b.x || a != b.y || a != b.z;}
 
 //==============================================================================
 //! Type representing a vector direction in Cartesian coordinates

include/openmc/settings.h

@@ -19,11 +19,9 @@ namespace openmc {
 extern "C" bool openmc_check_overlaps;
 extern "C" bool openmc_particle_restart_run;
 extern "C" bool openmc_restart_run;
-extern "C" bool openmc_trace;
 extern "C" int openmc_verbosity;
 extern "C" bool openmc_write_all_tracks;
 extern "C" double temperature_default;
-#pragma omp threadprivate(openmc_trace)
 
 // Defined in .cpp
 // TODO: Make strings instead of char* once Fortran is gone

include/openmc/simulation.h

@@ -7,7 +7,8 @@ extern "C" int openmc_current_batch;
 extern "C" int openmc_current_gen;
 extern "C" int64_t openmc_current_work;
 extern "C" int openmc_n_lost_particles;
+extern "C" bool openmc_trace;
 
-#pragma omp threadprivate(openmc_current_work)
+#pragma omp threadprivate(openmc_current_work, openmc_trace)
 
 #endif // OPENMC_SIMULATION_H

openmc/cell.py

@@ -57,7 +57,7 @@ class Cell(IDManagerMixin):
 
         .. math::
 
-           \left [ \begin{array}{ccc} \cos\theta \cos\psi & -\cos\theta \sin\psi
+           \left [ \begin{array}{ccc} \cos\theta \cos\psi & -\cos\phi \sin\psi
            + \sin\phi \sin\theta \cos\psi & \sin\phi \sin\psi + \cos\phi
            \sin\theta \cos\psi \\ \cos\theta \sin\psi & \cos\phi \cos\psi +
            \sin\phi \sin\theta \sin\psi & -\sin\phi \cos\psi + \cos\phi

src/cell.cpp

@@ -494,7 +494,7 @@ Cell::to_hdf5(hid_t cells_group) const
   } else if (type == FILL_UNIVERSE) {
     write_dataset(group, "fill_type", "universe");
     write_dataset(group, "fill", global_universes[fill]->id);
-    if (translation != 0) {
+    if (translation != Position(0, 0, 0)) {
       write_dataset(group, "translation", translation);
     }
     if (!rotation.empty()) {

src/geometry.cpp

@@ -7,7 +7,7 @@
 #include "openmc/constants.h"
 #include "openmc/error.h"
 #include "openmc/lattice.h"
-#include "openmc/settings.h"
+#include "openmc/simulation.h"
 #include "openmc/surface.h"
 
 
@@ -23,16 +23,14 @@ extern "C" bool
 check_cell_overlap(Particle* p) {
   int n_coord = p->n_coord;
 
-  // loop through each coordinate level
+  // Loop through each coordinate level
   for (int j = 0; j < n_coord; j++) {
     Universe& univ = *global_universes[p->coord[j].universe];
     int n = univ.cells.size();
 
-    // loop through each cell on this level
-    for (int i = 0; i < n; i++) {
-      int index_cell = univ.cells[i];
+    // Loop through each cell on this level
+    for (auto index_cell : univ.cells) {
       Cell& c = *global_cells[index_cell];
-
       if (c.contains(p->coord[j].xyz, p->coord[j].uvw, p->surface)) {
         if (index_cell != p->coord[j].cell) {
           std::stringstream err_msg;
@@ -65,26 +63,23 @@ find_cell(Particle* p, int search_surf) {
   }
 
   // If a surface was indicated, only search cells from the neighbor list of
-  // that surface.
-  int* search_cells;
-  int n_search_cells;
+  // that surface.  The surface index is signed, and the sign signifies whether
+  // the positive or negative side of the surface should be searched.
+  const std::vector<int>* search_cells;
   if (search_surf > 0) {
-    search_cells = global_surfaces[search_surf-1]->neighbor_pos.data();
-    n_search_cells = global_surfaces[search_surf-1]->neighbor_pos.size();
+    search_cells = &global_surfaces[search_surf-1]->neighbor_pos;
   } else if (search_surf < 0) {
-    search_cells = global_surfaces[-search_surf-1]->neighbor_neg.data();
-    n_search_cells = global_surfaces[-search_surf-1]->neighbor_neg.size();
+    search_cells = &global_surfaces[-search_surf-1]->neighbor_neg;
   } else {
     // No surface was indicated, search all cells in the universe.
-    search_cells = global_universes[i_universe]->cells.data();
-    n_search_cells = global_universes[i_universe]->cells.size();
+    search_cells = &global_universes[i_universe]->cells;
   }
 
   // Find which cell of this universe the particle is in.
   bool found = false;
   int32_t i_cell;
-  for (int i = 0; i < n_search_cells; i++) {
-    i_cell = search_cells[i];
+  for (int i = 0; i < search_cells->size(); i++) {
+    i_cell = (*search_cells)[i];
 
     // Make sure the search cell is in the same universe.
     if (global_cells[i_cell]->universe != i_universe) continue;
@@ -155,6 +150,8 @@ find_cell(Particle* p, int search_surf) {
         p->sqrtkT = c.sqrtkT[0];
       }
 
+      return true;
+
     } else if (c.type == FILL_UNIVERSE) {
       //========================================================================
       //! Found a lower universe, update this coord level then search the next.
@@ -198,7 +195,7 @@ find_cell(Particle* p, int search_surf) {
 
       // Update the coordinate level and recurse.
       ++p->n_coord;
-      find_cell(p, 0);
+      return find_cell(p, 0);
 
     } else if (c.type == FILL_LATTICE) {
       //========================================================================
@@ -245,11 +242,9 @@ find_cell(Particle* p, int search_surf) {
 
       // Update the coordinate level and recurse.
       ++p->n_coord;
-      find_cell(p, 0);
+      return find_cell(p, 0);
     }
   }
-
-  return found;
 }
 
 //==============================================================================

src/output.cpp

@@ -46,29 +46,29 @@ print_overlap_check() {
 #ifdef OPENMC_MPI
   std::vector<int64_t> temp(overlap_check_count);
   int err = MPI_Reduce(temp.data(), overlap_check_count.data(),
-                       overlap_check_count.size(), MPI_LONG, MPI_SUM, 0,
+                       overlap_check_count.size(), MPI_INT64_T, MPI_SUM, 0,
                        mpi::intracomm);
 #endif
 
   if (openmc_master) {
     header("cell overlap check summary", 1);
-    std::cout << " Cell ID      No. Overlap Checks" << std::endl;
+    std::cout << " Cell ID      No. Overlap Checks\n";
 
     std::vector<int32_t> sparse_cell_ids;
     for (int i = 0; i < n_cells; i++) {
       std::cout << " " << std::setw(8) << global_cells[i]->id << std::setw(17)
-                << overlap_check_count[i] << std::endl;;
+                << overlap_check_count[i] << "\n";
       if (overlap_check_count[i] < 10) {
         sparse_cell_ids.push_back(global_cells[i]->id);
       }
     }
 
-    std::cout << std::endl << " There were " << sparse_cell_ids.size()
-              << " cells with less than 10 overlap checks" << std::endl;
+    std::cout << "\n There were " << sparse_cell_ids.size()
+              << " cells with less than 10 overlap checks\n";
     for (auto id : sparse_cell_ids) {
       std::cout << " " << id;
     }
-    std::cout << std::endl;
+    std::cout << "\n";
   }
 }