From cb45c97b70d52ae6b8a4b9f74b3906560de15ac5 Mon Sep 17 00:00:00 2001 From: shriwise Date: Fri, 23 Feb 2018 11:21:53 -0600 Subject: [PATCH] Adding Surface population in CAD. Creating a derived CSGSurfce class for support of (currently) non-CAD operations related to hdf5 summary files. --- include/openmc/cell.h | 5 ++--- include/openmc/surface.h | 35 ++++++++++++++++++++++------------- src/cad.cpp | 23 +++++++++++++++++++++-- src/cad.h | 1 + src/cell.cpp | 2 +- src/surface.cpp | 35 +++++++++++++++++++++++++---------- 6 files changed, 72 insertions(+), 29 deletions(-) diff --git a/include/openmc/cell.h b/include/openmc/cell.h index 41209246e7..f819d41f44 100644 --- a/include/openmc/cell.h +++ b/include/openmc/cell.h @@ -155,10 +155,9 @@ protected: #ifdef CAD class CADCell : public Cell { - moab::DagMC *dagmc_ptr; - public: - explicit CADCell(); + moab::DagMC *dagmc_ptr; + explicit CADCell(); }; #endif diff --git a/include/openmc/surface.h b/include/openmc/surface.h index 61ae8ef53b..8a6b340ed7 100644 --- a/include/openmc/surface.h +++ b/include/openmc/surface.h @@ -105,23 +105,32 @@ public: //! \return Normal direction virtual Direction normal(Position r) const = 0; +}; + +class CSGSurface : public Surface +{ + public: + explicit CSGSurface(pugi::xml_node surf_node); + explicit CSGSurface(); + //! 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 void to_hdf5(hid_t group_id) const; -protected: - virtual void to_hdf5_inner(hid_t group_id) const = 0; + + protected: + virtual void to_hdf5_inner(hid_t group_id) const = 0; }; - + //============================================================================== //! A `Surface` representing a CAD-based surface in DAGMC. //============================================================================== #ifdef CAD class CADSurface : public Surface { - moab::DagMC* dagmc_ptr; public: + moab::DagMC* dagmc_ptr; explicit CADSurface(); double evaluate(const double xyz[3]) const; double distance(const double xyz[3], const double uvw[3], @@ -139,7 +148,7 @@ class CADSurface : public Surface //! `XPlane`-`YPlane` pairs. //============================================================================== -class PeriodicSurface : public Surface +class PeriodicSurface : public CSGSurface { public: int i_periodic_{C_NONE}; //!< Index of corresponding periodic surface @@ -248,7 +257,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: @@ -266,7 +275,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: @@ -284,7 +293,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: @@ -302,7 +311,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: @@ -320,7 +329,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: @@ -338,7 +347,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: @@ -356,7 +365,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: @@ -373,7 +382,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_; diff --git a/src/cad.cpp b/src/cad.cpp index f07c8cd39d..90b325866c 100644 --- a/src/cad.cpp +++ b/src/cad.cpp @@ -1,9 +1,13 @@ #include "cad.h" +moab::DagMC* DAGMC; + void load_cad_geometry_c() { - moab::DagMC* DAGMC = new moab::DagMC(); + if(!DAGMC) { + DAGMC = new moab::DagMC(); + } moab::ErrorCode rval = DAGMC->load_file("dagmc.h5m"); MB_CHK_ERR_CONT(rval); @@ -12,13 +16,28 @@ void load_cad_geometry_c() MB_CHK_ERR_CONT(rval); // initialize cell objects - for(int i = 0; i < DAGMC->num_entities(3); i++) + openmc::n_cells = DAGMC->num_entities(3); + for(int i = 0; i < openmc::n_cells; i++) { // set cell ids using global IDs openmc::CADCell c = openmc::CADCell(); c.id = DAGMC->id_by_index(3, i); + c.dagmc_ptr = DAGMC; openmc::cells_c.push_back(c); } + + // initialize surface objects + openmc::n_surfaces = DAGMC->num_entities(2); + openmc::surfaces_c = new openmc::Surface*[openmc::n_surfaces]; + for(int i = 0; i < openmc::n_surfaces; i++) + { + // set cell ids using global IDs + openmc::CADSurface* s = new openmc::CADSurface(); + s->id = DAGMC->id_by_index(2, i); + s->dagmc_ptr = DAGMC; + openmc::surfaces_c[i] = s; + } + return; } diff --git a/src/cad.h b/src/cad.h index adef8f4325..f88c4be3f6 100644 --- a/src/cad.h +++ b/src/cad.h @@ -3,6 +3,7 @@ #include "DagMC.hpp" #include "cell.h" +#include "surface.h" extern "C" void load_cad_geometry_c(); diff --git a/src/cell.cpp b/src/cell.cpp index b417b11995..35c69f371c 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -601,7 +601,7 @@ Cell::contains_complex(Position r, Direction u, int32_t on_surface) const #ifdef CAD CADCell::CADCell() : Cell{} {}; #endif - + //============================================================================== // Non-method functions extern "C" void diff --git a/src/surface.cpp b/src/surface.cpp index 0111d77bc5..53d9d2eebf 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -179,6 +179,7 @@ Surface::Surface(pugi::xml_node surf_node) } + bool Surface::sense(Position r, Direction u) const { @@ -209,8 +210,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_); @@ -267,7 +271,7 @@ BoundingBox CADSurface::bounding_box() const { return BoundingBox(); } //============================================================================== 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")); @@ -603,7 +607,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_); } @@ -637,7 +641,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_); } @@ -670,7 +674,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_); } @@ -703,7 +707,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_); } @@ -856,7 +860,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_); } @@ -889,7 +893,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_); } @@ -922,7 +926,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_); } @@ -955,7 +959,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_); } @@ -1253,6 +1257,17 @@ extern "C" { int surface_i_periodic(PeriodicSurface* surf) {return surf->i_periodic_;} + void surface_normal(Surface* surf, double xyz[3], double uvw[3]) + { + Position r {xyz}; + Direction u = surf->normal(r); + uvw[0] = u.x; + uvw[1] = u.y; + uvw[2] = u.z; + } + + void surface_to_hdf5(CSGSurface* surf, hid_t group) {surf->to_hdf5(group);} + bool surface_periodic(PeriodicSurface* surf, PeriodicSurface* other, double xyz[3], double uvw[3])