Some updates to get the CAD build working again.

include/openmc/settings.h

@@ -45,6 +45,7 @@ extern "C" bool ufs_on;                  //!< uniform fission site method on?
 extern "C" bool urr_ptables_on;          //!< use unresolved resonance prob. tables?
 extern "C" bool write_all_tracks;        //!< write track files for every particle?
 extern "C" bool write_initial_source;    //!< write out initial source file?
+extern "C" bool dagmc;                   //!< indicator of DAGMC geometry
 
 // Paths to various files
 extern std::string path_cross_sections;   //!< path to cross_sections.xml
@@ -86,7 +87,6 @@ extern "C" int trigger_batch_interval;   //!< Batch interval for triggers
 extern "C" int verbosity;                //!< How verbose to make output
 extern "C" double weight_cutoff;         //!< Weight cutoff for Russian roulette
 extern "C" double weight_survive;        //!< Survival weight after Russian roulette
-
 } // namespace settings
 
 //! Read settings from XML file

include/openmc/surface.h

@@ -147,6 +147,8 @@ class CADSurface : public Surface
   Direction normal(Position p) const;
   //! Get the bounding box of this surface.
   BoundingBox bounding_box() const;
+
+  void to_hdf5(hid_t group_id) const;
 };
 #endif 
 //==============================================================================

src/cad.cpp

@@ -2,7 +2,7 @@
 #ifdef CAD
 
 #include "cad.h"
-#include "error.h"
+#include "openmc/error.h"
 
 #include <string>
 #include <algorithm>
@@ -59,25 +59,25 @@ void load_cad_geometry_c()
       
       c->fill_ = openmc::C_NONE;
 
-      openmc::global_cells.push_back(c);
+      openmc::cells.push_back(c);
       openmc::cell_map[c->id_] = c->id_;
 
       // Populate the Universe vector and dict
       auto it = openmc::universe_map.find(cad_univ_id);
       if (it == openmc::universe_map.end()) {
-	openmc::global_universes.push_back(new openmc::Universe());
-	openmc::global_universes.back()-> id = cad_univ_id;
-	openmc::global_universes.back()->cells.push_back(i);
-	openmc::universe_map[cad_univ_id] = openmc::global_universes.size() - 1;
+	openmc::universes.push_back(new openmc::Universe());
+	openmc::universes.back()-> id_ = cad_univ_id;
+	openmc::universes.back()->cells_.push_back(i);
+	openmc::universe_map[cad_univ_id] = openmc::universes.size() - 1;
       }
       else {
-	openmc::global_universes[it->second]->cells.push_back(i);
+	openmc::universes[it->second]->cells_.push_back(i);
       }
 
       
       if(DAGMC->is_implicit_complement(vol_handle)) {
 	// assuming implicit complement is void for now
-        c->material.push_back(openmc::MATERIAL_VOID);	
+        c->material_.push_back(openmc::MATERIAL_VOID);	
 	continue;
       }
       
@@ -88,21 +88,21 @@ void load_cad_geometry_c()
 	TOLOWER(mat_value);
 	
 	if(mat_value == "void") {
-	  c->material.push_back(openmc::MATERIAL_VOID);
+	  c->material_.push_back(openmc::MATERIAL_VOID);
 	}
 	else {
-	  c->material.push_back(std::stoi(mat_value));
+	  c->material_.push_back(std::stoi(mat_value));
 	}
       }
       else {
 	std::cout << "Warning: volume without material found!" << std::endl;
-	c->material.push_back(openmc::MATERIAL_VOID);
+	c->material_.push_back(openmc::MATERIAL_VOID);
       }     
     }
 
   // initialize surface objects
   openmc::n_surfaces = DAGMC->num_entities(2);
-  openmc::global_surfaces.resize(openmc::n_surfaces);
+  openmc::surfaces.resize(openmc::n_surfaces);
   
   for(int i = 0; i < openmc::n_surfaces; i++)
     {
@@ -110,7 +110,7 @@ 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+1);
+      s->id_ = DAGMC->id_by_index(2, i+1);
       s->dagmc_ptr = DAGMC;
 
       if(DAGMC->has_prop(surf_handle, "boundary")) {
@@ -120,32 +120,32 @@ void load_cad_geometry_c()
 	TOLOWER(bc_value);
 
 	if(bc_value == "transmit" || bc_value == "transmission") {
-	  s->bc = openmc::BC_TRANSMIT;
+	  s->bc_ = openmc::BC_TRANSMIT;
 	}
 	else if(bc_value == "vacuum") {
-	  s->bc = openmc::BC_VACUUM;
+	  s->bc_ = openmc::BC_VACUUM;
 	}
 	else if(bc_value == "reflective" || bc_value == "reflect" || bc_value == "reflecting") {
-	  s->bc = openmc::BC_REFLECT;
+	  s->bc_ = openmc::BC_REFLECT;
 	}
 	else if(bc_value == "periodic") {
 	  openmc::fatal_error("Periodic boundary condition not supported in CAD.");
 	}
 	else {
 	  std::stringstream err_msg;
-	  err_msg << "Unknown boundary condition \"" << s->bc
-		  << "\" specified on surface " << s->id;
+	  err_msg << "Unknown boundary condition \"" << s->bc_
+		  << "\" specified on surface " << s->id_;
 	  openmc::fatal_error(err_msg);
 	}
       }
       // if no BC property is found, set to transmit
       else {
-	s->bc = openmc::BC_TRANSMIT;
+	s->bc_ = openmc::BC_TRANSMIT;
       }
 
       // add to global array and map
-      openmc::global_surfaces[i] = s;
-      openmc::surface_map[s->id] = s->id;
+      openmc::surfaces[i] = s;
+      openmc::surface_map[s->id_] = s->id_;
     }
 
   return;

src/cad.h

@@ -5,8 +5,8 @@
 #define CAD_H
 
 #include "DagMC.hpp"
-#include "cell.h"
-#include "surface.h"
+#include "openmc/cell.h"
+#include "openmc/surface.h"
 
 extern moab::DagMC* DAGMC;
 

src/cell.cpp

@@ -599,7 +599,7 @@ CADCell::CADCell() : Cell{} {};
 
 std::pair<double, int32_t> CADCell::distance(Position p, Direction u, int32_t on_surface) const {
 
-  moab::EntityHandle vol = dagmc_ptr->entity_by_id(3, id);
+  moab::EntityHandle vol = dagmc_ptr->entity_by_id(3, id_);
   moab::EntityHandle hit_surf;
   double dist;
   dagmc_ptr->ray_fire(vol, p.xyz, u.xyz, hit_surf, dist);
@@ -620,7 +620,7 @@ std::pair<double, int32_t> CADCell::distance(Position p, Direction u, int32_t on
   
 bool CADCell::contains(Position p, Direction u, int32_t on_surface) const {
 
-  moab::EntityHandle vol = dagmc_ptr->entity_by_id(3, id);
+  moab::EntityHandle vol = dagmc_ptr->entity_by_id(3, id_);
 
   int result = 0;
   dagmc_ptr->point_in_volume(vol, p.xyz, result, u.xyz);
@@ -772,8 +772,8 @@ extern "C" {
 #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 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);

src/input_xml.F90

@@ -219,17 +219,6 @@ contains
     ! Get proper XMLNode type given pointer
     root % ptr = root_ptr
 
-    ! Check for use of CAD geometry
-    if (check_for_node(root, "dagmc")) then
-#ifdef CAD
-      call get_node_value_bool(root, "dagmc", dagmc)
-#else
-      if (dagmc) then
-        call fatal_error("CAD mode unsupported for this build of OpenMC")
-      end if
-#endif
-    end if
-
     if (run_mode == MODE_EIGENVALUE) then
       ! Preallocate space for keff and entropy by generation
       call k_generation % reserve(n_max_batches*gen_per_batch)
@@ -597,10 +586,10 @@ contains
     do i = 1, n_cells
       c => cells(i)
       ! additional metadata spoofing
-      allocate(c % sqrtKT(1))
-      c % sqrtkT(1) = -1.0 ! rely on temps from elsewhere
+!      allocate(c % sqrtKT(1))
+!      c % sqrtkT(1) = -1.0 ! rely on temps from elsewhere
       univ_id = c % universe()
-      if (.not. allocated(c%region)) allocate(c%region(0))
+!      if (.not. allocated(c%region)) allocate(c%region(0))
       
       if (.not. cells_in_univ_dict % has(univ_id)) then
          n_universes = n_universes + 1
@@ -615,30 +604,30 @@ contains
 
     ! ==========================================================================
     ! SETUP UNIVERSES
-
+    
     ! Allocate universes, universe cell arrays, and assign base universe
-    allocate(universes(n_universes))
-    do i = 1, n_universes
-      associate (u => universes(i))
-        u % id = univ_ids % data(i)
-        ! Allocate cell list
-        n_cells_in_univ = cells_in_univ_dict % get(u % id)
-        allocate(u % cells(n_cells_in_univ))
-        u % cells(:) = 0
-      end associate
-    end do
+    ! allocate(universes(n_universes))
+    ! do i = 1, n_universes
+    !   associate (u => universes(i))
+    !     u % id = univ_ids % data(i)
+    !     ! Allocate cell list
+    !     n_cells_in_univ = cells_in_univ_dict % get(u % id)
+    !     allocate(u % cells(n_cells_in_univ))
+    !     u % cells(:) = 0
+    !   end associate
+    ! end do
 
-    root_universe = 0 + 1
+    root_universe = find_root_universe()
 
-    do i = 1, n_cells
-      ! Get index in universes array
-      j = universe_dict % get(cells(i) % universe())
-      ! Set the first zero entry in the universe cells array to the index in the
-      ! global cells array
-      associate (u => universes(j))
-        u % cells(find(u % cells, 0)) = i
-      end associate
-    end do
+    ! do i = 1, n_cells
+    !   ! Get index in universes array
+    !   j = universe_dict % get(cells(i) % universe())
+    !   ! Set the first zero entry in the universe cells array to the index in the
+    !   ! global cells array
+    !   associate (u => universes(j))
+    !     u % cells(find(u % cells, 0)) = i
+    !   end associate
+    ! end do
 
     ! Clear dictionary
     call cells_in_univ_dict%clear()

src/settings.F90

@@ -81,7 +81,7 @@ module settings
   integer(C_INT), bind(C) :: run_mode
 
   ! flag for use of CAD geometry
-  logical, bind(C) :: dagmc = .false.
+  logical(C_BOOL), bind(C) :: dagmc
 
   ! Restart run
   logical(C_BOOL), bind(C) :: restart_run

src/settings.cpp

@@ -56,7 +56,8 @@ bool ufs_on                  {false};
 bool urr_ptables_on          {true};
 bool write_all_tracks        {false};
 bool write_initial_source    {false};
-
+bool dagmc                   {false};
+  
 std::string path_cross_sections;
 std::string path_input;
 std::string path_multipole;
@@ -207,6 +208,17 @@ void read_settings_xml()
     verbosity = std::stoi(get_node_value(root, "verbosity"));
   }
 
+  // DAGMC geometry check
+  if (check_for_node(root, "dagmc")) {
+      dagmc = get_node_value_bool(root, "dagmc");
+  }
+
+#ifndef CAD
+  if (dagmc) {
+    fatal_error("CAD mode unsupported for this build of OpenMC");
+  }
+#endif
+  
   // To this point, we haven't displayed any output since we didn't know what
   // the verbosity is. Now that we checked for it, show the title if necessary
   if (openmc_master) {

src/summary.cpp

@@ -2,7 +2,7 @@
 #include "openmc/hdf5_interface.h"
 #include "openmc/lattice.h"
 #include "openmc/surface.h"
-
+#include "openmc/settings.h"
 
 namespace openmc {
 
@@ -10,7 +10,7 @@ extern "C" void
 write_geometry(hid_t file_id) {
 
 #ifdef CAD
-  if (dagmc) {
+  if (settings::dagmc) {
     return;
   }
 #endif

src/surface.cpp

@@ -257,7 +257,7 @@ double CADSurface::evaluate(Position r) const
 }
 
 double CADSurface::distance(Position p, Direction u, bool coincident) const {
-  moab::EntityHandle surf = dagmc_ptr->entity_by_id(2, id);
+  moab::EntityHandle surf = dagmc_ptr->entity_by_id(2, id_);
   moab::EntityHandle hit_surf;
   double dist;
   dagmc_ptr->ray_fire(surf, p.xyz, u.xyz, hit_surf, dist, NULL, 0, 0);
@@ -267,17 +267,19 @@ double CADSurface::distance(Position p, Direction u, bool coincident) const {
 
 Direction CADSurface::normal(Position p) const {
   Direction u;
-  moab::EntityHandle surf = dagmc_ptr->entity_by_id(2, id);
+  moab::EntityHandle surf = dagmc_ptr->entity_by_id(2, id_);
   dagmc_ptr->get_angle(surf, p.xyz, u.xyz);
   return u;
 }
 
 BoundingBox CADSurface::bounding_box() const {
-  moab::EntityHandle surf = dagmc_ptr->entity_by_id(2, id);
+  moab::EntityHandle surf = dagmc_ptr->entity_by_id(2, id_);
   double min[3], max[3];
   dagmc_ptr->getobb(surf, min, max);
   return BoundingBox(min[0], max[0], min[1], max[1], min[2], max[2]);
 }
+
+void CADSurface::to_hdf5(hid_t group_id) const {}
 #endif
 //==============================================================================
 // PeriodicSurface implementation