mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-28 06:05:58 -04:00
Move intersects method up to base class and combine the functionality using new get_index_in_direction function. Refs #1695
This commit is contained in:
aprilnovak
parent
45b7dfb5a0
commit
bcc072e0a8
2 changed files with 264 additions and 384 deletions
|
|
@ -126,9 +126,17 @@ public:
|
|||
|
||||
//! Get mesh index in a particular direction
|
||||
//!
|
||||
//! \param[in] r Position to get index for
|
||||
//! \param[in] r Coordinate to get index for
|
||||
//! \param[in] i Direction index
|
||||
virtual int get_index_in_direction(Position r, int i) const = 0;
|
||||
virtual int get_index_in_direction(double r, int i) const = 0;
|
||||
|
||||
//! Check where a line segment intersects the mesh and if it intersects at all
|
||||
//
|
||||
//! \param[in,out] r0 In: starting position, out: intersection point
|
||||
//! \param[in] r1 Ending position
|
||||
//! \param[out] ijk Indices of the mesh bin containing the intersection point
|
||||
//! \return Whether the line segment connecting r0 and r1 intersects mesh
|
||||
virtual bool intersects(Position& r0, Position r1, int* ijk) const;
|
||||
|
||||
//! Get a label for the mesh bin
|
||||
std::string bin_label(int bin) const override;
|
||||
|
|
@ -137,6 +145,11 @@ public:
|
|||
xt::xtensor<double, 1> lower_left_; //!< Lower-left coordinates of mesh
|
||||
xt::xtensor<double, 1> upper_right_; //!< Upper-right coordinates of mesh
|
||||
xt::xtensor<int, 1> shape_; //!< Number of mesh elements in each dimension
|
||||
|
||||
protected:
|
||||
virtual bool intersects_1d(Position& r0, Position r1, int* ijk) const;
|
||||
virtual bool intersects_2d(Position& r0, Position r1, int* ijk) const;
|
||||
virtual bool intersects_3d(Position& r0, Position r1, int* ijk) const;
|
||||
};
|
||||
|
||||
//==============================================================================
|
||||
|
|
@ -173,14 +186,6 @@ public:
|
|||
|
||||
// New methods
|
||||
|
||||
//! Check where a line segment intersects the mesh and if it intersects at all
|
||||
//
|
||||
//! \param[in,out] r0 In: starting position, out: intersection point
|
||||
//! \param[in] r1 Ending position
|
||||
//! \param[out] ijk Indices of the mesh bin containing the intersection point
|
||||
//! \return Whether the line segment connecting r0 and r1 intersects mesh
|
||||
bool intersects(Position& r0, Position r1, int* ijk) const;
|
||||
|
||||
//! Count number of bank sites in each mesh bin / energy bin
|
||||
//
|
||||
//! \param[in] bank Array of bank sites
|
||||
|
|
@ -193,11 +198,6 @@ public:
|
|||
|
||||
double volume_frac_; //!< Volume fraction of each mesh element
|
||||
xt::xtensor<double, 1> width_; //!< Width of each mesh element
|
||||
|
||||
private:
|
||||
bool intersects_1d(Position& r0, Position r1, int* ijk) const;
|
||||
bool intersects_2d(Position& r0, Position r1, int* ijk) const;
|
||||
bool intersects_3d(Position& r0, Position r1, int* ijk) const;
|
||||
};
|
||||
|
||||
|
||||
|
|
@ -228,16 +228,6 @@ public:
|
|||
|
||||
void to_hdf5(hid_t group) const override;
|
||||
|
||||
// New methods
|
||||
|
||||
//! Check where a line segment intersects the mesh and if it intersects at all
|
||||
//
|
||||
//! \param[in,out] r0 In: starting position, out: intersection point
|
||||
//! \param[in] r1 Ending position
|
||||
//! \param[out] ijk Indices of the mesh bin containing the intersection point
|
||||
//! \return Whether the line segment connecting r0 and r1 intersects mesh
|
||||
bool intersects(Position& r0, Position r1, int* ijk) const;
|
||||
|
||||
private:
|
||||
std::vector<std::vector<double>> grid_;
|
||||
};
|
||||
|
|
|
|||
608
src/mesh.cpp
608
src/mesh.cpp
|
|
@ -106,7 +106,7 @@ void StructuredMesh::get_indices(Position r, int* ijk, bool* in_mesh) const
|
|||
{
|
||||
*in_mesh = true;
|
||||
for (int i = 0; i < n_dimension_; ++i) {
|
||||
ijk[i] = get_index_in_direction(r, i);
|
||||
ijk[i] = get_index_in_direction(r[i], i);
|
||||
|
||||
if (ijk[i] < 1 || ijk[i] > shape_[i]) *in_mesh = false;
|
||||
}
|
||||
|
|
@ -140,6 +140,250 @@ void StructuredMesh::get_indices_from_bin(int bin, int* ijk) const
|
|||
}
|
||||
}
|
||||
|
||||
bool StructuredMesh::intersects(Position& r0, Position r1, int* ijk) const
|
||||
{
|
||||
switch(n_dimension_) {
|
||||
case 1:
|
||||
return intersects_1d(r0, r1, ijk);
|
||||
case 2:
|
||||
return intersects_2d(r0, r1, ijk);
|
||||
case 3:
|
||||
return intersects_3d(r0, r1, ijk);
|
||||
default:
|
||||
throw std::runtime_error{"Invalid number of mesh dimensions."};
|
||||
}
|
||||
}
|
||||
|
||||
bool StructuredMesh::intersects_1d(Position& r0, Position r1, int* ijk) const
|
||||
{
|
||||
// Copy coordinates of starting point
|
||||
double x0 = r0.x;
|
||||
double y0 = r0.y;
|
||||
double z0 = r0.z;
|
||||
|
||||
// Copy coordinates of ending point
|
||||
double x1 = r1.x;
|
||||
double y1 = r1.y;
|
||||
double z1 = r1.z;
|
||||
|
||||
// Copy coordinates of mesh lower_left and upper_right
|
||||
double xm0 = lower_left_[0];
|
||||
double xm1 = upper_right_[0];
|
||||
|
||||
double min_dist = INFTY;
|
||||
|
||||
// Check if line intersects left surface -- calculate the intersection point
|
||||
// (y,z)
|
||||
if ((x0 < xm0 && x1 > xm0) || (x0 > xm0 && x1 < xm0)) {
|
||||
double yi = y0 + (xm0 - x0) * (y1 - y0) / (x1 - x0);
|
||||
double zi = z0 + (xm0 - x0) * (z1 - z0) / (x1 - x0);
|
||||
if (check_intersection_point(xm0, x0, yi, yi, zi, zi, r0, min_dist)) {
|
||||
ijk[0] = 1;
|
||||
}
|
||||
}
|
||||
|
||||
// Check if line intersects right surface -- calculate the intersection point
|
||||
// (y,z)
|
||||
if ((x0 < xm1 && x1 > xm1) || (x0 > xm1 && x1 < xm1)) {
|
||||
double yi = y0 + (xm1 - x0) * (y1 - y0) / (x1 - x0);
|
||||
double zi = z0 + (xm1 - x0) * (z1 - z0) / (x1 - x0);
|
||||
if (check_intersection_point(xm1, x0, yi, yi, zi, zi, r0, min_dist)) {
|
||||
ijk[0] = shape_[0];
|
||||
}
|
||||
}
|
||||
|
||||
return min_dist < INFTY;
|
||||
}
|
||||
|
||||
bool StructuredMesh::intersects_2d(Position& r0, Position r1, int* ijk) const
|
||||
{
|
||||
// Copy coordinates of starting point
|
||||
double x0 = r0.x;
|
||||
double y0 = r0.y;
|
||||
double z0 = r0.z;
|
||||
|
||||
// Copy coordinates of ending point
|
||||
double x1 = r1.x;
|
||||
double y1 = r1.y;
|
||||
double z1 = r1.z;
|
||||
|
||||
// Copy coordinates of mesh lower_left
|
||||
double xm0 = lower_left_[0];
|
||||
double ym0 = lower_left_[1];
|
||||
|
||||
// Copy coordinates of mesh upper_right
|
||||
double xm1 = upper_right_[0];
|
||||
double ym1 = upper_right_[1];
|
||||
|
||||
double min_dist = INFTY;
|
||||
|
||||
// Check if line intersects left surface -- calculate the intersection point
|
||||
// (y,z)
|
||||
if ((x0 < xm0 && x1 > xm0) || (x0 > xm0 && x1 < xm0)) {
|
||||
double yi = y0 + (xm0 - x0) * (y1 - y0) / (x1 - x0);
|
||||
double zi = z0 + (xm0 - x0) * (z1 - z0) / (x1 - x0);
|
||||
if (yi >= ym0 && yi < ym1) {
|
||||
if (check_intersection_point(xm0, x0, yi, y0, zi, zi, r0, min_dist)) {
|
||||
ijk[0] = 1;
|
||||
ijk[1] = get_index_in_direction(yi, 1);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Check if line intersects back surface -- calculate the intersection point
|
||||
// (x,z)
|
||||
if ((y0 < ym0 && y1 > ym0) || (y0 > ym0 && y1 < ym0)) {
|
||||
double xi = x0 + (ym0 - y0) * (x1 - x0) / (y1 - y0);
|
||||
double zi = z0 + (ym0 - y0) * (z1 - z0) / (y1 - y0);
|
||||
if (xi >= xm0 && xi < xm1) {
|
||||
if (check_intersection_point(xi, x0, ym0, y0, zi, zi, r0, min_dist)) {
|
||||
ijk[0] = get_index_in_direction(xi, 0);
|
||||
ijk[1] = 1;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Check if line intersects right surface -- calculate the intersection point
|
||||
// (y,z)
|
||||
if ((x0 < xm1 && x1 > xm1) || (x0 > xm1 && x1 < xm1)) {
|
||||
double yi = y0 + (xm1 - x0) * (y1 - y0) / (x1 - x0);
|
||||
double zi = z0 + (xm1 - x0) * (z1 - z0) / (x1 - x0);
|
||||
if (yi >= ym0 && yi < ym1) {
|
||||
if (check_intersection_point(xm1, x0, yi, y0, zi, zi, r0, min_dist)) {
|
||||
ijk[0] = shape_[0];
|
||||
ijk[1] = get_index_in_direction(yi, 1);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Check if line intersects front surface -- calculate the intersection point
|
||||
// (x,z)
|
||||
if ((y0 < ym1 && y1 > ym1) || (y0 > ym1 && y1 < ym1)) {
|
||||
double xi = x0 + (ym1 - y0) * (x1 - x0) / (y1 - y0);
|
||||
double zi = z0 + (ym1 - y0) * (z1 - z0) / (y1 - y0);
|
||||
if (xi >= xm0 && xi < xm1) {
|
||||
if (check_intersection_point(xi, x0, ym1, y0, zi, zi, r0, min_dist)) {
|
||||
ijk[0] = get_index_in_direction(xi, 0);
|
||||
ijk[1] = shape_[1];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return min_dist < INFTY;
|
||||
}
|
||||
|
||||
bool StructuredMesh::intersects_3d(Position& r0, Position r1, int* ijk) const
|
||||
{
|
||||
// Copy coordinates of starting point
|
||||
double x0 = r0.x;
|
||||
double y0 = r0.y;
|
||||
double z0 = r0.z;
|
||||
|
||||
// Copy coordinates of ending point
|
||||
double x1 = r1.x;
|
||||
double y1 = r1.y;
|
||||
double z1 = r1.z;
|
||||
|
||||
// Copy coordinates of mesh lower_left
|
||||
double xm0 = lower_left_[0];
|
||||
double ym0 = lower_left_[1];
|
||||
double zm0 = lower_left_[2];
|
||||
|
||||
// Copy coordinates of mesh upper_right
|
||||
double xm1 = upper_right_[0];
|
||||
double ym1 = upper_right_[1];
|
||||
double zm1 = upper_right_[2];
|
||||
|
||||
double min_dist = INFTY;
|
||||
|
||||
// Check if line intersects left surface -- calculate the intersection point
|
||||
// (y,z)
|
||||
if ((x0 < xm0 && x1 > xm0) || (x0 > xm0 && x1 < xm0)) {
|
||||
double yi = y0 + (xm0 - x0) * (y1 - y0) / (x1 - x0);
|
||||
double zi = z0 + (xm0 - x0) * (z1 - z0) / (x1 - x0);
|
||||
if (yi >= ym0 && yi < ym1 && zi >= zm0 && zi < zm1) {
|
||||
if (check_intersection_point(xm0, x0, yi, y0, zi, z0, r0, min_dist)) {
|
||||
ijk[0] = 1;
|
||||
ijk[1] = get_index_in_direction(yi, 1);
|
||||
ijk[2] = get_index_in_direction(zi, 2);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Check if line intersects back surface -- calculate the intersection point
|
||||
// (x,z)
|
||||
if ((y0 < ym0 && y1 > ym0) || (y0 > ym0 && y1 < ym0)) {
|
||||
double xi = x0 + (ym0 - y0) * (x1 - x0) / (y1 - y0);
|
||||
double zi = z0 + (ym0 - y0) * (z1 - z0) / (y1 - y0);
|
||||
if (xi >= xm0 && xi < xm1 && zi >= zm0 && zi < zm1) {
|
||||
if (check_intersection_point(xi, x0, ym0, y0, zi, z0, r0, min_dist)) {
|
||||
ijk[0] = get_index_in_direction(xi, 0);
|
||||
ijk[1] = 1;
|
||||
ijk[2] = get_index_in_direction(zi, 2);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Check if line intersects bottom surface -- calculate the intersection
|
||||
// point (x,y)
|
||||
if ((z0 < zm0 && z1 > zm0) || (z0 > zm0 && z1 < zm0)) {
|
||||
double xi = x0 + (zm0 - z0) * (x1 - x0) / (z1 - z0);
|
||||
double yi = y0 + (zm0 - z0) * (y1 - y0) / (z1 - z0);
|
||||
if (xi >= xm0 && xi < xm1 && yi >= ym0 && yi < ym1) {
|
||||
if (check_intersection_point(xi, x0, yi, y0, zm0, z0, r0, min_dist)) {
|
||||
ijk[0] = get_index_in_direction(xi, 0);
|
||||
ijk[1] = get_index_in_direction(yi, 1);
|
||||
ijk[2] = 1;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Check if line intersects right surface -- calculate the intersection point
|
||||
// (y,z)
|
||||
if ((x0 < xm1 && x1 > xm1) || (x0 > xm1 && x1 < xm1)) {
|
||||
double yi = y0 + (xm1 - x0) * (y1 - y0) / (x1 - x0);
|
||||
double zi = z0 + (xm1 - x0) * (z1 - z0) / (x1 - x0);
|
||||
if (yi >= ym0 && yi < ym1 && zi >= zm0 && zi < zm1) {
|
||||
if (check_intersection_point(xm1, x0, yi, y0, zi, z0, r0, min_dist)) {
|
||||
ijk[0] = shape_[0];
|
||||
ijk[1] = get_index_in_direction(yi, 1);
|
||||
ijk[2] = get_index_in_direction(zi, 2);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Check if line intersects front surface -- calculate the intersection point
|
||||
// (x,z)
|
||||
if ((y0 < ym1 && y1 > ym1) || (y0 > ym1 && y1 < ym1)) {
|
||||
double xi = x0 + (ym1 - y0) * (x1 - x0) / (y1 - y0);
|
||||
double zi = z0 + (ym1 - y0) * (z1 - z0) / (y1 - y0);
|
||||
if (xi >= xm0 && xi < xm1 && zi >= zm0 && zi < zm1) {
|
||||
if (check_intersection_point(xi, x0, ym1, y0, zi, z0, r0, min_dist)) {
|
||||
ijk[0] = get_index_in_direction(xi, 0);
|
||||
ijk[1] = shape_[1];
|
||||
ijk[2] = get_index_in_direction(zi, 2);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Check if line intersects top surface -- calculate the intersection point
|
||||
// (x,y)
|
||||
if ((z0 < zm1 && z1 > zm1) || (z0 > zm1 && z1 < zm1)) {
|
||||
double xi = x0 + (zm1 - z0) * (x1 - x0) / (z1 - z0);
|
||||
double yi = y0 + (zm1 - z0) * (y1 - y0) / (z1 - z0);
|
||||
if (xi >= xm0 && xi < xm1 && yi >= ym0 && yi < ym1) {
|
||||
if (check_intersection_point(xi, x0, yi, y0, zm1, z0, r0, min_dist)) {
|
||||
ijk[0] = get_index_in_direction(xi, 0);
|
||||
ijk[1] = get_index_in_direction(yi, 1);
|
||||
ijk[2] = shape_[2];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return min_dist < INFTY;
|
||||
}
|
||||
|
||||
|
||||
//==============================================================================
|
||||
// RegularMesh implementation
|
||||
//==============================================================================
|
||||
|
|
@ -237,9 +481,9 @@ int RegularMesh::get_bin(Position r) const
|
|||
return get_bin_from_indices(ijk.data());
|
||||
}
|
||||
|
||||
int RegularMesh::get_index_in_direction(Position r, int i) const
|
||||
int RegularMesh::get_index_in_direction(double r, int i) const
|
||||
{
|
||||
return std::ceil((r[i] - lower_left_[i]) / width_[i]);
|
||||
return std::ceil((r - lower_left_[i]) / width_[i]);
|
||||
}
|
||||
|
||||
int RegularMesh::n_bins() const
|
||||
|
|
@ -254,249 +498,6 @@ int RegularMesh::n_surface_bins() const
|
|||
return 4 * n_dimension_ * n_bins();
|
||||
}
|
||||
|
||||
bool RegularMesh::intersects(Position& r0, Position r1, int* ijk) const
|
||||
{
|
||||
switch(n_dimension_) {
|
||||
case 1:
|
||||
return intersects_1d(r0, r1, ijk);
|
||||
case 2:
|
||||
return intersects_2d(r0, r1, ijk);
|
||||
case 3:
|
||||
return intersects_3d(r0, r1, ijk);
|
||||
default:
|
||||
throw std::runtime_error{"Invalid number of mesh dimensions."};
|
||||
}
|
||||
}
|
||||
|
||||
bool RegularMesh::intersects_1d(Position& r0, Position r1, int* ijk) const
|
||||
{
|
||||
// Copy coordinates of starting point
|
||||
double x0 = r0.x;
|
||||
double y0 = r0.y;
|
||||
double z0 = r0.z;
|
||||
|
||||
// Copy coordinates of ending point
|
||||
double x1 = r1.x;
|
||||
double y1 = r1.y;
|
||||
double z1 = r1.z;
|
||||
|
||||
// Copy coordinates of mesh lower_left and upper_right
|
||||
double xm0 = lower_left_[0];
|
||||
double xm1 = upper_right_[0];
|
||||
|
||||
double min_dist = INFTY;
|
||||
|
||||
// Check if line intersects left surface -- calculate the intersection point
|
||||
// (y,z)
|
||||
if ((x0 < xm0 && x1 > xm0) || (x0 > xm0 && x1 < xm0)) {
|
||||
double yi = y0 + (xm0 - x0) * (y1 - y0) / (x1 - x0);
|
||||
double zi = z0 + (xm0 - x0) * (z1 - z0) / (x1 - x0);
|
||||
if (check_intersection_point(xm0, x0, yi, yi, zi, zi, r0, min_dist)) {
|
||||
ijk[0] = 1;
|
||||
}
|
||||
}
|
||||
|
||||
// Check if line intersects right surface -- calculate the intersection point
|
||||
// (y,z)
|
||||
if ((x0 < xm1 && x1 > xm1) || (x0 > xm1 && x1 < xm1)) {
|
||||
double yi = y0 + (xm1 - x0) * (y1 - y0) / (x1 - x0);
|
||||
double zi = z0 + (xm1 - x0) * (z1 - z0) / (x1 - x0);
|
||||
if (check_intersection_point(xm1, x0, yi, yi, zi, zi, r0, min_dist)) {
|
||||
ijk[0] = shape_[0];
|
||||
}
|
||||
}
|
||||
|
||||
return min_dist < INFTY;
|
||||
}
|
||||
|
||||
bool RegularMesh::intersects_2d(Position& r0, Position r1, int* ijk) const
|
||||
{
|
||||
// Copy coordinates of starting point
|
||||
double x0 = r0.x;
|
||||
double y0 = r0.y;
|
||||
double z0 = r0.z;
|
||||
|
||||
// Copy coordinates of ending point
|
||||
double x1 = r1.x;
|
||||
double y1 = r1.y;
|
||||
double z1 = r1.z;
|
||||
|
||||
// Copy coordinates of mesh lower_left
|
||||
double xm0 = lower_left_[0];
|
||||
double ym0 = lower_left_[1];
|
||||
|
||||
// Copy coordinates of mesh upper_right
|
||||
double xm1 = upper_right_[0];
|
||||
double ym1 = upper_right_[1];
|
||||
|
||||
double min_dist = INFTY;
|
||||
|
||||
// Check if line intersects left surface -- calculate the intersection point
|
||||
// (y,z)
|
||||
if ((x0 < xm0 && x1 > xm0) || (x0 > xm0 && x1 < xm0)) {
|
||||
double yi = y0 + (xm0 - x0) * (y1 - y0) / (x1 - x0);
|
||||
double zi = z0 + (xm0 - x0) * (z1 - z0) / (x1 - x0);
|
||||
if (yi >= ym0 && yi < ym1) {
|
||||
if (check_intersection_point(xm0, x0, yi, y0, zi, zi, r0, min_dist)) {
|
||||
ijk[0] = 1;
|
||||
ijk[1] = std::ceil((yi - lower_left_[1]) / width_[1]);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Check if line intersects back surface -- calculate the intersection point
|
||||
// (x,z)
|
||||
if ((y0 < ym0 && y1 > ym0) || (y0 > ym0 && y1 < ym0)) {
|
||||
double xi = x0 + (ym0 - y0) * (x1 - x0) / (y1 - y0);
|
||||
double zi = z0 + (ym0 - y0) * (z1 - z0) / (y1 - y0);
|
||||
if (xi >= xm0 && xi < xm1) {
|
||||
if (check_intersection_point(xi, x0, ym0, y0, zi, zi, r0, min_dist)) {
|
||||
ijk[0] = std::ceil((xi - lower_left_[0]) / width_[0]);
|
||||
ijk[1] = 1;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Check if line intersects right surface -- calculate the intersection point
|
||||
// (y,z)
|
||||
if ((x0 < xm1 && x1 > xm1) || (x0 > xm1 && x1 < xm1)) {
|
||||
double yi = y0 + (xm1 - x0) * (y1 - y0) / (x1 - x0);
|
||||
double zi = z0 + (xm1 - x0) * (z1 - z0) / (x1 - x0);
|
||||
if (yi >= ym0 && yi < ym1) {
|
||||
if (check_intersection_point(xm1, x0, yi, y0, zi, zi, r0, min_dist)) {
|
||||
ijk[0] = shape_[0];
|
||||
ijk[1] = std::ceil((yi - lower_left_[1]) / width_[1]);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Check if line intersects front surface -- calculate the intersection point
|
||||
// (x,z)
|
||||
if ((y0 < ym1 && y1 > ym1) || (y0 > ym1 && y1 < ym1)) {
|
||||
double xi = x0 + (ym1 - y0) * (x1 - x0) / (y1 - y0);
|
||||
double zi = z0 + (ym1 - y0) * (z1 - z0) / (y1 - y0);
|
||||
if (xi >= xm0 && xi < xm1) {
|
||||
if (check_intersection_point(xi, x0, ym1, y0, zi, zi, r0, min_dist)) {
|
||||
ijk[0] = std::ceil((xi - lower_left_[0]) / width_[0]);
|
||||
ijk[1] = shape_[1];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return min_dist < INFTY;
|
||||
}
|
||||
|
||||
bool RegularMesh::intersects_3d(Position& r0, Position r1, int* ijk) const
|
||||
{
|
||||
// Copy coordinates of starting point
|
||||
double x0 = r0.x;
|
||||
double y0 = r0.y;
|
||||
double z0 = r0.z;
|
||||
|
||||
// Copy coordinates of ending point
|
||||
double x1 = r1.x;
|
||||
double y1 = r1.y;
|
||||
double z1 = r1.z;
|
||||
|
||||
// Copy coordinates of mesh lower_left
|
||||
double xm0 = lower_left_[0];
|
||||
double ym0 = lower_left_[1];
|
||||
double zm0 = lower_left_[2];
|
||||
|
||||
// Copy coordinates of mesh upper_right
|
||||
double xm1 = upper_right_[0];
|
||||
double ym1 = upper_right_[1];
|
||||
double zm1 = upper_right_[2];
|
||||
|
||||
double min_dist = INFTY;
|
||||
|
||||
// Check if line intersects left surface -- calculate the intersection point
|
||||
// (y,z)
|
||||
if ((x0 < xm0 && x1 > xm0) || (x0 > xm0 && x1 < xm0)) {
|
||||
double yi = y0 + (xm0 - x0) * (y1 - y0) / (x1 - x0);
|
||||
double zi = z0 + (xm0 - x0) * (z1 - z0) / (x1 - x0);
|
||||
if (yi >= ym0 && yi < ym1 && zi >= zm0 && zi < zm1) {
|
||||
if (check_intersection_point(xm0, x0, yi, y0, zi, z0, r0, min_dist)) {
|
||||
ijk[0] = 1;
|
||||
ijk[1] = std::ceil((yi - lower_left_[1]) / width_[1]);
|
||||
ijk[2] = std::ceil((zi - lower_left_[2]) / width_[2]);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Check if line intersects back surface -- calculate the intersection point
|
||||
// (x,z)
|
||||
if ((y0 < ym0 && y1 > ym0) || (y0 > ym0 && y1 < ym0)) {
|
||||
double xi = x0 + (ym0 - y0) * (x1 - x0) / (y1 - y0);
|
||||
double zi = z0 + (ym0 - y0) * (z1 - z0) / (y1 - y0);
|
||||
if (xi >= xm0 && xi < xm1 && zi >= zm0 && zi < zm1) {
|
||||
if (check_intersection_point(xi, x0, ym0, y0, zi, z0, r0, min_dist)) {
|
||||
ijk[0] = std::ceil((xi - lower_left_[0]) / width_[0]);
|
||||
ijk[1] = 1;
|
||||
ijk[2] = std::ceil((zi - lower_left_[2]) / width_[2]);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Check if line intersects bottom surface -- calculate the intersection
|
||||
// point (x,y)
|
||||
if ((z0 < zm0 && z1 > zm0) || (z0 > zm0 && z1 < zm0)) {
|
||||
double xi = x0 + (zm0 - z0) * (x1 - x0) / (z1 - z0);
|
||||
double yi = y0 + (zm0 - z0) * (y1 - y0) / (z1 - z0);
|
||||
if (xi >= xm0 && xi < xm1 && yi >= ym0 && yi < ym1) {
|
||||
if (check_intersection_point(xi, x0, yi, y0, zm0, z0, r0, min_dist)) {
|
||||
ijk[0] = std::ceil((xi - lower_left_[0]) / width_[0]);
|
||||
ijk[1] = std::ceil((yi - lower_left_[1]) / width_[1]);
|
||||
ijk[2] = 1;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Check if line intersects right surface -- calculate the intersection point
|
||||
// (y,z)
|
||||
if ((x0 < xm1 && x1 > xm1) || (x0 > xm1 && x1 < xm1)) {
|
||||
double yi = y0 + (xm1 - x0) * (y1 - y0) / (x1 - x0);
|
||||
double zi = z0 + (xm1 - x0) * (z1 - z0) / (x1 - x0);
|
||||
if (yi >= ym0 && yi < ym1 && zi >= zm0 && zi < zm1) {
|
||||
if (check_intersection_point(xm1, x0, yi, y0, zi, z0, r0, min_dist)) {
|
||||
ijk[0] = shape_[0];
|
||||
ijk[1] = std::ceil((yi - lower_left_[1]) / width_[1]);
|
||||
ijk[2] = std::ceil((zi - lower_left_[2]) / width_[2]);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Check if line intersects front surface -- calculate the intersection point
|
||||
// (x,z)
|
||||
if ((y0 < ym1 && y1 > ym1) || (y0 > ym1 && y1 < ym1)) {
|
||||
double xi = x0 + (ym1 - y0) * (x1 - x0) / (y1 - y0);
|
||||
double zi = z0 + (ym1 - y0) * (z1 - z0) / (y1 - y0);
|
||||
if (xi >= xm0 && xi < xm1 && zi >= zm0 && zi < zm1) {
|
||||
if (check_intersection_point(xi, x0, ym1, y0, zi, z0, r0, min_dist)) {
|
||||
ijk[0] = std::ceil((xi - lower_left_[0]) / width_[0]);
|
||||
ijk[1] = shape_[1];
|
||||
ijk[2] = std::ceil((zi - lower_left_[2]) / width_[2]);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Check if line intersects top surface -- calculate the intersection point
|
||||
// (x,y)
|
||||
if ((z0 < zm1 && z1 > zm1) || (z0 > zm1 && z1 < zm1)) {
|
||||
double xi = x0 + (zm1 - z0) * (x1 - x0) / (z1 - z0);
|
||||
double yi = y0 + (zm1 - z0) * (y1 - y0) / (z1 - z0);
|
||||
if (xi >= xm0 && xi < xm1 && yi >= ym0 && yi < ym1) {
|
||||
if (check_intersection_point(xi, x0, yi, y0, zm1, z0, r0, min_dist)) {
|
||||
ijk[0] = std::ceil((xi - lower_left_[0]) / width_[0]);
|
||||
ijk[1] = std::ceil((yi - lower_left_[1]) / width_[1]);
|
||||
ijk[2] = shape_[2];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return min_dist < INFTY;
|
||||
}
|
||||
|
||||
void RegularMesh::bins_crossed(const Particle& p, std::vector<int>& bins,
|
||||
std::vector<double>& lengths) const
|
||||
{
|
||||
|
|
@ -1173,9 +1174,9 @@ int RectilinearMesh::get_bin(Position r) const
|
|||
return get_bin_from_indices(ijk);
|
||||
}
|
||||
|
||||
int RectilinearMesh::get_index_in_direction(Position r, int i) const
|
||||
int RectilinearMesh::get_index_in_direction(double r, int i) const
|
||||
{
|
||||
return lower_bound_index(grid_[i].begin(), grid_[i].end(), r[i]) + 1;
|
||||
return lower_bound_index(grid_[i].begin(), grid_[i].end(), r) + 1;
|
||||
}
|
||||
|
||||
int RectilinearMesh::n_bins() const
|
||||
|
|
@ -1230,117 +1231,6 @@ void RectilinearMesh::to_hdf5(hid_t group) const
|
|||
close_group(mesh_group);
|
||||
}
|
||||
|
||||
bool RectilinearMesh::intersects(Position& r0, Position r1, int* ijk) const
|
||||
{
|
||||
// Copy coordinates of starting point
|
||||
double x0 = r0.x;
|
||||
double y0 = r0.y;
|
||||
double z0 = r0.z;
|
||||
|
||||
// Copy coordinates of ending point
|
||||
double x1 = r1.x;
|
||||
double y1 = r1.y;
|
||||
double z1 = r1.z;
|
||||
|
||||
// Copy coordinates of mesh lower_left
|
||||
double xm0 = grid_[0].front();
|
||||
double ym0 = grid_[1].front();
|
||||
double zm0 = grid_[2].front();
|
||||
|
||||
// Copy coordinates of mesh upper_right
|
||||
double xm1 = grid_[0].back();
|
||||
double ym1 = grid_[1].back();
|
||||
double zm1 = grid_[2].back();
|
||||
|
||||
double min_dist = INFTY;
|
||||
|
||||
// Check if line intersects left surface -- calculate the intersection point
|
||||
// (y,z)
|
||||
if ((x0 < xm0 && x1 > xm0) || (x0 > xm0 && x1 < xm0)) {
|
||||
double yi = y0 + (xm0 - x0) * (y1 - y0) / (x1 - x0);
|
||||
double zi = z0 + (xm0 - x0) * (z1 - z0) / (x1 - x0);
|
||||
if (yi >= ym0 && yi < ym1 && zi >= zm0 && zi < zm1) {
|
||||
if (check_intersection_point(xm0, x0, yi, y0, zi, z0, r0, min_dist)) {
|
||||
ijk[0] = 1;
|
||||
ijk[1] = lower_bound_index(grid_[1].begin(), grid_[1].end(), yi) + 1;
|
||||
ijk[2] = lower_bound_index(grid_[2].begin(), grid_[2].end(), zi) + 1;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Check if line intersects back surface -- calculate the intersection point
|
||||
// (x,z)
|
||||
if ((y0 < ym0 && y1 > ym0) || (y0 > ym0 && y1 < ym0)) {
|
||||
double xi = x0 + (ym0 - y0) * (x1 - x0) / (y1 - y0);
|
||||
double zi = z0 + (ym0 - y0) * (z1 - z0) / (y1 - y0);
|
||||
if (xi >= xm0 && xi < xm1 && zi >= zm0 && zi < zm1) {
|
||||
if (check_intersection_point(xi, x0, ym0, y0, zi, z0, r0, min_dist)) {
|
||||
ijk[0] = lower_bound_index(grid_[0].begin(), grid_[0].end(), xi) + 1;
|
||||
ijk[1] = 1;
|
||||
ijk[2] = lower_bound_index(grid_[2].begin(), grid_[2].end(), zi) + 1;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Check if line intersects bottom surface -- calculate the intersection
|
||||
// point (x,y)
|
||||
if ((z0 < zm0 && z1 > zm0) || (z0 > zm0 && z1 < zm0)) {
|
||||
double xi = x0 + (zm0 - z0) * (x1 - x0) / (z1 - z0);
|
||||
double yi = y0 + (zm0 - z0) * (y1 - y0) / (z1 - z0);
|
||||
if (xi >= xm0 && xi < xm1 && yi >= ym0 && yi < ym1) {
|
||||
if (check_intersection_point(xi, x0, yi, y0, zm0, z0, r0, min_dist)) {
|
||||
ijk[0] = lower_bound_index(grid_[0].begin(), grid_[0].end(), xi) + 1;
|
||||
ijk[1] = lower_bound_index(grid_[1].begin(), grid_[1].end(), yi) + 1;
|
||||
ijk[2] = 1;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Check if line intersects right surface -- calculate the intersection point
|
||||
// (y,z)
|
||||
if ((x0 < xm1 && x1 > xm1) || (x0 > xm1 && x1 < xm1)) {
|
||||
double yi = y0 + (xm1 - x0) * (y1 - y0) / (x1 - x0);
|
||||
double zi = z0 + (xm1 - x0) * (z1 - z0) / (x1 - x0);
|
||||
if (yi >= ym0 && yi < ym1 && zi >= zm0 && zi < zm1) {
|
||||
if (check_intersection_point(xm1, x0, yi, y0, zi, z0, r0, min_dist)) {
|
||||
ijk[0] = shape_[0];
|
||||
ijk[1] = lower_bound_index(grid_[1].begin(), grid_[1].end(), yi) + 1;
|
||||
ijk[2] = lower_bound_index(grid_[2].begin(), grid_[2].end(), zi) + 1;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Check if line intersects front surface -- calculate the intersection point
|
||||
// (x,z)
|
||||
if ((y0 < ym1 && y1 > ym1) || (y0 > ym1 && y1 < ym1)) {
|
||||
double xi = x0 + (ym1 - y0) * (x1 - x0) / (y1 - y0);
|
||||
double zi = z0 + (ym1 - y0) * (z1 - z0) / (y1 - y0);
|
||||
if (xi >= xm0 && xi < xm1 && zi >= zm0 && zi < zm1) {
|
||||
if (check_intersection_point(xi, x0, ym1, y0, zi, z0, r0, min_dist)) {
|
||||
ijk[0] = lower_bound_index(grid_[0].begin(), grid_[0].end(), xi) + 1;
|
||||
ijk[1] = shape_[1];
|
||||
ijk[2] = lower_bound_index(grid_[2].begin(), grid_[2].end(), zi) + 1;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Check if line intersects top surface -- calculate the intersection point
|
||||
// (x,y)
|
||||
if ((z0 < zm1 && z1 > zm1) || (z0 > zm1 && z1 < zm1)) {
|
||||
double xi = x0 + (zm1 - z0) * (x1 - x0) / (z1 - z0);
|
||||
double yi = y0 + (zm1 - z0) * (y1 - y0) / (z1 - z0);
|
||||
if (xi >= xm0 && xi < xm1 && yi >= ym0 && yi < ym1) {
|
||||
if (check_intersection_point(xi, x0, yi, y0, zm1, z0, r0, min_dist)) {
|
||||
ijk[0] = lower_bound_index(grid_[0].begin(), grid_[0].end(), xi) + 1;
|
||||
ijk[1] = lower_bound_index(grid_[1].begin(), grid_[1].end(), yi) + 1;
|
||||
ijk[2] = shape_[2];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return min_dist < INFTY;
|
||||
}
|
||||
|
||||
//==============================================================================
|
||||
// Helper functions for the C API
|
||||
//==============================================================================
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue