Skip to content

Let a mesh surface crossing change more than one index #4096

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
GuySten wants to merge 1 commit into
base: develop
Choose a base branch
from
+92 −49
Open

Let a mesh surface crossing change more than one index #4096

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
47 changes: 42 additions & 5 deletions include/openmc/mesh.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -295,21 +295,29 @@ class Mesh {
int n_dimension_ {-1}; //!< Number of dimensions
};

//! Maximum number of axes a structured mesh may have.
//
//! Mesh indices, shapes and the per-axis distance-to-surface arrays are all
//! sized by this, so it is the single place to change if a mesh type ever
//! needs more.
constexpr int MESH_MAX_AXES {3};

class StructuredMesh : public Mesh {
public:
StructuredMesh() = default;
StructuredMesh(pugi::xml_node node) : Mesh {node} {};
StructuredMesh(hid_t group) : Mesh {group} {};
virtual ~StructuredMesh() = default;

using MeshIndex = std::array<int, 3>;
using MeshIndex = std::array<int, MESH_MAX_AXES>;

struct MeshDistance {
MeshDistance() = default;
MeshDistance(int _index, bool _max_surface, double _distance)
: next_index {_index}, max_surface {_max_surface}, distance {_distance}
MeshDistance(MeshIndex _offset, bool _max_surface, double _distance)
: offset {_offset}, max_surface {_max_surface}, distance {_distance}
{}
int next_index {-1};
//! Change in each index when this surface is crossed
MeshIndex offset {0, 0, 0};
bool max_surface {true};
double distance {INFTY};
bool operator<(const MeshDistance& o) const
Expand All @@ -318,6 +326,12 @@ class StructuredMesh : public Mesh {
}
};

//! Fold an index back into range on any periodic axis.
//
//! Called after an index is advanced by a MeshDistance offset. Meshes with a
//! periodic axis override this; for the rest it is a no-op.
virtual void sanitize_index(MeshIndex& idx) const {}

Position sample_element(int32_t bin, uint64_t* seed) const override
{
return sample_element(get_indices_from_bin(bin), seed);
Expand Down Expand Up @@ -423,6 +437,17 @@ class StructuredMesh : public Mesh {
virtual MeshDistance distance_to_grid_boundary(const MeshIndex& ijk, int i,
const Position& r0, const Direction& u, double l) const = 0;

//! Distance to travel to re-enter the mesh from outside it
//
//! \param[in] ijk Current (out of range) indices
//! \param[in] distances Distance to the next surface on each axis
//! \param[in] traveled_distance Distance travelled so far
//! \param[out] k_max Axis whose surface is crossed last, or -1
//! \return Updated travelled distance
double distance_to_mesh(const MeshIndex& ijk,
const std::array<MeshDistance, MESH_MAX_AXES>& distances,
double traveled_distance, int& k_max) const;

//! Get a label for the mesh bin
std::string bin_label(int bin) const override;

Expand Down Expand Up @@ -458,7 +483,8 @@ class StructuredMesh : public Mesh {
virtual double volume(const MeshIndex& ijk) const = 0;

// Data members
std::array<int, 3> shape_; //!< Number of mesh elements in each dimension
std::array<int, MESH_MAX_AXES>
shape_; //!< Number of mesh elements in each dimension

protected:
};
Expand Down Expand Up @@ -587,6 +613,11 @@ class CylindricalMesh : public PeriodicStructuredMesh {
CylindricalMesh(hid_t group);

// Overridden methods
void sanitize_index(MeshIndex& idx) const override
{
idx[1] = sanitize_phi(idx[1]);
}

virtual MeshIndex get_indices(Position r, bool& in_mesh) const override;

int get_index_in_direction(double r, int i) const override;
Expand Down Expand Up @@ -654,6 +685,12 @@ class SphericalMesh : public PeriodicStructuredMesh {
SphericalMesh(hid_t group);

// Overridden methods
void sanitize_index(MeshIndex& idx) const override
{
idx[1] = sanitize_theta(idx[1]);
idx[2] = sanitize_phi(idx[2]);
}

virtual MeshIndex get_indices(Position r, bool& in_mesh) const override;

int get_index_in_direction(double r, int i) const override;
Expand Down
94 changes: 50 additions & 44 deletions src/mesh.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1208,7 +1208,7 @@ void StructuredMesh::raytrace_mesh(
}

// Calculate initial distances to next surfaces in all three dimensions
std::array<MeshDistance, 3> distances;
std::array<MeshDistance, MESH_MAX_AXES> distances;
for (int k = 0; k < n; ++k) {
distances[k] = distance_to_grid_boundary(ijk, k, local_r, u, 0.0);
}
Expand Down Expand Up @@ -1238,7 +1238,9 @@ void StructuredMesh::raytrace_mesh(

// Update cell and calculate distance to next surface in k-direction.
// The two other directions are still valid!
ijk[k] = distances[k].next_index;
for (int j = 0; j < MESH_MAX_AXES; ++j)
ijk[j] += distances[k].offset[j];
sanitize_index(ijk);
distances[k] =
distance_to_grid_boundary(ijk, k, local_r, u, traveled_distance);

Expand All @@ -1251,22 +1253,9 @@ void StructuredMesh::raytrace_mesh(
tally.surface(ijk, k, !distances[k].max_surface, true);

} else { // not inside mesh

// For all directions outside the mesh, find the distance that we need
// to travel to reach the next surface. Use the largest distance, as
// only this will cross all outer surfaces.
int k_max {-1};
for (int k = 0; k < n; ++k) {
if ((ijk[k] < 1 || ijk[k] > shape_[k]) &&
(distances[k].distance > traveled_distance)) {
traveled_distance = distances[k].distance;
k_max = k;
}
}
// Assure some distance is traveled
if (k_max == -1) {
traveled_distance += TINY_BIT;
}
int k_max;
traveled_distance =
distance_to_mesh(ijk, distances, traveled_distance, k_max);

// If r1 is not inside the mesh, exit here
if (traveled_distance >= total_distance)
Expand All @@ -1287,6 +1276,29 @@ void StructuredMesh::raytrace_mesh(
}
}

double StructuredMesh::distance_to_mesh(const MeshIndex& ijk,
const std::array<MeshDistance, MESH_MAX_AXES>& distances,
double traveled_distance, int& k_max) const
{
// For all directions outside the mesh, find the distance that we need
// to travel to reach the next surface. Use the largest distance, as
// only this will cross all outer surfaces.
const int n = n_dimension_;
k_max = -1;
for (int k = 0; k < n; ++k) {
if ((ijk[k] < 1 || ijk[k] > shape_[k]) &&
(distances[k].distance > traveled_distance)) {
traveled_distance = distances[k].distance;
k_max = k;
}
}
// Assure some distance is traveled
if (k_max == -1) {
traveled_distance += TINY_BIT;
}
return traveled_distance;
}

void StructuredMesh::bins_crossed(Position r0, Position r1, const Direction& u,
vector<int>& bins, vector<double>& lengths) const
{
Expand Down Expand Up @@ -1527,16 +1539,15 @@ StructuredMesh::MeshDistance RegularMesh::distance_to_grid_boundary(
double l) const
{
MeshDistance d;
d.next_index = ijk[i];
if (std::abs(u[i]) < FP_PRECISION)
return d;

d.max_surface = (u[i] > 0);
if (d.max_surface && (ijk[i] <= shape_[i])) {
d.next_index++;
++d.offset[i];
d.distance = (positive_grid_boundary(ijk, i) - r0[i]) / u[i];
} else if (!d.max_surface && (ijk[i] >= 1)) {
d.next_index--;
--d.offset[i];
d.distance = (negative_grid_boundary(ijk, i) - r0[i]) / u[i];
}

Expand Down Expand Up @@ -1699,16 +1710,15 @@ StructuredMesh::MeshDistance RectilinearMesh::distance_to_grid_boundary(
double l) const
{
MeshDistance d;
d.next_index = ijk[i];
if (std::abs(u[i]) < FP_PRECISION)
return d;

d.max_surface = (u[i] > 0);
if (d.max_surface && (ijk[i] <= shape_[i])) {
d.next_index++;
++d.offset[i];
d.distance = (positive_grid_boundary(ijk, i) - r0[i]) / u[i];
} else if (!d.max_surface && (ijk[i] > 0)) {
d.next_index--;
--d.offset[i];
d.distance = (negative_grid_boundary(ijk, i) - r0[i]) / u[i];
}
return d;
Expand Down Expand Up @@ -1966,18 +1976,16 @@ StructuredMesh::MeshDistance CylindricalMesh::find_z_crossing(
const Position& r, const Direction& u, double l, int shell) const
{
MeshDistance d;
d.next_index = shell;

// Direction of flight is within xy-plane. Will never intersect z.
if (std::abs(u.z) < FP_PRECISION)
return d;

d.max_surface = (u.z > 0.0);
if (d.max_surface && (shell <= shape_[2])) {
d.next_index += 1;
++d.offset[2];
d.distance = (grid_[2][shell] - r.z) / u.z;
} else if (!d.max_surface && (shell > 0)) {
d.next_index -= 1;
--d.offset[2];
d.distance = (grid_[2][shell - 1] - r.z) / u.z;
}
return d;
Expand All @@ -1990,15 +1998,14 @@ StructuredMesh::MeshDistance CylindricalMesh::distance_to_grid_boundary(
if (i == 0) {

return std::min(
MeshDistance(ijk[i] + 1, true, find_r_crossing(r0, u, l, ijk[i])),
MeshDistance(ijk[i] - 1, false, find_r_crossing(r0, u, l, ijk[i] - 1)));
MeshDistance({1, 0, 0}, true, find_r_crossing(r0, u, l, ijk[i])),
MeshDistance({-1, 0, 0}, false, find_r_crossing(r0, u, l, ijk[i] - 1)));

} else if (i == 1) {

return std::min(MeshDistance(sanitize_phi(ijk[i] + 1), true,
find_phi_crossing(r0, u, l, ijk[i])),
MeshDistance(sanitize_phi(ijk[i] - 1), false,
find_phi_crossing(r0, u, l, ijk[i] - 1)));
return std::min(
MeshDistance({0, 1, 0}, true, find_phi_crossing(r0, u, l, ijk[i])),
MeshDistance({0, -1, 0}, false, find_phi_crossing(r0, u, l, ijk[i] - 1)));

} else {
return find_z_crossing(r0, u, l, ijk[i]);
Expand Down Expand Up @@ -2322,20 +2329,19 @@ StructuredMesh::MeshDistance SphericalMesh::distance_to_grid_boundary(

if (i == 0) {
return std::min(
MeshDistance(ijk[i] + 1, true, find_r_crossing(r0, u, l, ijk[i])),
MeshDistance(ijk[i] - 1, false, find_r_crossing(r0, u, l, ijk[i] - 1)));
MeshDistance({1, 0, 0}, true, find_r_crossing(r0, u, l, ijk[i])),
MeshDistance({-1, 0, 0}, false, find_r_crossing(r0, u, l, ijk[i] - 1)));

} else if (i == 1) {
return std::min(MeshDistance(sanitize_theta(ijk[i] + 1), true,
find_theta_crossing(r0, u, l, ijk[i])),
MeshDistance(sanitize_theta(ijk[i] - 1), false,
find_theta_crossing(r0, u, l, ijk[i] - 1)));
return std::min(
MeshDistance({0, 1, 0}, true, find_theta_crossing(r0, u, l, ijk[i])),
MeshDistance(
{0, -1, 0}, false, find_theta_crossing(r0, u, l, ijk[i] - 1)));

} else {
return std::min(MeshDistance(sanitize_phi(ijk[i] + 1), true,
find_phi_crossing(r0, u, l, ijk[i])),
MeshDistance(sanitize_phi(ijk[i] - 1), false,
find_phi_crossing(r0, u, l, ijk[i] - 1)));
return std::min(
MeshDistance({0, 0, 1}, true, find_phi_crossing(r0, u, l, ijk[i])),
MeshDistance({0, 0, -1}, false, find_phi_crossing(r0, u, l, ijk[i] - 1)));
}
}

Expand Down
Loading