Add check_cell_overlap to C++

CMakeLists.txt

@@ -386,6 +386,7 @@ add_library(libopenmc SHARED
   src/cell.cpp
   src/initialize.cpp
   src/finalize.cpp
+  src/geometry.cpp
   src/geometry_aux.cpp
   src/hdf5_interface.cpp
   src/lattice.cpp

src/cell.cpp

@@ -7,8 +7,10 @@
 
 #include "constants.h"
 #include "error.h"
+#include "geometry.h"
 #include "hdf5_interface.h"
 #include "lattice.h"
+#include "settings.h"
 #include "surface.h"
 #include "xml_interface.h"
 
@@ -460,6 +462,9 @@ read_cells(pugi::xml_node *node)
       global_universes[it->second]->cells.push_back(i);
     }
   }
+
+  // Allocate the cell overlap count if necessary.
+  if (openmc_check_overlaps) overlap_check_count.resize(n_cells, 0);
 }
 
 //==============================================================================

src/geometry.F90

@@ -31,7 +31,14 @@ module geometry
       integer(C_INT32_T), intent(in), value :: search_univ
       integer(C_INT32_T), intent(in), value :: target_univ_id
       integer(C_INT)                        :: count
-    end function
+    end function count_universe_instances
+
+    function check_cell_overlap_c(p) bind(C, name="check_cell_overlap") &
+         result(is_overlapping)
+      import Particle, C_BOOL
+      type(Particle), intent(in) :: p
+      logical(C_BOOL)            :: is_overlapping
+    end function check_cell_overlap_c
   end interface
 
 contains
@@ -60,6 +67,9 @@ contains
     integer :: index_cell              ! index in cells array
     type(Cell),       pointer :: c     ! pointer to cell
     type(Universe),   pointer :: univ  ! universe to search in
+    logical :: overlap_c
+
+    overlap_c = check_cell_overlap_c(p)
 
     ! loop through each coordinate level
     n_coord = p % n_coord
@@ -76,6 +86,7 @@ contains
         if (cell_contains(c, p)) then
           ! the particle should only be contained in one cell per level
           if (index_cell /= p % coord(j) % cell) then
+            if (.not. overlap_c) call fatal_error("Cell overlap disagreement")
             call fatal_error("Overlapping cells detected: " &
                  &// trim(to_str(cells(index_cell) % id())) // ", " &
                  &// trim(to_str(cells(p % coord(j) % cell) % id())) &
@@ -88,6 +99,7 @@ contains
 
       end do
     end do
+    if (overlap_c) call fatal_error("Cell overlap disagreement")
 
   end subroutine check_cell_overlap
 

src/geometry.cpp

@@ -0,0 +1,52 @@
+#include "geometry.h"
+
+#include <sstream>
+
+#include "cell.h"
+#include "error.h"
+#include "particle.h"
+
+//TODO: remove this include
+#include <iostream>
+
+
+namespace openmc {
+
+std::vector<int64_t> overlap_check_count;
+
+extern "C" bool
+check_cell_overlap(Particle* p) {
+  int n_coord = p->n_coord;
+
+  // loop through each coordinate level
+  for (int j = 0; j < n_coord; j++) {
+    //p->n_coord = j + 1;
+    Universe& univ {*global_universes[p->coord[j].universe - 1]};
+    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];
+      Cell& c {*global_cells[index_cell]};
+
+      if (c.contains(p->coord[j].xyz, p->coord[j].uvw, p->surface)) {
+        //TODO: off-by-one indexing
+        if (index_cell != p->coord[j].cell - 1) {
+          std::stringstream err_msg;
+          err_msg << "Overlapping cells detected: " << c.id << ", "
+                  << global_cells[p->coord[j].cell-1]->id << " on universe "
+                  << univ.id;
+          fatal_error(err_msg);
+        }
+        ++overlap_check_count[index_cell];
+      }
+    }
+  }
+
+  return false;
+}
+
+extern "C" int64_t get_overlap_check_count(int i)
+{return overlap_check_count[i];}
+
+} // namespace openmc

src/geometry.h

@@ -1,10 +1,15 @@
 #ifndef OPENMC_GEOMETRY_H
 #define OPENMC_GEOMETRY_H
 
+#include <cstdint>
+#include <vector>
+
 namespace openmc {
 
 extern "C" int openmc_root_universe;
 
+extern std::vector<int64_t> overlap_check_count;
+
 } // namespace openmc
 
 #endif // OPENMC_GEOMETRY_H

src/output.F90

@@ -628,6 +628,14 @@ contains
 
   subroutine print_overlap_check
 
+    interface
+      function get_overlap_check_count(index_cell) result(count) bind(C)
+        import C_INT, C_INT64_T
+        integer(C_INT), intent(in), value :: index_cell
+        integer(C_INT64_T)                :: count
+      end function get_overlap_check_count
+    end interface
+
     integer :: i, j
     integer :: num_sparse = 0
 
@@ -638,6 +646,7 @@ contains
 
     do i = 1, n_cells
       write(ou,101) cells(i) % id(), overlap_check_cnt(i)
+      write(ou,101) cells(i) % id(), get_overlap_check_count(i-1)
       if (overlap_check_cnt(i) < 10) num_sparse = num_sparse + 1
     end do
     write(ou,*)