diff --git a/src/pe/ccd/HashGrids.cpp b/src/pe/ccd/HashGrids.cpp
index 3002eb6f62c4f6ca97b8fe7ef90a45d0523806e4..d5ca576cd453d015087956a11b9df628574cc0ab 100644
--- a/src/pe/ccd/HashGrids.cpp
+++ b/src/pe/ccd/HashGrids.cpp
@@ -363,12 +363,39 @@ size_t HashGrids::HashGrid::hashPoint(real_t x, real_t y, real_t z) const {
zHash = static_cast<size_t>( i ) & zHashMask_;
}
- WALBERLA_LOG_INFO("\t xH: " << xHash << ", yH: " << yHash << ", zH: " << zHash);
-
return xHash + yHash * xCellCount_ + zHash * xyCellCount_;
}
-void HashGrids::HashGrid::possibleRayIntersectingBodies(const raytracing::Ray& ray, const AABB& blockAABB) const {
+
+void HashGrids::HashGrid::insertRelatedCellIndicesForCenter(real_t x, real_t y, real_t z, std::unordered_set<size_t>& cellIndices, const AABB& blockAABB) const {
+ // cellIndices.insert(hashPoint(x, y, z)); // <- happens already in the loop down there -v
+
+ const Vec3& minCorner = blockAABB.minCorner();
+ const Vec3& maxCorner = blockAABB.maxCorner();
+ for (int i = -1; i <= 0; ++i) {
+ for (int j = -1; j <= 0; ++j) {
+ for (int k = -1; k <= 0; ++k) {
+ real_t x_shifted = x + i*cellSpan_;
+ real_t y_shifted = y + j*cellSpan_;
+ real_t z_shifted = z + k*cellSpan_;
+ if (x_shifted > minCorner[0] && y_shifted > minCorner[1] && z_shifted > minCorner[2] &&
+ x_shifted < maxCorner[0] && y_shifted < maxCorner[1] && z_shifted < maxCorner[2]) {
+ cellIndices.insert(hashPoint(x_shifted, y_shifted, z_shifted));
+ }
+ }
+ }
+ }
+}
+
+void HashGrids::HashGrid::possibleRayIntersectingBodies(const raytracing::Ray& ray, const AABB& blockAABB) const {
+ std::unordered_set<size_t> cellIndices;
+
+ const Vec3& direction = ray.getDirection();
+
+ real_t xCenterOffsetFactor = direction[0] < 0 ? -0.5 : 0.5;
+ real_t yCenterOffsetFactor = direction[1] < 0 ? -0.5 : 0.5;
+ real_t zCenterOffsetFactor = direction[2] < 0 ? -0.5 : 0.5;
+
for (int i = 0; i <= xCellCount_; i++) {
real_t xValue = blockAABB.minCorner()[0] + i*cellSpan_;
real_t lambda = (xValue - ray.getOrigin()[0]) * ray.getInvDirection()[0];
@@ -379,8 +406,15 @@ void HashGrids::HashGrid::possibleRayIntersectingBodies(const raytracing::Ray& r
lambda != lambda) {
WALBERLA_LOG_INFO("P_x" << i << " = (" << xValue << "/" << yValue << "/" << zValue << ") invalid");
} else {
- size_t arrayIndex = hashPoint(xValue, yValue, zValue);
- WALBERLA_LOG_INFO("P_x" << i << " = (" << xValue << "/" << yValue << "/" << zValue << ") maps to cell " << arrayIndex);
+ real_t centeredXValue = xValue + cellSpan_*xCenterOffsetFactor;
+ if (centeredXValue > blockAABB.xMax() || centeredXValue < blockAABB.xMin()) {
+ continue;
+ }
+ real_t centeredYValue = (size_t((yValue - blockAABB.yMin())*inverseCellSpan_) + real_t(0.5)) * cellSpan_ + blockAABB.yMin();
+ real_t centeredZValue = (size_t((zValue - blockAABB.zMin())*inverseCellSpan_) + real_t(0.5)) * cellSpan_ + blockAABB.zMin();
+ size_t arrayIndex = hashPoint(centeredXValue, centeredYValue, centeredZValue);
+ WALBERLA_LOG_INFO("P_x" << i << " = (" << xValue << "/" << yValue << "/" << zValue << ") maps to box center (" << centeredXValue << "/" << centeredYValue << "/" << centeredZValue << ") and cell " << arrayIndex);
+ insertRelatedCellIndicesForCenter(centeredXValue, centeredYValue, centeredZValue, cellIndices, blockAABB);
}
}
@@ -389,13 +423,20 @@ void HashGrids::HashGrid::possibleRayIntersectingBodies(const raytracing::Ray& r
real_t lambda = (yValue - ray.getOrigin()[1]) * ray.getInvDirection()[1];
real_t xValue = ray.getOrigin()[0] + lambda * ray.getDirection()[0];
real_t zValue = ray.getOrigin()[2] + lambda * ray.getDirection()[2];
- if (xValue > blockAABB.maxCorner()[0] || xValue < blockAABB.minCorner()[0] ||
- zValue > blockAABB.maxCorner()[2] || zValue < blockAABB.minCorner()[2] ||
+ if (xValue >= blockAABB.maxCorner()[0] || xValue < blockAABB.minCorner()[0] ||
+ zValue >= blockAABB.maxCorner()[2] || zValue < blockAABB.minCorner()[2] ||
lambda != lambda) {
WALBERLA_LOG_INFO("P_y" << i << " = (" << xValue << "/" << yValue << "/" << zValue << ") invalid");
} else {
- size_t arrayIndex = hashPoint(xValue, yValue, zValue);
- WALBERLA_LOG_INFO("P_y" << i << " = (" << xValue << "/" << yValue << "/" << zValue << ") maps to cell " << arrayIndex);
+ real_t centeredXValue = (size_t((xValue - blockAABB.xMin())*inverseCellSpan_) + real_t(0.5)) * cellSpan_ + blockAABB.xMin();
+ real_t centeredYValue = yValue + cellSpan_*yCenterOffsetFactor;
+ if (centeredYValue > blockAABB.yMax() || centeredYValue < blockAABB.yMin()) {
+ continue;
+ }
+ real_t centeredZValue = (size_t((zValue - blockAABB.zMin())*inverseCellSpan_) + real_t(0.5)) * cellSpan_ + blockAABB.zMin();
+ size_t arrayIndex = hashPoint(centeredXValue, centeredYValue, centeredZValue);
+ WALBERLA_LOG_INFO("P_y" << i << " = (" << xValue << "/" << yValue << "/" << zValue << ") maps to box center (" << centeredXValue << "/" << centeredYValue << "/" << centeredZValue << ") and cell " << arrayIndex);
+ insertRelatedCellIndicesForCenter(centeredXValue, centeredYValue, centeredZValue, cellIndices, blockAABB);
}
}
@@ -409,10 +450,23 @@ void HashGrids::HashGrid::possibleRayIntersectingBodies(const raytracing::Ray& r
lambda != lambda) {
WALBERLA_LOG_INFO("P_z" << i << " = (" << xValue << "/" << yValue << "/" << zValue << ") invalid");
} else {
- size_t arrayIndex = hashPoint(xValue, yValue, zValue);
- WALBERLA_LOG_INFO("P_z" << i << " = (" << xValue << "/" << yValue << "/" << zValue << ") maps to cell " << arrayIndex);
+ real_t centeredXValue = (size_t((xValue - blockAABB.xMin())*inverseCellSpan_) + real_t(0.5)) * cellSpan_ + blockAABB.xMin();
+ real_t centeredYValue = (size_t((yValue - blockAABB.yMin())*inverseCellSpan_) + real_t(0.5)) * cellSpan_ + blockAABB.yMin();
+ real_t centeredZValue = zValue + cellSpan_*zCenterOffsetFactor;
+ if (centeredZValue > blockAABB.zMax() || centeredZValue < blockAABB.zMin()) {
+ continue;
+ }
+ size_t arrayIndex = hashPoint(centeredXValue, centeredYValue, centeredZValue);
+ WALBERLA_LOG_INFO("P_z" << i << " = (" << xValue << "/" << yValue << "/" << zValue << ") maps to box center (" << centeredXValue << "/" << centeredYValue << "/" << centeredZValue << ") and cell " << arrayIndex);
+ insertRelatedCellIndicesForCenter(centeredXValue, centeredYValue, centeredZValue, cellIndices, blockAABB);
}
}
+
+ WALBERLA_LOG_INFO("Resulting indices:");
+ for (size_t const& index : cellIndices) {
+ std::cout << index << " ";
+ }
+ std::cout << std::endl;
}
diff --git a/src/pe/ccd/HashGrids.h b/src/pe/ccd/HashGrids.h
index 67bbcfe62ab92cf0c118a872b7e29328b0590ea5..5fc5ff175344b5545d2f709ef73cc75aa65b32dd 100644
--- a/src/pe/ccd/HashGrids.h
+++ b/src/pe/ccd/HashGrids.h
@@ -40,6 +40,7 @@
#include <sstream>
#include <vector>
+#include <unordered_set>
#include <pe/raytracing/Ray.h>
namespace walberla{
@@ -182,7 +183,7 @@ public: //ToDo fix to private again
//@}
//*******************************************************************************************
-
+ void insertRelatedCellIndicesForCenter(real_t x, real_t y, real_t z, std::unordered_set<size_t>& cellIndices, const AABB& blockAABB) const;
void possibleRayIntersectingBodies(const raytracing::Ray& ray, const AABB& blockAABB) const;
size_t hashPoint(real_t x, real_t y, real_t z) const;
diff --git a/tests/pe/Raytracing.cpp b/tests/pe/Raytracing.cpp
index c745d186ee2e15e691086f3ee041e8317f1bcf39..afc8cfc99e4421d02da796d060922dcad449b7ae 100644
--- a/tests/pe/Raytracing.cpp
+++ b/tests/pe/Raytracing.cpp
@@ -516,6 +516,10 @@ void hashgridsPlayground() {
Ray ray2(Vec3(-2.1, -2.1, 0), Vec3(1, 1, 0).getNormalized());
hashgrid.possibleRayIntersectingBodies(ray2, blockAABB);
+ WALBERLA_LOG_INFO("ray3:");
+ Ray ray3(Vec3(3, -1, 0), Vec3(-7, 2, 0).getNormalized());
+ hashgrid.possibleRayIntersectingBodies(ray3, blockAABB);
+
WALBERLA_LOG_INFO("");
WALBERLA_LOG_INFO(hashgrid.hashPoint(-0.5, -0.5, 0));
WALBERLA_LOG_INFO(hashgrid.hashPoint(-1.5, -1.5, 0));