Added distributors and test case to field module

src/field/all.h

@@ -45,6 +45,7 @@
 #include "allocation/all.h"
 #include "blockforest/all.h"
 #include "communication/all.h"
+#include "distributors/all.h"
 #include "interpolators/all.h"
 #include "iterators/all.h"
 #include "refinement/all.h"

src/field/distributors/DistributorCreators.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 DistributorCreators.h
+//! \ingroup field
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "core/DataTypes.h"
+#include "core/Set.h"
+
+#include "blockforest/BlockDataHandling.h"
+
+#include "domain_decomposition/StructuredBlockStorage.h"
+
+namespace walberla {
+namespace field {
+
+//**********************************************************************************************************************
+/*! DistributorCreators
+ *
+ * \ingroup field
+ *
+ * Distributor_T: A distributor that has a constructor
+ *  ( const shared_ptr<StructuredBlockStorage> & blockStorage, const IBlock & block, const BaseField_T & baseField,
+ *    const FlagField_T & flagField, const flag_t & evaluationMask )
+ * and distribution functions:
+ *  template< typename ForwardIterator_T > inline void distribute( const Vector3<real_t> & position, ForwardIterator_T distributeValueBegin )
+ *  template< typename ForwardIterator_T > inline void distribute( const real_t & x, const real_t & y, const real_t & z, ForwardIterator_T distributeValueBegin )
+ *
+ * See NearestNeighborDistributor for an example implementation.
+ *
+ * A distributor is aware of the flag field (FlagField_T) and distributes values only to cells flagged by a given mask.
+ *
+ * Distributors are used to spread a given value to the corresponding destination field.
+ * E.g. if a certain force has to be applied at some specific position onto the fluid, a distributor can be used
+ * to do so by distributing this force value (and conservation fo this force value is ensured) onto the force field.
+ *
+ */
+//**********************************************************************************************************************
+template< typename Distributor_T, typename FlagField_T >
+class DistributorHandling : public blockforest::AlwaysInitializeBlockDataHandling< Distributor_T >
+{
+public:
+
+   DistributorHandling( const shared_ptr<StructuredBlockStorage> & blockStorage,
+                        const BlockDataID & distributionDestinationFieldID,
+                        const ConstBlockDataID & flagFieldID,
+                        const Set< FlagUID > & cellsToEvaluate ) :
+   blockStorage_( blockStorage ), distributionDestinationFieldID_( distributionDestinationFieldID ), flagFieldID_( flagFieldID ), cellsToEvaluate_( cellsToEvaluate )
+   {}
+
+   Distributor_T * initialize( IBlock * const block )
+   {
+      typedef typename Distributor_T::BaseField_T DistributionDestinationField_T;
+      typedef typename FlagField_T::flag_t flag_t;
+      DistributionDestinationField_T * distributionDestinationField = block->getData< DistributionDestinationField_T >( distributionDestinationFieldID_ );
+      const FlagField_T * flagField = block->getData< FlagField_T >( flagFieldID_ );
+
+      WALBERLA_ASSERT_NOT_NULLPTR( distributionDestinationField );
+      WALBERLA_ASSERT_NOT_NULLPTR( flagField );
+
+      const flag_t evaluationMask = flagField->getMask( cellsToEvaluate_ );
+
+      return new Distributor_T( blockStorage_, *block, *distributionDestinationField, *flagField, evaluationMask );
+   }
+
+private:
+
+   shared_ptr<StructuredBlockStorage> blockStorage_;
+   BlockDataID distributionDestinationFieldID_;
+   ConstBlockDataID flagFieldID_;
+   Set< FlagUID > cellsToEvaluate_;
+   
+}; // class DistributorHandling
+
+
+
+template< typename Distributor_T, typename FlagField_T  >
+inline BlockDataID addDistributor( const shared_ptr< StructuredBlockStorage > & blocks,
+                                   const BlockDataID & distributionDestinationFieldID,
+                                   const ConstBlockDataID & flagFieldID,
+                                   const Set< FlagUID > & cellsToEvaluate,
+                                   const std::string & identifier = std::string(),
+                                   const Set<SUID> & requiredSelectors     = Set<SUID>::emptySet(),
+                                   const Set<SUID> & incompatibleSelectors = Set<SUID>::emptySet() )
+{
+   return blocks->addBlockData( make_shared< DistributorHandling< Distributor_T, FlagField_T > >( blocks, distributionDestinationFieldID, flagFieldID, cellsToEvaluate ), identifier, requiredSelectors, incompatibleSelectors );
+}
+
+
+} // namespace field
+} // namespace walberla

src/field/distributors/KernelDistributor.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 KernelDistributor.h
+//! \ingroup field
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "core/debug/Debug.h"
+#include "core/math/Vector3.h"
+
+#include "domain_decomposition/StructuredBlockStorage.h"
+
+#include "field/interpolators/KernelFieldInterpolator.h"
+#include "field/GhostLayerField.h"
+
+#include <vector>
+
+namespace walberla {
+namespace field {
+
+/*! Distributor for walberla::field::GhostLayerField with kernel strategy
+ *
+ * \ingroup field
+ *
+ * This distributor uses a smoothed dirac kernel function to distribute values to the field at the given position.
+ * The applied weights are given in the namespace field::kernelweights.
+ * Needs the full neighborhood of the containing cell, i.e. 27 cells.
+ * Never construct this distributor directly, but use the functionality from DistributorCreators.h instead.
+ */
+template< typename Field_T, typename FlagField_T >
+class KernelDistributor
+{
+public:
+
+   static const uint_t F_SIZE = Field_T::F_SIZE;
+
+   typedef Field_T                                BaseField_T;
+   typedef typename FlagField_T::flag_t           flag_t;
+   typedef KernelDistributor<Field_T,FlagField_T> OwnType;
+
+   KernelDistributor( const shared_ptr<StructuredBlockStorage> & blockStorage, const IBlock & block,
+                      BaseField_T & baseField, const FlagField_T & flagField,
+                      const flag_t & evaluationMask )
+   : blockStorage_( blockStorage ), block_( block ), baseField_( baseField ), flagField_( flagField ), evaluationMask_( evaluationMask )
+   {
+      WALBERLA_ASSERT(baseField.nrOfGhostLayers() > uint_t(0), "field for kernel distribution needs at least one ghost layer");
+   }
+
+   inline bool operator==( const OwnType & other ){ return baseField_ == other.baseField_; }
+
+   template< typename ForwardIterator_T >
+   inline void distribute( const Vector3<real_t> & position, ForwardIterator_T distributeValueBegin )
+   {
+      distribute( position[0], position[1], position[2], distributeValueBegin );
+   }
+
+   template< typename ForwardIterator_T >
+   inline void distribute( const real_t & x, const real_t & y, const real_t & z, ForwardIterator_T distributeValueBegin )
+   {
+      WALBERLA_ASSERT(block_.getAABB().contains(x,y,z),
+                      "Distribution position <" << x << ", " << y << ", " << z << "> is not contained inside the block of this distributor with AABB " << block_.getAABB() << " !");
+
+      Cell centerCell = blockStorage_->getBlockLocalCell( block_, x, y, z );
+
+      const real_t dx = blockStorage_->dx( blockStorage_->getLevel( block_ ) );
+      const real_t dy = blockStorage_->dy( blockStorage_->getLevel( block_ ) );
+      const real_t dz = blockStorage_->dz( blockStorage_->getLevel( block_ ) );
+      
+      const uint_t neighborhoodSize = cell_idx_t(1);
+
+      CellInterval cellNeighborhood( centerCell[0] - cell_idx_c(neighborhoodSize), centerCell[1] - cell_idx_c(neighborhoodSize), centerCell[2] - cell_idx_c(neighborhoodSize),
+                                     centerCell[0] + cell_idx_c(neighborhoodSize), centerCell[1] + cell_idx_c(neighborhoodSize), centerCell[2] + cell_idx_c(neighborhoodSize) );
+
+      const uint_t kernelSizeOneDirection = uint_t(2) * neighborhoodSize + uint_t(1);
+      std::vector<real_t> weights( kernelSizeOneDirection * kernelSizeOneDirection * kernelSizeOneDirection, real_t(0) );
+
+      uint_t counter = uint_t(0);
+      real_t sumOfWeights = real_t(0);
+      real_t sumOfWeightsUnavailable = real_t(0);
+
+      // get distribution weights and count available cells in surrounding cells
+      for( auto cellIt = cellNeighborhood.begin(); cellIt != cellNeighborhood.end(); ++cellIt )
+      {
+         Vector3<real_t> curCellCenter = blockStorage_->getBlockLocalCellCenter( block_, *cellIt );
+         if( flagField_.isPartOfMaskSet( *cellIt, evaluationMask_ ) )
+         {
+            weights[counter] = kernelweights::kernelWeightFunction( x, y, z, curCellCenter[0], curCellCenter[1], curCellCenter[2], dx, dy, dz );
+            sumOfWeights += weights[counter];
+         }
+         else
+         {
+            weights[counter] = real_t(0);
+            sumOfWeightsUnavailable += kernelweights::kernelWeightFunction( x, y, z, curCellCenter[0], curCellCenter[1], curCellCenter[2], dx, dy, dz );
+         }
+         ++counter;
+      }
+
+      // check if at least one cell was available, to prevent division by 0
+      if( sumOfWeights <= real_t(0) )
+         return;
+
+      // scale domain weights if some non-domain cells are in neighborhood
+      const real_t scalingFactor = real_t(1) + sumOfWeightsUnavailable / sumOfWeights;
+
+      // distribute the values to the neighboring domain cells with the corresponding (scaled) weighting
+      counter = uint_t(0);
+      for( auto cellIt = cellNeighborhood.begin(); cellIt != cellNeighborhood.end(); ++cellIt )
+      {
+         if ( weights[counter] > real_t(0) )
+         {
+            addWeightedCellValue( distributeValueBegin, *cellIt, scalingFactor * weights[counter] );
+         }
+         ++counter;
+      }
+   }
+
+private:
+
+   template< typename ForwardIterator_T >
+   void addWeightedCellValue( ForwardIterator_T distributeValueBegin, const Cell & curCell, const real_t & weighting )
+   {
+      for( uint_t f = uint_t(0); f < F_SIZE; ++f )
+      {
+         baseField_( curCell, f) += weighting * (*distributeValueBegin);
+         ++distributeValueBegin;
+      }
+   }
+
+   shared_ptr<StructuredBlockStorage> blockStorage_;
+   const IBlock & block_;
+   BaseField_T & baseField_;
+   const FlagField_T & flagField_;
+   flag_t evaluationMask_;
+
+};
+
+
+} // namespace field
+} // namespace walberla

src/field/distributors/NearestNeighborDistributor.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 NearestNeighborDistributor.h
+//! \ingroup field
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "core/math/Vector3.h"
+
+#include "domain_decomposition/StructuredBlockStorage.h"
+
+#include "field/GhostLayerField.h"
+
+#include <numeric>
+#include <vector>
+
+namespace walberla {
+namespace field {
+
+/*! Distributor for walberla::field::Field with nearest neighbor strategy
+ *
+ * \ingroup field
+ *
+ * This distributor distributes the given value to the nearest cell, flagged in the evaluation mask, of the given position.
+ * This is usually the cell that contains the distribution position.
+ * If this cell is a cell that is not within the evaluation mask, the direct neighborhood (8 cells) will be searched for an available cell.
+ * Never construct this distributor directly, but use the functionality from DistributorCreators.h instead.
+ */
+template< typename Field_T, typename FlagField_T >
+class NearestNeighborDistributor
+{
+public:
+
+   static const uint_t F_SIZE = Field_T::F_SIZE;
+
+   typedef Field_T                                         BaseField_T;
+   typedef typename FlagField_T::flag_t                    flag_t;
+   typedef NearestNeighborDistributor<Field_T,FlagField_T> OwnType;
+
+   NearestNeighborDistributor( const shared_ptr<StructuredBlockStorage> & blockStorage, const IBlock & block,
+                               BaseField_T & baseField, const FlagField_T & flagField,
+                               const flag_t & evaluationMask )
+   : blockStorage_( blockStorage ), block_( block ), baseField_( baseField ), flagField_( flagField ), evaluationMask_( evaluationMask )
+   {}
+
+   inline bool operator==( const OwnType & other ){ return baseField_ == other.baseField_; }
+
+   template< typename ForwardIterator_T >
+   inline void distribute( const Vector3<real_t> & position, ForwardIterator_T distributeValueBegin )
+   {
+      distribute( position[0], position[1], position[2], distributeValueBegin );
+   }
+
+   template< typename ForwardIterator_T >
+   inline void distribute( const real_t & x, const real_t & y, const real_t & z, ForwardIterator_T distributeValueBegin )
+   {
+      WALBERLA_ASSERT(block_.getAABB().contains(x,y,z),
+                      "Distribution position <" << x << ", " << y << ", " << z << "> is not contained inside the block of this distributor with AABB " << block_.getAABB() << " !");
+
+      Cell nearestCell = blockStorage_->getBlockLocalCell( block_, x, y, z );
+
+      if( flagField_.isPartOfMaskSet( nearestCell, evaluationMask_ ) )
+      {
+         for( uint_t f = uint_t(0); f < F_SIZE; ++f )
+         {
+            baseField_( nearestCell, f) += *distributeValueBegin;
+            ++distributeValueBegin;
+         }
+         return;
+      }
+      else
+      {
+         // look for available cell in direct neighborhood
+
+         CellInterval fieldXYZSize = baseField_.xyzSize();
+
+         Vector3<real_t> nearestCellCenter = blockStorage_->getBlockLocalCellCenter( block_, nearestCell );
+         const cell_idx_t xNeighbor = cell_idx_c( floor( x - nearestCellCenter[0] ) );
+         const cell_idx_t yNeighbor = cell_idx_c( floor( y - nearestCellCenter[1] ) );
+         const cell_idx_t zNeighbor = cell_idx_c( floor( z - nearestCellCenter[2] ) );
+
+         const cell_idx_t xMin = nearestCell.x() + xNeighbor;
+         const cell_idx_t yMin = nearestCell.y() + yNeighbor;
+         const cell_idx_t zMin = nearestCell.z() + zNeighbor;
+
+         for( cell_idx_t zC = zMin; zC <= zMin + cell_idx_t(1); ++zC)
+         {
+            for( cell_idx_t yC = yMin; yC <= yMin + cell_idx_t(1); ++yC)
+            {
+               for( cell_idx_t xC = xMin; xC <= xMin + cell_idx_t(1); ++xC)
+               {
+                  Cell curCell(xC,yC,zC);
+                  if( fieldXYZSize.contains(curCell) )
+                  {
+                     if (flagField_.isPartOfMaskSet(curCell, evaluationMask_))
+                     {
+                        for (uint_t f = uint_t(0); f < F_SIZE; ++f)
+                        {
+                           baseField_(curCell, f) += *distributeValueBegin;
+                           ++distributeValueBegin;
+                        }
+                        return;
+                     }
+                  }
+               }
+            }
+         }
+      }
+   }
+
+private:
+
+   shared_ptr<StructuredBlockStorage> blockStorage_;
+   const IBlock & block_;
+   BaseField_T & baseField_;
+   const FlagField_T & flagField_;
+   flag_t evaluationMask_;
+
+};
+
+} // namespace field
+} // namespace walberla

src/field/distributors/all.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 all.h
+//! \ingroup field
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "DistributorCreators.h"
+#include "KernelDistributor.h"
+#include "NearestNeighborDistributor.h"

tests/field/CMakeLists.txt

@@ -10,6 +10,9 @@ waLBerla_execute_test( NAME AccuracyEvaluationTest4 COMMAND $<TARGET_FILE:Accura
 waLBerla_compile_test( FILES communication/FieldPackInfoTest.cpp DEPENDS blockforest )
 waLBerla_execute_test( NAME  FieldPackInfoTest )
 
+waLBerla_compile_test( FILES distributors/DistributionTest.cpp)
+waLBerla_execute_test( NAME DistributionTest )
+
 waLBerla_compile_test( FILES FieldTest.cpp )
 waLBerla_execute_test( NAME FieldTest )