NearestNeighborFieldInterpolator.h 5.86 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
//======================================================================================================================
//
//  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 NearestNeighborFieldInterpolator.h
//! \ingroup field
//! \author Christoph Rettinger <christoph.rettinger@fau.de>
//
//======================================================================================================================

#pragma once

#include "core/debug/Debug.h"

#include "domain_decomposition/StructuredBlockStorage.h"

#include "field/FlagField.h"

namespace walberla {
namespace field {

/*! Interpolator for walberla::field::Field with nearest neighbor strategy
 *
 * \ingroup field
 *
 * This interpolator obtains the value from the nearest cell, flagged in the evaluation mask, of the interpolation position.
 * This is usually the cell that contains the interpolation 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 interpolator directly, but use the functionality from FieldInterpolatorCreators.h instead.
 */
template< typename Field_T, typename FlagField_T >
class NearestNeighborFieldInterpolator
{
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 NearestNeighborFieldInterpolator<Field_T,FlagField_T> OwnType;

53
   NearestNeighborFieldInterpolator( const weak_ptr<StructuredBlockStorage> & blockStorage, const IBlock & block,
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
                                     const 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 get( const Vector3<real_t> & position, ForwardIterator_T interpolationResultBegin )
   {
      get( position[0], position[1], position[2], interpolationResultBegin );
   }

   template< typename ForwardIterator_T >
   inline void get( const real_t & x, const real_t & y, const real_t & z, ForwardIterator_T interpolationResultBegin )
   {

      WALBERLA_ASSERT(block_.getAABB().contains(x,y,z),
                      "Interpolation position <" << x << ", " << y << ", " << z << "> is not contained inside the block of this interpolator with AABB " << block_.getAABB() << " !");

75
76
77
78
      WALBERLA_CHECK( !blockStorage_.expired() );
      auto blockStorage = blockStorage_.lock();

      Cell nearestCell = blockStorage->getBlockLocalCell( block_, x, y, z );
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94

      if( flagField_.isPartOfMaskSet( nearestCell, evaluationMask_ ) )
      {
         for( uint_t f = uint_t(0); f < F_SIZE; ++f )
         {
            WALBERLA_ASSERT( !math::isnan( baseField_( nearestCell, f) ), "NaN found in component " << f << " when interpolating from cell " << nearestCell );
            *interpolationResultBegin = baseField_( nearestCell, f);
            ++interpolationResultBegin;
         }
      }
      else
      {
         // look for available cell in direct neighborhood

         CellInterval fieldXYZSize = baseField_.xyzSize();

95
         Vector3<real_t> nearestCellCenter = blockStorage->getBlockLocalCellCenter( block_, nearestCell );
96
97
98
99
100
101
102
103
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
         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)
                        {
                           WALBERLA_ASSERT(!math::isnan(baseField_(curCell, f)),
                                           "NaN found in component " << f << " when interpolating from cell "
                                                                     << curCell);
                           *interpolationResultBegin = baseField_(curCell, f);
                           ++interpolationResultBegin;
                        }
                        break;
                     }
                  }
               }
            }
         }
      }
   }

private:

134
   weak_ptr<StructuredBlockStorage> blockStorage_;
135
136
137
138
139
140
141
142
   const IBlock & block_;
   const BaseField_T & baseField_;
   const FlagField_T & flagField_;
   flag_t evaluationMask_;
};

} // namespace field
} // namespace walberla