diff --git a/python/mesa_pd/templates/data/ShapeStorage.templ.h b/python/mesa_pd/templates/data/ShapeStorage.templ.h index 230f6a923ab392b90448961744e825a1b6905c8f..772628ea6ae98307622d83a2c373f6457316d287 100644 --- a/python/mesa_pd/templates/data/ShapeStorage.templ.h +++ b/python/mesa_pd/templates/data/ShapeStorage.templ.h @@ -58,8 +58,10 @@ struct ShapeStorage ReturnType doubleDispatch( ParticleStorage& ps, size_t idx, size_t idy, func& f ); }; //Make sure that no two different shapes have the same unique identifier! -{%- for shape in particle.shapes[1:] %} -static_assert( {{particle.shapes[0]}}::SHAPE_TYPE != {{shape}}::SHAPE_TYPE, "Shape types have to be different!" ); +{%- for shape1 in particle.shapes %} +{%- for shape2 in particle.shapes[loop.index:] %} +static_assert( {{shape1}}::SHAPE_TYPE != {{shape2}}::SHAPE_TYPE, "Shape types have to be different!" ); +{%- endfor %} {%- endfor %} template <typename ShapeT, typename... Args> diff --git a/src/mesa_pd/data/ShapeStorage.h b/src/mesa_pd/data/ShapeStorage.h index f7e158d0e40e9bc9c6f16e678b74001b7b27925c..a8baee5c713ddf991286a6e6ef6bc21055d4affb 100644 --- a/src/mesa_pd/data/ShapeStorage.h +++ b/src/mesa_pd/data/ShapeStorage.h @@ -64,6 +64,12 @@ static_assert( Sphere::SHAPE_TYPE != HalfSpace::SHAPE_TYPE, "Shape types have to static_assert( Sphere::SHAPE_TYPE != CylindricalBoundary::SHAPE_TYPE, "Shape types have to be different!" ); static_assert( Sphere::SHAPE_TYPE != Box::SHAPE_TYPE, "Shape types have to be different!" ); static_assert( Sphere::SHAPE_TYPE != Ellipsoid::SHAPE_TYPE, "Shape types have to be different!" ); +static_assert( HalfSpace::SHAPE_TYPE != CylindricalBoundary::SHAPE_TYPE, "Shape types have to be different!" ); +static_assert( HalfSpace::SHAPE_TYPE != Box::SHAPE_TYPE, "Shape types have to be different!" ); +static_assert( HalfSpace::SHAPE_TYPE != Ellipsoid::SHAPE_TYPE, "Shape types have to be different!" ); +static_assert( CylindricalBoundary::SHAPE_TYPE != Box::SHAPE_TYPE, "Shape types have to be different!" ); +static_assert( CylindricalBoundary::SHAPE_TYPE != Ellipsoid::SHAPE_TYPE, "Shape types have to be different!" ); +static_assert( Box::SHAPE_TYPE != Ellipsoid::SHAPE_TYPE, "Shape types have to be different!" ); template <typename ShapeT, typename... Args> size_t ShapeStorage::create(Args&&... args) diff --git a/src/mesa_pd/data/shape/Box.h b/src/mesa_pd/data/shape/Box.h index 7cc5606785586fea4c1156a3bb796b5fc34a1456..cb67d2d6306c9f09d18aa3f26ae893d4754233e5 100644 --- a/src/mesa_pd/data/shape/Box.h +++ b/src/mesa_pd/data/shape/Box.h @@ -61,7 +61,10 @@ void Box::updateMassAndInertia(const real_t density) edgeLength_[0]*edgeLength_[0] + edgeLength_[2]*edgeLength_[2] , edgeLength_[0]*edgeLength_[0] + edgeLength_[1]*edgeLength_[1] ) * (m / static_cast<real_t>( 12 )); + mass_ = m; invMass_ = real_t(1.0) / m; + + inertiaBF_ = I; invInertiaBF_ = I.getInverse(); } diff --git a/src/mesa_pd/data/shape/CylindricalBoundary.h b/src/mesa_pd/data/shape/CylindricalBoundary.h index abe3bca759184a7f8f167ebdc2d5cffb67a7634c..81a04d7e47daf073d6a09616f85369f209838d15 100644 --- a/src/mesa_pd/data/shape/CylindricalBoundary.h +++ b/src/mesa_pd/data/shape/CylindricalBoundary.h @@ -59,8 +59,11 @@ private: inline void CylindricalBoundary::updateMassAndInertia(const real_t /*density*/) { - invMass_ = real_t(0.0); - invInertiaBF_ = Mat3(real_t(0.0)); + mass_ = std::numeric_limits<real_t>::infinity(); + invMass_ = real_t(0); + + inertiaBF_ = Mat3(std::numeric_limits<real_t>::infinity()); + invInertiaBF_ = Mat3(real_t(0)); } inline diff --git a/src/mesa_pd/data/shape/Ellipsoid.h b/src/mesa_pd/data/shape/Ellipsoid.h index 8ff9411852f7f08de766a5409b31cc5e68a05d5d..7101209cbc2d29e7d9f5fb394686545af0169662 100644 --- a/src/mesa_pd/data/shape/Ellipsoid.h +++ b/src/mesa_pd/data/shape/Ellipsoid.h @@ -61,7 +61,9 @@ void Ellipsoid::updateMassAndInertia(const real_t density) real_c(0.2) * m * (semiAxes_[2] * semiAxes_[2] + semiAxes_[0] * semiAxes_[0]), real_c(0.2) * m * (semiAxes_[0] * semiAxes_[0] + semiAxes_[1] * semiAxes_[1])); + mass_ = m; invMass_ = real_t(1.0) / m; + inertiaBF_ = I; invInertiaBF_ = I.getInverse(); } diff --git a/tests/mesa_pd/mpi/ShapePackUnpack.cpp b/tests/mesa_pd/mpi/ShapePackUnpack.cpp index a26177f4896fc66680279ea4e197f72a3025fe34..3efe46404386cc2a9591c403e35825dfe4b1a157 100644 --- a/tests/mesa_pd/mpi/ShapePackUnpack.cpp +++ b/tests/mesa_pd/mpi/ShapePackUnpack.cpp @@ -37,6 +37,16 @@ namespace walberla { using namespace walberla::mesa_pd; +void checkIdenticalOrInf(const Mat3& mat0, const Mat3& mat1) { + for (uint_t i = 0; i <= 8; i++) { + WALBERLA_CHECK((math::isinf(mat0[i]) && math::isinf(mat1[i])) || realIsIdentical(mat0[i], mat1[i]), "Matrices don't match: " << mat0 << " vs. " << mat1); + } +} + +void checkIdenticalOrInf(real_t a, real_t b) { + WALBERLA_CHECK((math::isinf(a) && math::isinf(b)) || realIsIdentical(a, b), "Values don't match: " << a << " vs. " << b); +} + void checkBox() { using namespace walberla::mpi; @@ -57,8 +67,8 @@ void checkBox() WALBERLA_CHECK_EQUAL(bs1->getShapeType(), data::Box::SHAPE_TYPE); WALBERLA_CHECK_IDENTICAL(bs0->getMass(), bs1->getMass()); WALBERLA_CHECK_IDENTICAL(bs0->getInvMass(), bs1->getInvMass()); - WALBERLA_CHECK_IDENTICAL(bs0->getInertiaBF(), bs1->getInertiaBF()); - WALBERLA_CHECK_IDENTICAL(bs0->getInvInertiaBF(), bs1->getInvInertiaBF()); + checkIdenticalOrInf(bs0->getInertiaBF(), bs1->getInertiaBF()); + checkIdenticalOrInf(bs0->getInvInertiaBF(), bs1->getInvInertiaBF()); auto bx0 = static_cast<data::Box*> (bs0.get()); auto bx1 = static_cast<data::Box*> (bs1.get()); @@ -87,8 +97,8 @@ void checkCylindricalBoundary() WALBERLA_CHECK_EQUAL(bs1->getShapeType(), data::CylindricalBoundary::SHAPE_TYPE); WALBERLA_CHECK_IDENTICAL(bs0->getMass(), bs1->getMass()); WALBERLA_CHECK_IDENTICAL(bs0->getInvMass(), bs1->getInvMass()); - WALBERLA_CHECK_IDENTICAL(bs0->getInertiaBF(), bs1->getInertiaBF()); - WALBERLA_CHECK_IDENTICAL(bs0->getInvInertiaBF(), bs1->getInvInertiaBF()); + checkIdenticalOrInf(bs0->getInertiaBF(), bs1->getInertiaBF()); + checkIdenticalOrInf(bs0->getInvInertiaBF(), bs1->getInvInertiaBF()); auto cb0 = static_cast<data::CylindricalBoundary*> (bs0.get()); auto cb1 = static_cast<data::CylindricalBoundary*> (bs1.get()); @@ -117,8 +127,8 @@ void checkEllipsoid() WALBERLA_CHECK_EQUAL(bs1->getShapeType(), data::Ellipsoid::SHAPE_TYPE); WALBERLA_CHECK_IDENTICAL(bs0->getMass(), bs1->getMass()); WALBERLA_CHECK_IDENTICAL(bs0->getInvMass(), bs1->getInvMass()); - WALBERLA_CHECK_IDENTICAL(bs0->getInertiaBF(), bs1->getInertiaBF()); - WALBERLA_CHECK_IDENTICAL(bs0->getInvInertiaBF(), bs1->getInvInertiaBF()); + checkIdenticalOrInf(bs0->getInertiaBF(), bs1->getInertiaBF()); + checkIdenticalOrInf(bs0->getInvInertiaBF(), bs1->getInvInertiaBF()); auto el0 = static_cast<data::Ellipsoid*> (bs0.get()); auto el1 = static_cast<data::Ellipsoid*> (bs1.get()); @@ -144,9 +154,9 @@ void checkHalfSpace() WALBERLA_LOG_INFO("checking half space"); WALBERLA_CHECK_EQUAL(bs0->getShapeType(), data::HalfSpace::SHAPE_TYPE); WALBERLA_CHECK_EQUAL(bs1->getShapeType(), data::HalfSpace::SHAPE_TYPE); -// WALBERLA_CHECK_IDENTICAL(bs0->getMass(), bs1->getMass()); //INF + checkIdenticalOrInf(bs0->getMass(), bs1->getMass()); //INF WALBERLA_CHECK_IDENTICAL(bs0->getInvMass(), bs1->getInvMass()); -// WALBERLA_CHECK_IDENTICAL(bs0->getInertiaBF(), bs1->getInertiaBF()); //INF + checkIdenticalOrInf(bs0->getInertiaBF(), bs1->getInertiaBF()); //INF WALBERLA_CHECK_IDENTICAL(bs0->getInvInertiaBF(), bs1->getInvInertiaBF()); auto hs0 = static_cast<data::HalfSpace*> (bs0.get()); @@ -204,4 +214,4 @@ int main( int argc, char ** argv ) int main( int argc, char ** argv ) { return walberla::main(argc, argv); -} +} \ No newline at end of file