minor simplification of csg/dag geometry type checking

include/openmc/cell.h

@@ -349,12 +349,8 @@ class Cell {
 
   vector<int32_t> offset_; //!< Distribcell offset table
 
-  // Accessors
-  const GeometryType& geom_type() const { return geom_type_; }
-  GeometryType& geom_type() { return geom_type_; }
-
-private:
-  GeometryType geom_type_; //!< Geometric representation type (CSG, DAGMC)
+  // Right now, either CSG or DAGMC cells are used.
+  virtual GeometryType geom_type() const = 0;
 };
 
 struct CellInstanceItem {
@@ -368,7 +364,7 @@ class CSGCell : public Cell {
 public:
   //----------------------------------------------------------------------------
   // Constructors
-  CSGCell();
+  CSGCell() = default;
   explicit CSGCell(pugi::xml_node cell_node);
 
   //----------------------------------------------------------------------------
@@ -395,6 +391,8 @@ class CSGCell : public Cell {
 
   bool is_simple() const override { return region_.is_simple(); }
 
+  virtual GeometryType geom_type() const override { return GeometryType::CSG; }
+
 protected:
   //! Returns the beginning position of a parenthesis block (immediately before
   //! two surface tokens) in the RPN given a starting position at the end of

include/openmc/dagmc.h

@@ -53,6 +53,8 @@ class DAGSurface : public Surface {
 
   inline void to_hdf5_inner(hid_t group_id) const override {};
 
+  virtual GeometryType geom_type() const override { return GeometryType::DAG; }
+
   // Accessor methods
   moab::DagMC* dagmc_ptr() const { return dagmc_ptr_.get(); }
   int32_t dag_index() const { return dag_index_; }
@@ -77,6 +79,8 @@ class DAGCell : public Cell {
 
   void to_hdf5_inner(hid_t group_id) const override;
 
+  virtual GeometryType geom_type() const override { return GeometryType::DAG; }
+
   // Accessor methods
   moab::DagMC* dagmc_ptr() const { return dagmc_ptr_.get(); }
   int32_t dag_index() const { return dag_index_; }
@@ -164,6 +168,8 @@ class DAGUniverse : public Universe {
 
   void to_hdf5(hid_t universes_group) const override;
 
+  virtual GeometryType geom_type() const override { return GeometryType::DAG; }
+
   // Data Members
   std::shared_ptr<moab::DagMC>
     dagmc_instance_;        //!< DAGMC Instance for this universe

include/openmc/surface.h

@@ -90,30 +90,26 @@ class Surface {
   //! Get the BoundingBox for this surface.
   virtual BoundingBox bounding_box(bool /*pos_side*/) const { return {}; }
 
-  // Accessors
-  const GeometryType& geom_type() const { return geom_type_; }
-  GeometryType& geom_type() { return geom_type_; }
-
-private:
-  GeometryType geom_type_; //!< Geometry type indicator (CSG or DAGMC)
+  /* Must specify if this is a CSG or DAGMC-type surface. Only
+   * the DAGMC surface should return the DAG type geometry, so
+   * by default, this returns the CSG. The main difference is that
+   * if the geom_type is found to be DAG in the geometry handling code,
+   * some DAGMC-specific operations get carried out like resetting
+   * the particle's intersection history when necessary.
+   */
+  virtual GeometryType geom_type() const { return GeometryType::CSG; }
 
 protected:
   virtual void to_hdf5_inner(hid_t group_id) const = 0;
 };
 
-class CSGSurface : public Surface {
-public:
-  explicit CSGSurface(pugi::xml_node surf_node);
-  CSGSurface();
-};
-
 //==============================================================================
 //! A plane perpendicular to the x-axis.
 //
 //! The plane is described by the equation \f$x - x_0 = 0\f$
 //==============================================================================
 
-class SurfaceXPlane : public CSGSurface {
+class SurfaceXPlane : public Surface {
 public:
   explicit SurfaceXPlane(pugi::xml_node surf_node);
   double evaluate(Position r) const override;
@@ -131,7 +127,7 @@ class SurfaceXPlane : public CSGSurface {
 //! The plane is described by the equation \f$y - y_0 = 0\f$
 //==============================================================================
 
-class SurfaceYPlane : public CSGSurface {
+class SurfaceYPlane : public Surface {
 public:
   explicit SurfaceYPlane(pugi::xml_node surf_node);
   double evaluate(Position r) const override;
@@ -149,7 +145,7 @@ class SurfaceYPlane : public CSGSurface {
 //! The plane is described by the equation \f$z - z_0 = 0\f$
 //==============================================================================
 
-class SurfaceZPlane : public CSGSurface {
+class SurfaceZPlane : public Surface {
 public:
   explicit SurfaceZPlane(pugi::xml_node surf_node);
   double evaluate(Position r) const override;
@@ -167,7 +163,7 @@ class SurfaceZPlane : public CSGSurface {
 //! The plane is described by the equation \f$A x + B y + C z - D = 0\f$
 //==============================================================================
 
-class SurfacePlane : public CSGSurface {
+class SurfacePlane : public Surface {
 public:
   explicit SurfacePlane(pugi::xml_node surf_node);
   double evaluate(Position r) const override;
@@ -185,7 +181,7 @@ class SurfacePlane : public CSGSurface {
 //! \f$(y - y_0)^2 + (z - z_0)^2 - R^2 = 0\f$
 //==============================================================================
 
-class SurfaceXCylinder : public CSGSurface {
+class SurfaceXCylinder : public Surface {
 public:
   explicit SurfaceXCylinder(pugi::xml_node surf_node);
   double evaluate(Position r) const override;
@@ -204,7 +200,7 @@ class SurfaceXCylinder : public CSGSurface {
 //! \f$(x - x_0)^2 + (z - z_0)^2 - R^2 = 0\f$
 //==============================================================================
 
-class SurfaceYCylinder : public CSGSurface {
+class SurfaceYCylinder : public Surface {
 public:
   explicit SurfaceYCylinder(pugi::xml_node surf_node);
   double evaluate(Position r) const override;
@@ -223,7 +219,7 @@ class SurfaceYCylinder : public CSGSurface {
 //! \f$(x - x_0)^2 + (y - y_0)^2 - R^2 = 0\f$
 //==============================================================================
 
-class SurfaceZCylinder : public CSGSurface {
+class SurfaceZCylinder : public Surface {
 public:
   explicit SurfaceZCylinder(pugi::xml_node surf_node);
   double evaluate(Position r) const override;
@@ -242,7 +238,7 @@ class SurfaceZCylinder : public CSGSurface {
 //! \f$(x - x_0)^2 + (y - y_0)^2 + (z - z_0)^2 - R^2 = 0\f$
 //==============================================================================
 
-class SurfaceSphere : public CSGSurface {
+class SurfaceSphere : public Surface {
 public:
   explicit SurfaceSphere(pugi::xml_node surf_node);
   double evaluate(Position r) const override;
@@ -261,7 +257,7 @@ class SurfaceSphere : public CSGSurface {
 //! \f$(y - y_0)^2 + (z - z_0)^2 - R^2 (x - x_0)^2 = 0\f$
 //==============================================================================
 
-class SurfaceXCone : public CSGSurface {
+class SurfaceXCone : public Surface {
 public:
   explicit SurfaceXCone(pugi::xml_node surf_node);
   double evaluate(Position r) const override;
@@ -279,7 +275,7 @@ class SurfaceXCone : public CSGSurface {
 //! \f$(x - x_0)^2 + (z - z_0)^2 - R^2 (y - y_0)^2 = 0\f$
 //==============================================================================
 
-class SurfaceYCone : public CSGSurface {
+class SurfaceYCone : public Surface {
 public:
   explicit SurfaceYCone(pugi::xml_node surf_node);
   double evaluate(Position r) const override;
@@ -297,7 +293,7 @@ class SurfaceYCone : public CSGSurface {
 //! \f$(x - x_0)^2 + (y - y_0)^2 - R^2 (z - z_0)^2 = 0\f$
 //==============================================================================
 
-class SurfaceZCone : public CSGSurface {
+class SurfaceZCone : public Surface {
 public:
   explicit SurfaceZCone(pugi::xml_node surf_node);
   double evaluate(Position r) const override;
@@ -315,7 +311,7 @@ class SurfaceZCone : public CSGSurface {
 //! 0\f$
 //==============================================================================
 
-class SurfaceQuadric : public CSGSurface {
+class SurfaceQuadric : public Surface {
 public:
   explicit SurfaceQuadric(pugi::xml_node surf_node);
   double evaluate(Position r) const override;
@@ -333,7 +329,7 @@ class SurfaceQuadric : public CSGSurface {
 //! \f$(x-x_0)^2/B^2 + (\sqrt{(y-y_0)^2 + (z-z_0)^2} - A)^2/C^2 -1 \f$
 //==============================================================================
 
-class SurfaceXTorus : public CSGSurface {
+class SurfaceXTorus : public Surface {
 public:
   explicit SurfaceXTorus(pugi::xml_node surf_node);
   double evaluate(Position r) const override;
@@ -350,7 +346,7 @@ class SurfaceXTorus : public CSGSurface {
 //! \f$(y-y_0)^2/B^2 + (\sqrt{(x-x_0)^2 + (z-z_0)^2} - A)^2/C^2 -1 \f$
 //==============================================================================
 
-class SurfaceYTorus : public CSGSurface {
+class SurfaceYTorus : public Surface {
 public:
   explicit SurfaceYTorus(pugi::xml_node surf_node);
   double evaluate(Position r) const override;
@@ -367,7 +363,7 @@ class SurfaceYTorus : public CSGSurface {
 //! \f$(z-z_0)^2/B^2 + (\sqrt{(x-x_0)^2 + (y-y_0)^2} - A)^2/C^2 -1 \f$
 //==============================================================================
 
-class SurfaceZTorus : public CSGSurface {
+class SurfaceZTorus : public Surface {
 public:
   explicit SurfaceZTorus(pugi::xml_node surf_node);
   double evaluate(Position r) const override;

include/openmc/universe.h

@@ -38,13 +38,13 @@ class Universe {
 
   BoundingBox bounding_box() const;
 
-  const GeometryType& geom_type() const { return geom_type_; }
-  GeometryType& geom_type() { return geom_type_; }
+  /* By default, universes are CSG universes. The DAGMC
+   * universe overrides standard behaviors, and in the future,
+   * other things might too.
+   */
+  virtual GeometryType geom_type() const { return GeometryType::CSG; }
 
   unique_ptr<UniversePartitioner> partitioner_;
-
-private:
-  GeometryType geom_type_ = GeometryType::CSG;
 };
 
 //==============================================================================

src/cell.cpp

@@ -249,16 +249,8 @@ void Cell::to_hdf5(hid_t cell_group) const
 // CSGCell implementation
 //==============================================================================
 
-// default constructor
-CSGCell::CSGCell()
-{
-  geom_type() = GeometryType::CSG;
-}
-
 CSGCell::CSGCell(pugi::xml_node cell_node)
 {
-  geom_type() = GeometryType::CSG;
-
   if (check_for_node(cell_node, "id")) {
     id_ = std::stoi(get_node_value(cell_node, "id"));
   } else {

src/dagmc.cpp

@@ -129,8 +129,6 @@ void DAGUniverse::set_id()
 
 void DAGUniverse::initialize()
 {
-  geom_type() = GeometryType::DAG;
-
 #ifdef OPENMC_UWUW
   // read uwuw materials from the .h5m file if present
   read_uwuw_materials();
@@ -661,10 +659,7 @@ void DAGUniverse::override_assign_material(std::unique_ptr<DAGCell>& c) const
 //==============================================================================
 
 DAGCell::DAGCell(std::shared_ptr<moab::DagMC> dag_ptr, int32_t dag_idx)
-  : Cell {}, dagmc_ptr_(dag_ptr), dag_index_(dag_idx)
-{
-  geom_type() = GeometryType::DAG;
-};
+  : Cell {}, dagmc_ptr_(dag_ptr), dag_index_(dag_idx) {};
 
 std::pair<double, int32_t> DAGCell::distance(
   Position r, Direction u, int32_t on_surface, GeometryState* p) const
@@ -765,9 +760,7 @@ BoundingBox DAGCell::bounding_box() const
 
 DAGSurface::DAGSurface(std::shared_ptr<moab::DagMC> dag_ptr, int32_t dag_idx)
   : Surface {}, dagmc_ptr_(dag_ptr), dag_index_(dag_idx)
-{
-  geom_type() = GeometryType::DAG;
-} // empty constructor
+{} // empty constructor
 
 moab::EntityHandle DAGSurface::mesh_handle() const
 {