diff --git a/apps/benchmarks/FluidParticleCoupling/ForcesOnSphereNearPlane.cpp b/apps/benchmarks/FluidParticleCoupling/ForcesOnSphereNearPlane.cpp
index ee3c0f1431a0c73d137c1bd07943729bab44ba8b..1ef29cf39aa0157493dc4fd63a5a3a39c192c17e 100644
--- a/apps/benchmarks/FluidParticleCoupling/ForcesOnSphereNearPlane.cpp
+++ b/apps/benchmarks/FluidParticleCoupling/ForcesOnSphereNearPlane.cpp
@@ -309,6 +309,7 @@ void createPlaneSetup(const shared_ptr<mesa_pd::data::ParticleStorage> & ps, con
// create bounding planes
mesa_pd::data::Particle p0 = *ps->create(true);
p0.setPosition(simulationDomain.minCorner());
+ p0.setInteractionRadius(std::numeric_limits<real_t>::infinity());
p0.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(0,0,1) ));
p0.setOwner(mpi::MPIManager::instance()->rank());
p0.setType(0);
@@ -317,6 +318,7 @@ void createPlaneSetup(const shared_ptr<mesa_pd::data::ParticleStorage> & ps, con
mesa_pd::data::Particle p1 = *ps->create(true);
p1.setPosition(simulationDomain.maxCorner());
+ p1.setInteractionRadius(std::numeric_limits<real_t>::infinity());
p1.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(0,0,-1) ));
p1.setOwner(mpi::MPIManager::instance()->rank());
p1.setType(0);
diff --git a/apps/benchmarks/FluidParticleCoupling/LubricationForceEvaluation.cpp b/apps/benchmarks/FluidParticleCoupling/LubricationForceEvaluation.cpp
index e17092473fdbc3524286da109ec2b93552c71dd5..4d3b151b46ae8ef2556495dbf2a0d343b94ad570 100644
--- a/apps/benchmarks/FluidParticleCoupling/LubricationForceEvaluation.cpp
+++ b/apps/benchmarks/FluidParticleCoupling/LubricationForceEvaluation.cpp
@@ -443,6 +443,7 @@ int main( int argc, char **argv )
// create two planes
mesa_pd::data::Particle&& p0 = *ps->create(true);
p0.setPosition(referenceVector);
+ p0.setInteractionRadius(std::numeric_limits<real_t>::infinity());
p0.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(1,0,0) ));
p0.setOwner(mpi::MPIManager::instance()->rank());
mesa_pd::data::particle_flags::set(p0.getFlagsRef(), mesa_pd::data::particle_flags::INFINITE);
@@ -451,6 +452,7 @@ int main( int argc, char **argv )
mesa_pd::data::Particle&& p1 = *ps->create(true);
p1.setPosition(Vector3<real_t>(real_c(xSize),0,0));
+ p1.setInteractionRadius(std::numeric_limits<real_t>::infinity());
p1.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(-1,0,0) ));
p1.setOwner(mpi::MPIManager::instance()->rank());
mesa_pd::data::particle_flags::set(p1.getFlagsRef(), mesa_pd::data::particle_flags::INFINITE);
diff --git a/apps/benchmarks/FluidParticleCoupling/ObliqueDryCollision.cpp b/apps/benchmarks/FluidParticleCoupling/ObliqueDryCollision.cpp
index 9314699f4f67314fc0002c643e26c52d6971443d..784a10334a07e33021d984f336ed2617f4f5d3ce 100644
--- a/apps/benchmarks/FluidParticleCoupling/ObliqueDryCollision.cpp
+++ b/apps/benchmarks/FluidParticleCoupling/ObliqueDryCollision.cpp
@@ -119,11 +119,13 @@ int main( int argc, char ** argv )
data::Particle&& p = *ps->create();
p.setPosition(Vec3(0,0,2*radius));
p.setLinearVelocity(Vec3(uTin, 0., -uNin));
+ p.setInteractionRadius(radius);
p.setType(0);
// create plane
data::Particle&& p0 = *ps->create(true);
p0.setPosition(Vec3(0,0,0));
+ p0.setInteractionRadius(std::numeric_limits<real_t>::infinity());
p0.setShapeID(ss->create<data::HalfSpace>(Vector3<real_t>(0,0,1)));
p0.setType(0);
data::particle_flags::set(p0.getFlagsRef(), data::particle_flags::INFINITE);
diff --git a/apps/benchmarks/FluidParticleCoupling/ObliqueWetCollision.cpp b/apps/benchmarks/FluidParticleCoupling/ObliqueWetCollision.cpp
index 3bb0c57eef6b1138797cef785afe47869abb3f0f..2415cf72dbce81af49ef574b82774623bdbb1513 100644
--- a/apps/benchmarks/FluidParticleCoupling/ObliqueWetCollision.cpp
+++ b/apps/benchmarks/FluidParticleCoupling/ObliqueWetCollision.cpp
@@ -370,6 +370,7 @@ void createPlaneSetup(const shared_ptr<mesa_pd::data::ParticleStorage> & ps, con
// create bounding planes
mesa_pd::data::Particle&& p0 = *ps->create(true);
p0.setPosition(simulationDomain.minCorner());
+ p0.setInteractionRadius(std::numeric_limits<real_t>::infinity());
p0.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(0,0,1) ));
p0.setOwner(mpi::MPIManager::instance()->rank());
p0.setType(0);
@@ -381,6 +382,7 @@ void createPlaneSetup(const shared_ptr<mesa_pd::data::ParticleStorage> & ps, con
//only create top plane when no outflow BC should be set there
mesa_pd::data::Particle&& p1 = *ps->create(true);
p1.setPosition(simulationDomain.maxCorner());
+ p1.setInteractionRadius(std::numeric_limits<real_t>::infinity());
p1.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(0,0,-1) ));
p1.setOwner(mpi::MPIManager::instance()->rank());
p1.setType(0);
diff --git a/apps/benchmarks/FluidParticleCoupling/SettlingSphereInBox.cpp b/apps/benchmarks/FluidParticleCoupling/SettlingSphereInBox.cpp
index 3f2f613b252ece2fc72ad5c3df257e3c2f2f5223..899c27f674d8cdc9bd59914bc7e8c9a99700e287 100644
--- a/apps/benchmarks/FluidParticleCoupling/SettlingSphereInBox.cpp
+++ b/apps/benchmarks/FluidParticleCoupling/SettlingSphereInBox.cpp
@@ -290,6 +290,7 @@ void createPlaneSetup(const shared_ptr<mesa_pd::data::ParticleStorage> & ps, con
// create bounding planes
mesa_pd::data::Particle p0 = *ps->create(true);
p0.setPosition(simulationDomain.minCorner());
+ p0.setInteractionRadius(std::numeric_limits<real_t>::infinity());
p0.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(0,0,1) ));
p0.setOwner(mpi::MPIManager::instance()->rank());
p0.setType(0);
@@ -298,6 +299,7 @@ void createPlaneSetup(const shared_ptr<mesa_pd::data::ParticleStorage> & ps, con
mesa_pd::data::Particle p1 = *ps->create(true);
p1.setPosition(simulationDomain.maxCorner());
+ p1.setInteractionRadius(std::numeric_limits<real_t>::infinity());
p1.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(0,0,-1) ));
p1.setOwner(mpi::MPIManager::instance()->rank());
p1.setType(0);
@@ -306,6 +308,7 @@ void createPlaneSetup(const shared_ptr<mesa_pd::data::ParticleStorage> & ps, con
mesa_pd::data::Particle p2 = *ps->create(true);
p2.setPosition(simulationDomain.minCorner());
+ p2.setInteractionRadius(std::numeric_limits<real_t>::infinity());
p2.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(1,0,0) ));
p2.setOwner(mpi::MPIManager::instance()->rank());
p2.setType(0);
@@ -314,6 +317,7 @@ void createPlaneSetup(const shared_ptr<mesa_pd::data::ParticleStorage> & ps, con
mesa_pd::data::Particle p3 = *ps->create(true);
p3.setPosition(simulationDomain.maxCorner());
+ p3.setInteractionRadius(std::numeric_limits<real_t>::infinity());
p3.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(-1,0,0) ));
p3.setOwner(mpi::MPIManager::instance()->rank());
p3.setType(0);
@@ -322,6 +326,7 @@ void createPlaneSetup(const shared_ptr<mesa_pd::data::ParticleStorage> & ps, con
mesa_pd::data::Particle p4 = *ps->create(true);
p4.setPosition(simulationDomain.minCorner());
+ p4.setInteractionRadius(std::numeric_limits<real_t>::infinity());
p4.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(0,1,0) ));
p4.setOwner(mpi::MPIManager::instance()->rank());
p4.setType(0);
@@ -330,6 +335,7 @@ void createPlaneSetup(const shared_ptr<mesa_pd::data::ParticleStorage> & ps, con
mesa_pd::data::Particle p5 = *ps->create(true);
p5.setPosition(simulationDomain.maxCorner());
+ p5.setInteractionRadius(std::numeric_limits<real_t>::infinity());
p5.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(0,-1,0) ));
p5.setOwner(mpi::MPIManager::instance()->rank());
p5.setType(0);
diff --git a/apps/benchmarks/FluidParticleCoupling/SphereMovingWithPrescribedVelocity.cpp b/apps/benchmarks/FluidParticleCoupling/SphereMovingWithPrescribedVelocity.cpp
index 9c963fa9dcbc0580af7ce7476d92c0884bbaec27..0b9e6ec68b214d220b93ba75c56be6c8e35af7ed 100644
--- a/apps/benchmarks/FluidParticleCoupling/SphereMovingWithPrescribedVelocity.cpp
+++ b/apps/benchmarks/FluidParticleCoupling/SphereMovingWithPrescribedVelocity.cpp
@@ -288,6 +288,7 @@ void createPlaneSetup(const shared_ptr<mesa_pd::data::ParticleStorage> & ps, con
{
mesa_pd::data::Particle p2 = *ps->create(true);
p2.setPosition(simulationDomain.minCorner());
+ p2.setInteractionRadius(std::numeric_limits<real_t>::infinity());
p2.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(1,0,0) ));
p2.setOwner(mpi::MPIManager::instance()->rank());
p2.setType(0);
@@ -297,6 +298,7 @@ void createPlaneSetup(const shared_ptr<mesa_pd::data::ParticleStorage> & ps, con
mesa_pd::data::Particle p3 = *ps->create(true);
p3.setPosition(simulationDomain.maxCorner());
+ p3.setInteractionRadius(std::numeric_limits<real_t>::infinity());
p3.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(-1,0,0) ));
p3.setOwner(mpi::MPIManager::instance()->rank());
p3.setType(0);
@@ -306,6 +308,7 @@ void createPlaneSetup(const shared_ptr<mesa_pd::data::ParticleStorage> & ps, con
mesa_pd::data::Particle p4 = *ps->create(true);
p4.setPosition(simulationDomain.minCorner());
+ p4.setInteractionRadius(std::numeric_limits<real_t>::infinity());
p4.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(0,1,0) ));
p4.setOwner(mpi::MPIManager::instance()->rank());
p4.setType(0);
@@ -315,6 +318,7 @@ void createPlaneSetup(const shared_ptr<mesa_pd::data::ParticleStorage> & ps, con
mesa_pd::data::Particle p5 = *ps->create(true);
p5.setPosition(simulationDomain.maxCorner());
+ p5.setInteractionRadius(std::numeric_limits<real_t>::infinity());
p5.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(0,-1,0) ));
p5.setOwner(mpi::MPIManager::instance()->rank());
p5.setType(0);
diff --git a/apps/benchmarks/FluidParticleCoupling/SphereWallCollision.cpp b/apps/benchmarks/FluidParticleCoupling/SphereWallCollision.cpp
index 1ce306d76a368929dce44e0e52fb73e7753bab7f..13df3b0a875f416360629711ac01fb435ebd4f9c 100644
--- a/apps/benchmarks/FluidParticleCoupling/SphereWallCollision.cpp
+++ b/apps/benchmarks/FluidParticleCoupling/SphereWallCollision.cpp
@@ -367,6 +367,7 @@ void createPlaneSetup(const shared_ptr<mesa_pd::data::ParticleStorage> & ps, con
// create bounding planes
mesa_pd::data::Particle&& p0 = *ps->create(true);
p0.setPosition(simulationDomain.minCorner());
+ p0.setInteractionRadius(std::numeric_limits<real_t>::infinity());
p0.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(0,0,1) ));
p0.setOwner(mpi::MPIManager::instance()->rank());
p0.setType(0);
@@ -378,6 +379,7 @@ void createPlaneSetup(const shared_ptr<mesa_pd::data::ParticleStorage> & ps, con
//only create top plane when no outflow BC should be set there
mesa_pd::data::Particle&& p1 = *ps->create(true);
p1.setPosition(simulationDomain.maxCorner());
+ p1.setInteractionRadius(std::numeric_limits<real_t>::infinity());
p1.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(0,0,-1) ));
p1.setOwner(mpi::MPIManager::instance()->rank());
p1.setType(0);
diff --git a/apps/benchmarks/FluidParticleCouplingWithLoadBalancing/FluidParticleWorkloadDistribution.cpp b/apps/benchmarks/FluidParticleCouplingWithLoadBalancing/FluidParticleWorkloadDistribution.cpp
index 9f93af67b37cfec23a3543948c70da55c388542b..78a8570bcc36678bf30654eb32ac829c3b6469fe 100644
--- a/apps/benchmarks/FluidParticleCouplingWithLoadBalancing/FluidParticleWorkloadDistribution.cpp
+++ b/apps/benchmarks/FluidParticleCouplingWithLoadBalancing/FluidParticleWorkloadDistribution.cpp
@@ -389,6 +389,7 @@ void createPlane( const shared_ptr<mesa_pd::data::ParticleStorage> & ps, const s
{
mesa_pd::data::Particle&& p0 = *ps->create(true);
p0.setPosition(position);
+ p0.setInteractionRadius(std::numeric_limits<real_t>::infinity());
p0.setShapeID(ss->create<mesa_pd::data::HalfSpace>( normal.getNormalized() ));
p0.setOwner(mpi::MPIManager::instance()->rank());
p0.setType(0);
diff --git a/apps/benchmarks/FluidParticleCouplingWithLoadBalancing/FluidParticleWorkloadEvaluation.cpp b/apps/benchmarks/FluidParticleCouplingWithLoadBalancing/FluidParticleWorkloadEvaluation.cpp
index 062e3f9c95f25baa5ddadd3b6511d9695cb5e457..a78736fff7996ee079d6bc2534a0e6bd2722347f 100644
--- a/apps/benchmarks/FluidParticleCouplingWithLoadBalancing/FluidParticleWorkloadEvaluation.cpp
+++ b/apps/benchmarks/FluidParticleCouplingWithLoadBalancing/FluidParticleWorkloadEvaluation.cpp
@@ -460,6 +460,7 @@ int main( int argc, char **argv )
mesa_pd::data::Particle&& p0 = *ps->create(true);
p0.setPosition(generationDomain.minCorner());
+ p0.setInteractionRadius(std::numeric_limits<real_t>::infinity());
p0.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(0,0,1) ));
p0.setOwner(mpi::MPIManager::instance()->rank());
p0.setType(0);
@@ -468,6 +469,7 @@ int main( int argc, char **argv )
mesa_pd::data::Particle&& p1 = *ps->create(true);
p1.setPosition(generationDomain.maxCorner());
+ p1.setInteractionRadius(std::numeric_limits<real_t>::infinity());
p1.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(0,0,-1) ));
p1.setOwner(mpi::MPIManager::instance()->rank());
p1.setType(0);
diff --git a/apps/benchmarks/GranularGas/CreateParticles.cpp b/apps/benchmarks/GranularGas/CreateParticles.cpp
index 4c788c163066a68ce70f0448f458fbbb7877bfc7..a7a76190db22e8555f4154fc7bcf9a21bd441314 100644
--- a/apps/benchmarks/GranularGas/CreateParticles.cpp
+++ b/apps/benchmarks/GranularGas/CreateParticles.cpp
@@ -30,6 +30,7 @@ data::ParticleStorage::iterator createPlane( data::ParticleStorage& ps,
{
auto p0 = ps.create(true);
p0->getPositionRef() = pos;
+ p0->getInteractionRadiusRef() = std::numeric_limits<real_t>::infinity();
p0->getShapeIDRef() = ss.create<data::HalfSpace>( normal );
p0->getOwnerRef() = walberla::mpi::MPIManager::instance()->rank();
p0->getTypeRef() = 0;
diff --git a/python/mesa_pd/templates/kernel/InsertParticleIntoLinkedCells.templ.h b/python/mesa_pd/templates/kernel/InsertParticleIntoLinkedCells.templ.h
index 96bcd932df6bdcc44b2ae034f104f6ab59c2f7a1..9cf3952cd27e4e7346928be9187d64060d653d8b 100644
--- a/python/mesa_pd/templates/kernel/InsertParticleIntoLinkedCells.templ.h
+++ b/python/mesa_pd/templates/kernel/InsertParticleIntoLinkedCells.templ.h
@@ -77,7 +77,9 @@ inline void InsertParticleIntoLinkedCells::operator()(const size_t p_idx, Access
ac.setNextParticle(p_idx, lc.infiniteParticles_.exchange(int_c(p_idx)));
} else
{
+ WALBERLA_ASSERT_GREATER(ac.getInteractionRadius(p_idx), 0_r, "Did you forget to set the interaction radius?");
{%- for dim in range(3) %}
+ WALBERLA_ASSERT_LESS(2_r * ac.getInteractionRadius(p_idx), lc.cellDiameter_[0], "Interaction radius is to large for this cell size. Contacts might get lost.");
int hash{{dim}} = static_cast<int>(std::floor((ac.getPosition(p_idx)[{{dim}}] - minCorner[{{dim}}]) * lc.invCellDiameter_[{{dim}}]));
{%- endfor %}
{%- for dim in range(3) %}
diff --git a/python/mesa_pd/templates/kernel/InsertParticleIntoSparseLinkedCells.templ.h b/python/mesa_pd/templates/kernel/InsertParticleIntoSparseLinkedCells.templ.h
index ff10b7b6d3de00bca08b0a99ae90d996885d0284..3edaf2bb019e578f035f2927c6fe90f6e8fca08a 100644
--- a/python/mesa_pd/templates/kernel/InsertParticleIntoSparseLinkedCells.templ.h
+++ b/python/mesa_pd/templates/kernel/InsertParticleIntoSparseLinkedCells.templ.h
@@ -77,7 +77,9 @@ inline void InsertParticleIntoSparseLinkedCells::operator()(const size_t p_idx,
ac.setNextParticle(p_idx, lc.infiniteParticles_.exchange(int_c(p_idx)));
} else
{
+ WALBERLA_ASSERT_GREATER(ac.getInteractionRadius(p_idx), 0_r, "Did you forget to set the interaction radius?");
{%- for dim in range(3) %}
+ WALBERLA_ASSERT_LESS(2_r * ac.getInteractionRadius(p_idx), lc.cellDiameter_[0], "Interaction radius is to large for this cell size. Contacts might get lost.");
int hash{{dim}} = static_cast<int>(std::floor((ac.getPosition(p_idx)[{{dim}}] - minCorner[{{dim}}]) * lc.invCellDiameter_[{{dim}}]));
{%- endfor %}
{%- for dim in range(3) %}
diff --git a/python/mesa_pd/templates/mpi/SyncGhostOwners.templ.cpp b/python/mesa_pd/templates/mpi/SyncGhostOwners.templ.cpp
index 7336ec1e8699c5a738062d5e0689e8d13556c17e..d1b5865b86c5bf65eda79861ec05d55ee02f5e5f 100644
--- a/python/mesa_pd/templates/mpi/SyncGhostOwners.templ.cpp
+++ b/python/mesa_pd/templates/mpi/SyncGhostOwners.templ.cpp
@@ -72,6 +72,7 @@ void SyncGhostOwners::updateAndMigrate( data::ParticleStorage& ps,
std::set<walberla::mpi::MPIRank> recvRanks; // potential message senders
for( auto pIt = ps.begin(); pIt != ps.end(); ++pIt)
{
+ WALBERLA_ASSERT_GREATER(pIt->getInteractionRadius(), 0_r, "Did you forget to set the interaction radius?");
if (isSet( pIt->getFlags(), GHOST))
{
if (!isSet( pIt->getFlags(), NON_COMMUNICATING) || syncNonCommunicatingBodies)
diff --git a/python/mesa_pd/templates/mpi/SyncNextNeighbors.templ.cpp b/python/mesa_pd/templates/mpi/SyncNextNeighbors.templ.cpp
index 3dd7ca3c40988c74193df432df99dcb43498fbc7..bd1d1b284ed45e1dd71dd8c54422841898d0f382 100644
--- a/python/mesa_pd/templates/mpi/SyncNextNeighbors.templ.cpp
+++ b/python/mesa_pd/templates/mpi/SyncNextNeighbors.templ.cpp
@@ -87,6 +87,8 @@ void SyncNextNeighbors::generateSynchronizationMessages(walberla::mpi::BufferSys
// position update
for( auto pIt = ps.begin(); pIt != ps.end(); )
{
+ WALBERLA_ASSERT_GREATER(pIt->getInteractionRadius(), 0_r, "Did you forget to set the interaction radius?");
+
//skip all ghost particles
if (data::particle_flags::isSet( pIt->getFlags(), data::particle_flags::GHOST))
{
diff --git a/src/mesa_pd/common/AABBConversion.h b/src/mesa_pd/common/AABBConversion.h
index 33896f25f731c381b866adbde962f51bf490682d..c105d1ed2b336b91c671e203bba786601b613ffd 100644
--- a/src/mesa_pd/common/AABBConversion.h
+++ b/src/mesa_pd/common/AABBConversion.h
@@ -32,6 +32,7 @@ namespace mesa_pd {
math::AABB getAABBFromInteractionRadius(const Vector3<real_t> & pos, const real_t interactionRadius )
{
+ WALBERLA_ASSERT_GREATER(interactionRadius, 0_r, "Did you forget to set the interaction radius?");
return math::AABB( pos[0]-interactionRadius, pos[1]-interactionRadius, pos[2]-interactionRadius,
pos[0]+interactionRadius, pos[1]+interactionRadius, pos[2]+interactionRadius );
}
diff --git a/src/mesa_pd/kernel/InsertParticleIntoLinkedCells.h b/src/mesa_pd/kernel/InsertParticleIntoLinkedCells.h
index a2cbd5c578c89d9da0937d2f86b617b4818d5013..088e0577ff09e3899d5af11e4ca99fb3dcbb5fcf 100644
--- a/src/mesa_pd/kernel/InsertParticleIntoLinkedCells.h
+++ b/src/mesa_pd/kernel/InsertParticleIntoLinkedCells.h
@@ -72,8 +72,12 @@ inline void InsertParticleIntoLinkedCells::operator()(const size_t p_idx, Access
ac.setNextParticle(p_idx, lc.infiniteParticles_.exchange(int_c(p_idx)));
} else
{
+ WALBERLA_ASSERT_GREATER(ac.getInteractionRadius(p_idx), 0_r, "Did you forget to set the interaction radius?");
+ WALBERLA_ASSERT_LESS(2_r * ac.getInteractionRadius(p_idx), lc.cellDiameter_[0], "Interaction radius is to large for this cell size. Contacts might get lost.");
int hash0 = static_cast<int>(std::floor((ac.getPosition(p_idx)[0] - minCorner[0]) * lc.invCellDiameter_[0]));
+ WALBERLA_ASSERT_LESS(2_r * ac.getInteractionRadius(p_idx), lc.cellDiameter_[0], "Interaction radius is to large for this cell size. Contacts might get lost.");
int hash1 = static_cast<int>(std::floor((ac.getPosition(p_idx)[1] - minCorner[1]) * lc.invCellDiameter_[1]));
+ WALBERLA_ASSERT_LESS(2_r * ac.getInteractionRadius(p_idx), lc.cellDiameter_[0], "Interaction radius is to large for this cell size. Contacts might get lost.");
int hash2 = static_cast<int>(std::floor((ac.getPosition(p_idx)[2] - minCorner[2]) * lc.invCellDiameter_[2]));
if (hash0 < 0) hash0 = 0;
if (hash0 >= lc.numCellsPerDim_[0]) hash0 = lc.numCellsPerDim_[0] - 1;
diff --git a/src/mesa_pd/kernel/InsertParticleIntoSparseLinkedCells.h b/src/mesa_pd/kernel/InsertParticleIntoSparseLinkedCells.h
index 6975337c415013d365f7a90b6fe3d060f2a8cf48..44af0f3b6e93f7ac81636010ac04956671392790 100644
--- a/src/mesa_pd/kernel/InsertParticleIntoSparseLinkedCells.h
+++ b/src/mesa_pd/kernel/InsertParticleIntoSparseLinkedCells.h
@@ -72,8 +72,12 @@ inline void InsertParticleIntoSparseLinkedCells::operator()(const size_t p_idx,
ac.setNextParticle(p_idx, lc.infiniteParticles_.exchange(int_c(p_idx)));
} else
{
+ WALBERLA_ASSERT_GREATER(ac.getInteractionRadius(p_idx), 0_r, "Did you forget to set the interaction radius?");
+ WALBERLA_ASSERT_LESS(2_r * ac.getInteractionRadius(p_idx), lc.cellDiameter_[0], "Interaction radius is to large for this cell size. Contacts might get lost.");
int hash0 = static_cast<int>(std::floor((ac.getPosition(p_idx)[0] - minCorner[0]) * lc.invCellDiameter_[0]));
+ WALBERLA_ASSERT_LESS(2_r * ac.getInteractionRadius(p_idx), lc.cellDiameter_[0], "Interaction radius is to large for this cell size. Contacts might get lost.");
int hash1 = static_cast<int>(std::floor((ac.getPosition(p_idx)[1] - minCorner[1]) * lc.invCellDiameter_[1]));
+ WALBERLA_ASSERT_LESS(2_r * ac.getInteractionRadius(p_idx), lc.cellDiameter_[0], "Interaction radius is to large for this cell size. Contacts might get lost.");
int hash2 = static_cast<int>(std::floor((ac.getPosition(p_idx)[2] - minCorner[2]) * lc.invCellDiameter_[2]));
if (hash0 < 0) hash0 = 0;
if (hash0 >= lc.numCellsPerDim_[0]) hash0 = lc.numCellsPerDim_[0] - 1;
diff --git a/src/mesa_pd/mpi/SyncGhostOwners.cpp b/src/mesa_pd/mpi/SyncGhostOwners.cpp
index e9eceaf88138b404d93d466283d45f0c7f4d9910..7bec109a44fad5322dff9d0021c8230e2c4f2901 100644
--- a/src/mesa_pd/mpi/SyncGhostOwners.cpp
+++ b/src/mesa_pd/mpi/SyncGhostOwners.cpp
@@ -72,6 +72,7 @@ void SyncGhostOwners::updateAndMigrate( data::ParticleStorage& ps,
std::set<walberla::mpi::MPIRank> recvRanks; // potential message senders
for( auto pIt = ps.begin(); pIt != ps.end(); ++pIt)
{
+ WALBERLA_ASSERT_GREATER(pIt->getInteractionRadius(), 0_r, "Did you forget to set the interaction radius?");
if (isSet( pIt->getFlags(), GHOST))
{
if (!isSet( pIt->getFlags(), NON_COMMUNICATING) || syncNonCommunicatingBodies)
diff --git a/src/mesa_pd/mpi/SyncNextNeighbors.cpp b/src/mesa_pd/mpi/SyncNextNeighbors.cpp
index 44f3946e792828e8729dd2967427f07227b79bda..3f9e5a2db9104be08571256c99562dc03d851aa6 100644
--- a/src/mesa_pd/mpi/SyncNextNeighbors.cpp
+++ b/src/mesa_pd/mpi/SyncNextNeighbors.cpp
@@ -87,6 +87,8 @@ void SyncNextNeighbors::generateSynchronizationMessages(walberla::mpi::BufferSys
// position update
for( auto pIt = ps.begin(); pIt != ps.end(); )
{
+ WALBERLA_ASSERT_GREATER(pIt->getInteractionRadius(), 0_r, "Did you forget to set the interaction radius?");
+
//skip all ghost particles
if (data::particle_flags::isSet( pIt->getFlags(), data::particle_flags::GHOST))
{
diff --git a/src/mesa_pd/mpi/SyncNextNeighborsBlockForest.cpp b/src/mesa_pd/mpi/SyncNextNeighborsBlockForest.cpp
index f4adbf79fae95e17734d7184115f8a2eab988023..a8528aa1a43e400ab426010858ea586e208c4164 100644
--- a/src/mesa_pd/mpi/SyncNextNeighborsBlockForest.cpp
+++ b/src/mesa_pd/mpi/SyncNextNeighborsBlockForest.cpp
@@ -111,6 +111,8 @@ void SyncNextNeighborsBlockForest::generateSynchronizationMessages(walberla::mpi
// position update
for( auto pIt = ps.begin(); pIt != ps.end(); )
{
+ WALBERLA_ASSERT_GREATER(pIt->getInteractionRadius(), 0_r, "Did you forget to set the interaction radius?");
+
//skip all ghost particles
if (data::particle_flags::isSet( pIt->getFlags(), data::particle_flags::GHOST))
{
diff --git a/tests/lbm_mesapd_coupling/mapping/ParticleMapping.cpp b/tests/lbm_mesapd_coupling/mapping/ParticleMapping.cpp
index ea60887f3121a7feb4a89bec8cb961d13ed35e00..b2ecf7107f455b229cde21f9f6c22742c10aec00 100644
--- a/tests/lbm_mesapd_coupling/mapping/ParticleMapping.cpp
+++ b/tests/lbm_mesapd_coupling/mapping/ParticleMapping.cpp
@@ -769,6 +769,7 @@ int main( int argc, char **argv )
{
mesa_pd::data::Particle&& p = *ps->create(true);
p.setPosition(halfSpacePosition);
+ p.setInteractionRadius(std::numeric_limits<real_t>::infinity());
p.setOwner(mpi::MPIManager::instance()->rank());
p.setShapeID(halfSpaceShape);
mesa_pd::data::particle_flags::set(p.getFlagsRef(), mesa_pd::data::particle_flags::INFINITE);
diff --git a/tests/lbm_mesapd_coupling/momentum_exchange_method/SettlingSphere.cpp b/tests/lbm_mesapd_coupling/momentum_exchange_method/SettlingSphere.cpp
index 2cc78a22f4511557e7e9f7743f1c36db8502f3c7..202601669557a91b5b8e2b2f26c8d0e36ffbccc3 100644
--- a/tests/lbm_mesapd_coupling/momentum_exchange_method/SettlingSphere.cpp
+++ b/tests/lbm_mesapd_coupling/momentum_exchange_method/SettlingSphere.cpp
@@ -281,6 +281,7 @@ void createPlaneSetup(const shared_ptr<mesa_pd::data::ParticleStorage> & ps, con
// create bounding planes
mesa_pd::data::Particle p0 = *ps->create(true);
p0.setPosition(simulationDomain.minCorner());
+ p0.setInteractionRadius(std::numeric_limits<real_t>::infinity());
p0.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(0,0,1) ));
p0.setOwner(mpi::MPIManager::instance()->rank());
p0.setType(0);
@@ -289,6 +290,7 @@ void createPlaneSetup(const shared_ptr<mesa_pd::data::ParticleStorage> & ps, con
mesa_pd::data::Particle p1 = *ps->create(true);
p1.setPosition(simulationDomain.maxCorner());
+ p1.setInteractionRadius(std::numeric_limits<real_t>::infinity());
p1.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(0,0,-1) ));
p1.setOwner(mpi::MPIManager::instance()->rank());
p1.setType(0);
@@ -297,6 +299,7 @@ void createPlaneSetup(const shared_ptr<mesa_pd::data::ParticleStorage> & ps, con
mesa_pd::data::Particle p2 = *ps->create(true);
p2.setPosition(simulationDomain.minCorner());
+ p2.setInteractionRadius(std::numeric_limits<real_t>::infinity());
p2.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(1,0,0) ));
p2.setOwner(mpi::MPIManager::instance()->rank());
p2.setType(0);
@@ -305,6 +308,7 @@ void createPlaneSetup(const shared_ptr<mesa_pd::data::ParticleStorage> & ps, con
mesa_pd::data::Particle p3 = *ps->create(true);
p3.setPosition(simulationDomain.maxCorner());
+ p3.setInteractionRadius(std::numeric_limits<real_t>::infinity());
p3.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(-1,0,0) ));
p3.setOwner(mpi::MPIManager::instance()->rank());
p3.setType(0);
@@ -313,6 +317,7 @@ void createPlaneSetup(const shared_ptr<mesa_pd::data::ParticleStorage> & ps, con
mesa_pd::data::Particle p4 = *ps->create(true);
p4.setPosition(simulationDomain.minCorner());
+ p4.setInteractionRadius(std::numeric_limits<real_t>::infinity());
p4.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(0,1,0) ));
p4.setOwner(mpi::MPIManager::instance()->rank());
p4.setType(0);
@@ -321,6 +326,7 @@ void createPlaneSetup(const shared_ptr<mesa_pd::data::ParticleStorage> & ps, con
mesa_pd::data::Particle p5 = *ps->create(true);
p5.setPosition(simulationDomain.maxCorner());
+ p5.setInteractionRadius(std::numeric_limits<real_t>::infinity());
p5.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(0,-1,0) ));
p5.setOwner(mpi::MPIManager::instance()->rank());
p5.setType(0);
diff --git a/tests/lbm_mesapd_coupling/utility/LubricationCorrection.cpp b/tests/lbm_mesapd_coupling/utility/LubricationCorrection.cpp
index bef4e3e4644ca6efcbdbda122166d4f4b3f81c0c..b96ab7fb01ff40faedcb7ce4fb4ccd2b758067b2 100644
--- a/tests/lbm_mesapd_coupling/utility/LubricationCorrection.cpp
+++ b/tests/lbm_mesapd_coupling/utility/LubricationCorrection.cpp
@@ -147,24 +147,28 @@ int main( int argc, char **argv )
mesa_pd::data::Particle&& p1 = *ps->create();
p1.setPosition(position1);
+ p1.setInteractionRadius(sphereRadius);
p1.setShapeID(sphereShape);
p1.setLinearVelocity(Vector3<real_t>(relativeVelocity,-relativeVelocity,real_t(0.5)*relativeVelocity));
auto idxSph1 = p1.getIdx();
mesa_pd::data::Particle&& p2 = *ps->create();
p2.setPosition(position2);
+ p2.setInteractionRadius(sphereRadius);
p2.setShapeID(sphereShape);
p2.setLinearVelocity(Vector3<real_t>(real_t(0),real_t(2),real_t(0)));
auto idxSph2 = p2.getIdx();
mesa_pd::data::Particle&& p3 = *ps->create();
p3.setPosition(position3);
+ p3.setInteractionRadius(sphereRadius);
p3.setShapeID(sphereShape);
p3.setLinearVelocity(Vector3<real_t>(real_t(2),real_t(0),real_t(2)));
auto idxSph3 = p3.getIdx();
mesa_pd::data::Particle&& p4 = *ps->create();
p4.setPosition(position4);
+ p4.setInteractionRadius(sphereRadius);
p4.setShapeID(sphereShape);
p4.setLinearVelocity(Vector3<real_t>(real_t(0),real_t(0),real_t(-0.5)*relativeVelocity));
auto idxSph4 = p4.getIdx();
@@ -228,24 +232,28 @@ int main( int argc, char **argv )
mesa_pd::data::Particle&& p1 = *ps->create();
p1.setPosition(position1);
+ p1.setInteractionRadius(sphereRadius);
p1.setShapeID(sphereShape);
p1.setLinearVelocity(Vector3<real_t>(relativeVelocity,-relativeVelocity,real_t(0.5)*relativeVelocity));
auto idxSph1 = p1.getIdx();
mesa_pd::data::Particle&& p2 = *ps->create();
p2.setPosition(position2);
+ p2.setInteractionRadius(sphereRadius);
p2.setShapeID(sphereShape);
p2.setLinearVelocity(Vector3<real_t>(real_t(0),real_t(2),real_t(0)));
auto idxSph2 = p2.getIdx();
mesa_pd::data::Particle&& p3 = *ps->create();
p3.setPosition(position3);
+ p3.setInteractionRadius(sphereRadius);
p3.setShapeID(sphereShape);
p3.setLinearVelocity(Vector3<real_t>(real_t(2),real_t(0),real_t(2)));
auto idxSph3 = p3.getIdx();
mesa_pd::data::Particle&& p4 = *ps->create();
p4.setPosition(position4);
+ p4.setInteractionRadius(sphereRadius);
p4.setShapeID(sphereShape);
p4.setLinearVelocity(Vector3<real_t>(real_t(0),real_t(0),real_t(-0.5)*relativeVelocity));
auto idxSph4 = p4.getIdx();
@@ -328,6 +336,7 @@ int main( int argc, char **argv )
mesa_pd::data::Particle&& p1 = *ps->create();
p1.setPosition(position1);
+ p1.setInteractionRadius(sphereRadius);
p1.setShapeID(sphereShape);
p1.setLinearVelocity(Vector3<real_t>(real_t(0),real_t(0),relativeVelocity));
auto idxSph1 = p1.getIdx();
@@ -335,29 +344,34 @@ int main( int argc, char **argv )
// mix up order to test Sph - HSp and HSp - Sph variants
mesa_pd::data::Particle&& pW = *ps->create(true);
pW.setPosition(wallPosition);
+ pW.setInteractionRadius(std::numeric_limits<real_t>::infinity());
pW.setShapeID(planeShape);
auto idxWall = pW.getIdx();
mesa_pd::data::Particle&& p2 = *ps->create();
p2.setPosition(position2);
+ p2.setInteractionRadius(sphereRadius);
p2.setShapeID(sphereShape);
p2.setLinearVelocity(Vector3<real_t>(real_t(0),real_t(0),-relativeVelocity));
auto idxSph2 = p2.getIdx();
mesa_pd::data::Particle&& p3 = *ps->create();
p3.setPosition(position3);
+ p3.setInteractionRadius(sphereRadius);
p3.setShapeID(sphereShape);
p3.setLinearVelocity(Vector3<real_t>(real_t(0),real_t(0),relativeVelocity));
auto idxSph3 = p3.getIdx();
mesa_pd::data::Particle&& p4 = *ps->create();
p4.setPosition(position4);
+ p4.setInteractionRadius(sphereRadius);
p4.setShapeID(sphereShape);
p4.setLinearVelocity(Vector3<real_t>(real_t(0),real_t(0),relativeVelocity));
auto idxSph4 = p4.getIdx();
mesa_pd::data::Particle&& p5 = *ps->create();
p5.setPosition(position5);
+ p5.setInteractionRadius(sphereRadius);
p5.setShapeID(sphereShape);
p5.setLinearVelocity(Vector3<real_t>(real_t(0),real_t(0),relativeVelocity));
auto idxSph5 = p5.getIdx();
diff --git a/tests/mesa_pd/ContactDetection.cpp b/tests/mesa_pd/ContactDetection.cpp
index 41e62d75f7a486ed08a3099c24c35576181ec09c..94d6c757514feb3adcf3fc10d77d06387c4e4f8d 100644
--- a/tests/mesa_pd/ContactDetection.cpp
+++ b/tests/mesa_pd/ContactDetection.cpp
@@ -15,13 +15,11 @@
//
//! \file
//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
//
//======================================================================================================================
-#include <mesa_pd/vtk/ParticleVtkOutput.h>
-
#include <mesa_pd/collision_detection/AnalyticContactDetection.h>
-#include <mesa_pd/collision_detection/GeneralContactDetection.h>
#include <mesa_pd/data/LinkedCells.h>
#include <mesa_pd/data/ParticleAccessor.h>
#include <mesa_pd/data/ParticleStorage.h>
@@ -31,10 +29,8 @@
#include <mesa_pd/kernel/InsertParticleIntoLinkedCells.h>
#include <mesa_pd/kernel/ParticleSelector.h>
#include <mesa_pd/mpi/ContactFilter.h>
-#include <mesa_pd/mpi/notifications/ForceTorqueNotification.h>
#include <mesa_pd/mpi/SyncNextNeighbors.h>
-#include <blockforest/BlockForest.h>
#include <blockforest/Initialization.h>
#include <core/Abort.h>
#include <core/Environment.h>
@@ -42,9 +38,7 @@
#include <core/mpi/Reduce.h>
#include <core/grid_generator/SCIterator.h>
#include <core/logging/Logging.h>
-#include <core/waLBerlaBuildInfo.h>
-#include <functional>
#include <memory>
#include <string>
@@ -68,7 +62,13 @@ private:
std::shared_ptr<data::ShapeStorage> ss_;
};
-int main( const int particlesPerAxis = 2, const real_t radius = real_t(0.9) )
+/*
+ * Tests linked cells
+ * Fully periodic, simple cubic sphere array is created
+ * Together with linked cells of size 2 * sphere spacing.
+ * This allows to analytically calculate how many particles should reside in each linked cell and how many contacts should be found and treated.
+ */
+int main( const int particlesPerAxisPerProcess = 2 )
{
using namespace walberla::timing;
@@ -80,15 +80,18 @@ int main( const int particlesPerAxis = 2, const real_t radius = real_t(0.9) )
math::seedRandomGenerator( static_cast<unsigned int>(1337 * walberla::mpi::MPIManager::instance()->worldRank()) );
- const real_t spacing(1.0);
+ const real_t radius = real_t(0.9);
+ const real_t generationSpacing = 1_r;
+ const int linkedCellMultipleOfGenerationSpacing = 2;
+ WALBERLA_CHECK_EQUAL(particlesPerAxisPerProcess % linkedCellMultipleOfGenerationSpacing, 0, "Only a multiple of " << linkedCellMultipleOfGenerationSpacing << " is allowed");
WALBERLA_LOG_INFO_ON_ROOT("*** BLOCKFOREST ***");
- const int centerParticles = particlesPerAxis * particlesPerAxis * particlesPerAxis;
- const int faceParticles = particlesPerAxis * particlesPerAxis * 6;
- const int edgeParticles = particlesPerAxis * 4 * 3;
+ const int centerParticles = particlesPerAxisPerProcess * particlesPerAxisPerProcess * particlesPerAxisPerProcess;
+ const int faceParticles = particlesPerAxisPerProcess * particlesPerAxisPerProcess * 6;
+ const int edgeParticles = particlesPerAxisPerProcess * 4 * 3;
const int cornerParticles = 8;
- WALBERLA_LOG_DEVEL_VAR_ON_ROOT(particlesPerAxis);
+ WALBERLA_LOG_DEVEL_VAR_ON_ROOT(particlesPerAxisPerProcess);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(centerParticles);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(faceParticles);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(edgeParticles);
@@ -96,14 +99,16 @@ int main( const int particlesPerAxis = 2, const real_t radius = real_t(0.9) )
// create forest
const real_t eps = real_c(0.01); //shift to make contact points unambigous
+ Vector3<uint_t>blocksPerDirection (2,2,2);
+ int numBlocks = int(blocksPerDirection[0] * blocksPerDirection[1] * blocksPerDirection[2]);
+ Vector3<bool> periodicity(true, true, true);
auto forest = blockforest::createBlockForest( math::AABB(real_t(0),
real_t(0),
real_t(0),
- real_c(particlesPerAxis) * real_t(2),
- real_c(particlesPerAxis) * real_t(2),
- real_c(particlesPerAxis) * real_t(2)).getTranslated(Vec3(eps,eps,eps)),
- Vector3<uint_t>(2,2,2),
- Vector3<bool>(true, true, true) );
+ real_c(particlesPerAxisPerProcess) * generationSpacing * real_c(blocksPerDirection[0]),
+ real_c(particlesPerAxisPerProcess) * generationSpacing * real_c(blocksPerDirection[1]),
+ real_c(particlesPerAxisPerProcess) * generationSpacing * real_c(blocksPerDirection[2])).getTranslated(Vec3(eps,eps,eps)),
+ blocksPerDirection, periodicity);
domain::BlockForestDomain domain(forest);
auto localDomain = forest->begin()->getAABB();
@@ -115,13 +120,14 @@ int main( const int particlesPerAxis = 2, const real_t radius = real_t(0.9) )
auto ps = std::make_shared<data::ParticleStorage>(100);
auto ss = std::make_shared<data::ShapeStorage>();
ParticleAccessorWithShape accessor(ps, ss);
- data::LinkedCells lc(localDomain.getExtended(spacing), spacing );
+ const real_t linkedCellSize = real_c(linkedCellMultipleOfGenerationSpacing) * generationSpacing;
+ data::LinkedCells lc(localDomain.getExtended(linkedCellSize), linkedCellSize );
auto smallSphere = ss->create<data::Sphere>( radius );
ss->shapes[smallSphere]->updateMassAndInertia(real_t(2707));
for (auto& iBlk : *forest)
{
- for (auto pt : grid_generator::SCGrid(iBlk.getAABB(), Vector3<real_t>(spacing, spacing, spacing) * real_c(0.5), spacing))
+ for (auto pt : grid_generator::SCGrid(iBlk.getAABB(), Vector3<real_t>(generationSpacing, generationSpacing, generationSpacing) * real_c(0.5), generationSpacing))
{
WALBERLA_CHECK(iBlk.getAABB().contains(pt));
@@ -133,7 +139,7 @@ int main( const int particlesPerAxis = 2, const real_t radius = real_t(0.9) )
}
}
int64_t numParticles = int64_c(ps->size());
- walberla::mpi::reduceInplace(numParticles, walberla::mpi::SUM);
+ walberla::mpi::allReduceInplace(numParticles, walberla::mpi::SUM);
WALBERLA_LOG_INFO_ON_ROOT("#particles created: " << numParticles);
WALBERLA_LOG_INFO_ON_ROOT("*** SETUP - END ***");
@@ -146,14 +152,32 @@ int main( const int particlesPerAxis = 2, const real_t radius = real_t(0.9) )
// initial sync
WALBERLA_CHECK_EQUAL(ps->size(), centerParticles);
SNN(*ps, domain);
- WALBERLA_CHECK_EQUAL(ps->size(), centerParticles + faceParticles + edgeParticles + cornerParticles);
+ int numKnownParticles = centerParticles + faceParticles + edgeParticles + cornerParticles;
+ WALBERLA_CHECK_EQUAL(ps->size(), numKnownParticles);
lc.clear();
ps->forEachParticle(true, kernel::SelectAll(), accessor, ipilc, accessor, lc);
- WALBERLA_CHECK_EQUAL(lc.cells_.size(), centerParticles + faceParticles + edgeParticles + cornerParticles);
- int cell0 = 0;
- int cell1 = 0;
+ int interiorCellsPerDirection = particlesPerAxisPerProcess / linkedCellMultipleOfGenerationSpacing;
+
+ int numLinkedCellsPerProcess = (interiorCellsPerDirection+2) * (interiorCellsPerDirection+2) * (interiorCellsPerDirection+2);
+ WALBERLA_CHECK_EQUAL(lc.cells_.size(), numLinkedCellsPerProcess);
+
+ const int numLinkedCellsCenter = interiorCellsPerDirection * interiorCellsPerDirection * interiorCellsPerDirection;
+ const int numLinkedCellsFace = interiorCellsPerDirection * interiorCellsPerDirection * 6;
+ const int numLinkedCellsEdge = interiorCellsPerDirection * 12;
+ const int numLinkedCellsCorner = 8;
+
+ const int numParticlesPerCenterCell = centerParticles / numLinkedCellsCenter;
+ const int numParticlesPerFaceCell = faceParticles / numLinkedCellsFace;
+ const int numParticlesPerEdgeCell = edgeParticles / numLinkedCellsEdge;
+ const int numParticlesPerCornerCell = cornerParticles / numLinkedCellsCorner;
+
+ int cellsCenter = 0;
+ int cellsFace = 0;
+ int cellsEdge = 0;
+ int cellsCorner = 0;
+
for (const auto& idx : lc.cells_)
{
int particleCounter = 0;
@@ -163,18 +187,18 @@ int main( const int particlesPerAxis = 2, const real_t radius = real_t(0.9) )
++particleCounter;
p_idx = ps->getNextParticle(uint_c(p_idx));
}
- if (particleCounter==0)
- {
- ++cell0;
- } else if (particleCounter==1)
- {
- ++cell1;
- } else {
- WALBERLA_CHECK(false);
- }
+
+ if (particleCounter == numParticlesPerCenterCell) ++cellsCenter;
+ else if (particleCounter == numParticlesPerFaceCell) ++cellsFace;
+ else if (particleCounter == numParticlesPerEdgeCell) ++cellsEdge;
+ else if (particleCounter == numParticlesPerCornerCell) ++cellsCorner;
+ else WALBERLA_CHECK(false, "Linked cell found with unexpected number of particles " << particleCounter);
+
}
- WALBERLA_CHECK_EQUAL(cell0, 0);
- WALBERLA_CHECK_EQUAL(cell1, centerParticles + faceParticles + edgeParticles + cornerParticles); //ghost particles only at face, no ghost particle at edges and corners
+ WALBERLA_CHECK_EQUAL(cellsCenter, numLinkedCellsCenter);
+ WALBERLA_CHECK_EQUAL(cellsFace, numLinkedCellsFace);
+ WALBERLA_CHECK_EQUAL(cellsEdge, numLinkedCellsEdge);
+ WALBERLA_CHECK_EQUAL(cellsCorner, numLinkedCellsCorner);
std::atomic<int64_t> contactsChecked (0);
std::atomic<int64_t> contactsDetected(0);
@@ -201,9 +225,10 @@ int main( const int particlesPerAxis = 2, const real_t radius = real_t(0.9) )
WALBERLA_LOG_DEVEL_ON_ROOT( "contacts checked/detected/treated: " << contactsChecked << " / " << contactsDetected << " / " << contactsTreated );
- WALBERLA_CHECK_EQUAL(contactsChecked, (centerParticles*26 + faceParticles*17 + edgeParticles*11 + cornerParticles*7) / 2 );
+ WALBERLA_CHECK_LESS_EQUAL(contactsChecked, numKnownParticles * numKnownParticles / 2 ); // has to be better than n^2 checking
WALBERLA_CHECK_EQUAL(contactsDetected, (centerParticles*26 + faceParticles*17 + edgeParticles*11 + cornerParticles*7) / 2 );
- //TODO: ContactsTreated
+ const int totalNumberPairWiseContactsInDomain = 26 * int(numParticles) / 2; //fully-periodic simple cubic grid
+ WALBERLA_CHECK_EQUAL(contactsTreated, totalNumberPairWiseContactsInDomain / numBlocks);
return EXIT_SUCCESS;
}
@@ -216,8 +241,6 @@ int main( int argc, char* argv[] )
walberla::Environment env(argc, argv);
WALBERLA_UNUSED(env);
walberla::mesa_pd::main( 2 );
- walberla::mesa_pd::main( 3 );
walberla::mesa_pd::main( 4 );
- walberla::mesa_pd::main( 5 );
return EXIT_SUCCESS;
}
diff --git a/tests/mesa_pd/kernel/GenerateAnalyticContacts.cpp b/tests/mesa_pd/kernel/GenerateAnalyticContacts.cpp
index 4762aed06e159d322033faace72d48c807fedc2d..fcd2b82c8ebaef80843686cada57d5d96a764a5a 100644
--- a/tests/mesa_pd/kernel/GenerateAnalyticContacts.cpp
+++ b/tests/mesa_pd/kernel/GenerateAnalyticContacts.cpp
@@ -79,7 +79,7 @@ int main( int argc, char ** argv )
auto ps = std::make_shared<data::ParticleStorage>(100);
auto ss = std::make_shared<data::ShapeStorage>();
ParticleAccessorWithShape accessor(ps, ss);
- data::LinkedCells lc(math::AABB(-1,-1,-1,4,4,4), real_t(1));
+ data::LinkedCells lc(math::AABB(-1,-1,-1,4,4,4), real_t(1.3));
//initialize particles
const real_t radius = real_t(0.6);
diff --git a/tests/mesa_pd/kernel/GenerateLinkedCells.cpp b/tests/mesa_pd/kernel/GenerateLinkedCells.cpp
index 3ee5183cfcce026112cb192929db1de10bfff5aa..2aaf342aeec27dc3b5699bfa3844e6c2a2bb2e2b 100644
--- a/tests/mesa_pd/kernel/GenerateLinkedCells.cpp
+++ b/tests/mesa_pd/kernel/GenerateLinkedCells.cpp
@@ -66,6 +66,7 @@ int main( int argc, char ** argv )
{
data::Particle&& p = *storage->create();
p.getPositionRef() = (*it);
+ p.setInteractionRadius(0.1_r);
p.setShapeID(smallSphere);
}
diff --git a/tests/mesa_pd/kernel/LinkedCellsVsBruteForce.cpp b/tests/mesa_pd/kernel/LinkedCellsVsBruteForce.cpp
index c52e39327afae37a42460cdcb2510bd5f7195002..8743d6ac3bcf21b8a2d9cf82c6d3defed3405c01 100644
--- a/tests/mesa_pd/kernel/LinkedCellsVsBruteForce.cpp
+++ b/tests/mesa_pd/kernel/LinkedCellsVsBruteForce.cpp
@@ -100,7 +100,7 @@ int main( int argc, char ** argv )
//init data structures
auto ps = std::make_shared<data::ParticleStorage>(100);
auto ss = std::make_shared<data::ShapeStorage>();
- data::LinkedCells lc(blk.getAABB(), real_t(1));
+ data::LinkedCells lc(blk.getAABB(), real_t(1.1));
std::vector<collision_detection::AnalyticContactDetection> cs1(100);
std::vector<collision_detection::AnalyticContactDetection> cs2(100);