Relying on LibMesh's BoundingBox and Sphere.

This commit is contained in:
Patrick Shriwise 2020-01-13 09:31:02 -06:00
parent 57041cee5d
commit 80d35e2bb8
2 changed files with 13 additions and 11 deletions

View file

@ -24,6 +24,7 @@
#endif
#ifdef LIBMESH
#include "libmesh/bounding_box.h"
#include "libmesh/libmesh.h"
#include "libmesh/elem.h"
#include "libmesh/equation_systems.h"
@ -32,6 +33,7 @@
#include "libmesh/dof_map.h"
#include "libmesh/mesh.h"
#include "libmesh/point.h"
#include "libmesh/sphere.h"
#endif
namespace openmc {
@ -578,7 +580,8 @@ private:
std::vector<std::unique_ptr<libMesh::PointLocatorBase>> point_locators_; //!< pointers to locators for each thread
std::unique_ptr<libMesh::EquationSystems> equation_systems_; //!< pointer to the equation systems of the mesh (for result output)
std::map<std::string, unsigned int> variable_map_; //!< mappint of variable names (scores) to their numbers on the mesh
BoundingBox bbox_; //!< bounding box of the mesh
libMesh::BoundingBox bbox_; //!< bounding box of the mesh
libMesh::Sphere bsphere_; //<! bounding sphere of the mesh
std::string eq_system_name_; //!< name of the equation system holding OpenMC results
libMesh::Elem* first_element_; //!< pointer to the first element in the mesh (maybe should be a key?)
std::set<libMesh::Elem*> boundary_elements_; //<! all boundary elements in the mesh

View file

@ -28,6 +28,10 @@
#include "openmc/tallies/filter.h"
#include "openmc/xml_interface.h"
#ifdef LIBMESH
#include "libmesh/mesh_tools.h"
#endif
namespace openmc {
//==============================================================================
@ -2083,19 +2087,13 @@ LibMesh::LibMesh(pugi::xml_node node) : UnstructuredMeshBase(node) {
first_element_ = *m_->elements_begin();
// bounding box for the mesh
bbox_ = {INFTY, -INFTY, INFTY, -INFTY, INFTY, -INFTY};
bbox_ = libMesh::MeshTools::create_bounding_box(*m_);
bsphere_ = libMesh::MeshTools::bounding_sphere(*m_);
// determine boundary elements and create bounding box
for (int i = 0; i < m_->n_elem(); i++) {
auto e = m_->elem_ptr(i);
// update bounding box for each node
for (int j = 0; j < e->n_nodes(); j++) {
auto n = e->node_ref(j);
Position r(n(0), n(1), n(2));
bbox_.update(r);
}
for (int k = 0; k < e->n_neighbors(); k++) {
if (!e->neighbor_ptr(k)) {
boundary_elements_.insert(e);
@ -2210,11 +2208,12 @@ LibMesh::bins_crossed(const Particle* p,
int
LibMesh::get_bin(Position r) const
{
if (!bbox_.contains(r)) { return -1; }
// look-up a tet using the point locator
libMesh::Point p(r.x, r.y, r.z);
// quick rejection check
// if (bsphere_.above_surface(p) || !bbox_.contains_point(p)) { return -1; }
if (!bbox_.contains_point(p)) { return -1; }
int thread = omp_get_thread_num();
auto e = (*point_locators_[thread])(p);