Merge pull request #1075 from pshriwise/dagmc

Initial CAD-based geometry implementation

.travis.yml

@@ -13,6 +13,8 @@ addons:
       - libmpich-dev
       - libhdf5-serial-dev
       - libhdf5-mpich-dev
+      - libblas-dev
+      - liblapack-dev
 cache:
   directories:
   - $HOME/nndc_hdf5
@@ -27,6 +29,7 @@ env:
     - OPENMC_CROSS_SECTIONS=$HOME/nndc_hdf5/cross_sections.xml
     - OPENMC_ENDF_DATA=$HOME/endf-b-vii.1
     - OPENMC_MULTIPOLE_LIBRARY=$HOME/WMP_Library
+    - LD_LIBRARY_PATH=$HOME/MOAB/lib:$HOME/DAGMC/lib
     - PATH=$PATH:$HOME/NJOY2016/build
     - DISPLAY=:99.0
     - COVERALLS_PARALLEL=true
@@ -35,6 +38,8 @@ env:
     - OMP=y MPI=n PHDF5=n
     - OMP=n MPI=y PHDF5=n
     - OMP=n MPI=y PHDF5=y
+    - OMP=n MPI=y PHDF5=y DAGMC=y
+    - OMP=y MPI=y PHDF5=y DAGMC=y        
 notifications:
   webhooks: https://coveralls.io/webhook?repo_token=$COVERALLS_REPO_TOKEN
 install:

CMakeLists.txt

@@ -19,6 +19,9 @@ option(debug    "Compile with debug flags"                       OFF)
 option(optimize "Turn on all compiler optimization flags"        OFF)
 option(coverage "Compile with coverage analysis flags"           OFF)
 option(mpif08   "Use Fortran 2008 MPI interface"                 OFF)
+option(dagmc    "Enable support for DAGMC (CAD) geometry"        OFF)
+
+# Maximum number of nested coordinates levels
 set(maxcoord 10 CACHE STRING "Maximum number of nested coordinate levels")
 
 #===============================================================================
@@ -36,6 +39,17 @@ if(MPI_ENABLED AND mpif08)
   message("-- Using Fortran 2008 MPI bindings")
 endif()
 
+#===============================================================================
+# DAGMC Geometry Support - need DAGMC/MOAB
+#===============================================================================
+if(dagmc)
+  find_package(DAGMC REQUIRED)
+  if(NOT DAGMC_FOUND)
+    message(FATAL_ERROR "Could not find DAGMC installation")
+  endif()
+  link_directories(${DAGMC_LIBRARY_DIRS})
+endif()
+
 #===============================================================================
 # HDF5 for binary output
 #===============================================================================
@@ -295,6 +309,7 @@ add_library(libopenmc SHARED
   src/algorithm.F90
   src/bank_header.F90
   src/api.F90
+  src/dagmc_header.F90
   src/cmfd_data.F90
   src/cmfd_execute.F90
   src/cmfd_header.F90
@@ -383,6 +398,7 @@ add_library(libopenmc SHARED
   src/tallies/tally_header.F90
   src/tallies/trigger.F90
   src/tallies/trigger_header.F90
+  src/dagmc.cpp
   src/cell.cpp
   src/cmfd_execute.cpp
   src/distribution.cpp
@@ -428,6 +444,7 @@ add_library(libopenmc SHARED
   src/thermal.cpp
   src/xml_interface.cpp
   src/xsdata.cpp)
+
 set_target_properties(libopenmc PROPERTIES
   OUTPUT_NAME openmc
   LINKER_LANGUAGE Fortran)
@@ -473,10 +490,15 @@ endif()
 target_link_libraries(libopenmc ${ldflags} ${HDF5_LIBRARIES} pugixml
                       faddeeva xtensor)
 
+if(dagmc)
+  target_compile_definitions(libopenmc PRIVATE DAGMC)
+  target_link_libraries(libopenmc ${DAGMC_LIBRARIES})
+  target_include_directories(libopenmc PRIVATE ${DAGMC_INCLUDE_DIRS})
+endif()
+
 #===============================================================================
 # openmc executable
 #===============================================================================
-
 add_executable(openmc src/main.cpp)
 target_compile_options(openmc PRIVATE ${cxxflags})
 target_link_libraries(openmc libopenmc)

cmake/Modules/FindDAGMC.cmake

@@ -0,0 +1,18 @@
+# Try to find DAGMC
+#
+# Once done this will define
+#
+#  DAGMC_FOUND - system has DAGMC
+#  DAGMC_INCLUDE_DIRS - the DAGMC include directory
+#  DAGMC_LIBRARIES - Link these to use DAGMC
+#  DAGMC_DEFINITIONS - Compiler switches required for using DAGMC
+
+find_path(DAGMC_CMAKE_CONFIG NAMES DAGMCConfig.cmake
+          HINTS ${DAGMC_ROOT} $ENV{DAGMC_ROOT}
+          PATHS ENV LD_LIBRARY_PATH
+          PATH_SUFFIXES lib Lib cmake lib/cmake
+          NO_DEFAULT_PATH)
+
+message(STATUS "Found DAGMC in ${DAGMC_CMAKE_CONFIG}")
+
+include(${DAGMC_CMAKE_CONFIG}/DAGMCConfig.cmake)

docs/source/io_formats/summary.rst

@@ -24,6 +24,8 @@ The current version of the summary file format is 6.0.
              - **n_universes** (*int*) -- Number of unique universes in the
                problem.
              - **n_lattices** (*int*) -- Number of lattices in the problem.
+             - **dagmc** (*int*) -- Indicates that a DAGMC geometry was used
+               if present.
 
 **/geometry/cells/cell <uid>/**
 

include/openmc/cell.h

@@ -13,6 +13,9 @@
 #include "openmc/constants.h"
 #include "openmc/position.h"
 
+#ifdef DAGMC
+#include "DagMC.hpp"
+#endif
 
 namespace openmc {
 
@@ -107,9 +110,8 @@ public:
 
   std::vector<int32_t> offset_;  //!< Distribcell offset table
 
-  Cell() {};
-
   explicit Cell(pugi::xml_node cell_node);
+  Cell() {};
 
   //! \brief Determine if a cell contains the particle at a given location.
   //!
@@ -130,21 +132,57 @@ public:
   //! \param on_surface The signed index of a surface that the coordinate is
   //!   known to be on.  This index takes precedence over surface sense
   //!   calculations.
+  virtual bool
+  contains(Position r, Direction u, int32_t on_surface) const = 0;
+
+  //! Find the oncoming boundary of this cell.
+  virtual std::pair<double, int32_t>
+  distance(Position r, Direction u, int32_t on_surface) const = 0;
+
+  //! Write all information needed to reconstruct the cell to an HDF5 group.
+  //! @param group_id An HDF5 group id.
+  virtual void to_hdf5(hid_t group_id) const = 0;
+
+  virtual ~Cell() {}
+};
+ 
+class CSGCell : public Cell
+{
+public:
+
+  CSGCell();
+  
+  explicit CSGCell(pugi::xml_node cell_node);
+
   bool
   contains(Position r, Direction u, int32_t on_surface) const;
 
-  //! Find the oncoming boundary of this cell.
   std::pair<double, int32_t>
   distance(Position r, Direction u, int32_t on_surface) const;
 
-  //! \brief Write cell information to an HDF5 group.
-  //! \param group_id An HDF5 group id.
   void to_hdf5(hid_t group_id) const;
 
+
+  
 protected:
   bool contains_simple(Position r, Direction u, int32_t on_surface) const;
   bool contains_complex(Position r, Direction u, int32_t on_surface) const;
 };
 
+#ifdef DAGMC
+class DAGCell : public Cell
+{
+public:
+  moab::DagMC* dagmc_ptr_;
+  DAGCell();
+
+  std::pair<double, int32_t> distance(Position r, Direction u, int32_t on_surface) const;
+  bool contains(Position r, Direction u, int32_t on_surface) const;
+
+  void to_hdf5(hid_t group_id) const;
+
+};
+#endif
+ 
 } // namespace openmc
 #endif // OPENMC_CELL_H

include/openmc/dagmc.h

@@ -0,0 +1,28 @@
+
+#ifndef OPENMC_DAGMC_H
+#define OPENMC_DAGMC_H
+
+#ifdef DAGMC
+
+#include "DagMC.hpp"
+#include "openmc/cell.h"
+#include "openmc/surface.h"
+
+namespace openmc {
+
+extern moab::DagMC* DAG;
+
+extern "C" void load_dagmc_geometry();
+extern "C" void free_memory_dagmc();
+
+}
+
+#endif // DAGMC
+
+#endif // OPENMC_DAGMC_H
+
+#ifdef DAGMC
+extern "C" constexpr bool dagmc_enabled = true;
+#else
+extern "C" constexpr bool dagmc_enabled = false;
+#endif

include/openmc/position.h

@@ -25,6 +25,7 @@ struct Position {
   Position& operator-=(double);
   Position& operator*=(Position);
   Position& operator*=(double);
+
   const double& operator[](int i) const {
     switch (i) {
       case 0: return x;

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

@@ -12,6 +12,9 @@
 #include "openmc/constants.h"
 #include "openmc/position.h"
 
+#ifdef DAGMC
+#include "DagMC.hpp"
+#endif
 
 namespace openmc {
 
@@ -65,6 +68,7 @@ public:
   std::vector<int> neighbor_neg_; //!< List of cells on negative side
 
   explicit Surface(pugi::xml_node surf_node);
+  Surface();
 
   virtual ~Surface() {}
 
@@ -104,12 +108,40 @@ public:
   //! Write all information needed to reconstruct the surface to an HDF5 group.
   //! \param group_id An HDF5 group id.
   //TODO: this probably needs to include i_periodic for PeriodicSurface
+  virtual void to_hdf5(hid_t group_id) const = 0;
+
+};
+
+class CSGSurface : public Surface
+{
+public:
+  explicit CSGSurface(pugi::xml_node surf_node);
+  CSGSurface();
+
   void to_hdf5(hid_t group_id) const;
 
 protected:
   virtual void to_hdf5_inner(hid_t group_id) const = 0;
 };
 
+//==============================================================================
+//! A `Surface` representing a DAGMC-based surface in DAGMC.
+//==============================================================================
+#ifdef DAGMC
+class DAGSurface : public Surface
+{
+public:
+  moab::DagMC* dagmc_ptr_;
+  DAGSurface();
+  double evaluate(Position r) const;
+  double distance(Position r, Direction u, bool coincident) const;
+  Direction normal(Position r) const;
+  //! Get the bounding box of this surface.
+  BoundingBox bounding_box() const;
+
+  void to_hdf5(hid_t group_id) const;
+};
+#endif
 //==============================================================================
 //! A `Surface` that supports periodic boundary conditions.
 //!
@@ -118,7 +150,7 @@ protected:
 //! `XPlane`-`YPlane` pairs.
 //==============================================================================
 
-class PeriodicSurface : public Surface
+class PeriodicSurface : public CSGSurface
 {
 public:
   int i_periodic_{C_NONE};    //!< Index of corresponding periodic surface
@@ -227,7 +259,7 @@ public:
 //! \f$(y - y_0)^2 + (z - z_0)^2 - R^2 = 0\f$
 //==============================================================================
 
-class SurfaceXCylinder : public Surface
+class SurfaceXCylinder : public CSGSurface
 {
   double y0_, z0_, radius_;
 public:
@@ -245,7 +277,7 @@ public:
 //! \f$(x - x_0)^2 + (z - z_0)^2 - R^2 = 0\f$
 //==============================================================================
 
-class SurfaceYCylinder : public Surface
+class SurfaceYCylinder : public CSGSurface
 {
   double x0_, z0_, radius_;
 public:
@@ -263,7 +295,7 @@ public:
 //! \f$(x - x_0)^2 + (y - y_0)^2 - R^2 = 0\f$
 //==============================================================================
 
-class SurfaceZCylinder : public Surface
+class SurfaceZCylinder : public CSGSurface
 {
   double x0_, y0_, radius_;
 public:
@@ -281,7 +313,7 @@ public:
 //! \f$(x - x_0)^2 + (y - y_0)^2 + (z - z_0)^2 - R^2 = 0\f$
 //==============================================================================
 
-class SurfaceSphere : public Surface
+class SurfaceSphere : public CSGSurface
 {
   double x0_, y0_, z0_, radius_;
 public:
@@ -299,7 +331,7 @@ public:
 //! \f$(y - y_0)^2 + (z - z_0)^2 - R^2 (x - x_0)^2 = 0\f$
 //==============================================================================
 
-class SurfaceXCone : public Surface
+class SurfaceXCone : public CSGSurface
 {
   double x0_, y0_, z0_, radius_sq_;
 public:
@@ -317,7 +349,7 @@ public:
 //! \f$(x - x_0)^2 + (z - z_0)^2 - R^2 (y - y_0)^2 = 0\f$
 //==============================================================================
 
-class SurfaceYCone : public Surface
+class SurfaceYCone : public CSGSurface
 {
   double x0_, y0_, z0_, radius_sq_;
 public:
@@ -335,7 +367,7 @@ public:
 //! \f$(x - x_0)^2 + (y - y_0)^2 - R^2 (z - z_0)^2 = 0\f$
 //==============================================================================
 
-class SurfaceZCone : public Surface
+class SurfaceZCone : public CSGSurface
 {
   double x0_, y0_, z0_, radius_sq_;
 public:
@@ -352,7 +384,7 @@ public:
 //! \f$A x^2 + B y^2 + C z^2 + D x y + E y z + F x z + G x + H y + J z + K = 0\f$
 //==============================================================================
 
-class SurfaceQuadric : public Surface
+class SurfaceQuadric : public CSGSurface
 {
   // Ax^2 + By^2 + Cz^2 + Dxy + Eyz + Fxz + Gx + Hy + Jz + K = 0
   double A_, B_, C_, D_, E_, F_, G_, H_, J_, K_;

openmc/capi/__init__.py

@@ -12,7 +12,7 @@ objects in the :mod:`openmc.capi` subpackage, for example:
 
 """
 
-from ctypes import CDLL
+from ctypes import CDLL, c_bool
 import os
 import sys
 
@@ -38,6 +38,8 @@ else:
     from unittest.mock import Mock
     _dll = Mock()
 
+dagmc_enabled = bool(c_bool.in_dll(_dll, "dagmc_enabled"))
+    
 from .error import *
 from .core import *
 from .nuclide import *

openmc/geometry.py

@@ -266,7 +266,10 @@ class Geometry(object):
             Dictionary mapping cell IDs to :class:`openmc.Cell` instances
 
         """
-        return self.root_universe.get_all_cells()
+        if self.root_universe is not None:
+            return self.root_universe.get_all_cells()
+        else:
+            return []
 
     def get_all_universes(self):
         """Return all universes in the geometry.

openmc/settings.py

@@ -37,6 +37,8 @@ class Settings(object):
         weight assigned to particles that are not killed after Russian
         roulette. Value of energy should be a float indicating energy in eV
         below which particle type will be killed.
+    dagmc : bool
+        Indicate that a CAD-based DAGMC geometry will be used.
     electron_treatment : {'led', 'ttb'}
         Whether to deposit all energy from electrons locally ('led') or create
         secondary bremsstrahlung photons ('ttb').
@@ -224,6 +226,8 @@ class Settings(object):
         self._create_fission_neutrons = None
         self._log_grid_bins = None
 
+        self._dagmc = None
+        
     @property
     def run_mode(self):
         return self._run_mode
@@ -368,6 +372,10 @@ class Settings(object):
     def log_grid_bins(self):
         return self._log_grid_bins
 
+    @property
+    def dagmc(self):
+        return self._dagmc
+    
     @run_mode.setter
     def run_mode(self, run_mode):
         cv.check_value('run mode', run_mode, _RUN_MODES)
@@ -511,6 +519,11 @@ class Settings(object):
         cv.check_type('photon transport', photon_transport, bool)
         self._photon_transport = photon_transport
 
+    @dagmc.setter
+    def dagmc(self, dagmc):
+        cv.check_type('dagmc geometry', dagmc, bool)
+        self._dagmc = dagmc
+        
     @ptables.setter
     def ptables(self, ptables):
         cv.check_type('probability tables', ptables, bool)
@@ -945,6 +958,11 @@ class Settings(object):
             elem = ET.SubElement(root, "log_grid_bins")
             elem.text = str(self._log_grid_bins)
 
+    def _create_dagmc_subelement(self, root):
+        if self._dagmc is not None:
+            elem = ET.SubElement(root, "dagmc")
+            element.text = str(self._dagmc).lower()
+        
     def export_to_xml(self, path='settings.xml'):
         """Export simulation settings to an XML file.
 
@@ -992,7 +1010,8 @@ class Settings(object):
         self._create_volume_calcs_subelement(root_element)
         self._create_create_fission_neutrons_subelement(root_element)
         self._create_log_grid_bins_subelement(root_element)
-
+        self._create_dagmc_subelement(root_element)
+        
         # Clean the indentation in the file to be user-readable
         clean_indentation(root_element)
 

openmc/summary.py

@@ -97,6 +97,9 @@ class Summary(object):
                 self._macroscopics = name.decode()
 
     def _read_geometry(self):
+        if "dagmc" in self._f['geometry'].attrs.keys():
+            return
+
         # Read in and initialize the Materials and Geometry
         self._read_materials()
         self._read_surfaces()

src/api.F90

@@ -29,6 +29,10 @@ module openmc_api
   use timer_header
   use volume_calc,     only: openmc_calculate_volumes
 
+#ifdef DAGMC
+  use dagmc_header,      only: free_memory_dagmc
+#endif
+
   implicit none
 
   private
@@ -149,6 +153,7 @@ contains
     root_universe = -1
     run_CE = .true.
     run_mode = -1
+    dagmc = .false.
     satisfy_triggers = .false.
     call openmc_set_seed(DEFAULT_SEED)
     source_latest = .false.
@@ -325,6 +330,9 @@ contains
     call free_memory_tally_filter()
     call free_memory_tally_derivative()
     call free_memory_bank()
+#ifdef DAGMC
+    call free_memory_dagmc()
+#endif
 
     ! Deallocate CMFD
     call deallocate_cmfd(cmfd)

src/cell.cpp

@@ -15,7 +15,6 @@
 #include "openmc/surface.h"
 #include "openmc/xml_interface.h"
 
-
 namespace openmc {
 
 //==============================================================================
@@ -209,7 +208,9 @@ Universe::to_hdf5(hid_t universes_group) const
 // Cell implementation
 //==============================================================================
 
-Cell::Cell(pugi::xml_node cell_node)
+CSGCell::CSGCell() {} // empty constructor
+  
+CSGCell::CSGCell(pugi::xml_node cell_node)
 {
   if (check_for_node(cell_node, "id")) {
     id_ = std::stoi(get_node_value(cell_node, "id"));
@@ -393,7 +394,7 @@ Cell::Cell(pugi::xml_node cell_node)
 //==============================================================================
 
 bool
-Cell::contains(Position r, Direction u, int32_t on_surface) const
+CSGCell::contains(Position r, Direction u, int32_t on_surface) const
 {
   if (simple_) {
     return contains_simple(r, u, on_surface);
@@ -405,7 +406,7 @@ Cell::contains(Position r, Direction u, int32_t on_surface) const
 //==============================================================================
 
 std::pair<double, int32_t>
-Cell::distance(Position r, Direction u, int32_t on_surface) const
+CSGCell::distance(Position r, Direction u, int32_t on_surface) const
 {
   double min_dist {INFTY};
   int32_t i_surf {std::numeric_limits<int32_t>::max()};
@@ -434,12 +435,12 @@ Cell::distance(Position r, Direction u, int32_t on_surface) const
 //==============================================================================
 
 void
-Cell::to_hdf5(hid_t cells_group) const
+CSGCell::to_hdf5(hid_t cell_group) const
 {
   // Create a group for this cell.
   std::stringstream group_name;
   group_name << "cell " << id_;
-  auto group = create_group(cells_group, group_name);
+  auto group = create_group(cell_group, group_name);
 
   if (!name_.empty()) {
     write_string(group, "name", name_, false);
@@ -513,7 +514,7 @@ Cell::to_hdf5(hid_t cells_group) const
 //==============================================================================
 
 bool
-Cell::contains_simple(Position r, Direction u, int32_t on_surface) const
+CSGCell::contains_simple(Position r, Direction u, int32_t on_surface) const
 {
   for (int32_t token : rpn_) {
     if (token < OP_UNION) {
@@ -537,7 +538,7 @@ Cell::contains_simple(Position r, Direction u, int32_t on_surface) const
 //==============================================================================
 
 bool
-Cell::contains_complex(Position r, Direction u, int32_t on_surface) const
+CSGCell::contains_complex(Position r, Direction u, int32_t on_surface) const
 {
   // Make a stack of booleans.  We don't know how big it needs to be, but we do
   // know that rpn.size() is an upper-bound.
@@ -585,6 +586,50 @@ Cell::contains_complex(Position r, Direction u, int32_t on_surface) const
   }
 }
 
+//==============================================================================
+// DAGMC Cell implementation
+//==============================================================================
+#ifdef DAGMC
+DAGCell::DAGCell() : Cell{} {};
+
+std::pair<double, int32_t>
+DAGCell::distance(Position r, Direction u, int32_t on_surface) const
+{
+  moab::ErrorCode rval;
+  moab::EntityHandle vol = dagmc_ptr_->entity_by_id(3, id_);
+  moab::EntityHandle hit_surf;
+  double dist;
+  double pnt[3] = {r.x, r.y, r.z};
+  double dir[3] = {u.x, u.y, u.z};
+  rval = dagmc_ptr_->ray_fire(vol, pnt, dir, hit_surf, dist);
+  MB_CHK_ERR_CONT(rval);
+  int surf_idx;
+  if (hit_surf != 0) {
+    surf_idx = dagmc_ptr_->index_by_handle(hit_surf);
+  } else {  // indicate that particle is lost
+    surf_idx = -1;
+  }
+  
+  return {dist, surf_idx};
+}
+  
+bool DAGCell::contains(Position r, Direction u, int32_t on_surface) const
+{
+  moab::ErrorCode rval;
+  moab::EntityHandle vol = dagmc_ptr_->entity_by_id(3, id_);
+
+  int result = 0;
+  double pnt[3] = {r.x, r.y, r.z};
+  double dir[3] = {u.x, u.y, u.z};
+  rval = dagmc_ptr_->point_in_volume(vol, pnt, result, dir);
+  MB_CHK_ERR_CONT(rval);  
+  return result;
+}
+
+void DAGCell::to_hdf5(hid_t group_id) const { return; }
+
+#endif
+
 //==============================================================================
 // Non-method functions
 //==============================================================================
@@ -601,7 +646,7 @@ read_cells(pugi::xml_node* node)
   // Loop over XML cell elements and populate the array.
   cells.reserve(n_cells);
   for (pugi::xml_node cell_node: node->children("cell")) {
-    cells.push_back(new Cell(cell_node));
+    cells.push_back(new CSGCell(cell_node));
   }
 
   // Populate the Universe vector and map.
@@ -725,6 +770,21 @@ extern "C" {
 
   int cell_type(Cell* c) {return c->type_;}
 
+#ifdef DAGMC
+
+  int32_t next_cell(DAGCell* cur_cell, DAGSurface* 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);
+  }
+
+#endif
+
   int32_t cell_universe(Cell* c) {return c->universe_;}
 
   int32_t cell_fill(Cell* c) {return c->fill_;}
@@ -753,7 +813,7 @@ extern "C" {
   {
     cells.reserve(cells.size() + n);
     for (int32_t i = 0; i < n; i++) {
-      cells.push_back(new Cell());
+      cells.push_back(new CSGCell());
     }
     n_cells = cells.size();
   }

src/dagmc.cpp

@@ -0,0 +1,144 @@
+
+#include "openmc/dagmc.h"
+#include "openmc/error.h"
+#include "openmc/string_functions.h"
+#include "openmc/settings.h"
+#include "openmc/geometry.h"
+
+#include <string>
+#include <sstream>
+#include <algorithm>
+
+#ifdef DAGMC
+
+namespace openmc {
+
+moab::DagMC* DAG;
+
+void load_dagmc_geometry()
+{
+  if (!DAG) {
+    DAG = new moab::DagMC();
+  }
+
+  int32_t dagmc_univ_id = 0; // universe is always 0 for DAGMC
+
+  moab::ErrorCode rval = DAG->load_file("dagmc.h5m");
+  MB_CHK_ERR_CONT(rval);
+
+  rval = DAG->init_OBBTree();
+  MB_CHK_ERR_CONT(rval);
+
+  std::vector<std::string> prop_keywords;
+  prop_keywords.push_back("mat");
+  prop_keywords.push_back("boundary");
+
+  std::map<std::string, std::string> ph;
+  DAG->parse_properties(prop_keywords, ph, ":");
+  MB_CHK_ERR_CONT(rval);
+
+  // initialize cell objects
+  n_cells = DAG->num_entities(3);
+
+  // Allocate the cell overlap count if necessary.
+  if (settings::check_overlaps) overlap_check_count.resize(n_cells, 0);
+
+  for (int i = 0; i < n_cells; i++) {
+    moab::EntityHandle vol_handle = DAG->entity_by_index(3, i+1);
+
+    // set cell ids using global IDs
+    DAGCell* c = new DAGCell();
+    c->id_ = DAG->id_by_index(3, i+1);
+    c->dagmc_ptr_ = DAG;
+    c->universe_ = dagmc_univ_id; // set to zero for now
+    c->fill_ = C_NONE; // no fill, single universe
+
+    cells.push_back(c);
+    cell_map[c->id_] = c->id_;
+
+    // Populate the Universe vector and dict
+    auto it = universe_map.find(dagmc_univ_id);
+    if (it == universe_map.end()) {
+      universes.push_back(new Universe());
+      universes.back()-> id_ = dagmc_univ_id;
+      universes.back()->cells_.push_back(i);
+      universe_map[dagmc_univ_id] = universes.size() - 1;
+    } else {
+      universes[it->second]->cells_.push_back(i);
+    }
+
+    if (DAG->is_implicit_complement(vol_handle)) {
+      // assuming implicit complement is void for now
+      c->material_.push_back(MATERIAL_VOID);
+      continue;
+    }
+
+    if (DAG->has_prop(vol_handle, "mat")){
+      std::string mat_value;
+      rval = DAG->prop_value(vol_handle, "mat", mat_value);
+      MB_CHK_ERR_CONT(rval);
+      to_lower(mat_value);
+
+      if (mat_value == "void" || mat_value == "vacuum") {
+        c->material_.push_back(MATERIAL_VOID);
+      } else {
+        c->material_.push_back(std::stoi(mat_value));
+      }
+    } else {
+      std::stringstream err_msg;
+      err_msg << "Volume " << c->id_ << " has no material assignment.";
+      fatal_error(err_msg.str());
+    }
+  }
+
+  // initialize surface objects
+  n_surfaces = DAG->num_entities(2);
+  surfaces.resize(n_surfaces);
+
+  for (int i = 0; i < n_surfaces; i++) {
+    moab::EntityHandle surf_handle = DAG->entity_by_index(2, i+1);
+
+    // set cell ids using global IDs
+    DAGSurface* s = new DAGSurface();
+    s->id_ = DAG->id_by_index(2, i+1);
+    s->dagmc_ptr_ = DAG;
+
+    if (DAG->has_prop(surf_handle, "boundary")) {
+      std::string bc_value;
+      rval = DAG->prop_value(surf_handle, "boundary", bc_value);
+      MB_CHK_ERR_CONT(rval);
+      to_lower(bc_value);
+
+      if (bc_value == "transmit" || bc_value == "transmission") {
+        s->bc_ = BC_TRANSMIT;
+      } else if (bc_value == "vacuum") {
+        s->bc_ = BC_VACUUM;
+      } else if (bc_value == "reflective" || bc_value == "reflect" || bc_value == "reflecting") {
+        s->bc_ = BC_REFLECT;
+      } else if (bc_value == "periodic") {
+        fatal_error("Periodic boundary condition not supported in DAGMC.");
+      } else {
+        std::stringstream err_msg;
+        err_msg << "Unknown boundary condition \"" << s->bc_
+                << "\" specified on surface " << s->id_;
+        fatal_error(err_msg);
+      }
+    } else {   // if no BC property is found, set to transmit
+      s->bc_ = BC_TRANSMIT;
+    }
+
+    // add to global array and map
+    surfaces[i] = s;
+    surface_map[s->id_] = s->id_;
+  }
+
+  return;
+}
+
+void free_memory_dagmc()
+{
+  delete DAG;
+}
+
+}
+#endif

src/dagmc_header.F90

@@ -0,0 +1,20 @@
+#ifdef DAGMC
+
+module dagmc_header
+
+  use, intrinsic :: ISO_C_BINDING
+
+  implicit none
+
+  interface
+    subroutine load_dagmc_geometry() bind(C)
+    end subroutine load_dagmc_geometry
+
+    subroutine free_memory_dagmc() bind(C)
+    end subroutine free_memory_dagmc
+
+  end interface
+
+end module dagmc_header
+
+#endif

src/geometry.F90

@@ -55,10 +55,34 @@ module geometry
 
     subroutine neighbor_lists() bind(C)
     end subroutine neighbor_lists
+
+#ifdef DAGMC
+
+    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
+
+#endif
+
   end interface
 
 contains
 
+#ifdef DAGMC
+
+  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
+
+#endif
+
 !===============================================================================
 ! 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/input_xml.F90

@@ -11,6 +11,9 @@ module input_xml
   use error,            only: fatal_error, warning, write_message, openmc_err_msg
   use geometry,         only: neighbor_lists
   use geometry_header
+#ifdef DAGMC
+  use dagmc_header
+#endif
   use hdf5_interface
   use list_header,      only: ListChar, ListInt, ListReal
   use material_header
@@ -176,7 +179,9 @@ contains
 
     ! After reading input and basic geometry setup is complete, build lists of
     ! neighboring cells for efficient tracking
-    call neighbor_lists()
+    if (.not. dagmc) then
+      call neighbor_lists()
+    end if
 
     ! Assign temperatures to cells that don't have temperatures already assigned
     call assign_temperatures()
@@ -347,6 +352,63 @@ contains
 
   end subroutine read_settings_xml_f
 
+
+#ifdef DAGMC
+
+!===============================================================================
+! READ_GEOMETRY_DAGMC reads data from a DAGMC .h5m file, checking
+! for material properties and surface boundary conditions
+! some universe information is spoofed for now
+!===============================================================================
+
+  subroutine read_geometry_dagmc()
+
+    integer :: i, j
+    integer :: univ_id
+    integer :: n_cells_in_univ
+    logical :: file_exists
+    character(MAX_LINE_LEN) :: filename
+    type(Cell),     pointer :: c
+    type(VectorInt) :: univ_ids      ! List of all universe IDs
+    type(DictIntInt) :: cells_in_univ_dict ! Used to count how many cells each
+                                           ! universe contains
+
+    ! Check if dagmc.h5m exists
+    filename = trim(path_input) // "dagmc.h5m"
+    inquire(FILE=filename, EXIST=file_exists)
+    if (.not. file_exists) then
+      call fatal_error("Geometry DAGMC file '" // trim(filename) // "' does not &
+           &exist!")
+    end if
+
+    call write_message("Reading DAGMC geometry...", 5)
+    call load_dagmc_geometry()
+    call allocate_surfaces()
+    call allocate_cells()
+
+    ! setup universe data structs
+    do i = 1, n_cells
+      c => cells(i)
+      ! additional metadata spoofing
+      univ_id = c % universe()
+
+      if (.not. cells_in_univ_dict % has(univ_id)) then
+        n_universes = n_universes + 1
+        n_cells_in_univ = 1
+        call universe_dict % set(univ_id, n_universes)
+        call univ_ids % push_back(univ_id)
+      else
+        n_cells_in_univ = 1 + cells_in_univ_dict % get(univ_id)
+      end if
+      call cells_in_univ_dict % set(univ_id, n_cells_in_univ)
+    end do
+
+    root_universe = find_root_universe()
+
+  end subroutine read_geometry_dagmc
+
+#endif
+
 !===============================================================================
 ! READ_GEOMETRY_XML reads data from a geometry.xml file and parses it, checking
 ! for errors and placing properly-formatted data in the right data structures
@@ -374,6 +436,12 @@ contains
     type(VectorInt) :: univ_ids      ! List of all universe IDs
     type(DictIntInt) :: cells_in_univ_dict ! Used to count how many cells each
                                            ! universe contains
+#ifdef DAGMC
+    if (dagmc) then
+      call read_geometry_dagmc()
+      return
+    end if
+#endif
 
     ! Display output message
     call write_message("Reading geometry XML file...", 5)
@@ -401,7 +469,6 @@ contains
 
     ! Allocate surfaces array
     allocate(surfaces(n_surfaces))
-
     do i = 1, n_surfaces
       surfaces(i) % ptr = surface_pointer(i - 1);
 
@@ -559,6 +626,40 @@ contains
 
   end subroutine read_geometry_xml
 
+  subroutine allocate_surfaces()
+    integer :: i
+
+    ! Allocate surfaces array
+    allocate(surfaces(n_surfaces))
+
+    do i = 1, n_surfaces
+      surfaces(i) % ptr = surface_pointer(i - 1);
+      ! Add surface to dictionary
+      call surface_dict % set(surfaces(i) % id(), i)
+    end do
+
+  end subroutine allocate_surfaces
+
+  subroutine allocate_cells()
+    integer :: i
+    type(Cell), pointer :: c
+
+    ! Allocate cells array
+    allocate(cells(n_cells))
+
+    do i = 1, n_cells
+      c => cells(i)
+      c % ptr = cell_pointer(i - 1)
+      ! Check to make sure 'id' hasn't been used
+      if (cell_dict % has(c % id())) then
+        call fatal_error("Two or more cells use the same unique ID: " &
+               // to_str(c % id()))
+      end if
+      ! Add cell to dictionary
+      call cell_dict % set(c % id(), i)
+    end do
+  end subroutine allocate_cells
+
 !===============================================================================
 ! READ_MATERIAL_XML reads data from a materials.xml file and parses it, checking
 ! for errors and placing properly-formatted data in the right data structures

src/relaxng/settings.rnc

@@ -57,6 +57,8 @@ element settings {
 
   element ptables { xsd:boolean }? &
 
+  element dagmc { xsd:boolean }? &
+
   element run_cmfd { xsd:boolean }? &
 
   element run_mode { xsd:string }? &

src/relaxng/settings.rng

@@ -254,6 +254,11 @@
         <data type="boolean"/>
       </element>
     </optional>
+    <optional>
+      <element name="dagmc">
+        <data type="boolean"/>
+      </element>
+    </optional>
     <optional>
       <element name="run_cmfd">
         <data type="boolean"/>

src/settings.F90

@@ -80,6 +80,9 @@ module settings
   ! Mode to run in (fixed source, eigenvalue, plotting, etc)
   integer(C_INT), bind(C) :: run_mode
 
+  ! flag for use of DAGMC geometry
+  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 DAGMC
+  if (dagmc) {
+    fatal_error("DAGMC 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) {
@@ -396,6 +408,13 @@ void read_settings_xml()
 #endif
   }
 
+#ifdef _OPENMP
+  if (dagmc && omp_get_max_threads() > 1) {
+    warning("Forcing number of threads to 1 for DAGMC simulation.");
+    omp_set_num_threads(1);
+  }
+#endif
+  
   // ==========================================================================
   // EXTERNAL SOURCE
 

src/summary.cpp

@@ -2,13 +2,22 @@
 #include "openmc/hdf5_interface.h"
 #include "openmc/lattice.h"
 #include "openmc/surface.h"
-
+#include "openmc/settings.h"
 
 namespace openmc {
 
 extern "C" void
 write_geometry(hid_t file_id) {
+
   auto geom_group = create_group(file_id, "geometry");
+  
+#ifdef DAGMC
+  if (settings::dagmc) {
+    write_attribute(geom_group, "dagmc", 1);
+    return;
+  }
+#endif
+
   write_attribute(geom_group, "n_cells", cells.size());
   write_attribute(geom_group, "n_surfaces", surfaces.size());
   write_attribute(geom_group, "n_universes", universes.size());

src/surface.cpp

@@ -138,6 +138,8 @@ void read_coeffs(pugi::xml_node surf_node, int surf_id, double &c1, double &c2,
 // Surface implementation
 //==============================================================================
 
+Surface::Surface() {} // empty constructor
+
 Surface::Surface(pugi::xml_node surf_node)
 {
   if (check_for_node(surf_node, "id")) {
@@ -207,8 +209,11 @@ Surface::reflect(Position r, Direction u) const
   return u -= (2.0 * projection / magnitude) * n;
 }
 
+CSGSurface::CSGSurface() : Surface{} {};
+CSGSurface::CSGSurface(pugi::xml_node surf_node) : Surface{surf_node} {};
+
 void
-Surface::to_hdf5(hid_t group_id) const
+CSGSurface::to_hdf5(hid_t group_id) const
 {
   std::string group_name {"surface "};
   group_name += std::to_string(id_);
@@ -239,12 +244,62 @@ Surface::to_hdf5(hid_t group_id) const
   close_group(surf_group);
 }
 
+//==============================================================================
+// DAGSurface implementation
+//==============================================================================
+#ifdef DAGMC
+DAGSurface::DAGSurface() : Surface{} {} // empty constructor
+
+double DAGSurface::evaluate(Position r) const
+{
+  return 0.0;
+}
+
+double
+DAGSurface::distance(Position r, Direction u, bool coincident) const
+{
+  moab::ErrorCode rval;
+  moab::EntityHandle surf = dagmc_ptr_->entity_by_id(2, id_);
+  moab::EntityHandle hit_surf;
+  double dist;
+  double pnt[3] = {r.x, r.y, r.z};
+  double dir[3] = {u.x, u.y, u.z};
+  rval = dagmc_ptr_->ray_fire(surf, pnt, dir, hit_surf, dist, NULL, 0, 0);
+  MB_CHK_ERR_CONT(rval);
+  if (dist < 0.0) dist = INFTY;
+  return dist;
+}
+
+Direction DAGSurface::normal(Position r) const
+{
+  moab::ErrorCode rval;
+  Direction u;
+  moab::EntityHandle surf = dagmc_ptr_->entity_by_id(2, id_);
+  double pnt[3] = {r.x, r.y, r.z};
+  double dir[3] = {u.x, u.y, u.z};
+  rval = dagmc_ptr_->get_angle(surf, pnt, dir);
+  MB_CHK_ERR_CONT(rval);
+  return u;
+}
+
+BoundingBox DAGSurface::bounding_box() const
+{
+  moab::ErrorCode rval;
+  moab::EntityHandle surf = dagmc_ptr_->entity_by_id(2, id_);
+  double min[3], max[3];
+  rval = dagmc_ptr_->getobb(surf, min, max);
+  MB_CHK_ERR_CONT(rval);
+  return {min[0], max[0], min[1], max[1], min[2], max[2]};
+}
+
+void DAGSurface::to_hdf5(hid_t group_id) const {}
+#endif
 //==============================================================================
 // PeriodicSurface implementation
 //==============================================================================
 
 PeriodicSurface::PeriodicSurface(pugi::xml_node surf_node)
-  : Surface {surf_node}
+  : CSGSurface {surf_node}
 {
   if (check_for_node(surf_node, "periodic_surface_id")) {
     i_periodic_ = std::stoi(get_node_value(surf_node, "periodic_surface_id"));
@@ -580,7 +635,7 @@ axis_aligned_cylinder_normal(Position r, double offset1, double offset2)
 //==============================================================================
 
 SurfaceXCylinder::SurfaceXCylinder(pugi::xml_node surf_node)
-  : Surface(surf_node)
+  : CSGSurface(surf_node)
 {
   read_coeffs(surf_node, id_, y0_, z0_, radius_);
 }
@@ -614,7 +669,7 @@ void SurfaceXCylinder::to_hdf5_inner(hid_t group_id) const
 //==============================================================================
 
 SurfaceYCylinder::SurfaceYCylinder(pugi::xml_node surf_node)
-  : Surface(surf_node)
+  : CSGSurface(surf_node)
 {
   read_coeffs(surf_node, id_, x0_, z0_, radius_);
 }
@@ -647,7 +702,7 @@ void SurfaceYCylinder::to_hdf5_inner(hid_t group_id) const
 //==============================================================================
 
 SurfaceZCylinder::SurfaceZCylinder(pugi::xml_node surf_node)
-  : Surface(surf_node)
+  : CSGSurface(surf_node)
 {
   read_coeffs(surf_node, id_, x0_, y0_, radius_);
 }
@@ -680,7 +735,7 @@ void SurfaceZCylinder::to_hdf5_inner(hid_t group_id) const
 //==============================================================================
 
 SurfaceSphere::SurfaceSphere(pugi::xml_node surf_node)
-  : Surface(surf_node)
+  : CSGSurface(surf_node)
 {
   read_coeffs(surf_node, id_, x0_, y0_, z0_, radius_);
 }
@@ -833,7 +888,7 @@ axis_aligned_cone_normal(Position r, double offset1, double offset2,
 //==============================================================================
 
 SurfaceXCone::SurfaceXCone(pugi::xml_node surf_node)
-  : Surface(surf_node)
+  : CSGSurface(surf_node)
 {
   read_coeffs(surf_node, id_, x0_, y0_, z0_, radius_sq_);
 }
@@ -866,7 +921,7 @@ void SurfaceXCone::to_hdf5_inner(hid_t group_id) const
 //==============================================================================
 
 SurfaceYCone::SurfaceYCone(pugi::xml_node surf_node)
-  : Surface(surf_node)
+  : CSGSurface(surf_node)
 {
   read_coeffs(surf_node, id_, x0_, y0_, z0_, radius_sq_);
 }
@@ -899,7 +954,7 @@ void SurfaceYCone::to_hdf5_inner(hid_t group_id) const
 //==============================================================================
 
 SurfaceZCone::SurfaceZCone(pugi::xml_node surf_node)
-  : Surface(surf_node)
+  : CSGSurface(surf_node)
 {
   read_coeffs(surf_node, id_, x0_, y0_, z0_, radius_sq_);
 }
@@ -932,7 +987,7 @@ void SurfaceZCone::to_hdf5_inner(hid_t group_id) const
 //==============================================================================
 
 SurfaceQuadric::SurfaceQuadric(pugi::xml_node surf_node)
-  : Surface(surf_node)
+  : CSGSurface(surf_node)
 {
   read_coeffs(surf_node, id_, A_, B_, C_, D_, E_, F_, G_, H_, J_, K_);
 }

src/tracking.F90

@@ -7,6 +7,10 @@ module tracking
   use geometry_header,    only: cells
   use geometry,           only: find_cell, distance_to_boundary, cross_lattice,&
                                 check_cell_overlap
+#ifdef DAGMC
+  use geometry,           only: next_cell
+#endif
+
   use material_header,    only: materials, Material
   use message_passing
   use mgxs_interface
@@ -309,6 +313,7 @@ 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
@@ -467,6 +472,21 @@ contains
     ! ==========================================================================
     ! SEARCH NEIGHBOR LISTS FOR NEXT CELL
 
+#ifdef DAGMC
+    if (dagmc) then
+      i_cell = next_cell(cells(p % last_cell(1) + 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-1 ! decrement for C++ indexing
+      p % cell_instance = 1
+      p % material = cells(i_cell) % material(1)
+      p % sqrtKT = cells(i_cell) % sqrtKT(1)
+      return
+    end if
+#endif
+
     call find_cell(p, found, p % surface)
     if (found) return
 

src/xml_interface.F90

@@ -12,6 +12,7 @@ module xml_interface
   public :: check_for_node
   public :: get_node_list
   public :: get_node_value
+  public :: get_node_value_bool
   public :: get_node_array
   public :: node_value_string
   public :: node_word_count

tests/regression_tests/dagmc/materials.xml

@@ -0,0 +1,16 @@
+<?xml version="1.0"?>
+<materials>
+
+  <material id="40">
+    <density value="11" units="g/cc" />
+    <nuclide name="U235" ao="1.0" />
+  </material>
+
+  <material id="41">
+    <density value="1.0" units="g/cc" />
+    <nuclide name="H1" ao="2.0" />
+    <nuclide name="O16" ao="1.0" />
+    <sab name="c_H_in_H2O"/>
+  </material>
+  
+</materials>

tests/regression_tests/dagmc/results_true.dat

@@ -0,0 +1,5 @@
+k-combined:
+1.115067E+00 5.423808E-02
+tally 1:
+8.543144E+00
+1.530584E+01

tests/regression_tests/dagmc/settings.xml

@@ -0,0 +1,16 @@
+<?xml version="1.0"?>
+<settings>
+
+  <run_mode>eigenvalue</run_mode>
+  <dagmc>true</dagmc>
+  <batches>5</batches>
+  <inactive>0</inactive>
+  <particles>100</particles>
+  <!-- Starting source -->
+  <source>
+    <space type="box">
+      <parameters>-4 -4 -4  4  4  4</parameters>
+    </space>
+  </source>
+
+</settings>

tests/regression_tests/dagmc/tallies.xml

@@ -0,0 +1,13 @@
+<?xml version="1.0"?>
+<tallies>
+
+  <filter id="1" type="cell">
+    <bins>1</bins>
+  </filter>
+
+  <tally id="1">
+    <filters>1</filters>
+    <scores>total </scores>
+  </tally>
+
+</tallies>

tests/regression_tests/dagmc/test.py

@@ -0,0 +1,12 @@
+from tests.testing_harness import TestHarness
+import os
+import pytest
+import openmc
+
+pytestmark = pytest.mark.skipif(
+    not openmc.capi.dagmc_enabled,
+    reason="DAGMC CAD geometry is not enabled.")
+
+def test_dagmc():
+    harness = TestHarness('statepoint.5.h5')
+    harness.main()    

tools/ci/travis-install-dagmc.sh

@@ -0,0 +1,36 @@
+
+#!/bin/bash
+set -ex
+
+# MOAB Variables
+MOAB_BRANCH='Version5.0'
+MOAB_REPO='https://bitbucket.org/fathomteam/moab/'
+MOAB_INSTALL_DIR=$HOME/MOAB/
+
+# DAGMC Variables
+DAGMC_BRANCH='develop'
+DAGMC_REPO='https://github.com/svalinn/dagmc'
+DAGMC_INSTALL_DIR=$HOME/DAGMC/
+
+CURRENT_DIR=$(pwd)
+
+# MOAB Install
+cd $HOME
+mkdir MOAB && cd MOAB
+git clone -b $MOAB_BRANCH $MOAB_REPO
+mkdir build && cd build
+cmake ../moab -DENABLE_HDF5=ON -DCMAKE_INSTALL_PREFIX=$MOAB_INSTALL_DIR
+make -j && make -j test install
+rm -rf $HOME/MOAB/moab
+export LD_LIBRARY_PATH=$MOAB_INSTALL_DIR/lib:$LD_LIBRARY_PATH
+
+# DAGMC Install
+mkdir DAGMC && cd DAGMC
+git clone -b $DAGMC_BRANCH $DAGMC_REPO
+mkdir build && cd build
+cmake ../dagmc -DBUILD_TALLY=ON -DCMAKE_INSTALL_PREFIX=$DAGMC_INSTALL_DIR
+make -j install
+rm -rf $HOME/DAGMC/dagmc
+export LD_LIBRARY_PATH=$DAGMC_INSTALL_DIR/lib:$LD_LIBRARY_PATH
+
+cd $CURRENT_DIR

tools/ci/travis-install.py

@@ -2,7 +2,6 @@ import os
 import shutil
 import subprocess
 
-
 def which(program):
     def is_exe(fpath):
         return os.path.isfile(fpath) and os.access(fpath, os.X_OK)
@@ -20,7 +19,7 @@ def which(program):
     return None
 
 
-def install(omp=False, mpi=False, phdf5=False):
+def install(omp=False, mpi=False, phdf5=False, dagmc=False):
     # Create build directory and change to it
     shutil.rmtree('build', ignore_errors=True)
     os.mkdir('build')
@@ -48,23 +47,26 @@ def install(omp=False, mpi=False, phdf5=False):
     else:
         cmake_cmd.append('-DHDF5_PREFER_PARALLEL=OFF')
 
+    if dagmc:
+        cmake_cmd.append('-Ddagmc=ON')
+        
     # Build and install
     cmake_cmd.append('..')
     print(' '.join(cmake_cmd))
     subprocess.check_call(cmake_cmd)
-    subprocess.check_call(['make', '-j'])
+    subprocess.check_call(['make', '-j4'])
     subprocess.check_call(['sudo', 'make', 'install'])
 
-
 def main():
     # Convert Travis matrix environment variables into arguments for install()
     omp = (os.environ.get('OMP') == 'y')
     mpi = (os.environ.get('MPI') == 'y')
     phdf5 = (os.environ.get('PHDF5') == 'y')
 
-    # Build and install
-    install(omp, mpi, phdf5)
+    dagmc = (os.environ.get('DAGMC') == 'y')
 
+    # Build and install
+    install(omp, mpi, phdf5, dagmc)
 
 if __name__ == '__main__':
     main()

tools/ci/travis-install.sh

@@ -4,6 +4,11 @@ set -ex
 # Install NJOY 2016
 ./tools/ci/travis-install-njoy.sh
 
+# Install DAGMC if needed
+if [[ $DAGMC = 'y' ]]; then
+    ./tools/ci/travis-install-dagmc.sh
+fi
+
 # Upgrade pip before doing anything else
 pip install --upgrade pip