Corrections after rebase.

This commit is contained in:
Patrick Shriwise 2019-11-26 04:57:16 -06:00
parent 66346e8356
commit ff0af7d594
2 changed files with 133 additions and 92 deletions

View file

@ -37,62 +37,9 @@ extern std::unordered_map<int32_t, int32_t> mesh_map;
} // namespace model
class Mesh {
public:
// Constructor
Mesh() {}; // empty constructor
Mesh(pugi::xml_node node);
//! Determine which bins were crossed by a particle
//
//! \param[in] p Particle to check
//! \param[out] bins Bins that were crossed
//! \param[out] lengths Fraction of tracklength in each bin
virtual void bins_crossed(const Particle* p, std::vector<int>& bins,
std::vector<double>& lengths) const = 0;
//! Write mesh data to an HDF5 group
//
//! \param[in] group HDF5 group
virtual void to_hdf5(hid_t group) const = 0;
//! Determine which surface bins were crossed by a particle
//
//! \param[in] p Particle to check
//! \param[out] bins Surface bins that were crossed
virtual void
surface_bins_crossed(const Particle* p, std::vector<int>& bins) const = 0;
//! Get bin at a given position in space
//
//! \param[in] r Position to get bin for
//! \return Mesh bin
virtual int get_bin(Position r) const = 0;
virtual std::string get_label_for_bin(int bin) const = 0;
//! Count number of bank sites in each mesh bin / energy bin
//
//! \param[in] bank Array of bank sites
//! \param[out] Whether any bank sites are outside the mesh
//! \return Array indicating number of sites in each mesh/energy bin
virtual xt::xarray<double>
count_sites(const std::vector<Particle::Bank>& bank, bool* outside) const;
virtual double get_volume_frac(int bin = -1) const = 0;
virtual int num_bins() const = 0;
int id_ {-1}; //!< User-specified ID
int n_dimension_; //!< Number of dimensions
};
//==============================================================================
//! Tessellation of n-dimensional Euclidean space by congruent squares or cubes
//==============================================================================
class RegularMesh : public Mesh {
class Mesh
{
public:
// Constructors and destructor
Mesh() = default;
@ -171,6 +118,60 @@ public:
xt::xtensor<double, 1> upper_right_; //!< Upper-right coordinates of mesh
};
// class Mesh {
// public:
// // Constructors
// Mesh() = default;
// Mesh(pugi::xml_node node);
// // Destructor
// virtual ~Mesh() = default;
// //! Determine which bins were crossed by a particle
// //
// //! \param[in] p Particle to check
// //! \param[out] bins Bins that were crossed
// //! \param[out] lengths Fraction of tracklength in each bin
// virtual void bins_crossed(const Particle* p, std::vector<int>& bins,
// std::vector<double>& lengths) const = 0;
// //! Write mesh data to an HDF5 group
// //
// //! \param[in] group HDF5 group
// virtual void to_hdf5(hid_t group) const = 0;
// //! Determine which surface bins were crossed by a particle
// //
// //! \param[in] p Particle to check
// //! \param[out] bins Surface bins that were crossed
// virtual void
// surface_bins_crossed(const Particle* p, std::vector<int>& bins) const = 0;
// //! Get bin at a given position in space
// //
// //! \param[in] r Position to get bin for
// //! \return Mesh bin
// virtual int get_bin(Position r) const = 0;
// virtual std::string get_label_for_bin(int bin) const = 0;
// //! Count number of bank sites in each mesh bin / energy bin
// //
// //! \param[in] bank Array of bank sites
// //! \param[out] Whether any bank sites are outside the mesh
// //! \return Array indicating number of sites in each mesh/energy bin
// virtual xt::xarray<double>
// count_sites(const std::vector<Particle::Bank>& bank, bool* outside) const;
// virtual double get_volume_frac(int bin = -1) const = 0;
// virtual int num_bins() const = 0;
// int id_ {-1}; //!< User-specified ID
// int n_dimension_; //!< Number of dimensions
// };
//==============================================================================
//! Tessellation of n-dimensional Euclidean space by congruent squares or cubes
//==============================================================================
@ -393,6 +394,15 @@ intersect_track(const moab::CartVect& start,
//! \return MOAB EntityHandle of tet
moab::EntityHandle get_ent_handle_from_bin(int bin) const;
int get_bin_from_indices(const int* ijk) const override;
void get_indices(Position r, int* ijk, bool* in_mesh) const override;
void get_indices_from_bin(int bin, int* ijk) const override;
std::pair<std::vector<double>, std::vector<double>>
plot(Position plot_ll, Position plot_ur) const override;
//! Builds a KDTree for all tetrahedra in the mesh. All
//! triangles representing 2D faces of the mesh are
//! added to the tree as well.
@ -420,9 +430,11 @@ public:
std::string get_label_for_bin(int bin) const;
double get_volume_frac(int bin = -1) const;
// double get_volume_frac(int bin = -1) const;
int num_bins() const;
int n_bins() const override;
int n_surface_bins() const override;
std::string filename_; //<! Path to unstructured mesh file

View file

@ -69,8 +69,6 @@ inline bool check_intersection_point(double x1, double x0, double y1,
// Mesh implementation
//==============================================================================
Mesh::Mesh(pugi::xml_node node)
{
Mesh::Mesh(pugi::xml_node node) {
// Copy mesh id
if (check_for_node(node, "id")) {
@ -800,7 +798,7 @@ RegularMesh::count_sites(const Particle::Bank* bank, int64_t length,
bool* outside) const
{
// Determine shape of array for counts
std::size_t m = num_bins();
std::size_t m = n_bins();
std::vector<std::size_t> shape = {m};
// Create array of zeros
@ -850,28 +848,22 @@ RegularMesh::count_sites(const Particle::Bank* bank, int64_t length,
return counts;
}
std::string RegularMesh::get_label_for_bin(int bin) const {
int ijk[n_dimension_];
get_indices_from_bin(bin, ijk);
// std::string RegularMesh::get_label_for_bin(int bin) const {
// int ijk[n_dimension_];
// get_indices_from_bin(bin, ijk);
std::stringstream out;
out << "Mesh Index (" << ijk[0];
if (n_dimension_ > 1) out << ", " << ijk[1];
if (n_dimension_ > 2) out << ", " << ijk[2];
out << ")";
// std::stringstream out;
// out << "Mesh Index (" << ijk[0];
// if (n_dimension_ > 1) out << ", " << ijk[1];
// if (n_dimension_ > 2) out << ", " << ijk[2];
// out << ")";
return out.str();
}
// return out.str();
// }
double RegularMesh::get_volume_frac(int bin) const {
return volume_frac_;
}
int RegularMesh::num_bins() const {
int n_bins = 1;
for (auto v : shape_) n_bins *= v;
return n_bins;
}
// double RegularMesh::get_volume_frac(int bin) const {
// return volume_frac_;
// }
//==============================================================================
// RectilinearMesh implementation
@ -1936,10 +1928,35 @@ UnstructuredMesh::point_in_tet(const moab::CartVect& r, moab::EntityHandle tet)
return in_tet;
}
int
UnstructuredMesh::get_bin_from_indices(const int* ijk) const {
if (ijk[0] >= n_bins()) {
std::stringstream s;
s << "Invalid bin: " << ijk[0];
fatal_error(s);
}
int bin = ehs_[ijk[0]] - ehs_[0];
return bin;
}
void
UnstructuredMesh::get_indices(Position r, int* ijk, bool* in_mesh) const {
int bin = get_bin(r);
ijk[0]= bin;
*in_mesh = bin != -1;
}
void UnstructuredMesh::get_indices_from_bin(int bin, int* ijk) const {
ijk[0] = bin;
}
std::pair<std::vector<double>, std::vector<double>>
UnstructuredMesh::plot(Position plot_ll, Position plot_ur) const { return {}; }
int
UnstructuredMesh::get_bin_from_ent_handle(moab::EntityHandle eh) const {
int bin = eh - ehs_[0];
if (bin >= num_bins()) {
if (bin >= n_bins()) {
std::stringstream s;
s << "Invalid bin: " << bin;
fatal_error(s);
@ -1949,7 +1966,7 @@ UnstructuredMesh::get_bin_from_ent_handle(moab::EntityHandle eh) const {
moab::EntityHandle
UnstructuredMesh::get_ent_handle_from_bin(int bin) const {
if (bin >= num_bins()) {
if (bin >= n_bins()) {
std::stringstream s;
s << "Invalid bin index: " << bin;
fatal_error(s);
@ -1957,19 +1974,31 @@ UnstructuredMesh::get_ent_handle_from_bin(int bin) const {
return ehs_[bin];
}
double UnstructuredMesh::get_volume_frac(int bin) const {
if (bin == -1) { return 0.0; }
if (bin > ehs_.size() || bin < -1) {
std::stringstream msg;
msg << "Invalid bin " << bin << " for umesh with id " << id_;
fatal_error(msg);
}
moab::EntityHandle tet = get_ent_handle_from_bin(bin);
return tet_volume(tet);
// double UnstructuredMesh::get_volume_frac(int bin) const {
// if (bin == -1) { return 0.0; }
// if (bin > ehs_.size() || bin < -1) {
// std::stringstream msg;
// msg << "Invalid bin " << bin << " for umesh with id " << id_;
// fatal_error(msg);
// }
// moab::EntityHandle tet = get_ent_handle_from_bin(bin);
// return tet_volume(tet);
// }
int UnstructuredMesh::n_bins() const {
return ehs_.size();
}
int UnstructuredMesh::num_bins() const {
return ehs_.size();
int UnstructuredMesh::n_surface_bins() const {
// collect all triangles in the set of tets for this mesh
moab::Range tris;
moab::ErrorCode rval;
rval = mbi_->get_entities_by_type(0, moab::MBTRI, tris);
if (rval != moab::MB_SUCCESS) {
warning("Failed to get all triangles in the mesh instance");
return -1;
}
return 2 * tris.size();
}
#endif