Some hacking to get a simple geometry running. The simulations agree though 🎉

src/cad.cpp

@@ -23,12 +23,12 @@ void load_cad_geometry_c()
     {
       // set cell ids using global IDs
       openmc::CADCell* c = new openmc::CADCell();
-      c->id_ = DAGMC->id_by_index(3, i);
+      c->id_ = DAGMC->id_by_index(3, i+1);
       c->dagmc_ptr = DAGMC;
       c->universe_ = cad_univ_id; // set to zero for now
       c->material_.push_back(40); // TEMPORARY
       openmc::cells_c.push_back(c);
-      openmc::cell_dict[c->id_] = i;
+      openmc::cell_dict[c->id_] = c->id;
 
       // Populate the Universe vector and dict
       auto it = openmc::universe_dict.find(cad_univ_id);
@@ -51,11 +51,11 @@ void load_cad_geometry_c()
     {
       // set cell ids using global IDs
       openmc::CADSurface* s = new openmc::CADSurface();
-      s->id = DAGMC->id_by_index(2, i);
+      s->id = DAGMC->id_by_index(2, i+1);
       s->dagmc_ptr = DAGMC;
       s->bc = openmc::BC_TRANSMIT;
       openmc::surfaces_c[i] = s;
-      openmc::surface_dict[s->id] = i;
+      openmc::surface_dict[s->id] = s->id;
     }
 
   return;

src/cell.cpp

@@ -598,14 +598,26 @@ CSGCell::contains_complex(Position r, Direction u, int32_t on_surface) const
 CADCell::CADCell() : Cell{} {};
 
 std::pair<double, int32_t> CADCell::distance(const double xyz[3], const double uvw[3], int32_t on_surface) const {
-  std::cout << "Cell Distance" << std::endl;
-  std::pair<double, int> result(1.0, 0);
+
+  moab::EntityHandle vol = dagmc_ptr->entity_by_id(3, id);
+  moab::EntityHandle hit_surf;
+  double dist;
+  dagmc_ptr->ray_fire(vol, xyz, uvw, hit_surf, dist);
+
+  int surf_idx = dagmc_ptr->index_by_handle(hit_surf);
+  
+  std::pair<double, int> result(dist, surf_idx);
+  
   return result;
 }
   
 bool CADCell::contains(const double xyz[3], const double uvw[3], int32_t on_surface) const {
-  std::cout << "Cell Contains" << std::endl;
-  return true;
+  moab::EntityHandle vol = dagmc_ptr->entity_by_id(3, id);
+
+  int result = 0;
+  dagmc_ptr->point_in_volume(vol, xyz, result, uvw);
+  
+  return bool(result);
 }
 
 void CADCell::to_hdf5(hid_t group_id) const { return; }
@@ -749,6 +761,23 @@ extern "C" {
 
   int cell_type(Cell* c) {return c->type_;}
 
+#ifdef CAD
+  int32_t next_cell(CADCell* cur_cell, CADSurface *surf_xed ) {
+    moab::EntityHandle surf = surf_xed->dagmc_ptr->entity_by_id(2,surf_xed->id);
+    moab::EntityHandle vol = cur_cell->dagmc_ptr->entity_by_id(3,cur_cell->id);
+
+    moab::EntityHandle new_vol;
+    cur_cell->dagmc_ptr->next_vol(surf, vol, new_vol);
+
+    return cur_cell->dagmc_ptr->index_by_handle(new_vol);
+  }
+
+  bool is_implicit_complement(CADCell *c) {
+    moab::EntityHandle handle = c->dagmc_ptr->entity_by_id(3,c->id);
+    return c->dagmc_ptr->is_implicit_complement(handle);
+  }
+#endif
+
   int32_t cell_universe(Cell* c) {return c->universe_;}
 
   int32_t cell_fill(Cell* c) {return c->fill_;}

src/geometry.F90

@@ -55,10 +55,38 @@ module geometry
 
     subroutine neighbor_lists() bind(C)
     end subroutine neighbor_lists
-  end interface
+
+    function next_cell_c(current_cell, surface_crossed) &
+      bind(C, name="next_cell") result(new_cell)
+      import C_PTR, C_INT32_T
+      type(C_PTR), intent(in), value :: current_cell
+      type(C_PTR), intent(in), value :: surface_crossed
+      integer(C_INT32_T)             :: new_cell
+    end function next_cell_c
+
+    function is_implicit_complement_C(cell) &
+         bind(C, name="is_implicit_complement") result(res)
+      import C_PTR, C_BOOL
+      type(C_PTR), intent(in), value :: cell
+      logical(C_BOOL)                :: res
+    end function is_implicit_complement_C
+ end interface
 
 contains
 
+  function next_cell(c, s) result(new_cell)
+    type(Cell), intent(in) :: c
+    type(Surface), intent(in) :: s
+    integer :: new_cell
+    new_cell = next_cell_c(c%ptr, s%ptr)
+  end function next_cell
+
+  function is_implicit_complement(c) result(res)
+    type(Cell), intent(in) :: c
+    logical:: res
+    res = is_implicit_complement_c(c%ptr)
+  end function is_implicit_complement
+  
 !===============================================================================
 ! FIND_CELL determines what cell a source particle is in within a particular
 ! universe. If the base universe is passed, the particle should be found as long

src/tracking.F90

@@ -6,7 +6,7 @@ module tracking
   use error,              only: warning, write_message
   use geometry_header,    only: cells
   use geometry,           only: find_cell, distance_to_boundary, cross_lattice,&
-                                check_cell_overlap
+                                check_cell_overlap, next_cell, is_implicit_complement
   use material_header,    only: materials, Material
   use message_passing
   use mgxs_interface
@@ -309,11 +309,34 @@ contains
     real(8) :: norm       ! "norm" of surface normal
     real(8) :: xyz(3)     ! Saved global coordinate
     integer :: i_surface  ! index in surfaces
+    integer :: i_cell     ! index of new cell
     logical :: rotational ! if rotational periodic BC applied
     logical :: found      ! particle found in universe?
     class(Surface), pointer :: surf
     class(Surface), pointer :: surf2 ! periodic partner surface
+    
+#ifdef CAD
+    if (dagmc) then
+       i_cell = next_cell(cells(p % last_cell(1)), surfaces(ABS(p % surface)))
+       ! save material and temp
+       p % last_material = p % material
+       p % last_sqrtkT = p % sqrtKT
+       ! set new cell value
+       p % coord(1) % cell = i_cell
+       p % cell_instance = 1
+       p % material = cells(i_cell) % material(1)
+       p % sqrtKT = cells(i_cell) % sqrtKT(1)
 
+       if (is_implicit_complement(cells(i_cell))) then
+          p % alive = .false.
+       end if
+       
+       return
+       
+    end if
+#endif
+    
+    
     i_surface = abs(p % surface)
     surf => surfaces(i_surface)
     if (verbosity >= 10 .or. trace) then