Merge branch 'mesapd_box_and_shapetype_fixes' into 'master'

Fixes for MESA-PD: Shape type duplication checks and mass/inertia on some shapes

See merge request walberla/walberla!263

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>

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)

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();
 }
 

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

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();
 }
 

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