From cddf39b201eb152f89a3ebed536fd36ba1797c39 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Mon, 20 Aug 2018 19:52:07 -0400 Subject: [PATCH] Move FILL_UNIVERSE chunk of find_cell to C++ --- src/cell.cpp | 71 +++++++++++++++++++++++++++++++++++++++++++++++ src/cell.h | 10 +++++++ src/geometry.F90 | 22 --------------- src/geometry.cpp | 48 ++++++++++++++++++++++++++++++-- src/input_xml.F90 | 22 --------------- src/position.h | 19 ++++++++++++- src/summary.F90 | 11 -------- 7 files changed, 144 insertions(+), 59 deletions(-) diff --git a/src/cell.cpp b/src/cell.cpp index a68ab9d0b..9786a89fc 100644 --- a/src/cell.cpp +++ b/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(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(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 rot {rotation[0], rotation[1], rotation[2]}; + write_dataset(cell_group, "rotation", rot); + } + } } //============================================================================== diff --git a/src/cell.h b/src/cell.h index 8b7c3a78b..4f3a8978f 100644 --- a/src/cell.h +++ b/src/cell.h @@ -85,6 +85,16 @@ public: std::vector 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 rotation; + std::vector offset; //!< Distribcell offset table Cell() {}; diff --git a/src/geometry.F90 b/src/geometry.F90 index 29e03bbb5..817da18de 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -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 diff --git a/src/geometry.cpp b/src/geometry.cpp index 7131a4513..695afe0be 100644 --- a/src/geometry.cpp +++ b/src/geometry.cpp @@ -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; + } } } diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 7d8a8393d..27b8be946 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -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) diff --git a/src/position.h b/src/position.h index 808e0bcf6..62af3a122 100644 --- a/src/position.h +++ b/src/position.h @@ -1,6 +1,8 @@ #ifndef OPENMC_POSITION_H #define OPENMC_POSITION_H +#include + 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 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 \ No newline at end of file +#endif // OPENMC_POSITION_H diff --git a/src/summary.F90 b/src/summary.F90 index 6c311e982..e2cc5f786 100644 --- a/src/summary.F90 +++ b/src/summary.F90 @@ -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