mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-26 21:25:36 -04:00
Move FILL_UNIVERSE chunk of find_cell to C++
This commit is contained in:
Sterling Harper
parent
2ed0fceb7d
commit
cddf39b201
7 changed files with 144 additions and 59 deletions
71
src/cell.cpp
71
src/cell.cpp
|
|
@ -317,6 +317,65 @@ Cell::Cell(pugi::xml_node cell_node)
|
|||
break;
|
||||
}
|
||||
}
|
||||
|
||||
// Read the translation vector.
|
||||
if (check_for_node(cell_node, "translation")) {
|
||||
if (fill == C_NONE) {
|
||||
std::stringstream err_msg;
|
||||
err_msg << "Cannot apply a translation to cell " << id
|
||||
<< " because it is not filled with another universe";
|
||||
fatal_error(err_msg);
|
||||
}
|
||||
|
||||
auto xyz {get_node_array<double>(cell_node, "translation")};
|
||||
if (xyz.size() != 3) {
|
||||
std::stringstream err_msg;
|
||||
err_msg << "Non-3D translation vector applied to cell " << id;
|
||||
fatal_error(err_msg);
|
||||
}
|
||||
translation = xyz;
|
||||
}
|
||||
|
||||
// Read the rotation transform.
|
||||
if (check_for_node(cell_node, "rotation")) {
|
||||
if (fill == C_NONE) {
|
||||
std::stringstream err_msg;
|
||||
err_msg << "Cannot apply a rotation to cell " << id
|
||||
<< " because it is not filled with another universe";
|
||||
fatal_error(err_msg);
|
||||
}
|
||||
|
||||
auto rot {get_node_array<double>(cell_node, "rotation")};
|
||||
if (rot.size() != 3) {
|
||||
std::stringstream err_msg;
|
||||
err_msg << "Non-3D rotation vector applied to cell " << id;
|
||||
fatal_error(err_msg);
|
||||
}
|
||||
|
||||
// Store the rotation angles.
|
||||
rotation.reserve(12);
|
||||
rotation.push_back(rot[0]);
|
||||
rotation.push_back(rot[1]);
|
||||
rotation.push_back(rot[2]);
|
||||
|
||||
// Compute and store the rotation matrix.
|
||||
auto phi = -rot[0] * PI / 180.0;
|
||||
auto theta = -rot[1] * PI / 180.0;
|
||||
auto psi = -rot[2] * PI / 180.0;
|
||||
rotation.push_back(std::cos(theta) * std::cos(psi));
|
||||
rotation.push_back(-std::cos(phi) * std::sin(psi)
|
||||
+ std::sin(phi) * std::sin(theta) * std::cos(psi));
|
||||
rotation.push_back(std::sin(phi) * std::sin(psi)
|
||||
+ std::cos(phi) * std::sin(theta) * std::cos(psi));
|
||||
rotation.push_back(std::cos(theta) * std::sin(psi));
|
||||
rotation.push_back(std::cos(phi) * std::cos(psi)
|
||||
+ std::sin(phi) * std::sin(theta) * std::sin(psi));
|
||||
rotation.push_back(-std::sin(phi) * std::cos(psi)
|
||||
+ std::cos(phi) * std::sin(theta) * std::sin(psi));
|
||||
rotation.push_back(-std::sin(theta));
|
||||
rotation.push_back(std::sin(phi) * std::cos(theta));
|
||||
rotation.push_back(std::cos(phi) * std::cos(theta));
|
||||
}
|
||||
}
|
||||
|
||||
//==============================================================================
|
||||
|
|
@ -393,6 +452,18 @@ Cell::to_hdf5(hid_t cell_group) const
|
|||
}
|
||||
write_string(cell_group, "region", region_spec.str(), false);
|
||||
}
|
||||
|
||||
if (type == FILL_UNIVERSE) {
|
||||
write_dataset(cell_group, "fill_type", "universe");
|
||||
write_dataset(cell_group, "fill", global_universes[fill]->id);
|
||||
if (translation != 0) {
|
||||
write_dataset(cell_group, "translation", translation);
|
||||
}
|
||||
if (!rotation.empty()) {
|
||||
std::array<double, 3> rot {rotation[0], rotation[1], rotation[2]};
|
||||
write_dataset(cell_group, "rotation", rot);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
//==============================================================================
|
||||
|
|
|
|||
10
src/cell.h
10
src/cell.h
|
|
@ -85,6 +85,16 @@ public:
|
|||
std::vector<std::int32_t> rpn;
|
||||
bool simple; //!< Does the region contain only intersections?
|
||||
|
||||
Position translation {0, 0, 0}; //!< Translation vector for filled universe
|
||||
|
||||
//! \brief Rotational tranfsormation of the filled universe.
|
||||
//
|
||||
//! The vector is empty if there is no rotation. Otherwise, the first three
|
||||
//! values are the rotation angles respectively about the x-, y-, and z-, axes
|
||||
//! in degrees. The next 9 values give the rotation matrix in row-major
|
||||
//! order.
|
||||
std::vector<double> rotation;
|
||||
|
||||
std::vector<int32_t> offset; //!< Distribcell offset table
|
||||
|
||||
Cell() {};
|
||||
|
|
|
|||
|
|
@ -84,30 +84,8 @@ contains
|
|||
if (found) then
|
||||
associate(c => cells(i_cell))
|
||||
CELL_TYPE: if (c % type() == FILL_UNIVERSE) then
|
||||
! ======================================================================
|
||||
! CELL CONTAINS LOWER UNIVERSE, RECURSIVELY FIND CELL
|
||||
|
||||
! Store lower level coordinates
|
||||
p % coord(j + 1) % xyz = p % coord(j) % xyz
|
||||
p % coord(j + 1) % uvw = p % coord(j) % uvw
|
||||
|
||||
! Move particle to next level and set universe
|
||||
j = j + 1
|
||||
p % n_coord = j
|
||||
p % coord(j) % universe = c % fill() + 1
|
||||
|
||||
! Apply translation
|
||||
if (allocated(c % translation)) then
|
||||
p % coord(j) % xyz = p % coord(j) % xyz - c % translation
|
||||
end if
|
||||
|
||||
! Apply rotation
|
||||
if (allocated(c % rotation_matrix)) then
|
||||
p % coord(j) % xyz = matmul(c % rotation_matrix, p % coord(j) % xyz)
|
||||
p % coord(j) % uvw = matmul(c % rotation_matrix, p % coord(j) % uvw)
|
||||
p % coord(j) % rotated = .true.
|
||||
end if
|
||||
|
||||
call find_cell(p, found)
|
||||
j = p % n_coord
|
||||
|
||||
|
|
|
|||
|
|
@ -102,11 +102,11 @@ find_cell(Particle* p, int n_search_cells, int* search_cells) {
|
|||
if (found) {
|
||||
Cell& c {*global_cells[i_cell]};
|
||||
if (c.type == FILL_MATERIAL) {
|
||||
//=======================================================================
|
||||
//! Found a material cell which means this is the lowest coord level.
|
||||
|
||||
// Find the distribcell instance number.
|
||||
if (c.material.size() > 1 || c.sqrtkT.size() > 1) {
|
||||
//=====================================================================
|
||||
//! Found a material cell which means this is the lowest coord level.
|
||||
|
||||
//TODO: off-by-one indexing
|
||||
int distribcell_index = c.distribcell_index - 1;
|
||||
int offset = 0;
|
||||
|
|
@ -148,6 +148,48 @@ find_cell(Particle* p, int n_search_cells, int* search_cells) {
|
|||
} else {
|
||||
p->sqrtkT = c.sqrtkT[0];
|
||||
}
|
||||
|
||||
} else if (c.type == FILL_UNIVERSE) {
|
||||
//========================================================================
|
||||
//! Found a lower universe, update this coord level then search the next.
|
||||
|
||||
// Add another coordinate level.
|
||||
//++p->n_coord;
|
||||
p->coord[p->n_coord].universe = c.fill + 1;
|
||||
|
||||
// Set the position and direction.
|
||||
for (int i = 0; i < 3; i++) {
|
||||
p->coord[p->n_coord].xyz[i] = p->coord[p->n_coord-1].xyz[i];
|
||||
p->coord[p->n_coord].uvw[i] = p->coord[p->n_coord-1].uvw[i];
|
||||
}
|
||||
|
||||
// Apply translation.
|
||||
p->coord[p->n_coord].xyz[0] -= c.translation.x;
|
||||
p->coord[p->n_coord].xyz[1] -= c.translation.y;
|
||||
p->coord[p->n_coord].xyz[2] -= c.translation.z;
|
||||
|
||||
// Apply rotation.
|
||||
if (!c.rotation.empty()) {
|
||||
auto x = p->coord[p->n_coord].xyz[0];
|
||||
auto y = p->coord[p->n_coord].xyz[1];
|
||||
auto z = p->coord[p->n_coord].xyz[2];
|
||||
p->coord[p->n_coord].xyz[0] = x*c.rotation[3] + y*c.rotation[4]
|
||||
+ z*c.rotation[5];
|
||||
p->coord[p->n_coord].xyz[1] = x*c.rotation[6] + y*c.rotation[7]
|
||||
+ z*c.rotation[8];
|
||||
p->coord[p->n_coord].xyz[2] = x*c.rotation[9] + y*c.rotation[10]
|
||||
+ z*c.rotation[11];
|
||||
auto u = p->coord[p->n_coord].uvw[0];
|
||||
auto v = p->coord[p->n_coord].uvw[1];
|
||||
auto w = p->coord[p->n_coord].uvw[2];
|
||||
p->coord[p->n_coord].uvw[0] = u*c.rotation[3] + v*c.rotation[4]
|
||||
+ w*c.rotation[5];
|
||||
p->coord[p->n_coord].uvw[1] = u*c.rotation[6] + v*c.rotation[7]
|
||||
+ w*c.rotation[8];
|
||||
p->coord[p->n_coord].uvw[2] = u*c.rotation[9] + v*c.rotation[10]
|
||||
+ w*c.rotation[11];
|
||||
p->coord[p->n_coord].rotated = true;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -1176,28 +1176,6 @@ contains
|
|||
cos(phi)*cos(theta) /), (/ 3,3 /))
|
||||
end if
|
||||
|
||||
! Translation vector
|
||||
if (check_for_node(node_cell, "translation")) then
|
||||
! Translations can only be applied to cells that are being filled with
|
||||
! another universe
|
||||
if (c % fill() == C_NONE) then
|
||||
call fatal_error("Cannot apply a translation to cell " &
|
||||
// trim(to_str(c % id())) // " because it is not filled with &
|
||||
&another universe")
|
||||
end if
|
||||
|
||||
! Read number of translation parameters
|
||||
n = node_word_count(node_cell, "translation")
|
||||
if (n /= 3) then
|
||||
call fatal_error("Incorrect number of translation parameters on &
|
||||
&cell " // to_str(c % id()))
|
||||
end if
|
||||
|
||||
! Copy translation vector
|
||||
allocate(c % translation(3))
|
||||
call get_node_array(node_cell, "translation", c % translation)
|
||||
end if
|
||||
|
||||
! Add cell to dictionary
|
||||
call cell_dict % set(c % id(), i)
|
||||
|
||||
|
|
|
|||
|
|
@ -1,6 +1,8 @@
|
|||
#ifndef OPENMC_POSITION_H
|
||||
#define OPENMC_POSITION_H
|
||||
|
||||
#include <vector>
|
||||
|
||||
namespace openmc {
|
||||
|
||||
//==============================================================================
|
||||
|
|
@ -12,6 +14,7 @@ struct Position {
|
|||
Position() = default;
|
||||
Position(double x_, double y_, double z_) : x{x_}, y{y_}, z{z_} { };
|
||||
Position(const double xyz[]) : x{xyz[0]}, y{xyz[1]}, z{xyz[2]} { };
|
||||
Position(const std::vector<double> xyz) : x{xyz[0]}, y{xyz[1]}, z{xyz[2]} { };
|
||||
|
||||
// Unary operators
|
||||
Position& operator+=(Position);
|
||||
|
|
@ -63,6 +66,20 @@ inline Position operator*(Position a, Position b) { return a *= b; }
|
|||
inline Position operator*(Position a, double b) { return a *= b; }
|
||||
inline Position operator*(double a, Position b) { return b *= a; }
|
||||
|
||||
inline bool operator==(Position a, Position b)
|
||||
{return a.x == b.x && a.y == b.y && a.z == b.z;}
|
||||
inline bool operator==(Position a, double b)
|
||||
{return a.x == b && a.y == b && a.z == b;}
|
||||
inline bool operator==(double a, Position b)
|
||||
{return a == b.x && a == b.y && a == b.z;}
|
||||
|
||||
inline bool operator!=(Position a, Position b)
|
||||
{return a.x != b.x || a.y != b.y || a.z != b.z;}
|
||||
inline bool operator!=(Position a, double b)
|
||||
{return a.x != b || a.y != b || a.z != b;}
|
||||
inline bool operator!=(double a, Position b)
|
||||
{return a != b.x || a != b.y || a != b.z;}
|
||||
|
||||
//==============================================================================
|
||||
//! Type representing a vector direction in Cartesian coordinates
|
||||
//==============================================================================
|
||||
|
|
@ -71,4 +88,4 @@ using Direction = Position;
|
|||
|
||||
} // namespace openmc
|
||||
|
||||
#endif // OPENMC_POSITION_H
|
||||
#endif // OPENMC_POSITION_H
|
||||
|
|
|
|||
|
|
@ -225,17 +225,6 @@ contains
|
|||
call write_dataset(cell_group, "temperature", cell_temperatures)
|
||||
deallocate(cell_temperatures)
|
||||
|
||||
case (FILL_UNIVERSE)
|
||||
call write_dataset(cell_group, "fill_type", "universe")
|
||||
call write_dataset(cell_group, "fill", universes(c%fill()+1)%id)
|
||||
|
||||
if (allocated(c%translation)) then
|
||||
call write_dataset(cell_group, "translation", c%translation)
|
||||
end if
|
||||
if (allocated(c%rotation)) then
|
||||
call write_dataset(cell_group, "rotation", c%rotation)
|
||||
end if
|
||||
|
||||
case (FILL_LATTICE)
|
||||
call write_dataset(cell_group, "fill_type", "lattice")
|
||||
! Do not access the 'lattices' array with 'c % fill() + 1' directly; it
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue