Move DAGCell implementation into the dagmc files.

This commit is contained in:
Patrick Shriwise 2021-04-01 23:24:06 -05:00
parent e0aed3d2b3
commit 230e054ec3
4 changed files with 111 additions and 114 deletions

View file

@ -260,28 +260,6 @@ protected:
vector<int32_t>::iterator start, const vector<int32_t>& rpn);
};
//==============================================================================
#ifdef DAGMC
class DAGCell : public Cell
{
public:
DAGCell();
bool contains(Position r, Direction u, int32_t on_surface) const;
std::pair<double, int32_t>
distance(Position r, Direction u, int32_t on_surface, Particle* p) const;
BoundingBox bounding_box() const;
void to_hdf5_inner(hid_t group_id) const;
std::shared_ptr<moab::DagMC> dagmc_ptr_; //!< Pointer to DagMC instance
int32_t dag_index_; //!< DagMC index of cell
};
#endif
//==============================================================================
//! Speeds up geometry searches by grouping cells in a search tree.
//
@ -351,7 +329,6 @@ void read_cells(pugi::xml_node node);
#ifdef DAGMC
class DAGUniverse;
int32_t next_cell(DAGUniverse* dag_univ, DAGCell* cur_cell, DAGSurface* surf_xed);
#endif
} // namespace openmc

View file

@ -10,13 +10,32 @@ extern "C" const bool DAGMC_ENABLED;
#include "DagMC.hpp"
#include "openmc/cell.h"
#include "openmc/particle.h"
#include "openmc/position.h"
#include "openmc/surface.h"
#include "openmc/xml_interface.h"
class UWUW;
namespace openmc {
class DAGCell : public Cell {
public:
DAGCell();
bool contains(Position r, Direction u, int32_t on_surface) const;
std::pair<double, int32_t>
distance(Position r, Direction u, int32_t on_surface, Particle* p) const;
BoundingBox bounding_box() const;
void to_hdf5_inner(hid_t group_id) const;
std::shared_ptr<moab::DagMC> dagmc_ptr_; //!< Pointer to DagMC instance
int32_t dag_index_; //!< DagMC index of cell
};
class DAGUniverse : public Universe {
public:
@ -26,7 +45,6 @@ public:
void initialize(); //!< Sets up the DAGMC instance and OpenMC internals
void read_uwuw_materials(); //!< Reads UWUW materials and returns an ID map
bool uses_uwuw() const;
int32_t implicit_complement_idx() const;
@ -52,7 +70,9 @@ public:
bool adjust_material_ids_;
};
// DAGMC Functions
void read_dagmc_universes(pugi::xml_node node);
int32_t next_cell(DAGUniverse* dag_univ, DAGCell* cur_cell, DAGSurface* surf_xed);
} // namespace openmc

View file

@ -825,81 +825,6 @@ CSGCell::contains_complex(Position r, Direction u, int32_t on_surface) const
}
}
//==============================================================================
// DAGMC Cell implementation
//==============================================================================
#ifdef DAGMC
DAGCell::DAGCell() : Cell{} {
geom_type_ = GeometryType::DAG;
simple_ = true;
};
std::pair<double, int32_t>
DAGCell::distance(Position r, Direction u, int32_t on_surface, Particle* p) const
{
Expects(p);
// if we've changed direction or we're not on a surface,
// reset the history and update last direction
if (u != p->last_dir()) { p->last_dir() = u; p->history().reset(); }
if (on_surface == 0) { p->history().reset(); }
const auto& univ = model::universes[p->coord_[p->n_coord_ - 1].universe];
DAGUniverse* dag_univ = static_cast<DAGUniverse*>(univ.get());
if (!dag_univ) fatal_error("DAGMC call made for particle in a non-DAGMC universe");
moab::ErrorCode rval;
moab::EntityHandle vol = dagmc_ptr_->entity_by_index(3, dag_index_);
moab::EntityHandle hit_surf;
double dist;
double pnt[3] = {r.x, r.y, r.z};
double dir[3] = {u.x, u.y, u.z};
rval = dagmc_ptr_->ray_fire(vol, pnt, dir, hit_surf, dist, &p->history());
MB_CHK_ERR_CONT(rval);
int surf_idx;
if (hit_surf != 0) {
surf_idx = dag_univ->surf_idx_offset_ + dagmc_ptr_->index_by_handle(hit_surf);
} else {
// indicate that particle is lost
surf_idx = -1;
dist = INFINITY;
if (!dagmc_ptr_->is_implicit_complement(vol) || model::universe_map[dag_univ->id_] == model::root_universe) {
p->mark_as_lost(fmt::format("No intersection found with DAGMC cell {}", id_));
}
}
return {dist, surf_idx};
}
bool DAGCell::contains(Position r, Direction u, int32_t on_surface) const
{
moab::ErrorCode rval;
moab::EntityHandle vol = dagmc_ptr_->entity_by_index(3, dag_index_);
int result = 0;
double pnt[3] = {r.x, r.y, r.z};
double dir[3] = {u.x, u.y, u.z};
rval = dagmc_ptr_->point_in_volume(vol, pnt, result, dir);
MB_CHK_ERR_CONT(rval);
return result;
}
void DAGCell::to_hdf5_inner(hid_t group_id) const {
write_string(group_id, "geom_type", "dagmc", false);
}
BoundingBox DAGCell::bounding_box() const
{
moab::ErrorCode rval;
moab::EntityHandle vol = dagmc_ptr_->entity_by_index(3, dag_index_);
double min[3], max[3];
rval = dagmc_ptr_->getobb(vol, min, max);
MB_CHK_ERR_CONT(rval);
return {min[0], max[0], min[1], max[1], min[2], max[2]};
}
#endif
//==============================================================================
// UniversePartitioner implementation
//==============================================================================
@ -1356,21 +1281,6 @@ openmc_extend_cells(int32_t n, int32_t* index_start, int32_t* index_end)
return 0;
}
#ifdef DAGMC
int32_t next_cell(DAGUniverse* dag_univ, DAGCell* cur_cell, DAGSurface* surf_xed)
{
moab::EntityHandle surf =
surf_xed->dagmc_ptr_->entity_by_index(2, surf_xed->dag_index_);
moab::EntityHandle vol =
cur_cell->dagmc_ptr_->entity_by_index(3, cur_cell->dag_index_);
moab::EntityHandle new_vol;
cur_cell->dagmc_ptr_->next_vol(surf, vol, new_vol);
return cur_cell->dagmc_ptr_->index_by_handle(new_vol) + dag_univ->cell_idx_offset_;
}
#endif
extern "C" int cells_size() { return model::cells.size(); }
} // namespace openmc

View file

@ -7,6 +7,7 @@
#include "openmc/file_utils.h"
#include "openmc/geometry.h"
#include "openmc/geometry_aux.h"
#include "openmc/hdf5_interface.h"
#include "openmc/material.h"
#include "openmc/string_utils.h"
#include "openmc/settings.h"
@ -489,5 +490,94 @@ DAGUniverse::read_uwuw_materials() {
}
}
//==============================================================================
// DAGMC Cell implementation
//==============================================================================
DAGCell::DAGCell() : Cell{} {
geom_type_ = GeometryType::DAG;
simple_ = true;
};
std::pair<double, int32_t>
DAGCell::distance(Position r, Direction u, int32_t on_surface, Particle* p) const
{
Expects(p);
// if we've changed direction or we're not on a surface,
// reset the history and update last direction
if (u != p->last_dir_) { p->last_dir_ = u; p->history_.reset(); }
if (on_surface == 0) { p->history_.reset(); }
const auto& univ = model::universes[p->coord_[p->n_coord_ - 1].universe];
DAGUniverse* dag_univ = static_cast<DAGUniverse*>(univ.get());
if (!dag_univ) fatal_error("DAGMC call made for particle in a non-DAGMC universe");
moab::ErrorCode rval;
moab::EntityHandle vol = dagmc_ptr_->entity_by_index(3, dag_index_);
moab::EntityHandle hit_surf;
double dist;
double pnt[3] = {r.x, r.y, r.z};
double dir[3] = {u.x, u.y, u.z};
rval = dagmc_ptr_->ray_fire(vol, pnt, dir, hit_surf, dist, &p->history_);
MB_CHK_ERR_CONT(rval);
int surf_idx;
if (hit_surf != 0) {
surf_idx = dag_univ->surf_idx_offset_ + dagmc_ptr_->index_by_handle(hit_surf);
} else {
// indicate that particle is lost
surf_idx = -1;
dist = INFINITY;
if (!dagmc_ptr_->is_implicit_complement(vol) || model::universe_map[dag_univ->id_] == model::root_universe) {
p->mark_as_lost(fmt::format("No intersection found with DAGMC cell {}", id_));
}
}
return {dist, surf_idx};
}
bool DAGCell::contains(Position r, Direction u, int32_t on_surface) const
{
moab::ErrorCode rval;
moab::EntityHandle vol = dagmc_ptr_->entity_by_index(3, dag_index_);
int result = 0;
double pnt[3] = {r.x, r.y, r.z};
double dir[3] = {u.x, u.y, u.z};
rval = dagmc_ptr_->point_in_volume(vol, pnt, result, dir);
MB_CHK_ERR_CONT(rval);
return result;
}
void DAGCell::to_hdf5_inner(hid_t group_id) const {
write_string(group_id, "geom_type", "dagmc", false);
}
BoundingBox DAGCell::bounding_box() const
{
moab::ErrorCode rval;
moab::EntityHandle vol = dagmc_ptr_->entity_by_index(3, dag_index_);
double min[3], max[3];
rval = dagmc_ptr_->getobb(vol, min, max);
MB_CHK_ERR_CONT(rval);
return {min[0], max[0], min[1], max[1], min[2], max[2]};
}
// DAGMC Functions
int32_t next_cell(DAGUniverse* dag_univ, DAGCell* cur_cell, DAGSurface* surf_xed)
{
moab::EntityHandle surf =
surf_xed->dagmc_ptr_->entity_by_index(2, surf_xed->dag_index_);
moab::EntityHandle vol =
cur_cell->dagmc_ptr_->entity_by_index(3, cur_cell->dag_index_);
moab::EntityHandle new_vol;
cur_cell->dagmc_ptr_->next_vol(surf, vol, new_vol);
return cur_cell->dagmc_ptr_->index_by_handle(new_vol) + dag_univ->cell_idx_offset_;
}
}
#endif