Skip to content

Commit 8e6842d

Browse files
smharpersmharper
committed
Move rotational periodic BCs to the new class
Remove the now unnecessary Surface::periodic_translate function.
1 parent 16d4c58 commit 8e6842d

4 files changed

Lines changed: 51 additions & 104 deletions

File tree

‎include/openmc/boundary_condition.h‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -77,6 +77,9 @@ class RotationalPeriodicBC : public PeriodicBC
7777

7878
void
7979
handle_particle(Particle& p, const Surface& surf) const override;
80+
81+
protected:
82+
double angle_;
8083
};
8184

8285
} // namespace openmc

‎include/openmc/surface.h‎

Lines changed: 0 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -196,18 +196,6 @@ class PeriodicSurface : public CSGSurface
196196
int i_periodic_{C_NONE}; //!< Index of corresponding periodic surface
197197

198198
explicit PeriodicSurface(pugi::xml_node surf_node);
199-
200-
//! Translate a particle onto this surface from a periodic partner surface.
201-
//! \param other A pointer to the partner surface in this periodic BC.
202-
//! \param r A point on the partner surface that will be translated onto
203-
//! this surface.
204-
//! \param u A direction that will be rotated for systems with rotational
205-
//! periodicity.
206-
//! \return true if this surface and its partner make a rotationally-periodic
207-
//! boundary condition.
208-
virtual bool periodic_translate(const PeriodicSurface* other, Position& r,
209-
Direction& u) const = 0;
210-
211199
};
212200

213201
//==============================================================================

‎src/boundary_condition.cpp‎

Lines changed: 48 additions & 21 deletions
Original file line numberDiff line numberDiff line change
@@ -4,6 +4,7 @@
44

55
#include <fmt/core.h>
66

7+
#include "openmc/constants.h"
78
#include "openmc/surface.h"
89

910
namespace openmc {
@@ -104,32 +105,58 @@ TranslationalPeriodicBC::handle_particle(Particle& p, const Surface& surf) const
104105
RotationalPeriodicBC::RotationalPeriodicBC(int i_surf, int j_surf)
105106
: PeriodicBC(i_surf, j_surf)
106107
{
108+
Surface& surf1 {*model::surfaces[i_surf_]};
109+
Surface& surf2 {*model::surfaces[j_surf_]};
110+
111+
// The following blocks will resolve the type of each surface and compute the
112+
// appropriate rotation angle
113+
114+
if (const auto* xplane = dynamic_cast<const SurfaceXPlane*>(&surf1)) {
115+
if (const auto* yplane = dynamic_cast<const SurfaceYPlane*>(&surf2)) {
116+
angle_ = 0.5 * PI;
117+
}
118+
} else if (const auto* yplane = dynamic_cast<const SurfaceYPlane*>(&surf1)) {
119+
if (const auto* xplane = dynamic_cast<const SurfaceXPlane*>(&surf2)) {
120+
angle_ = -0.5 * PI;
121+
}
122+
}
107123
}
108124

109125
void
110126
RotationalPeriodicBC::handle_particle(Particle& p, const Surface& surf) const
111127
{
112-
// Get a pointer of the first surface downcast to PeriodicSurface
113-
auto surf_p = dynamic_cast<const PeriodicSurface*>(&surf);
114-
115-
// Get a pointer to the partner PeriodicSurface
116-
auto other =
117-
surf.id_ == model::surfaces[i_surf_]->id_ ?
118-
dynamic_cast<const PeriodicSurface*>(model::surfaces[j_surf_].get()) :
119-
dynamic_cast<const PeriodicSurface*>(model::surfaces[i_surf_].get());
120-
121-
// Compute the new particle location and direction
122-
Position r {p.r()};
123-
Direction u {p.u()};
124-
bool rotational = other->periodic_translate(surf_p, r, u);
125-
126-
// Pick the particle's new surface
127-
// TODO: off-by-one
128-
int new_surface = rotational ?
129-
surf_p->i_periodic_ + 1 :
130-
((p.surface_ > 0) ? surf_p->i_periodic_ + 1 : -(surf_p->i_periodic_ + 1));
131-
132-
p.cross_periodic_bc(surf, r, u, new_surface);
128+
// TODO: off-by-one on surface indices throughout this function.
129+
int i_particle_surf = std::abs(p.surface_) - 1;
130+
131+
// Figure out which of the two BC surfaces were struck to figure out if a
132+
// forward or backward rotation is required. Also specify the new surface.
133+
double theta;
134+
int new_surface;
135+
if (i_particle_surf == i_surf_) {
136+
theta = -angle_;
137+
new_surface = p.surface_ > 0 ? -(j_surf_ + 1) : j_surf_ + 1;
138+
} else if (i_particle_surf == j_surf_) {
139+
theta = angle_;
140+
new_surface = p.surface_ > 0 ? -(i_surf_ + 1) : i_surf_ + 1;
141+
} else {
142+
throw std::runtime_error("Called BoundaryCondition::handle_particle after "
143+
"hitting a surface, but that surface is not recognized by the BC.");
144+
}
145+
146+
// Rotate the particle's position and direction about the z-axis.
147+
Position r = p.r();
148+
Direction u = p.u();
149+
Position new_r = {
150+
std::cos(theta)*r[0] - std::sin(theta)*r[1],
151+
std::sin(theta)*r[0] + std::cos(theta)*r[1],
152+
r[2]};
153+
Direction new_u = {
154+
std::cos(theta)*u[0] - std::sin(theta)*u[1],
155+
std::sin(theta)*u[0] + std::cos(theta)*u[1],
156+
u[2]};
157+
158+
// Pass the new location, direction, and surface to the particle.
159+
p.cross_periodic_bc(surf, new_r, new_u, new_surface);
133160
}
134161

135162
} // namespace openmc

‎src/surface.cpp‎

Lines changed: 0 additions & 71 deletions
Original file line numberDiff line numberDiff line change
@@ -352,29 +352,6 @@ void SurfaceXPlane::to_hdf5_inner(hid_t group_id) const
352352
write_dataset(group_id, "coefficients", coeffs);
353353
}
354354

355-
bool SurfaceXPlane::periodic_translate(const PeriodicSurface* other,
356-
Position& r, Direction& u) const
357-
{
358-
Direction other_n = other->normal(r);
359-
if (other_n.x == 1 && other_n.y == 0 && other_n.z == 0) {
360-
r.x = x0_;
361-
return false;
362-
} else {
363-
// Assume the partner is an YPlane (the only supported partner). Use the
364-
// evaluate function to find y0, then adjust position/Direction for
365-
// rotational symmetry.
366-
double y0_ = -other->evaluate({0., 0., 0.});
367-
r.y = r.x - x0_ + y0_;
368-
r.x = x0_;
369-
370-
double ux = u.x;
371-
u.x = -u.y;
372-
u.y = ux;
373-
374-
return true;
375-
}
376-
}
377-
378355
BoundingBox
379356
SurfaceXPlane::bounding_box(bool pos_side) const
380357
{
@@ -417,30 +394,6 @@ void SurfaceYPlane::to_hdf5_inner(hid_t group_id) const
417394
write_dataset(group_id, "coefficients", coeffs);
418395
}
419396

420-
bool SurfaceYPlane::periodic_translate(const PeriodicSurface* other,
421-
Position& r, Direction& u) const
422-
{
423-
Direction other_n = other->normal(r);
424-
if (other_n.x == 0 && other_n.y == 1 && other_n.z == 0) {
425-
// The periodic partner is also aligned along y. Just change the y coord.
426-
r.y = y0_;
427-
return false;
428-
} else {
429-
// Assume the partner is an XPlane (the only supported partner). Use the
430-
// evaluate function to find x0, then adjust position/Direction for rotational
431-
// symmetry.
432-
double x0_ = -other->evaluate({0., 0., 0.});
433-
r.x = r.y - y0_ + x0_;
434-
r.y = y0_;
435-
436-
double ux = u.x;
437-
u.x = u.y;
438-
u.y = -ux;
439-
440-
return true;
441-
}
442-
}
443-
444397
BoundingBox
445398
SurfaceYPlane::bounding_box(bool pos_side) const
446399
{
@@ -483,14 +436,6 @@ void SurfaceZPlane::to_hdf5_inner(hid_t group_id) const
483436
write_dataset(group_id, "coefficients", coeffs);
484437
}
485438

486-
bool SurfaceZPlane::periodic_translate(const PeriodicSurface* other,
487-
Position& r, Direction& u) const
488-
{
489-
// Assume the other plane is aligned along z. Just change the z coord.
490-
r.z = z0_;
491-
return false;
492-
}
493-
494439
BoundingBox
495440
SurfaceZPlane::bounding_box(bool pos_side) const
496441
{
@@ -544,22 +489,6 @@ void SurfacePlane::to_hdf5_inner(hid_t group_id) const
544489
write_dataset(group_id, "coefficients", coeffs);
545490
}
546491

547-
bool SurfacePlane::periodic_translate(const PeriodicSurface* other, Position& r,
548-
Direction& u) const
549-
{
550-
// This function assumes the other plane shares this plane's normal direction.
551-
552-
// Determine the distance to intersection.
553-
double d = evaluate(r) / (A_*A_ + B_*B_ + C_*C_);
554-
555-
// Move the particle that distance along the normal vector.
556-
r.x -= d * A_;
557-
r.y -= d * B_;
558-
r.z -= d * C_;
559-
560-
return false;
561-
}
562-
563492
//==============================================================================
564493
// Generic functions for x-, y-, and z-, cylinders
565494
//==============================================================================

0 commit comments

Comments
 (0)