diff --git a/src/pe/ccd/HashGrids.cpp b/src/pe/ccd/HashGrids.cpp
index b9d170c2e6d0c3bd1cc920357caa6965729fd8ca..21f63787f9eea0ced53c4a673d71dbe1da7578e2 100644
--- a/src/pe/ccd/HashGrids.cpp
+++ b/src/pe/ccd/HashGrids.cpp
@@ -503,6 +503,7 @@ void HashGrids::HashGrid::remove( BodyID body, Cell* cell )
*/
void HashGrids::HashGrid::enlarge()
{
+ WALBERLA_LOG_INFO("enlarging grid"); // ToDo remove again
BodyID* bodies = new BodyID[ bodyCount_ ];
BodyID* body = bodies;
@@ -1202,7 +1203,7 @@ size_t HashGrids::intersectionTestCount = 0; //ToDo remove again
* note that the initial number of cells does not necessarily have to be equal for all three
* coordinate directions.
*/
-const size_t HashGrids::xCellCount = 16;
+const size_t HashGrids::xCellCount = 4; //ToDo change to 16 again
//*************************************************************************************************
@@ -1219,7 +1220,7 @@ const size_t HashGrids::xCellCount = 16;
* note that the initial number of cells does not necessarily have to be equal for all three
* coordinate directions.
*/
-const size_t HashGrids::yCellCount = 16;
+const size_t HashGrids::yCellCount = 4; //ToDo change to 16 again
//*************************************************************************************************
@@ -1236,7 +1237,7 @@ const size_t HashGrids::yCellCount = 16;
* note that the initial number of cells does not necessarily have to be equal for all three
* coordinate directions.
*/
-const size_t HashGrids::zCellCount = 16;
+const size_t HashGrids::zCellCount = 4; //ToDo change to 16 again
//*************************************************************************************************
@@ -1301,7 +1302,7 @@ const size_t HashGrids::minimalGridDensity = 8;
*
* Possible settings: any integral value greater-or-equal to 0.
*/
-const size_t HashGrids::gridActivationThreshold = 32;
+const size_t HashGrids::gridActivationThreshold = 0; // ToDo change to 32 again
//*************************************************************************************************
diff --git a/src/pe/ccd/HashGrids.h b/src/pe/ccd/HashGrids.h
index f7a3d3346477361b75bdc7712f5446ecec1cd78f..b56cdf71c66585852235464b3f73c45b75eab5be 100644
--- a/src/pe/ccd/HashGrids.h
+++ b/src/pe/ccd/HashGrids.h
@@ -186,6 +186,15 @@ public:
template<typename BodyTuple>
BodyID getRayIntersectingBody(const raytracing::Ray& ray, const AABB& blockAABB, real_t& t, Vec3& n) const;
+ template<typename BodyTuple>
+ BodyID getRayIntersectingBodyOLD(const raytracing::Ray& ray, const AABB& blockAABB,
+ real_t& t_closest, Vec3& n_closest) const;
+
+ template<typename BodyTuple>
+ bool getBodyIntersectionForBlockCell(const Vector3<int32_t>& blockCell, const AABB& blockAABB,
+ const raytracing::Ray& ray,
+ real_t& t, Vec3& n, BodyID& body) const;
+
void clear();
//@}
//*******************************************************************************************
@@ -209,6 +218,11 @@ public:
bool getBodyIntersectionForCellCenter(real_t x, real_t y, real_t z, const AABB& blockAABB,
const raytracing::Ray& ray, size_t c,
real_t& t, Vec3& n, BodyID& body) const;
+ template<typename BodyTuple>
+ bool getBodyIntersectionForCellBetweenPoints(const Vec3& firstPoint, const Vec3& secondPoint,
+ const AABB& blockAABB,
+ const raytracing::Ray& ray, size_t c,
+ real_t& t, Vec3& n, BodyID& body) const;
//@}
//*******************************************************************************************
@@ -302,7 +316,7 @@ public:
BodyID getClosestBodyIntersectingWithRay(const raytracing::Ray& ray, const AABB& blockAABB,
real_t& t, Vec3& n);
-protected:
+//protected: // ToDo uncomment to change to protected again
//**Utility functions***************************************************************************
/*!\name Utility functions */
//@{
@@ -311,7 +325,7 @@ protected:
//@}
//**********************************************************************************************
-private:
+//private: // ToDo uncomment to change to private again
//**Add functions*******************************************************************************
/*!\name Add functions */
//@{
@@ -511,60 +525,84 @@ template<typename BodyTuple>
bool HashGrids::HashGrid::getBodyIntersectionForCellCenter(real_t x, real_t y, real_t z, const AABB& blockAABB,
const raytracing::Ray& ray, size_t c,
real_t& t, Vec3& n, BodyID& body) const {
- //const Vec3& minCorner = blockAABB.minCorner();
- //const Vec3& maxCorner = blockAABB.maxCorner();
+ real_t t_local;
+ Vec3 n_local;
+ bool intersected = false;
+
+ raytracing::IntersectsFunctor intersectsFunc(ray, t_local, n_local);
+
+ const Cell& centerCell = cell_[hashPoint(x, y, z)];
+
+ for (unsigned int i = 0; i < 27; ++i) {
+ const Cell* nbCell = ¢erCell + centerCell.neighborOffset_[i];
+ const BodyVector* nbBodies = nbCell->bodies_;
+
+ if (nbBodies != NULL) {
+ for (const BodyID& cellBody: *nbBodies) {
+ HashGrids::intersectionTestCount++;
+ bool intersects = SingleCast<BodyTuple, raytracing::IntersectsFunctor, bool>::execute(cellBody, intersectsFunc);
+ if (intersects && t_local < t) {
+ body = cellBody;
+ t = t_local;
+ n = n_local;
+ intersected = true;
+ }
+ }
+ }
+ }
+
+ return intersected;
+}
+
+/*!\brief Computes closest ray-body intersection of cell with center at point x,y,z and neighboring ones.
+ *
+ * \param blockCell index of cell within block grid.
+ * \param blockAABB AABB of the block this grid corresponds to.
+ * \param ray Ray being casted trough grid.
+ * \param t Reference for calculated distance from ray origin.
+ * \param n Reference for normal of intersection point.
+ * \param body Reference for closest body.
+ *
+ * Inserts bodies of specified cell into bodiesContainer and additionally considers bodies in neighboring cells
+ * in all negative coordinate direction to take bodies in consideration, which protrude from their
+ * cell into the intersected cell (and thus, possibly intersect with the ray as well).
+ */
+template<typename BodyTuple>
+bool HashGrids::HashGrid::getBodyIntersectionForBlockCell(const Vector3<int32_t>& blockCell, const AABB& blockAABB,
+ const raytracing::Ray& ray,
+ real_t& t, Vec3& n, BodyID& body) const {
+ const real_t inf = std::numeric_limits<real_t>::max();
+
+ real_t t_closest = inf;
real_t t_local;
Vec3 n_local;
bool intersected = false;
raytracing::IntersectsFunctor intersectsFunc(ray, t_local, n_local);
+
+ real_t x = (real_c(blockCell[0]) + real_t(0.5)) * cellSpan_;
+ real_t y = (real_c(blockCell[1]) + real_t(0.5)) * cellSpan_;
+ real_t z = (real_c(blockCell[2]) + real_t(0.5)) * cellSpan_;
- //const Vec3& dir = ray.getDirection();
-
- int minX = -1, minY = -1, minZ = -1;
- int maxX = 0, maxY = 0, maxZ = 0;
-
- /*if (c == 0) {
- minX = -2;
- maxY = 1;
- maxZ = 1;
- } else if (c == 1) {
- minY = -2;
- maxX = 1;
- maxZ = 1;
- } else if (c == 2) {
- minZ = -2;
- maxY = 1;
- maxX = 1;
- }*/
-
- for (int i = minX; i <= maxX; ++i) {
- const real_t x_shifted = x + i*cellSpan_;
- for (int j = minY; j <= maxY; ++j) {
- const real_t y_shifted = y + j*cellSpan_;
- for (int k = minZ; k <= maxZ; ++k) {
- const real_t z_shifted = z + k*cellSpan_;
- // commenting upper line of condition fixes a bug where objects with their minCorner
- // out of the blocks bounds would not get considered for intersections.
- //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]*/ true) {
- size_t hash = hashPoint(x_shifted, y_shifted, z_shifted);
-
- const Cell& cell = cell_[hash];
- if (cell.bodies_ != NULL) {
- for (const BodyID& cellBody: *cell.bodies_) {
- HashGrids::intersectionTestCount++;
- bool intersects = SingleCast<BodyTuple, raytracing::IntersectsFunctor, bool>::execute(cellBody, intersectsFunc);
- if (intersects && t_local < t) {
- body = cellBody;
- t = t_local;
- n = n_local;
- intersected = true;
- }
- }
- }
- //}
+ const Cell& centerCell = cell_[hashPoint(x, y, z)];
+
+ //WALBERLA_LOG_INFO("checking cell with index " << hashPoint(x, y, z));
+
+ for (unsigned int i = 0; i < 27; ++i) {
+ const Cell* nbCell = ¢erCell + centerCell.neighborOffset_[i];
+ const BodyVector* nbBodies = nbCell->bodies_;
+
+ if (nbBodies != NULL) {
+ for (const BodyID& cellBody: *nbBodies) {
+ HashGrids::intersectionTestCount++;
+ bool intersects = SingleCast<BodyTuple, raytracing::IntersectsFunctor, bool>::execute(cellBody, intersectsFunc);
+ if (intersects && t_local < t_closest) {
+ body = cellBody;
+ t = t_closest = t_local;
+ n = n_local;
+ intersected = true;
+ }
}
}
}
@@ -572,7 +610,7 @@ bool HashGrids::HashGrid::getBodyIntersectionForCellCenter(real_t x, real_t y, r
return intersected;
}
-/*!\brief Calculates ray-cell intersections and determines the closest body from the ray origin.
+/*!\brief Calculates ray-cell intersections and determines the closest body to the ray origin.
*
* \param ray Ray getting shot through this grid.
* \param blockAABB AABB of the block this grid corresponds to.
@@ -587,155 +625,120 @@ bool HashGrids::HashGrid::getBodyIntersectionForCellCenter(real_t x, real_t y, r
template<typename BodyTuple>
BodyID HashGrids::HashGrid::getRayIntersectingBody(const raytracing::Ray& ray, const AABB& blockAABB,
real_t& t_closest, Vec3& n_closest) const {
- //real_t inf = std::numeric_limits<real_t>::max();
+ const real_t inf = std::numeric_limits<real_t>::max();
+
+ BodyID body_closest = NULL;
+ real_t t_local_closest = inf;
+ Vec3 n_local_closest;
+ BodyID body = NULL;
+ real_t t;
+ Vec3 n;
- const Vec3& rayDirection = ray.getDirection();
- const Vec3& minCorner = blockAABB.minCorner();
- const Vec3& maxCorner = blockAABB.maxCorner();
- const Vec3& rayOrigin = ray.getOrigin();
- const Vec3& rayInvDirection = ray.getInvDirection();
+ int32_t blockXCellCountMin = int32_c(blockAABB.xMin() * inverseCellSpan_) - 1;
+ int32_t blockXCellCountMax = int32_c(std::ceil(blockAABB.xMax() * inverseCellSpan_)) + 1;
+ int32_t blockYCellCountMin = int32_c(blockAABB.yMin() * inverseCellSpan_) - 1;
+ int32_t blockYCellCountMax = int32_c(std::ceil(blockAABB.yMax() * inverseCellSpan_)) + 1;
+ int32_t blockZCellCountMin = int32_c(blockAABB.zMin() * inverseCellSpan_) - 1;
+ int32_t blockZCellCountMax = int32_c(std::ceil(blockAABB.zMax() * inverseCellSpan_)) + 1;
+
+ //WALBERLA_LOG_INFO("cellspan: " << cellSpan_ << " lims: " << blockXCellCountMin << " " << blockXCellCountMax);
- real_t xCenterOffsetFactor = rayDirection[0] < 0 ? -0.5 : 0.5;
- real_t yCenterOffsetFactor = rayDirection[1] < 0 ? -0.5 : 0.5;
- real_t zCenterOffsetFactor = rayDirection[2] < 0 ? -0.5 : 0.5;
- //real_t t_closest = inf;
- //Vec3 n_closest;
- BodyID body_closest = NULL;
-
- int i_first_intersection = -1;
- for (size_t i = 0; i <= xCellCount_; i++) {
- size_t i_dir = i;
- if (rayDirection[0] < 0) {
- i_dir = xCellCount_ - i;
- }
- real_t xValue = minCorner[0] + i_dir*cellSpan_; // calculate x value
- // get lambda in ray equation for which ray intersects with yz plane with this x value
- real_t lambda = (xValue - rayOrigin[0]) * rayInvDirection[0];
- // lambda is the distance from the ray origin to the cell intersection point
- if (lambda < 0) {
- // intersection is in rays negative direction
- continue;
- }
- if (rayDirection[0] >= 0 && lambda > t_closest) {
- // cell is too far away to generate useful body intersections
- // this condition only works for checks in positive direction
- break;
- }
- if (rayDirection[0] < 0 && lambda > t_closest+cellSpan_) {
- break;
- }
- real_t yValue = rayOrigin[1] + lambda * rayDirection[1]; // get y and z values for this intersection
- real_t zValue = rayOrigin[2] + lambda * rayDirection[2];
- if (yValue > maxCorner[1]+cellSpan_ || yValue < minCorner[1]-cellSpan_ ||
- zValue > maxCorner[2]+cellSpan_ || zValue < minCorner[2]-cellSpan_ ||
- lambda != lambda) {
- // calculated intersection with yz plane is out of the blocks bounds
- continue;
- }
- real_t centeredXValue = xValue + cellSpan_*xCenterOffsetFactor;
- if (centeredXValue > maxCorner[0] || centeredXValue < minCorner[0]) {
- // intersection was with one of the outer bounds of the block and in this direction no more cells exist
- continue;
- }
- // calculate the y and z values of the center of the cell
- real_t centeredYValue = (int((yValue - minCorner[1])*inverseCellSpan_) + real_t(0.5)) * cellSpan_ + minCorner[1];
- real_t centeredZValue = (int((zValue - minCorner[2])*inverseCellSpan_) + real_t(0.5)) * cellSpan_ + minCorner[2];
- bool intersected = getBodyIntersectionForCellCenter<BodyTuple>(centeredXValue, centeredYValue, centeredZValue, blockAABB,
- ray, 0,
- t_closest, n_closest, body_closest);
- if (intersected && i_first_intersection == -1) {
- i_first_intersection = int(i);
- }
- if (i_first_intersection != -1 && i_first_intersection == i-2) {
- break;
+ Vec3 firstPoint;
+ if (blockAABB.contains(ray.getOrigin(), cellSpan_)) {
+ firstPoint = ray.getOrigin();
+ } else {
+ real_t t_start;
+ if (intersects(blockAABB, ray, t_start, cellSpan_)) {
+ firstPoint = ray.getOrigin() + ray.getDirection()*t_start;
+ } else {
+ //WALBERLA_LOG_INFO("ray does not intersect block")
+ return NULL;
}
}
- i_first_intersection = -1;
- for (size_t i = 0; i < yCellCount_; i++) {
- size_t i_dir = i;
- if (rayDirection[1] < 0) {
- i_dir = yCellCount_ - i;
- }
- real_t yValue = minCorner[1] + i_dir*cellSpan_;
- real_t lambda = (yValue - rayOrigin[1]) * rayInvDirection[1];
- if (lambda < 0) {
- continue;
- }
- if (rayDirection[1] >= 0 && lambda > t_closest) {
- break;
- }
- if (rayDirection[1] < 0 && lambda > t_closest+cellSpan_) {
- break;
- }
- real_t xValue = rayOrigin[0] + lambda * rayDirection[0];
- real_t zValue = rayOrigin[2] + lambda * rayDirection[2];
- if (xValue >= maxCorner[0]+cellSpan_ || xValue < minCorner[0]-cellSpan_ ||
- zValue >= maxCorner[2]+cellSpan_ || zValue < minCorner[2]-cellSpan_ ||
- lambda != lambda) {
- continue;
- }
- real_t centeredXValue = (int((xValue - minCorner[0])*inverseCellSpan_) + real_t(0.5)) * cellSpan_ + minCorner[0];
- real_t centeredYValue = yValue + cellSpan_*yCenterOffsetFactor;
- if (centeredYValue > maxCorner[1] || centeredYValue < minCorner[1]) {
- continue;
- }
- real_t centeredZValue = (int((zValue - minCorner[2])*inverseCellSpan_) + real_t(0.5)) * cellSpan_ + minCorner[2];
- bool intersected = getBodyIntersectionForCellCenter<BodyTuple>(centeredXValue, centeredYValue, centeredZValue, blockAABB,
- ray, 1,
- t_closest, n_closest, body_closest);
- if (intersected && i_first_intersection == -1) {
- i_first_intersection = int(i);
- }
- if (i_first_intersection != -1 && i_first_intersection == i-2) {
- break;
+ Vector3<int32_t> firstCell(int32_c(std::floor(firstPoint[0]*inverseCellSpan_)),
+ int32_c(std::floor(firstPoint[1]*inverseCellSpan_)),
+ int32_c(std::floor(firstPoint[2]*inverseCellSpan_)));
+
+ const int8_t stepX = ray.xDir() >= 0 ? 1 : -1;
+ const int8_t stepY = ray.yDir() >= 0 ? 1 : -1;
+ const int8_t stepZ = ray.zDir() >= 0 ? 1 : -1;
+
+ Vec3 near((stepX >= 0) ? (firstCell[0]+1)*cellSpan_-firstPoint[0] : firstPoint[0]-firstCell[0]*cellSpan_,
+ (stepY >= 0) ? (firstCell[1]+1)*cellSpan_-firstPoint[1] : firstPoint[1]-firstCell[1]*cellSpan_,
+ (stepZ >= 0) ? (firstCell[2]+1)*cellSpan_-firstPoint[2] : firstPoint[2]-firstCell[2]*cellSpan_);
+
+ real_t tMaxX = (ray.xDir() != 0) ? std::abs(near[0]*ray.xInvDir()) : inf;
+ real_t tMaxY = (ray.yDir() != 0) ? std::abs(near[1]*ray.yInvDir()) : inf;
+ real_t tMaxZ = (ray.zDir() != 0) ? std::abs(near[2]*ray.zInvDir()) : inf;
+
+ real_t tDeltaX = (ray.xDir() != 0) ? std::abs(cellSpan_*ray.xInvDir()) : inf;
+ real_t tDeltaY = (ray.yDir() != 0) ? std::abs(cellSpan_*ray.yInvDir()) : inf;
+ real_t tDeltaZ = (ray.zDir() != 0) ? std::abs(cellSpan_*ray.zInvDir()) : inf;
+
+ Vector3<int32_t> currentCell = firstCell;
+
+ // First cell needs extra treatment, as it might lay out of the blocks upper bounds
+ // due to the nature of how it is calculated: If the first point lies on a upper border
+ // it maps to the cell "above" the grid.
+ if (currentCell[0] < blockXCellCountMax &&
+ currentCell[1] < blockYCellCountMax &&
+ currentCell[2] < blockZCellCountMax) {
+ //WALBERLA_LOG_INFO("found block cell at " << currentCell);
+ bool intersected = getBodyIntersectionForBlockCell<BodyTuple>(currentCell, blockAABB, ray,
+ t, n, body);
+ if (intersected && t < t_local_closest) {
+ body_closest = body;
+ t_local_closest = t;
+ n_local_closest = n;
}
}
- i_first_intersection = -1;
- for (size_t i = 0; i < zCellCount_; i++) {
- size_t i_dir = i;
- if (rayDirection[2] < 0) {
- i_dir = zCellCount_ - i;
- }
- real_t zValue = minCorner[2] + i_dir*cellSpan_;
- real_t lambda = (zValue - rayOrigin[2]) * rayInvDirection[2];
- if (lambda < 0) {
- continue;
- }
- if (rayDirection[2] >= 0 && lambda > t_closest) {
- break;
- }
- if (rayDirection[2] < 0 && lambda > t_closest+cellSpan_) {
- break;
- }
- real_t xValue = rayOrigin[0] + lambda * rayDirection[0];
- real_t yValue = rayOrigin[1] + lambda * rayDirection[1];
- if (xValue > maxCorner[0]+cellSpan_ || xValue < minCorner[0]-cellSpan_ ||
- yValue > maxCorner[1]+cellSpan_ || yValue < minCorner[1]-cellSpan_ ||
- lambda != lambda) {
- continue;
- }
- real_t centeredXValue = (int((xValue - minCorner[0])*inverseCellSpan_) + real_t(0.5)) * cellSpan_ + minCorner[0];
- real_t centeredYValue = (int((yValue - minCorner[1])*inverseCellSpan_) + real_t(0.5)) * cellSpan_ + minCorner[1];
- real_t centeredZValue = zValue + cellSpan_*zCenterOffsetFactor;
- if (centeredZValue > maxCorner[2] || centeredZValue < minCorner[2]) {
- continue;
- }
- bool intersected = getBodyIntersectionForCellCenter<BodyTuple>(centeredXValue, centeredYValue, centeredZValue, blockAABB,
- ray, 2,
- t_closest, n_closest, body_closest);
- if (intersected && i_first_intersection == -1) {
- i_first_intersection = int(i);
+ while (true) {
+ if (tMaxX < tMaxY) {
+ if (tMaxX < tMaxZ) {
+ tMaxX += tDeltaX;
+ currentCell[0] += stepX;
+ if (currentCell[0] >= blockXCellCountMax || currentCell[0] <= blockXCellCountMin) {
+ break;
+ }
+ } else {
+ tMaxZ += tDeltaZ;
+ currentCell[2] += stepZ;
+ if (currentCell[2] >= blockZCellCountMax || currentCell[2] <= blockZCellCountMin) {
+ break;
+ }
+ }
+ } else {
+ if (tMaxY < tMaxZ) {
+ tMaxY += tDeltaY;
+ currentCell[1] += stepY;
+ if (currentCell[1] >= blockYCellCountMax || currentCell[1] <= blockYCellCountMin) {
+ break;
+ }
+ } else {
+ tMaxZ += tDeltaZ;
+ currentCell[2] += stepZ;
+ if (currentCell[2] >= blockZCellCountMax || currentCell[2] <= blockZCellCountMin) {
+ break;
+ }
+ }
}
- if (i_first_intersection != -1 && i_first_intersection == i-2) {
- break;
+
+ //WALBERLA_LOG_INFO("found block cell at " << currentCell);
+
+ bool intersected = getBodyIntersectionForBlockCell<BodyTuple>(currentCell, blockAABB, ray,
+ t, n, body);
+ if (intersected && t < t_local_closest) {
+ body_closest = body;
+ t_local_closest = t;
+ n_local_closest = n;
}
}
- //n = n_closest;
- //t = t_closest;
+ t_closest = t_local_closest;
+ n_closest = n_local_closest;
return body_closest;
}
@@ -760,10 +763,16 @@ BodyID HashGrids::getClosestBodyIntersectingWithRay(const raytracing::Ray& ray,
real_t t_local;
Vec3 n_local;
+ //int i = 0;
for(auto grid: gridList_) {
- BodyID body = grid->getRayIntersectingBody<BodyTuple>(ray, blockAABB, t_closest, n_closest);
- if (body != NULL) {
+ //i++; //Todo
+ //if (i != 2) continue;
+
+ BodyID body = grid->getRayIntersectingBody<BodyTuple>(ray, blockAABB, t_local, n_local);
+ if (body != NULL && t_local < t_closest) {
body_closest = body;
+ t_closest = t_local;
+ n_closest = n_local;
}
}
@@ -777,9 +786,9 @@ BodyID HashGrids::getClosestBodyIntersectingWithRay(const raytracing::Ray& ray,
}
}
- n = n_closest;
t = t_closest;
-
+ n = n_closest;
+
return body_closest;
}
diff --git a/src/pe/raytracing/Intersects.h b/src/pe/raytracing/Intersects.h
index 217b65cb23120d863a8f40bd666623f8f83efa50..5b1c754489628544ed6df499d4ae87becc3aea3b 100644
--- a/src/pe/raytracing/Intersects.h
+++ b/src/pe/raytracing/Intersects.h
@@ -120,7 +120,7 @@ inline bool intersects(const BoxID box, const Ray& ray, real_t& t, Vec3& n) {
const Vec3& invDirection = transformedRay.getInvDirection();
const Vec3& origin = transformedRay.getOrigin();
- real_t inf = std::numeric_limits<real_t>::max();
+ const real_t inf = std::numeric_limits<real_t>::max();
size_t tminAxis = 0, tmaxAxis = 0;
real_t txmin, txmax;
@@ -181,9 +181,9 @@ inline bool intersects(const BoxID box, const Ray& ray, real_t& t, Vec3& n) {
}
inline bool intersects(const CapsuleID capsule, const Ray& ray, real_t& t, Vec3& n) {
- Ray transformedRay = ray.transformedToBF(capsule);
- Vec3 direction = transformedRay.getDirection();
- Vec3 origin = transformedRay.getOrigin();
+ const Ray transformedRay = ray.transformedToBF(capsule);
+ const Vec3& direction = transformedRay.getDirection();
+ const Vec3& origin = transformedRay.getOrigin();
real_t halfLength = capsule->getLength()/real_t(2);
real_t inf = std::numeric_limits<real_t>::max();
diff --git a/src/pe/raytracing/Ray.h b/src/pe/raytracing/Ray.h
index 3e248444459cd78fd6c96f266cea7098d2acdabb..343a6d0327ce42d55c26dacac12b3a20aa70b491 100644
--- a/src/pe/raytracing/Ray.h
+++ b/src/pe/raytracing/Ray.h
@@ -73,11 +73,61 @@ public:
return direction_;
}
+ /*!\brief Returns the normalized direction vector of the ray for a given axis.
+ */
+ inline real_t getDirection (size_t axis) const {
+ WALBERLA_ASSERT(axis >= 0 && axis <= 2, "No valid axis index passed.");
+ return direction_[axis];
+ }
+
+ /*!\brief Returns the x component of the ray direction.
+ */
+ inline real_t xDir () const {
+ return direction_[0];
+ }
+
+ /*!\brief Returns the y component of the ray direction.
+ */
+ inline real_t yDir () const {
+ return direction_[1];
+ }
+
+ /*!\brief Returns the z component of the ray direction.
+ */
+ inline real_t zDir () const {
+ return direction_[2];
+ }
+
/*!\brief Returns the inverse of the direction vector of the ray.
*/
inline const Vec3& getInvDirection () const {
return inv_direction_;
}
+
+ /*!\brief Returns the inverse of the direction vector of the ray for a given axis.
+ */
+ inline const real_t getInvDirection (size_t axis) const {
+ WALBERLA_ASSERT(axis >= 0 && axis <= 2, "No valid axis index passed.");
+ return inv_direction_[axis];
+ }
+
+ /*!\brief Returns the x component of the inverse ray direction.
+ */
+ inline real_t xInvDir () const {
+ return inv_direction_[0];
+ }
+
+ /*!\brief Returns the y component of the inverse ray direction.
+ */
+ inline real_t yInvDir () const {
+ return inv_direction_[1];
+ }
+
+ /*!\brief Returns the z component of the inverse ray direction.
+ */
+ inline real_t zInvDir () const {
+ return inv_direction_[2];
+ }
/*!\brief Returns the signs of the inverted direction vector of the ray.
*
diff --git a/src/pe/raytracing/Raytracer.h b/src/pe/raytracing/Raytracer.h
index 8a58415479ebdca57c7603272588ae2e1c9173ac..5a62e69e04eaf1a26c70d4d1b644409cb9293e4d 100644
--- a/src/pe/raytracing/Raytracer.h
+++ b/src/pe/raytracing/Raytracer.h
@@ -445,13 +445,11 @@ void Raytracer::rayTrace(const size_t timestep, WcTimingTree* tt) {
BodyID body = hashgrids->getClosestBodyIntersectingWithRay<BodyTypeTuple>(ray, blockAabb, t, n);
#if defined(COMPARE_NAIVE_AND_HASHGRIDS_RAYTRACING)
- if (body != NULL && t < t_hashgrids_closest) {
-#ifdef ignore_this_just_for_correct_indentation
- }
-#endif
+ if (body != NULL && t < t_hashgrids_closest)
#else
- if (body != NULL && t < t_closest) {
+ if (body != NULL && t < t_closest)
#endif
+ {
#if defined(COMPARE_NAIVE_AND_HASHGRIDS_RAYTRACING)
t_hashgrids_closest = t;
body_hashgrids_closest = body;
@@ -602,7 +600,8 @@ void Raytracer::rayTrace(const size_t timestep, WcTimingTree* tt) {
WALBERLA_LOG_INFO("Performed " << naiveIntersectionTests << " naive intersection tests");
#endif
#if defined(COMPARE_NAIVE_AND_HASHGRIDS_RAYTRACING)
- WALBERLA_LOG_INFO("Performed " << HashGrids::inter << " naive intersection tests");
+ WALBERLA_LOG_INFO("Performed " << ccd::HashGrids::intersectionTestCount << " intersection tests in hashgrids");
+ WALBERLA_LOG_INFO("Saved " << (int(naiveIntersectionTests)-int(ccd::HashGrids::intersectionTestCount)) << " tests (" << ((real_t(1)-real_c(ccd::HashGrids::intersectionTestCount)/real_c(naiveIntersectionTests))*100) << "%).")
#endif
if (tt != NULL) tt->start("Reduction");
diff --git a/tests/pe/Raytracing.cpp b/tests/pe/Raytracing.cpp
index 1d1edece6fff6943fc81116fc828bc586562a6f6..2b7f750f938299ea511f61d8a98fa18775262126 100644
--- a/tests/pe/Raytracing.cpp
+++ b/tests/pe/Raytracing.cpp
@@ -267,21 +267,58 @@ void RaytracerTest() {
createPlane(*globalBodyStorage, 0, Vec3(-1,1,1), Vec3(8,2,2), iron); // tilted plane in right bottom back corner
- createSphere(*globalBodyStorage, *forest, storageID, 2, Vec3(6,9.5,9.5), real_t(0.5));
+ //createSphere(*globalBodyStorage, *forest, storageID, 2, Vec3(6,9.5,9.5), real_t(0.5));
createSphere(*globalBodyStorage, *forest, storageID, 3, Vec3(4,5.5,5), real_t(1));
- createSphere(*globalBodyStorage, *forest, storageID, 6, Vec3(3,8.5,5), real_t(1));
+ //createSphere(*globalBodyStorage, *forest, storageID, 6, Vec3(3,8.5,5), real_t(1));
BoxID box = createBox(*globalBodyStorage, *forest, storageID, 7, Vec3(5,6.5,5), Vec3(2,4,3));
if (box != NULL) box->rotate(0,math::M_PI/4,math::M_PI/4);
- createBox(*globalBodyStorage, *forest, storageID, 8, Vec3(5,1,8), Vec3(2,2,2));
+ //createBox(*globalBodyStorage, *forest, storageID, 8, Vec3(5,1,8), Vec3(2,2,2));
// Test scene v1 end
// Test scene v2 additions start
- createBox(*globalBodyStorage, *forest, storageID, 9, Vec3(9,9,5), Vec3(1,1,10));
+ //createBox(*globalBodyStorage, *forest, storageID, 9, Vec3(9,9,5), Vec3(1,1,10));
createCapsule(*globalBodyStorage, *forest, storageID, 10, Vec3(3, 9, 1), real_t(0.5), real_t(7), iron);
CapsuleID capsule = createCapsule(*globalBodyStorage, *forest, storageID, 11, Vec3(7, 3.5, 7.5), real_t(1), real_t(2), iron);
if (capsule != NULL) capsule->rotate(0,math::M_PI/3,math::M_PI/4-math::M_PI/8);
// Test scene v2 end
+ for (auto blockIt = forest->begin(); blockIt != forest->end(); ++blockIt) {
+ ccd::HashGrids* hashgrids = blockIt->getData<ccd::HashGrids>(ccdID);
+ hashgrids->update();
+ }
+
+ WALBERLA_LOG_INFO("--- CAPSULE");
+ WALBERLA_LOG_INFO(" hash: " << capsule->getHash());
+ WALBERLA_LOG_INFO(" grid: " << capsule->getGrid());
+ WALBERLA_LOG_INFO(" AABB: " << capsule->getAABB());
+ WALBERLA_LOG_INFO("--- ROTATED BOX");
+ WALBERLA_LOG_INFO(" hash: " << box->getHash());
+ WALBERLA_LOG_INFO(" grid: " << box->getGrid());
+ WALBERLA_LOG_INFO(" AABB: " << box->getAABB());
+
+ Ray ray(Vec3(-5, 6.5, 5), Vec3(1, 0, 0));
+ real_t t = INFINITY;
+ Vec3 n;
+ BodyID body = NULL;
+
+ // output info about grids
+ for (auto blockIt = forest->begin(); blockIt != forest->end(); ++blockIt) {
+ ccd::HashGrids* hashgrids = blockIt->getData<ccd::HashGrids>(ccdID);
+ for (auto grid: hashgrids->gridList_) {
+ WALBERLA_LOG_INFO("--- GRID " << grid << " ---");
+ WALBERLA_LOG_INFO(" cellSpan: " << grid->getCellSpan());
+ WALBERLA_LOG_INFO(" dims: " << grid->xCellCount_ << "/" << grid->yCellCount_ << "/" << grid->zCellCount_);
+ WALBERLA_LOG_INFO(" items: " << grid->bodyCount_);
+ WALBERLA_LOG_INFO(" enlargement threshold: " << grid->enlargementThreshold_);
+ }
+ body = hashgrids->getClosestBodyIntersectingWithRay<BodyTuple>(ray, blockIt->getAABB(), t, n);
+ if (body != NULL) {
+ WALBERLA_LOG_INFO("body found");
+ } else {
+ WALBERLA_LOG_INFO("body not found");
+ }
+ }
+
//raytracer.setTBufferOutputDirectory("tbuffer");
//raytracer.setTBufferOutputEnabled(true);
raytracer.setImageOutputDirectory("image");
@@ -379,7 +416,7 @@ void RaytracerSpheresTest() {
raytracer.setImageOutputDirectory("image");
raytracer.setImageOutputEnabled(true);
- raytracer.rayTrace<BodyTuple>(0);
+ raytracer.rayTrace<BodyTuple>(1);
}
void hashgridsPlayground() {
@@ -502,7 +539,7 @@ void HashGridsTest(size_t boxes, size_t capsules, size_t spheres) {
tt.start("Setup");
shared_ptr<BodyStorage> globalBodyStorage = make_shared<BodyStorage>();
- shared_ptr<BlockForest> forest = createBlockForest(AABB(0,0,0,4,4,4), Vec3(2,2,1), Vec3(false, false, false));
+ shared_ptr<BlockForest> forest = createBlockForest(AABB(0,0,0,4,4,4), Vec3(1,1,1), Vec3(false, false, false));
auto storageID = forest->addBlockData(createStorageDataHandling<BodyTuple>(), "Storage");
auto ccdID = forest->addBlockData(ccd::createHashGridsDataHandling(globalBodyStorage, storageID), "CCD");
@@ -526,10 +563,10 @@ void HashGridsTest(size_t boxes, size_t capsules, size_t spheres) {
// generate bodies for test
std::vector<BodyID> bodies;
for (size_t i = 0; i < boxes; i++) {
- real_t len = math::realRandom(real_t(0.2), real_t(0.5));
- real_t x_min = math::realRandom(forestAABB.xMin()+len, forestAABB.xMax()-len);
- real_t y_min = math::realRandom(forestAABB.yMin()+len, forestAABB.yMax()-len);
- real_t z_min = math::realRandom(forestAABB.zMin()+len, forestAABB.zMax()-len);
+ real_t len = math::realRandom(real_t(0.1), real_t(0.2)); //0.2 0.5
+ real_t x_min = math::realRandom(forestAABB.xMin()+len, forestAABB.xMax());
+ real_t y_min = math::realRandom(forestAABB.yMin()+len, forestAABB.yMax());
+ real_t z_min = math::realRandom(forestAABB.zMin()+len, forestAABB.zMax());
//real_t z_min = len+0.1;
walberla::id_t id = walberla::id_t(i);
BoxID box_ = createBox(*globalBodyStorage, *forest, storageID, id, Vec3(x_min, y_min, z_min), Vec3(len, len, len));
@@ -540,12 +577,12 @@ void HashGridsTest(size_t boxes, size_t capsules, size_t spheres) {
}
for (size_t i = 0; i < capsules; i++) {
- real_t len = math::realRandom(real_t(0.2), real_t(0.5));
+ real_t len = math::realRandom(real_t(0.1), real_t(0.3)); // 0.2 0.5
real_t radius = real_t(0.1);
real_t maxlen = len + 2*radius;
- real_t x = math::realRandom(forestAABB.xMin()+maxlen, forestAABB.xMax()-maxlen);
- real_t y = math::realRandom(forestAABB.yMin()+maxlen, forestAABB.yMax()-maxlen);
- real_t z = math::realRandom(forestAABB.zMin()+maxlen, forestAABB.zMax()-maxlen);
+ real_t x = math::realRandom(forestAABB.xMin()+maxlen, forestAABB.xMax());
+ real_t y = math::realRandom(forestAABB.yMin()+maxlen, forestAABB.yMax());
+ real_t z = math::realRandom(forestAABB.zMin()+maxlen, forestAABB.zMax());
walberla::id_t id = walberla::id_t(boxes+i);
CapsuleID capsule = createCapsule(*globalBodyStorage, *forest, storageID, id, Vec3(x, y, z), radius, len);
WALBERLA_CHECK(capsule != NULL);
@@ -555,11 +592,11 @@ void HashGridsTest(size_t boxes, size_t capsules, size_t spheres) {
}
for (size_t i = 0; i < spheres; i++) {
- real_t radius = math::realRandom(real_t(0.2), real_t(0.3));
+ real_t radius = math::realRandom(real_t(0.1), real_t(0.2)); // 0.2 0.3
// forestAABB.xMax()-radius gerechtfertigt?
- real_t x = math::realRandom(forestAABB.xMin()+radius, forestAABB.xMax()-radius);
- real_t y = math::realRandom(forestAABB.yMin()+radius, forestAABB.yMax()-radius);
- real_t z = math::realRandom(forestAABB.zMin()+radius, forestAABB.zMax()-radius);
+ real_t x = math::realRandom(forestAABB.xMin()+radius, forestAABB.xMax());
+ real_t y = math::realRandom(forestAABB.yMin()+radius, forestAABB.yMax());
+ real_t z = math::realRandom(forestAABB.zMin()+radius, forestAABB.zMax());
walberla::id_t id = walberla::id_t(boxes+capsules+i);
SphereID sphere = createSphere(*globalBodyStorage, *forest, storageID, id, Vec3(x, y, z), radius);
WALBERLA_CHECK(sphere != NULL);
@@ -590,7 +627,7 @@ void HashGridsTest(size_t boxes, size_t capsules, size_t spheres) {
- Lighting lighting(Vec3((forestAABB.xMin()+forestAABB.xMax())/real_t(2)+1, (forestAABB.yMin()+forestAABB.yMax())/real_t(2),
+ Lighting lighting(Vec3(forestAABB.xSize()/real_t(2)+1, forestAABB.ySize()/real_t(2),
real_t(2)*forestAABB.zMax()+2), // 8, 5, 9.5 gut für ebenen, 0,5,8
Color(1, 1, 1), //diffuse
Color(1, 1, 1), //specular
@@ -598,21 +635,29 @@ void HashGridsTest(size_t boxes, size_t capsules, size_t spheres) {
Raytracer raytracer(forest, storageID, globalBodyStorage, ccdID,
size_t(640), size_t(480),
49.13,
- Vec3((forestAABB.xMin()+forestAABB.xMax())/real_t(2), (forestAABB.yMin()+forestAABB.yMax())/real_t(2),
+ Vec3(forestAABB.xSize()/real_t(2), forestAABB.ySize()/real_t(2),
real_t(2)*forestAABB.zMax()),
- Vec3((forestAABB.xMin()+forestAABB.xMax())/real_t(2), (forestAABB.yMin()+forestAABB.yMax())/real_t(2),
+ Vec3(forestAABB.xSize()/real_t(2), forestAABB.ySize()/real_t(2),
0),
Vec3(0,1,0), //-5,5,5; -1,5,5
lighting,
Color(0.2,0.2,0.2),
real_t(2),
customHashgridsBodyToShadingParams);
- raytracer.setImageOutputDirectory("image");
+#if defined(USE_NAIVE_INTERSECTION_FINDING)
+ raytracer.setImageOutputDirectory("image/naive");
+#endif
+#if !defined(USE_NAIVE_INTERSECTION_FINDING) && !defined(COMPARE_NAIVE_AND_HASHGRIDS_RAYTRACING)
+ raytracer.setImageOutputDirectory("image/hashgrids");
+#endif
+#if defined(COMPARE_NAIVE_AND_HASHGRIDS_RAYTRACING)
+ raytracer.setImageOutputDirectory("image/comparison");
+#endif
raytracer.setImageOutputEnabled(true);
- raytracer.setFilenameTimestepWidth(9);
+ raytracer.setFilenameTimestepWidth(12);
tt.stop("Setup");
- WALBERLA_LOG_INFO("output to: " << (boxes*1000000 + capsules*1000 + spheres));
- raytracer.rayTrace<BodyTuple>(boxes*1000000 + capsules*1000 + spheres, &tt);
+ WALBERLA_LOG_INFO("output to: " << (boxes*100000000 + capsules*10000 + spheres));
+ raytracer.rayTrace<BodyTuple>(boxes*100000000 + capsules*10000 + spheres, &tt);
/*HashGrids::HashGrid* grid;
for (auto bodyId: problematicBodyIDs) {
@@ -632,26 +677,254 @@ void HashGridsTest(size_t boxes, size_t capsules, size_t spheres) {
WALBERLA_ROOT_SECTION() {
std::cout << temp;
}
+}
+
+Vec3 minCornerToGpos(const Vec3& minCorner, real_t lengths) {
+ return minCorner + Vec3(lengths/2, lengths/2, lengths/2);
+}
+
+void HashGridsTestScene() {
+#if defined(USE_NAIVE_INTERSECTION_FINDING)
+ WALBERLA_LOG_INFO("Using naive method for intersection testing");
+#endif
+#if !defined(USE_NAIVE_INTERSECTION_FINDING) && !defined(COMPARE_NAIVE_AND_HASHGRIDS_RAYTRACING)
+ WALBERLA_LOG_INFO("Using hashgrids for intersection testing");
+#endif
+#if defined(COMPARE_NAIVE_AND_HASHGRIDS_RAYTRACING)
+ WALBERLA_LOG_INFO("Comparing hashgrids and naive method for intersection testing");
+#endif
+ using namespace walberla::pe::ccd;
+ WcTimingTree tt;
+ tt.start("Setup");
+
+ shared_ptr<BodyStorage> globalBodyStorage = make_shared<BodyStorage>();
+ shared_ptr<BlockForest> forest = createBlockForest(AABB(0,0,0,8,8,8), Vec3(1,1,1), Vec3(false, false, false));
+ auto storageID = forest->addBlockData(createStorageDataHandling<BodyTuple>(), "Storage");
+ auto ccdID = forest->addBlockData(ccd::createHashGridsDataHandling(globalBodyStorage, storageID), "CCD");
+
+ const AABB& forestAABB = forest->getDomain();
+
+ std::vector<BodyID> bodies;
+
+ // create bodies
+ size_t id = 0;
+ real_t len = 2;
+ //bodies.push_back(createBox(*globalBodyStorage, *forest, storageID, ++id, minCornerToGpos(Vec3(1.1, 1.1, 1.1), len), Vec3(len, len, len)));
+ len = 0.6;
+ //bodies.push_back(createBox(*globalBodyStorage, *forest, storageID, ++id, minCornerToGpos(Vec3(0.6, 0.5, 0.5), len), Vec3(len, len, len)));
+ //bodies.push_back(createBox(*globalBodyStorage, *forest, storageID, ++id, minCornerToGpos(Vec3(1.7, 0.5, 0.5), len), Vec3(len, len, len)));
+ //bodies.push_back(createBox(*globalBodyStorage, *forest, storageID, ++id, minCornerToGpos(Vec3(7.2, 0.5, 0.5), len), Vec3(len, len, len)));
+ //bodies.push_back(createBox(*globalBodyStorage, *forest, storageID, ++id, minCornerToGpos(Vec3(5, 0.5, 0.5), len), Vec3(len, len, len)));
+
+ real_t x_min = 0, y_min = 0;
+ len = 1.2;
+ real_t gap = 0.4;
+ for (int i = 0; ; ++i) {
+ x_min = forestAABB.xMin() + i*(gap+len);
+ if (x_min > forestAABB.max(0)) {
+ break;
+ }
+ for (int j = 0; ; ++j) {
+ y_min = forestAABB.yMin() + j*(gap+len);
+ if (y_min > forestAABB.max(1)) {
+ break;
+ }
+
+ bodies.push_back(createBox(*globalBodyStorage, *forest, storageID, ++id, minCornerToGpos(Vec3(x_min, y_min, 0), len), Vec3(len, len, len)));
+ }
+ }
+
+ // update hashgrids
+ for (auto blockIt = forest->begin(); blockIt != forest->end(); ++blockIt) {
+ ccd::HashGrids* hashgrids = blockIt->getData<ccd::HashGrids>(ccdID);
+ hashgrids->update();
+ }
+
+ // output info about bodies
+ /*for (auto body: bodies) {
+ WALBERLA_LOG_INFO("--- BODY " << body->getID() << " ---");
+ WALBERLA_LOG_INFO(" corners: " << body->getAABB().minCorner() << ", " << body->getAABB().maxCorner());
+ WALBERLA_LOG_INFO(" hash: " << body->getHash());
+ //WALBERLA_LOG_INFO(" body index: " << body->getCellId());
+ WALBERLA_LOG_INFO(" grid: " << body->getGrid());
+ }*/
+
+ // output info about grids
+ for (auto blockIt = forest->begin(); blockIt != forest->end(); ++blockIt) {
+ ccd::HashGrids* hashgrids = blockIt->getData<ccd::HashGrids>(ccdID);
+ for (auto grid: hashgrids->gridList_) {
+ WALBERLA_LOG_INFO("--- GRID " << grid << " ---");
+ WALBERLA_LOG_INFO(" cellSpan: " << grid->getCellSpan());
+ WALBERLA_LOG_INFO(" dims: " << grid->xCellCount_ << "/" << grid->yCellCount_ << "/" << grid->zCellCount_);
+ WALBERLA_LOG_INFO(" items: " << grid->bodyCount_);
+ WALBERLA_LOG_INFO(" enlargement threshold: " << grid->enlargementThreshold_);
+ }
+ }
+
+
+ /*Ray ray(Vec3(4,4,10), Vec3(0,0,-1));
+ real_t t;
+ Vec3 n;
+
+ for (auto blockIt = forest->begin(); blockIt != forest->end(); ++blockIt) {
+ ccd::HashGrids* hashgrids = blockIt->getData<ccd::HashGrids>(ccdID);
+ BodyID body = hashgrids->getClosestBodyIntersectingWithRay<BodyTuple>(ray, blockIt->getAABB(), t, n);
+ if (body != NULL) {
+ WALBERLA_LOG_INFO("found body " << body->getID() << ", t = " << t);
+ WALBERLA_LOG_INFO(" hash: " << body->getHash());
+ WALBERLA_LOG_INFO(" grid: " << body->getGrid());
+ } else {
+ WALBERLA_LOG_INFO("did not find body");
+ }
+ }*/
+
+
+ // setup raytracer
+ Lighting lighting(Vec3(1,2,15),
+ Color(1, 1, 1), //diffuse
+ Color(1, 1, 1), //specular
+ Color(0.4, 0.4, 0.4)); //ambient
+ Raytracer raytracer(forest, storageID, globalBodyStorage, ccdID,
+ size_t(480), size_t(480),
+ 49.13,
+ Vec3(4,4,18), //6,6,8
+ Vec3(4,4,10), //6,6,4
+ Vec3(0,1,0),
+ lighting,
+ Color(0.2,0.2,0.2),
+ real_t(2));
+
+#if defined(USE_NAIVE_INTERSECTION_FINDING)
+ raytracer.setImageOutputDirectory("image/naive");
+#endif
+#if !defined(USE_NAIVE_INTERSECTION_FINDING) && !defined(COMPARE_NAIVE_AND_HASHGRIDS_RAYTRACING)
+ raytracer.setImageOutputDirectory("image/hashgrids");
+#endif
+#if defined(COMPARE_NAIVE_AND_HASHGRIDS_RAYTRACING)
+ raytracer.setImageOutputDirectory("image/comparison");
+#endif
+ raytracer.setImageOutputEnabled(true);
+ raytracer.setFilenameTimestepWidth(12);
+ tt.stop("Setup");
+ WALBERLA_LOG_INFO("output to: 1");
+ raytracer.rayTrace<BodyTuple>(1, &tt);
+}
+
+void GridCellIntersectionPlayground() {
+ real_t cellSpan = 2;
+
+ AABB blockAABB(Vec3(-2.5,-2.5,0), Vec3(3.5,5.5,12));
+
+ size_t blockXCellCount = uint_c(std::ceil(blockAABB.xSize() / cellSpan));
+ size_t blockYCellCount = uint_c(std::ceil(blockAABB.ySize() / cellSpan));
+ size_t blockZCellCount = uint_c(std::ceil(blockAABB.zSize() / cellSpan));
+
+ //Ray ray(Vec3(2,-0.5,0), Vec3(2,13,0).getNormalized());
+ //Ray ray(Vec3(1.5,6,0), Vec3(-4,-9,0).getNormalized());
+ //Ray ray(Vec3(2,-2,0), Vec3(5,8,0).getNormalized());
+ //Ray ray(Vec3(2,-2,0), Vec3(-5,8,0).getNormalized());
+ Ray ray(Vec3(0.5,-3,0), Vec3(-4,7,0).getNormalized());
+
+ const real_t inf = std::numeric_limits<real_t>::max();
+
+ real_t invCellSpan = real_t(1)/cellSpan;
+
+ std::vector<Vec3> intersectionPoints;
+ std::vector<real_t> distances;
+ std::vector<size_t> axii;
+
+ Vec3 firstPoint;
+ if (blockAABB.contains(ray.getOrigin())) {
+ firstPoint = ray.getOrigin();
+ } else {
+ real_t t_start;
+ if (intersects(blockAABB, ray, t_start)) {
+ firstPoint = ray.getOrigin() + ray.getDirection()*t_start;
+ } else {
+ WALBERLA_LOG_INFO("Ray does not intersect block.");
+ return;
+ }
+ }
+
+ Vector3<int32_t> firstCell(int32_c(std::floor(std::max(firstPoint[0]-blockAABB.xMin()*invCellSpan, real_t(0)))),
+ int32_c(std::floor(std::max(firstPoint[1]-blockAABB.yMin()*invCellSpan, real_t(0)))),
+ int32_c(std::floor(std::max(firstPoint[2]-blockAABB.zMin()*invCellSpan, real_t(0)))));
+
+ const int8_t stepX = ray.xDir() >= 0 ? 1 : -1;
+ const int8_t stepY = ray.yDir() >= 0 ? 1 : -1;
+ const int8_t stepZ = ray.zDir() >= 0 ? 1 : -1;
+
+ Vec3 near((stepX >= 0) ? (firstCell[0]+1)*cellSpan-firstPoint[0] : firstPoint[0]-firstCell[0]*cellSpan,
+ (stepY >= 0) ? (firstCell[1]+1)*cellSpan-firstPoint[1] : firstPoint[1]-firstCell[1]*cellSpan,
+ (stepZ >= 0) ? (firstCell[2]+1)*cellSpan-firstPoint[2] : firstPoint[2]-firstCell[2]*cellSpan);
+
+ real_t tMaxX = (ray.xDir() != 0) ? std::abs(near[0]*ray.xInvDir()) : inf;
+ real_t tMaxY = (ray.yDir() != 0) ? std::abs(near[1]*ray.yInvDir()) : inf;
+ real_t tMaxZ = (ray.zDir() != 0) ? std::abs(near[2]*ray.zInvDir()) : inf;
+
+ real_t tDeltaX = (ray.xDir() != 0) ? std::abs(cellSpan*ray.xInvDir()) : inf;
+ real_t tDeltaY = (ray.yDir() != 0) ? std::abs(cellSpan*ray.yInvDir()) : inf;
+ real_t tDeltaZ = (ray.zDir() != 0) ? std::abs(cellSpan*ray.zInvDir()) : inf;
- WALBERLA_LOG_INFO("Performed " << HashGrids::intersectionTestCount << " intersection tests in hashgrids");
+ Vector3<int32_t> currentCell = firstCell;
+
+ // First cell needs extra treatment, as it might lay out of the blocks upper bounds
+ // due to the nature of how it is calculated: If the first point lies on a upper border
+ // it maps to the cell "above" the grid.
+ if (currentCell[0] < blockXCellCount && currentCell[1] < blockYCellCount && currentCell[2] < blockZCellCount) {
+ WALBERLA_LOG_INFO("found valid cell at " << currentCell);
+ }
+
+ while (true) {
+ if (tMaxX < tMaxY) {
+ if (tMaxX < tMaxZ) {
+ tMaxX += tDeltaX;
+ currentCell[0] += stepX;
+ if (currentCell[0] == blockXCellCount || currentCell[0] == -1) {
+ break;
+ }
+ } else {
+ tMaxZ += tDeltaZ;
+ currentCell[2] += stepZ;
+ if (currentCell[2] == blockZCellCount || currentCell[2] == -1) {
+ break;
+ }
+ }
+ } else {
+ if (tMaxY < tMaxZ) {
+ tMaxY += tDeltaY;
+ currentCell[1] += stepY;
+ if (currentCell[1] == blockYCellCount || currentCell[1] == -1) {
+ break;
+ }
+ } else {
+ tMaxZ += tDeltaZ;
+ currentCell[2] += stepZ;
+ if (currentCell[2] == blockZCellCount || currentCell[2] == -1) {
+ break;
+ }
+ }
+ }
+
+ WALBERLA_LOG_INFO("found cell at " << currentCell);
+ }
}
int main( int argc, char** argv )
{
walberla::debug::enterTestMode();
walberla::MPIManager::instance()->initializeMPI( &argc, &argv );
-
SetBodyTypeIDs<BodyTuple>::execute();
+ math::seedRandomGenerator( static_cast<unsigned int>(1337 * mpi::MPIManager::instance()->worldRank()) );
//SphereIntersectsTest();
//PlaneIntersectsTest();
//BoxIntersectsTest();
//AABBIntersectsTest();
//CapsuleIntersectsTest();
- //RaytracerTest();
+ RaytracerTest();
//RaytracerSpheresTest();
- math::seedRandomGenerator( static_cast<unsigned int>(1337 * mpi::MPIManager::instance()->worldRank()) );
std::vector<size_t> boxes = {127, 70, 20, 150};
std::vector<size_t> capsules = {127, 60, 140, 100};
@@ -660,8 +933,11 @@ int main( int argc, char** argv )
//for (size_t i = 0; i < boxes.size(); ++i) {
// HashGridsTest(boxes[i], capsules[i], spheres[i]);
//}
- HashGridsTest(boxes[2], capsules[2], spheres[2]);
-
+ //HashGridsTest(boxes[1], boxes[1], boxes[1]);
+
+ //HashGridsTestScene();
+
+ //GridCellIntersectionPlayground();
//HashGridsTest();