diff --git a/include/openmc/boundary_condition.h b/include/openmc/boundary_condition.h new file mode 100644 index 0000000000..cbc1575aac --- /dev/null +++ b/include/openmc/boundary_condition.h @@ -0,0 +1,83 @@ +#ifndef OPENMC_BOUNDARY_CONDITION_H +#define OPENMC_BOUNDARY_CONDITION_H + +#include "openmc/position.h" + +namespace openmc { + +//============================================================================== +// Global variables +//============================================================================== + +//class BoundaryCondition; +// +//namespace model { +// extern std::vector> boundary_conditions; +//} // namespace model + +class Particle; +class Surface; + +class BoundaryCondition +{ +public: + virtual void + handle_particle(Particle& p, const Surface& surf) const = 0; +}; + +class VacuumBC : public BoundaryCondition +{ +public: + void + handle_particle(Particle& p, const Surface& surf) const override; +}; + +class ReflectiveBC : public BoundaryCondition +{ +public: + void + handle_particle(Particle& p, const Surface& surf) const override; +}; + +class WhiteBC : public BoundaryCondition +{ +public: + void + handle_particle(Particle& p, const Surface& surf) const override; +}; + +class PeriodicBC : public BoundaryCondition +{ +public: + PeriodicBC(int i_surf, int j_surf) + : i_surf_(i_surf), j_surf_(j_surf) + {}; + +protected: + int i_surf_; + int j_surf_; +}; + +class TranslationalPeriodicBC : public PeriodicBC +{ +public: + TranslationalPeriodicBC(int i_surf, int j_surf); + + void + handle_particle(Particle& p, const Surface& surf) const override; + +protected: + Position translation_; +}; + +class RotationalPeriodicBC : public PeriodicBC +{ +public: + RotationalPeriodicBC(int i_surf, int j_surf); + + void + handle_particle(Particle& p, const Surface& surf) const override; +}; + +} // namespace openmc +#endif // OPENMC_BOUNDARY_CONDITION_H diff --git a/src/boundary_condition.cpp b/src/boundary_condition.cpp index d577f6bbf9..1c19eba595 100644 --- a/src/boundary_condition.cpp +++ b/src/boundary_condition.cpp @@ -1,4 +1,9 @@ #include "openmc/boundary_condition.h" + +#include + +#include + #include "openmc/surface.h" namespace openmc { @@ -27,8 +32,82 @@ WhiteBC::handle_particle(Particle& p, const Surface& surf) const p.cross_reflective_bc(surf, u); } +TranslationalPeriodicBC::TranslationalPeriodicBC(int i_surf, int j_surf) + : PeriodicBC(i_surf, j_surf) +{ + Surface& surf1 {*model::surfaces[i_surf_]}; + Surface& surf2 {*model::surfaces[j_surf_]}; + + // The following blocks will resolve the type of each surface and compute the + // appropriate translation vector + + // Check for a pair of x-planes + if (const auto* xplane1 = dynamic_cast(&surf1)) { + if (const auto* xplane2 = dynamic_cast(&surf2)) { + translation_ = {xplane2->x0_ - xplane1->x0_, 0, 0}; + } else { + throw std::invalid_argument(fmt::format("Invalid pair of periodic " + "surfaces ({} and {}). For a translational periodic BC, both surfaces " + "must be of the same type (e.g. both x-plane).", surf1.id_, surf2.id_)); + } + + // Check for a pair of y-planes + } else if (const auto* yplane1 = dynamic_cast(&surf1)) { + if (const auto* yplane2 = dynamic_cast(&surf2)) { + translation_ = {0, yplane2->y0_ - yplane1->y0_, 0}; + } else { + throw std::invalid_argument(fmt::format("Invalid pair of periodic " + "surfaces ({} and {}). For a translational periodic BC, both surfaces " + "must be of the same type (e.g. both x-plane).", surf1.id_, surf2.id_)); + } + + // Check for a pair of z-planes + } else if (const auto* zplane1 = dynamic_cast(&surf1)) { + if (const auto* zplane2 = dynamic_cast(&surf2)) { + translation_ = {0, 0, zplane2->z0_ - zplane1->z0_}; + } else { + throw std::invalid_argument(fmt::format("Invalid pair of periodic " + "surfaces ({} and {}). For a translational periodic BC, both surfaces " + "must be of the same type (e.g. both x-plane).", surf1.id_, surf2.id_)); + } + } + + // TODO: Check for a pair of general planes + +} + void -PeriodicBC::handle_particle(Particle& p, const Surface& surf) const +TranslationalPeriodicBC::handle_particle(Particle& p, const Surface& surf) const +{ + // TODO: off-by-one on surface indices throughout this function. + int i_particle_surf = std::abs(p.surface_) - 1; + + // Figure out which of the two BC surfaces were struck then find the + // particle's new location and surface. + Position new_r; + int new_surface; + if (i_particle_surf == i_surf_) { + new_r = p.r() + translation_; + new_surface = p.surface_ > 0 ? j_surf_ + 1 : -(j_surf_ + 1); + } else if (i_particle_surf == j_surf_) { + new_r = p.r() - translation_; + new_surface = p.surface_ > 0 ? i_surf_ + 1 : -(i_surf_ + 1); + } else { + throw std::runtime_error("Called BoundaryCondition::handle_particle after " + "hitting a surface, but that surface is not recognized by the BC."); + } + + // Pass the new location and surface to the particle. + p.cross_periodic_bc(surf, new_r, p.u(), new_surface); +} + +RotationalPeriodicBC::RotationalPeriodicBC(int i_surf, int j_surf) + : PeriodicBC(i_surf, j_surf) +{ +} + +void +RotationalPeriodicBC::handle_particle(Particle& p, const Surface& surf) const { // Get a pointer of the first surface downcast to PeriodicSurface auto surf_p = dynamic_cast(&surf); @@ -38,8 +117,6 @@ PeriodicBC::handle_particle(Particle& p, const Surface& surf) const surf.id_ == model::surfaces[i_surf_]->id_ ? dynamic_cast(model::surfaces[j_surf_].get()) : dynamic_cast(model::surfaces[i_surf_].get()); - //auto other = dynamic_cast( - // model::surfaces[surf_p->i_periodic_].get()); // Compute the new particle location and direction Position r {p.r()}; diff --git a/src/surface.cpp b/src/surface.cpp index 3006ea34ba..5540cf860a 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -1249,8 +1249,19 @@ void read_surfaces(pugi::xml_node node) int j_surf = model::surface_map[periodic_pair.second]; Surface& surf1 {*model::surfaces[i_surf]}; Surface& surf2 {*model::surfaces[j_surf]}; - surf1.new_bc_ = std::make_shared(i_surf, j_surf); - surf2.new_bc_ = surf1.new_bc_; + + Direction norm1 = surf1.normal({0, 0, 0}); + Direction norm2 = surf2.normal({0, 0, 0}); + double dot_prod = norm1.dot(norm2); + + if (std::abs(1.0 - dot_prod) < FP_PRECISION) { + surf1.new_bc_ = std::make_shared(i_surf, j_surf); + surf2.new_bc_ = surf1.new_bc_; + } else { + surf1.new_bc_ = std::make_shared(i_surf, j_surf); + surf2.new_bc_ = surf1.new_bc_; + } + dynamic_cast(&surf1)->i_periodic_ = j_surf; dynamic_cast(&surf2)->i_periodic_ = i_surf; }