pe_coupling: New functionality, adapted test case

src/pe_coupling/utility/BodiesForceAndTorqueContainer.h → src/pe_coupling/utility/BodiesForceTorqueContainer.h

@@ -13,7 +13,7 @@
 //  You should have received a copy of the GNU General Public License along
 //  with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
 //
-//! \file BodiesForceAndTorqueContainer.h
+//! \file BodiesForceTorqueContainer.h
 //! \ingroup pe_coupling
 //! \author Christoph Rettinger <christoph.rettinger@fau.de>
 //
@@ -32,11 +32,11 @@
 namespace walberla {
 namespace pe_coupling {
 
-class BodiesForceAndTorqueContainer
+class BodiesForceTorqueContainer
 {  
 public:
 
-   BodiesForceAndTorqueContainer( const shared_ptr<StructuredBlockStorage> & blockStorage, const BlockDataID & bodyStorageID )
+   BodiesForceTorqueContainer( const shared_ptr<StructuredBlockStorage> & blockStorage, const BlockDataID & bodyStorageID )
    : blockStorage_( blockStorage ), bodyStorageID_( bodyStorageID )
      { }
 
@@ -60,7 +60,7 @@ public:
       {
          for( auto bodyIt = pe::BodyIterator::begin(*blockIt, bodyStorageID_); bodyIt != pe::BodyIterator::end(); ++bodyIt )
          {
-            auto & f = bodyForceAndTorqueMap_[ bodyIt->getSystemID() ];
+            auto & f = bodyForceTorqueMap_[ bodyIt->getSystemID() ];
 
             // only add if body has not been added already before (from another block)
             if( f.empty() )
@@ -94,7 +94,8 @@ public:
       {
          for( auto bodyIt = pe::LocalBodyIterator::begin(*blockIt, bodyStorageID_); bodyIt != pe::LocalBodyIterator::end(); ++bodyIt )
          {
-            const auto &f = bodyForceAndTorqueMap_[bodyIt->getSystemID()];
+            const auto &f = bodyForceTorqueMap_[bodyIt->getSystemID()];
+            WALBERLA_ASSERT( !f.empty(), "When attempting to set force/torque on local body " << bodyIt->getSystemID() << " at position " << bodyIt->getPosition() << ", body was not found in map!");
             bodyIt->addForce ( f[0], f[1], f[2] );
             bodyIt->addTorque( f[3], f[4], f[5] );
          }
@@ -103,14 +104,38 @@ public:
 
    void clear()
    {
-      bodyForceAndTorqueMap_.clear();
+      bodyForceTorqueMap_.clear();
+   }
+
+   void swap( BodiesForceTorqueContainer & other )
+   {
+      std::swap( this->bodyForceTorqueMap_, other.bodyForceTorqueMap_ );
    }
 
 private:
 
    shared_ptr<StructuredBlockStorage> blockStorage_;
    const BlockDataID bodyStorageID_;
-   std::map< walberla::id_t, std::vector<real_t> > bodyForceAndTorqueMap_;
+   std::map< walberla::id_t, std::vector<real_t> > bodyForceTorqueMap_;
+};
+
+
+class BodyContainerSwapper
+{
+public:
+   BodyContainerSwapper( const shared_ptr<BodiesForceTorqueContainer> & cont1,
+                         const shared_ptr<BodiesForceTorqueContainer> & cont2 )
+   : cont1_( cont1 ), cont2_( cont2 )
+   { }
+
+   void operator()()
+   {
+      cont1_->swap( *cont2_ );
+   }
+
+private:
+   shared_ptr<BodiesForceTorqueContainer> cont1_;
+   shared_ptr<BodiesForceTorqueContainer> cont2_;
 };
 
 } // namespace pe_coupling

src/pe_coupling/utility/ForceOnBodiesAdder.h

+//======================================================================================================================
+//
+//  This file is part of waLBerla. waLBerla is free software: you can 
+//  redistribute it and/or modify it under the terms of the GNU General Public
+//  License as published by the Free Software Foundation, either version 3 of 
+//  the License, or (at your option) any later version.
+//  
+//  waLBerla is distributed in the hope that it will be useful, but WITHOUT 
+//  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 
+//  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License 
+//  for more details.
+//  
+//  You should have received a copy of the GNU General Public License along
+//  with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file ForceOnBodiesAdder.h
+//! \ingroup pe_coupling
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "core/math/Vector3.h"
+#include "domain_decomposition/StructuredBlockStorage.h"
+#include "pe/rigidbody/BodyIterators.h"
+
+
+namespace walberla {
+namespace pe_coupling {
+
+class ForceOnBodiesAdder
+{  
+public:
+
+   ForceOnBodiesAdder( const shared_ptr<StructuredBlockStorage> & blockStorage, const BlockDataID & bodyStorageID,
+                       const Vector3<real_t> & force )
+   : blockStorage_( blockStorage ), bodyStorageID_( bodyStorageID ), force_( force )
+     { }
+
+   // set a constant force on all (only local, to avoid force duplication) bodies
+   void operator()()
+   {
+      for( auto blockIt = blockStorage_->begin(); blockIt != blockStorage_->end(); ++blockIt )
+      {
+         for( auto bodyIt = pe::LocalBodyIterator::begin( *blockIt, bodyStorageID_); bodyIt != pe::LocalBodyIterator::end(); ++bodyIt )
+         {
+            bodyIt->addForce ( force_ );
+         }
+      }
+   }
+
+   void updateForce( const Vector3<real_t> & newForce )
+   {
+      force_ = newForce;
+   }
+
+private:
+
+   shared_ptr<StructuredBlockStorage> blockStorage_;
+   const BlockDataID bodyStorageID_;
+   Vector3<real_t> force_;
+};
+
+} // namespace pe_coupling
+} // namespace walberla

src/pe_coupling/utility/ForceTorqueOnBodiesResetter.h

+//======================================================================================================================
+//
+//  This file is part of waLBerla. waLBerla is free software: you can 
+//  redistribute it and/or modify it under the terms of the GNU General Public
+//  License as published by the Free Software Foundation, either version 3 of 
+//  the License, or (at your option) any later version.
+//  
+//  waLBerla is distributed in the hope that it will be useful, but WITHOUT 
+//  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 
+//  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License 
+//  for more details.
+//  
+//  You should have received a copy of the GNU General Public License along
+//  with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file ForceTorqueOnBodiesResetter.h
+//! \ingroup pe_coupling
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "core/math/Vector3.h"
+#include "domain_decomposition/StructuredBlockStorage.h"
+#include "pe/rigidbody/BodyIterators.h"
+
+
+namespace walberla {
+namespace pe_coupling {
+
+class ForceTorqueOnBodiesResetter
+{  
+public:
+
+   ForceTorqueOnBodiesResetter( const shared_ptr<StructuredBlockStorage> & blockStorage, const BlockDataID & bodyStorageID )
+   : blockStorage_( blockStorage ), bodyStorageID_( bodyStorageID )
+     { }
+
+   // resets forces and torques on all (local and remote) bodies
+   void operator()()
+   {
+      for( auto blockIt = blockStorage_->begin(); blockIt != blockStorage_->end(); ++blockIt )
+      {
+         for( auto bodyIt = pe::BodyIterator::begin( *blockIt, bodyStorageID_); bodyIt != pe::BodyIterator::end(); ++bodyIt )
+         {
+            bodyIt->resetForceAndTorque();
+         }
+      }
+   }
+
+
+private:
+
+   shared_ptr<StructuredBlockStorage> blockStorage_;
+   const BlockDataID bodyStorageID_;
+};
+
+} // namespace pe_coupling
+} // namespace walberla

src/pe_coupling/utility/ForceTorqueOnBodiesScaler.h

+//======================================================================================================================
+//
+//  This file is part of waLBerla. waLBerla is free software: you can 
+//  redistribute it and/or modify it under the terms of the GNU General Public
+//  License as published by the Free Software Foundation, either version 3 of 
+//  the License, or (at your option) any later version.
+//  
+//  waLBerla is distributed in the hope that it will be useful, but WITHOUT 
+//  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 
+//  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License 
+//  for more details.
+//  
+//  You should have received a copy of the GNU General Public License along
+//  with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file ForceTorqueOnBodiesScaler.h
+//! \ingroup pe_coupling
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "core/math/Vector3.h"
+#include "domain_decomposition/StructuredBlockStorage.h"
+#include "pe/rigidbody/BodyIterators.h"
+
+
+namespace walberla {
+namespace pe_coupling {
+
+// scales force/torquew on all bodies (local and remote) by a constant scalar value
+// can e.g. be used to average the force/torque over two time steps
+class ForceTorqueOnBodiesScaler
+{  
+public:
+
+   ForceTorqueOnBodiesScaler( const shared_ptr<StructuredBlockStorage> & blockStorage, const BlockDataID & bodyStorageID,
+                              const real_t & scalingFactor )
+   : blockStorage_( blockStorage ), bodyStorageID_( bodyStorageID ), scalingFactor_( scalingFactor )
+     { }
+
+   // resets forces and torques on all (local and remote) bodies
+   void operator()()
+   {
+      Vector3<real_t> force(real_t(0));
+      Vector3<real_t> torque(real_t(0));
+      for( auto blockIt = blockStorage_->begin(); blockIt != blockStorage_->end(); ++blockIt )
+      {
+         for( auto bodyIt = pe::BodyIterator::begin( *blockIt, bodyStorageID_); bodyIt != pe::BodyIterator::end(); ++bodyIt )
+         {
+            force  = scalingFactor_ * bodyIt->getForce();
+            torque = scalingFactor_ * bodyIt->getTorque();
+
+            bodyIt->resetForceAndTorque();
+
+            bodyIt->setForce ( force );
+            bodyIt->setTorque( torque );
+         }
+      }
+   }
+
+
+private:
+
+   shared_ptr<StructuredBlockStorage> blockStorage_;
+   const BlockDataID bodyStorageID_;
+   const real_t scalingFactor_;
+};
+
+} // namespace pe_coupling
+} // namespace walberla

src/pe_coupling/utility/TimeStep.h

@@ -124,6 +124,7 @@ public:
                   for( auto bodyIt = pe::LocalBodyIterator::begin(*blockIt, bodyStorageID_); bodyIt != pe::LocalBodyIterator::end(); ++bodyIt )
                   {
                      const auto & f = forceMap[ bodyIt->getSystemID() ];
+                     WALBERLA_ASSERT( !f.empty(), "When attempting to set force/torque on local body " << bodyIt->getSystemID() << " at position " << bodyIt->getPosition() << ", body was not found in map!");
                      bodyIt->addForce ( f[0], f[1], f[2] );
                      bodyIt->addTorque( f[3], f[4], f[5] );
                   }

src/pe_coupling/utility/all.h

@@ -22,6 +22,9 @@
 
 #pragma once
 
-#include "BodiesForceAndTorqueContainer.h"
+#include "BodiesForceTorqueContainer.h"
+#include "ForceOnBodiesAdder.h"
+#include "ForceTorqueOnBodiesResetter.h"
+#include "ForceTorqueOnBodiesScaler.h"
 #include "LubricationCorrection.h"
 #include "TimeStep.h"

tests/pe_coupling/momentum_exchange_method/DragForceSphereMEM.cpp

@@ -58,6 +58,7 @@
 
 #include "pe_coupling/mapping/all.h"
 #include "pe_coupling/momentum_exchange_method/all.h"
+#include "pe_coupling/utility/all.h"
 
 #include <vector>
 #include <iomanip>
@@ -331,30 +332,6 @@ private:
 
 };
 
-
-class ResetForce
-{
-public:
-   ResetForce( const shared_ptr< StructuredBlockStorage > & blocks,
-               const BlockDataID & bodyStorageID )
-   : blocks_( blocks ), bodyStorageID_( bodyStorageID )
-   { }
-
-   void operator()()
-   {
-      for( auto blockIt = blocks_->begin(); blockIt != blocks_->end(); ++blockIt )
-      {
-         for( auto bodyIt = pe::BodyIterator::begin( *blockIt, bodyStorageID_); bodyIt != pe::BodyIterator::end(); ++bodyIt )
-         {
-            bodyIt->resetForceAndTorque();
-         }
-      }
-   }
-private:
-   shared_ptr< StructuredBlockStorage > blocks_;
-   const BlockDataID bodyStorageID_;
-};
-
 class PDFFieldCopy
 {
 public:
@@ -583,7 +560,7 @@ int main( int argc, char **argv )
                   << AfterFunction( SharedFunctor< DragForceEvaluator >( forceEval ), "drag force evaluation" );
 
    // resetting force
-   timeloop.addFuncAfterTimeStep( ResetForce( blocks, bodyStorageID ), "reset force on sphere");
+   timeloop.addFuncAfterTimeStep( pe_coupling::ForceTorqueOnBodiesResetter( blocks, bodyStorageID ), "reset force on sphere");
 
    timeloop.addFuncAfterTimeStep( RemainingTimeLogger( timeloop.getNrOfTimeSteps() ), "Remaining Time Logger" );
 

tests/pe_coupling/momentum_exchange_method/DragForceSphereMEMRefinement.cpp

@@ -61,6 +61,7 @@
 
 #include "pe_coupling/mapping/all.h"
 #include "pe_coupling/momentum_exchange_method/all.h"
+#include "pe_coupling/utility/all.h"
 
 #include "vtk/all.h"
 #include "field/vtk/all.h"
@@ -405,30 +406,6 @@ class ForceEval
 };
 
 
-class ResetForce
-{
-   public:
-      ResetForce( const shared_ptr< StructuredBlockStorage > & blocks,
-                  const BlockDataID & bodyStorageID )
-      : blocks_( blocks ), bodyStorageID_( bodyStorageID )
-      { }
-
-      void operator()()
-      {
-         for( auto blockIt = blocks_->begin(); blockIt != blocks_->end(); ++blockIt )
-         {
-            for( auto bodyIt = pe::BodyIterator::begin( *blockIt, bodyStorageID_); bodyIt != pe::BodyIterator::end(); ++bodyIt )
-            {
-                bodyIt->resetForceAndTorque();
-            }
-         }
-      }
-private:
-      shared_ptr< StructuredBlockStorage > blocks_;
-      const BlockDataID bodyStorageID_;
-};
-
-
 //////////
 // MAIN //
 //////////
@@ -596,7 +573,7 @@ int main( int argc, char **argv )
    timeloop.addFuncAfterTimeStep( makeSharedFunctor( forceEval ), "drag force evaluation" );
 
    // resetting force
-   timeloop.addFuncAfterTimeStep( ResetForce( blocks, bodyStorageID ), "reset force on sphere");
+   timeloop.addFuncAfterTimeStep( pe_coupling::ForceTorqueOnBodiesResetter( blocks, bodyStorageID ), "reset force on sphere");
 
    if( vtkIO )
    {

tests/pe_coupling/momentum_exchange_method/TorqueSphereMEM.cpp

@@ -58,6 +58,7 @@
 
 #include "pe_coupling/mapping/all.h"
 #include "pe_coupling/momentum_exchange_method/all.h"
+#include "pe_coupling/utility/all.h"
 
 #include <vector>
 #include <iomanip>
@@ -276,30 +277,6 @@ class TorqueEval
 
 };
 
-
-class ResetForce
-{
-   public:
-      ResetForce( const shared_ptr< StructuredBlockStorage > & blocks,
-                  const BlockDataID & bodyStorageID )
-      : blocks_( blocks ), bodyStorageID_( bodyStorageID )
-      { }
-
-      void operator()()
-      {
-         for( auto blockIt = blocks_->begin(); blockIt != blocks_->end(); ++blockIt )
-         {
-            for( auto bodyIt = pe::BodyIterator::begin( *blockIt, bodyStorageID_); bodyIt != pe::BodyIterator::end(); ++bodyIt )
-            {
-                bodyIt->resetForceAndTorque();
-            }
-         }
-      }
-private:
-      shared_ptr< StructuredBlockStorage > blocks_;
-      const BlockDataID bodyStorageID_;
-};
-
 class PDFFieldCopy
 {
 public:
@@ -539,7 +516,7 @@ int main( int argc, char **argv )
                      << AfterFunction( SharedFunctor< TorqueEval >( torqueEval ), "torque evaluation" );
    }
    // resetting force
-   timeloop.addFuncAfterTimeStep( ResetForce( blocks, bodyStorageID ), "reset force on sphere");
+   timeloop.addFuncAfterTimeStep( pe_coupling::ForceTorqueOnBodiesResetter( blocks, bodyStorageID ), "reset force on sphere");
 
    timeloop.addFuncAfterTimeStep( RemainingTimeLogger( timeloop.getNrOfTimeSteps() ), "Remaining Time Logger" );