KernelDistributor.h 6.5 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
//======================================================================================================================
//
//  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;

57
   KernelDistributor( const weak_ptr<StructuredBlockStorage> & blockStorage, const IBlock & block,
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
                      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() << " !");

79
80
      WALBERLA_CHECK( !blockStorage_.expired() );
      auto blockStorage = blockStorage_.lock();
81

82
83
84
85
86
      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_ ) );
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
      
      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 )
      {
103
         Vector3<real_t> curCellCenter = blockStorage->getBlockLocalCellCenter( block_, *cellIt );
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
         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;
      }
   }

148
   weak_ptr<StructuredBlockStorage> blockStorage_;
149
150
151
152
153
154
155
156
157
158
   const IBlock & block_;
   BaseField_T & baseField_;
   const FlagField_T & flagField_;
   flag_t evaluationMask_;

};


} // namespace field
} // namespace walberla