1 #ifndef OPENMC_SURFACE_H 2 #define OPENMC_SURFACE_H 3 4 #include <limits> // For numeric_limits 5 #include <string> 6 #include <unordered_map> 7 8 #include "hdf5.h" 9 #include "pugixml.hpp" 10 11 #include "dagmc.h" 12 #include "openmc/boundary_condition.h" 13 #include "openmc/constants.h" 14 #include "openmc/memory.h" // for unique_ptr 15 #include "openmc/particle.h" 16 #include "openmc/position.h" 17 #include "openmc/vector.h" 18 19 namespace openmc { 20 21 //============================================================================== 22 // Global variables 23 //============================================================================== 24 25 class Surface; 26 27 namespace model { 28 extern std::unordered_map<int, int> surface_map; 29 extern vector<unique_ptr<Surface>> surfaces; 30 } // namespace model 31 32 //============================================================================== 33 //! Coordinates for an axis-aligned cuboid bounding a geometric object. 34 //============================================================================== 35 36 struct BoundingBox 37 { 38 double xmin = -INFTY; 39 double xmax = INFTY; 40 double ymin = -INFTY; 41 double ymax = INFTY; 42 double zmin = -INFTY; 43 double zmax = INFTY; 44 45 46 inline BoundingBox operator &(const BoundingBox& other) { 47 BoundingBox result = *this; 48 return result &= other; 49 } 50 51 inline BoundingBox operator |(const BoundingBox& other) { 52 BoundingBox result = *this; 53 return result |= other; 54 } 55 56 // intersect operator 57 inline BoundingBox& operator &=(const BoundingBox& other) { 58 xmin = std::max(xmin, other.xmin); 59 xmax = std::min(xmax, other.xmax); 60 ymin = std::max(ymin, other.ymin); 61 ymax = std::min(ymax, other.ymax); 62 zmin = std::max(zmin, other.zmin); 63 zmax = std::min(zmax, other.zmax); 64 return *this; 65 } 66 67 // union operator 68 inline BoundingBox& operator |=(const BoundingBox& other) { 69 xmin = std::min(xmin, other.xmin); 70 xmax = std::max(xmax, other.xmax); 71 ymin = std::min(ymin, other.ymin); 72 ymax = std::max(ymax, other.ymax); 73 zmin = std::min(zmin, other.zmin); 74 zmax = std::max(zmax, other.zmax); 75 return *this; 76 } 77 78 }; 79 80 //============================================================================== 81 //! A geometry primitive used to define regions of 3D space. 82 //============================================================================== 83 84 class Surface 85 { 86 public: 87 88 int id_; //!< Unique ID 89 std::string name_; //!< User-defined name 90 std::shared_ptr<BoundaryCondition> bc_ {nullptr}; //!< Boundary condition 91 bool surf_source_ {false}; //!< Activate source banking for the surface? 92 93 explicit Surface(pugi::xml_node surf_node); 94 Surface(); 95 ~Surface()96 virtual ~Surface() {} 97 98 //! Determine which side of a surface a point lies on. 99 //! \param r The 3D Cartesian coordinate of a point. 100 //! \param u A direction used to "break ties" and pick a sense when the 101 //! point is very close to the surface. 102 //! \return true if the point is on the "positive" side of the surface and 103 //! false otherwise. 104 bool sense(Position r, Direction u) const; 105 106 //! Determine the direction of a ray reflected from the surface. 107 //! \param[in] r The point at which the ray is incident. 108 //! \param[in] u Incident direction of the ray 109 //! \param[inout] p Pointer to the particle 110 //! \return Outgoing direction of the ray 111 virtual Direction reflect(Position r, Direction u, Particle* p) const; 112 113 virtual Direction diffuse_reflect(Position r, Direction u, 114 uint64_t* seed) const; 115 116 //! Evaluate the equation describing the surface. 117 //! 118 //! Surfaces can be described by some function f(x, y, z) = 0. This member 119 //! function evaluates that mathematical function. 120 //! \param r A 3D Cartesian coordinate. 121 virtual double evaluate(Position r) const = 0; 122 123 //! Compute the distance between a point and the surface along a ray. 124 //! \param r A 3D Cartesian coordinate. 125 //! \param u The direction of the ray. 126 //! \param coincident A hint to the code that the given point should lie 127 //! exactly on the surface. 128 virtual double distance(Position r, Direction u, bool coincident) const = 0; 129 130 //! Compute the local outward normal direction of the surface. 131 //! \param r A 3D Cartesian coordinate. 132 //! \return Normal direction 133 virtual Direction normal(Position r) const = 0; 134 135 //! Write all information needed to reconstruct the surface to an HDF5 group. 136 //! \param group_id An HDF5 group id. 137 virtual void to_hdf5(hid_t group_id) const = 0; 138 139 //! Get the BoundingBox for this surface. bounding_box(bool)140 virtual BoundingBox bounding_box(bool /*pos_side*/) const { return {}; } 141 }; 142 143 class CSGSurface : public Surface 144 { 145 public: 146 explicit CSGSurface(pugi::xml_node surf_node); 147 CSGSurface(); 148 149 void to_hdf5(hid_t group_id) const; 150 151 protected: 152 virtual void to_hdf5_inner(hid_t group_id) const = 0; 153 }; 154 155 //============================================================================== 156 //! A `Surface` representing a DAGMC-based surface in DAGMC. 157 //============================================================================== 158 #ifdef DAGMC 159 class DAGSurface : public Surface 160 { 161 public: 162 DAGSurface(); 163 164 double evaluate(Position r) const; 165 double distance(Position r, Direction u, bool coincident) const; 166 Direction normal(Position r) const; 167 Direction reflect(Position r, Direction u, Particle* p) const; 168 169 void to_hdf5(hid_t group_id) const; 170 171 moab::DagMC* dagmc_ptr_; //!< Pointer to DagMC instance 172 int32_t dag_index_; //!< DagMC index of surface 173 }; 174 #endif 175 176 //============================================================================== 177 //! A plane perpendicular to the x-axis. 178 // 179 //! The plane is described by the equation \f$x - x_0 = 0\f$ 180 //============================================================================== 181 182 class SurfaceXPlane : public CSGSurface 183 { 184 public: 185 explicit SurfaceXPlane(pugi::xml_node surf_node); 186 double evaluate(Position r) const; 187 double distance(Position r, Direction u, bool coincident) const; 188 Direction normal(Position r) const; 189 void to_hdf5_inner(hid_t group_id) const; 190 BoundingBox bounding_box(bool pos_side) const; 191 192 double x0_; 193 }; 194 195 //============================================================================== 196 //! A plane perpendicular to the y-axis. 197 // 198 //! The plane is described by the equation \f$y - y_0 = 0\f$ 199 //============================================================================== 200 201 class SurfaceYPlane : public CSGSurface 202 { 203 public: 204 explicit SurfaceYPlane(pugi::xml_node surf_node); 205 double evaluate(Position r) const; 206 double distance(Position r, Direction u, bool coincident) const; 207 Direction normal(Position r) const; 208 void to_hdf5_inner(hid_t group_id) const; 209 BoundingBox bounding_box(bool pos_side) const; 210 211 double y0_; 212 }; 213 214 //============================================================================== 215 //! A plane perpendicular to the z-axis. 216 // 217 //! The plane is described by the equation \f$z - z_0 = 0\f$ 218 //============================================================================== 219 220 class SurfaceZPlane : public CSGSurface 221 { 222 public: 223 explicit SurfaceZPlane(pugi::xml_node surf_node); 224 double evaluate(Position r) const; 225 double distance(Position r, Direction u, bool coincident) const; 226 Direction normal(Position r) const; 227 void to_hdf5_inner(hid_t group_id) const; 228 BoundingBox bounding_box(bool pos_side) const; 229 230 double z0_; 231 }; 232 233 //============================================================================== 234 //! A general plane. 235 // 236 //! The plane is described by the equation \f$A x + B y + C z - D = 0\f$ 237 //============================================================================== 238 239 class SurfacePlane : public CSGSurface 240 { 241 public: 242 explicit SurfacePlane(pugi::xml_node surf_node); 243 double evaluate(Position r) const; 244 double distance(Position r, Direction u, bool coincident) const; 245 Direction normal(Position r) const; 246 void to_hdf5_inner(hid_t group_id) const; 247 248 double A_, B_, C_, D_; 249 }; 250 251 //============================================================================== 252 //! A cylinder aligned along the x-axis. 253 // 254 //! The cylinder is described by the equation 255 //! \f$(y - y_0)^2 + (z - z_0)^2 - R^2 = 0\f$ 256 //============================================================================== 257 258 class SurfaceXCylinder : public CSGSurface 259 { 260 public: 261 explicit SurfaceXCylinder(pugi::xml_node surf_node); 262 double evaluate(Position r) const; 263 double distance(Position r, Direction u, bool coincident) const; 264 Direction normal(Position r) const; 265 void to_hdf5_inner(hid_t group_id) const; 266 BoundingBox bounding_box(bool pos_side) const; 267 268 double y0_, z0_, radius_; 269 }; 270 271 //============================================================================== 272 //! A cylinder aligned along the y-axis. 273 // 274 //! The cylinder is described by the equation 275 //! \f$(x - x_0)^2 + (z - z_0)^2 - R^2 = 0\f$ 276 //============================================================================== 277 278 class SurfaceYCylinder : public CSGSurface 279 { 280 public: 281 explicit SurfaceYCylinder(pugi::xml_node surf_node); 282 double evaluate(Position r) const; 283 double distance(Position r, Direction u, bool coincident) const; 284 Direction normal(Position r) const; 285 void to_hdf5_inner(hid_t group_id) const; 286 BoundingBox bounding_box(bool pos_side) const; 287 288 double x0_, z0_, radius_; 289 }; 290 291 //============================================================================== 292 //! A cylinder aligned along the z-axis. 293 // 294 //! The cylinder is described by the equation 295 //! \f$(x - x_0)^2 + (y - y_0)^2 - R^2 = 0\f$ 296 //============================================================================== 297 298 class SurfaceZCylinder : public CSGSurface 299 { 300 public: 301 explicit SurfaceZCylinder(pugi::xml_node surf_node); 302 double evaluate(Position r) const; 303 double distance(Position r, Direction u, bool coincident) const; 304 Direction normal(Position r) const; 305 void to_hdf5_inner(hid_t group_id) const; 306 BoundingBox bounding_box(bool pos_side) const; 307 308 double x0_, y0_, radius_; 309 }; 310 311 //============================================================================== 312 //! A sphere. 313 // 314 //! The cylinder is described by the equation 315 //! \f$(x - x_0)^2 + (y - y_0)^2 + (z - z_0)^2 - R^2 = 0\f$ 316 //============================================================================== 317 318 class SurfaceSphere : public CSGSurface 319 { 320 public: 321 explicit SurfaceSphere(pugi::xml_node surf_node); 322 double evaluate(Position r) const; 323 double distance(Position r, Direction u, bool coincident) const; 324 Direction normal(Position r) const; 325 void to_hdf5_inner(hid_t group_id) const; 326 BoundingBox bounding_box(bool pos_side) const; 327 328 double x0_, y0_, z0_, radius_; 329 }; 330 331 //============================================================================== 332 //! A cone aligned along the x-axis. 333 // 334 //! The cylinder is described by the equation 335 //! \f$(y - y_0)^2 + (z - z_0)^2 - R^2 (x - x_0)^2 = 0\f$ 336 //============================================================================== 337 338 class SurfaceXCone : public CSGSurface 339 { 340 public: 341 explicit SurfaceXCone(pugi::xml_node surf_node); 342 double evaluate(Position r) const; 343 double distance(Position r, Direction u, bool coincident) const; 344 Direction normal(Position r) const; 345 void to_hdf5_inner(hid_t group_id) const; 346 347 double x0_, y0_, z0_, radius_sq_; 348 }; 349 350 //============================================================================== 351 //! A cone aligned along the y-axis. 352 // 353 //! The cylinder is described by the equation 354 //! \f$(x - x_0)^2 + (z - z_0)^2 - R^2 (y - y_0)^2 = 0\f$ 355 //============================================================================== 356 357 class SurfaceYCone : public CSGSurface 358 { 359 public: 360 explicit SurfaceYCone(pugi::xml_node surf_node); 361 double evaluate(Position r) const; 362 double distance(Position r, Direction u, bool coincident) const; 363 Direction normal(Position r) const; 364 void to_hdf5_inner(hid_t group_id) const; 365 366 double x0_, y0_, z0_, radius_sq_; 367 }; 368 369 //============================================================================== 370 //! A cone aligned along the z-axis. 371 // 372 //! The cylinder is described by the equation 373 //! \f$(x - x_0)^2 + (y - y_0)^2 - R^2 (z - z_0)^2 = 0\f$ 374 //============================================================================== 375 376 class SurfaceZCone : public CSGSurface 377 { 378 public: 379 explicit SurfaceZCone(pugi::xml_node surf_node); 380 double evaluate(Position r) const; 381 double distance(Position r, Direction u, bool coincident) const; 382 Direction normal(Position r) const; 383 void to_hdf5_inner(hid_t group_id) const; 384 385 double x0_, y0_, z0_, radius_sq_; 386 }; 387 388 //============================================================================== 389 //! A general surface described by a quadratic equation. 390 // 391 //! \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$ 392 //============================================================================== 393 394 class SurfaceQuadric : public CSGSurface 395 { 396 public: 397 explicit SurfaceQuadric(pugi::xml_node surf_node); 398 double evaluate(Position r) const; 399 double distance(Position r, Direction u, bool coincident) const; 400 Direction normal(Position r) const; 401 void to_hdf5_inner(hid_t group_id) const; 402 403 // Ax^2 + By^2 + Cz^2 + Dxy + Eyz + Fxz + Gx + Hy + Jz + K = 0 404 double A_, B_, C_, D_, E_, F_, G_, H_, J_, K_; 405 }; 406 407 //============================================================================== 408 // Non-member functions 409 //============================================================================== 410 411 void read_surfaces(pugi::xml_node node); 412 413 void free_memory_surfaces(); 414 415 } // namespace openmc 416 #endif // OPENMC_SURFACE_H 417