diff --git a/src/field/interpolators/FieldInterpolatorCreators.h b/src/field/interpolators/FieldInterpolatorCreators.h
new file mode 100644
index 0000000000000000000000000000000000000000..9b7c4f656c7cf158856945c0a7f09741f9d8475b
--- /dev/null
+++ b/src/field/interpolators/FieldInterpolatorCreators.h
@@ -0,0 +1,107 @@
+//======================================================================================================================
+//
+// 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 FieldInterpolatorCreators.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 {
+
+//**********************************************************************************************************************
+/*! FieldInterpolatorCreators
+ *
+ * \ingroup field
+ *
+ * Interpolator_T: A field interpolator 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 getter functions:
+ * template< typename ForwardIterator_T > inline void get( const Vector3<real_t> & position, ForwardIterator_T interpolationResultBegin )
+ * template< typename ForwardIterator_T > inline void get( const real_t & x, const real_t & y, const real_t & z, ForwardIterator_T interpolationResultBegin )
+ *
+ * See TrilinearFieldInterpolator for an example implementation.
+ *
+ * A field interpolator is aware of the flag field (FlagField_T) and uses only values that contain flags from a given mask.
+ *
+ * Field Interpolators can be used to sample the underlying field at certain positions.
+ * E.g. the fluid velocity can be interpolated to a desired global position from a velocity field.
+ *
+ */
+//**********************************************************************************************************************
+template< typename Interpolator_T, typename FlagField_T >
+class InterpolatorHandling : public blockforest::AlwaysInitializeBlockDataHandling< Interpolator_T >
+{
+public:
+
+ InterpolatorHandling( const shared_ptr<StructuredBlockStorage> & blockStorage,
+ const ConstBlockDataID & interpolatedFieldID,
+ const ConstBlockDataID & flagFieldID,
+ const Set< FlagUID > & cellsToEvaluate ) :
+ blockStorage_( blockStorage ), interpolatedFieldID_( interpolatedFieldID ), flagFieldID_( flagFieldID ), cellsToEvaluate_( cellsToEvaluate )
+ {}
+
+ Interpolator_T * initialize( IBlock * const block )
+ {
+ typedef typename Interpolator_T::BaseField_T InterpolatedField_T;
+ typedef typename FlagField_T::flag_t flag_t;
+ const InterpolatedField_T * interpolatedField = block->getData< InterpolatedField_T >( interpolatedFieldID_ );
+ const FlagField_T * flagField = block->getData< FlagField_T >( flagFieldID_ );
+
+ WALBERLA_ASSERT_NOT_NULLPTR( interpolatedField );
+ WALBERLA_ASSERT_NOT_NULLPTR( flagField );
+
+ const flag_t evaluationMask = flagField->getMask( cellsToEvaluate_ );
+
+ return new Interpolator_T( blockStorage_, *block, *interpolatedField, *flagField, evaluationMask );
+ }
+
+private:
+
+ shared_ptr<StructuredBlockStorage> blockStorage_;
+ ConstBlockDataID interpolatedFieldID_;
+ ConstBlockDataID flagFieldID_;
+ Set< FlagUID > cellsToEvaluate_;
+
+}; // class InterpolatorHandling
+
+
+
+template< typename Interpolator_T, typename FlagField_T >
+inline BlockDataID addFieldInterpolator( const shared_ptr< StructuredBlockStorage > & blocks,
+ const ConstBlockDataID & interpolatedFieldID,
+ 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< InterpolatorHandling< Interpolator_T, FlagField_T > >( blocks, interpolatedFieldID, flagFieldID, cellsToEvaluate ), identifier, requiredSelectors, incompatibleSelectors );
+}
+
+
+} // namespace field
+} // namespace walberla
diff --git a/src/field/interpolators/KernelFieldInterpolator.h b/src/field/interpolators/KernelFieldInterpolator.h
new file mode 100644
index 0000000000000000000000000000000000000000..8e22b16609890258147fdea20d35b946192f674d
--- /dev/null
+++ b/src/field/interpolators/KernelFieldInterpolator.h
@@ -0,0 +1,269 @@
+//======================================================================================================================
+//
+// 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 KernelFieldInterpolator.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/FlagField.h"
+
+#include "stencil/D3Q27.h"
+
+#include <numeric>
+
+namespace walberla {
+namespace field {
+
+namespace kernelweights
+{
+
+// corresponds to the smoothed dirac delta function from Roma et al - An Adaptive Version of the Immersed Boundary Method
+// f(r) != 0 for abs(r) < 1.5 -> requires three neighborhood cells
+real_t smoothedDeltaFunction( const real_t & r )
+{
+ real_t rAbs = std::fabs(r);
+ if( rAbs <= real_t(0.5) )
+ {
+ return (real_t(1) + std::sqrt( - real_t(3) * r * r + real_t(1) ) ) * (real_t(1) / real_t(3));
+ } else if ( rAbs < real_t(1.5) )
+ {
+ return (real_t(5) - real_t(3) * rAbs - std::sqrt( - real_t(3) * ( real_t(1) - rAbs ) * ( real_t(1) - rAbs ) + real_t(1) ) ) * ( real_t(1) / real_t(6) );
+ } else
+ {
+ return real_t(0);
+ }
+
+}
+
+// X: Lagrangian position, x: Eulerian position (usually cell center), global coordinates
+// dx, dy, dz: mesh spacing
+real_t kernelWeightFunction( const real_t & X, const real_t & Y, const real_t & Z,
+ const real_t & x, const real_t & y, const real_t & z,
+ const real_t & dx = real_t(1), const real_t & dy = real_t(1), const real_t & dz = real_t(1) )
+{
+ return smoothedDeltaFunction( ( X - x ) / dx ) * smoothedDeltaFunction( ( Y - y ) / dy ) * smoothedDeltaFunction( ( Z - z ) / dz );
+}
+
+// X: Lagrangian position, x: Eulerian position (usually cell center), global coordinates
+// dx, dy, dz: mesh spacing
+real_t kernelWeightFunction( const Vector3<real_t> & X, const Vector3<real_t> & x,
+ const real_t & dx = real_t(1), const real_t & dy = real_t(1), const real_t & dz = real_t(1) )
+{
+ return smoothedDeltaFunction( ( X[0] - x[0] ) / dx ) * smoothedDeltaFunction( ( X[1] - x[1] ) / dy ) * smoothedDeltaFunction( ( X[2] - x[2] ) / dz );
+}
+
+} // namespace kernelweights
+
+
+/*! Interpolator for walberla::field::GhostLayerField with kernel strategy
+ *
+ * \ingroup field
+ *
+ * This interpolator uses a smoothed dirac kernel function to interpolate values.
+ * 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 interpolator directly, but use the functionality from the FieldInterpolatorCreator.h instead.
+ */
+template< typename Field_T, typename FlagField_T >
+class KernelFieldInterpolator
+{
+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 KernelFieldInterpolator<Field_T,FlagField_T> OwnType;
+
+ KernelFieldInterpolator( const shared_ptr<StructuredBlockStorage> & blockStorage, const IBlock & block,
+ const 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 interpolation needs at least one ghost layer");
+ }
+
+
+ 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() << " !");
+
+ 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 = uint_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);
+ // store the calculated weighting factors of the available cells, i.e. cells flagged by the evaluation mask
+ std::vector<real_t> weights( kernelSizeOneDirection * kernelSizeOneDirection * kernelSizeOneDirection, real_t(0) );
+
+ // store the calculated weighting factors of the UNavailable cells, i.e. cells not included in the evaluation mask
+ std::vector<real_t> weightsUnavailable( kernelSizeOneDirection * kernelSizeOneDirection * kernelSizeOneDirection, real_t(0) );
+
+ cell_idx_t cellIdx0x = cellNeighborhood.xMin();
+ cell_idx_t cellIdx0y = cellNeighborhood.yMin();
+ cell_idx_t cellIdx0z = cellNeighborhood.zMin();
+ uint_t nx = kernelSizeOneDirection;
+ uint_t nxy = kernelSizeOneDirection * kernelSizeOneDirection;
+ Vector3<real_t> cellCenter0 = blockStorage_->getBlockLocalCellCenter( block_, Cell(cellIdx0x, cellIdx0y, cellIdx0z) ); // = cell in neighborhood with smallest x-, y-, and z-indices
+
+ // calculate kernel weights of all cells in neighborhood
+ for( uint_t k = uint_t(0); k < kernelSizeOneDirection; ++k)
+ {
+ for( uint_t j = uint_t(0); j < kernelSizeOneDirection; ++j)
+ {
+ for( uint_t i = uint_t(0); i < kernelSizeOneDirection; ++i)
+ {
+ Vector3<real_t> curCellCenter( cellCenter0[0] + real_c(i) * dx, cellCenter0[1] + real_c(j) * dy, cellCenter0[2] + real_c(k) * dz);
+ if( flagField_.isPartOfMaskSet( cellIdx0x + cell_idx_c(i), cellIdx0y + cell_idx_c(j), cellIdx0z + cell_idx_c(k), evaluationMask_ ) )
+ {
+ weights[k*nxy+j*nx+i] = kernelweights::kernelWeightFunction( x, y, z, curCellCenter[0], curCellCenter[1], curCellCenter[2], dx, dy, dz );
+ }
+ else
+ {
+ weightsUnavailable[k*nxy+j*nx+i] = kernelweights::kernelWeightFunction( x, y, z, curCellCenter[0], curCellCenter[1], curCellCenter[2], dx, dy, dz );
+ }
+ }
+ }
+ }
+
+/*
+ // zero entries in weights array (and non-zero entries in weightsUnavailable) indicate non-fluid nodes
+ // weight correction to incorporate extrapolation for kernel interpolation, center element can not be redistributed
+ // without this part, the kernel interpolator can not extrapolate values
+ // this is based on experimental findings by the author and general correctness is not ensured, thus it is commented out
+ for( auto stenIt = stencil::D3Q27::beginNoCenter(); stenIt != stencil::D3Q27::end(); ++stenIt )
+ {
+ int cx = stenIt.cx();
+ int cy = stenIt.cy();
+ int cz = stenIt.cz();
+ int i = cx + 1;
+ int j = cy + 1;
+ int k = cz + 1;
+
+ if( weightsUnavailable[k*nxy+j*nx+i] > real_t(0) )
+ {
+ // check if we can distribute the non-fluid weight to nearby fluid weights that lie in one line inside the neighborhood
+ // it should generally not matter in which direction it is distributed
+ // check x direction
+ if( weights[k*nxy+j*nx+i-cx] > real_t(0) && weights[k*nxy+j*nx+i-2*cx] > real_t(0) )
+ {
+ weights[k*nxy+j*nx+i-cx] += real_t(2)*weightsUnavailable[k*nxy+j*nx+i];
+ weights[k*nxy+j*nx+i-2*cx] -= weightsUnavailable[k*nxy+j*nx+i];
+ weightsUnavailable[k*nxy+j*nx+i] = real_t(0);
+ continue;
+ }
+ // check y direction
+ if( weights[k*nxy+(j-cy)*nx+i] > real_t(0) && weights[k*nxy+(j-2*cy)*nx+i] > real_t(0) )
+ {
+ weights[k*nxy+(j-cy)*nx+i] += real_t(2)*weightsUnavailable[k*nxy+j*nx+i];
+ weights[k*nxy+(j-2*cy)*nx+i] -= weightsUnavailable[k*nxy+j*nx+i];
+ weightsUnavailable[k*nxy+j*nx+i] = real_t(0);
+ continue;
+ }
+ // check z direction
+ if( weights[(k-cz)*nxy+j*nx+i] > real_t(0) && weights[(k-2*cz)*nxy+j*nx+i] > real_t(0) )
+ {
+ weights[(k-cz)*nxy+j*nx+i] += real_t(2)*weightsUnavailable[k*nxy+j*nx+i];
+ weights[(k-2*cz)*nxy+j*nx+i] -= weightsUnavailable[k*nxy+j*nx+i];
+ weightsUnavailable[k*nxy+j*nx+i] = real_t(0);
+ continue;
+ }
+ }
+ }
+*/
+ // scale available weights by the total amount of unavailable weights such that afterwards sum over all weights is 1
+ const real_t sumOfWeightsUnavailable = std::accumulate(weightsUnavailable.begin(), weightsUnavailable.end(), real_t(0) );
+ const real_t sumOfWeights = std::accumulate(weights.begin(), weights.end(), real_t(0) );
+
+ // check if at least one cell was available, to prevent division by 0
+ if( sumOfWeights <= real_t(0) )
+ return;
+
+ const real_t scalingFactor = real_t(1) + sumOfWeightsUnavailable / sumOfWeights;
+
+ for( auto weightIt = weights.begin(); weightIt != weights.end(); ++weightIt )
+ {
+ *weightIt *= scalingFactor;
+ }
+
+ // interpolate value to interpolation position using the previously calculated weights
+ // values to not be included have a weight of 0
+ for( uint_t k = uint_t(0); k < kernelSizeOneDirection; ++k)
+ {
+ for( uint_t j = uint_t(0); j < kernelSizeOneDirection; ++j)
+ {
+ for( uint_t i = uint_t(0); i < kernelSizeOneDirection; ++i)
+ {
+ if ( weights[k*nxy+j*nx+i] > real_t(0) )
+ {
+ Cell curCell( cellIdx0x + cell_idx_c(i), cellIdx0y + cell_idx_c(j), cellIdx0z + cell_idx_c(k) );
+ addWeightedCellValue( interpolationResultBegin, curCell, weights[k*nxy+j*nx+i] );
+ }
+ }
+ }
+ }
+ }
+
+private:
+
+ template< typename ForwardIterator_T >
+ void addWeightedCellValue( ForwardIterator_T interpolationResultBegin, const Cell & curCell, const real_t & weighting )
+ {
+ 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 += weighting * baseField_( curCell, f);
+ ++interpolationResultBegin;
+ }
+ }
+
+ shared_ptr<StructuredBlockStorage> blockStorage_;
+ const IBlock & block_;
+ const BaseField_T & baseField_;
+ const FlagField_T & flagField_;
+ flag_t evaluationMask_;
+
+};
+
+
+} // namespace field
+} // namespace walberla
diff --git a/src/field/interpolators/NearestNeighborFieldInterpolator.h b/src/field/interpolators/NearestNeighborFieldInterpolator.h
new file mode 100644
index 0000000000000000000000000000000000000000..677310feb3fcd00c6da69fb8d6c1b8df62a97086
--- /dev/null
+++ b/src/field/interpolators/NearestNeighborFieldInterpolator.h
@@ -0,0 +1,139 @@
+//======================================================================================================================
+//
+// 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;
+
+ NearestNeighborFieldInterpolator( const shared_ptr<StructuredBlockStorage> & blockStorage, const IBlock & block,
+ 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() << " !");
+
+ Cell nearestCell = blockStorage_->getBlockLocalCell( block_, x, y, z );
+
+ 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();
+
+ 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)
+ {
+ WALBERLA_ASSERT(!math::isnan(baseField_(curCell, f)),
+ "NaN found in component " << f << " when interpolating from cell "
+ << curCell);
+ *interpolationResultBegin = baseField_(curCell, f);
+ ++interpolationResultBegin;
+ }
+ break;
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+
+private:
+
+ shared_ptr<StructuredBlockStorage> blockStorage_;
+ const IBlock & block_;
+ const BaseField_T & baseField_;
+ const FlagField_T & flagField_;
+ flag_t evaluationMask_;
+};
+
+} // namespace field
+} // namespace walberla
diff --git a/src/field/interpolators/TrilinearFieldInterpolator.h b/src/field/interpolators/TrilinearFieldInterpolator.h
new file mode 100644
index 0000000000000000000000000000000000000000..4a878381cb9622e55fae5d8cfee5105bddf40def
--- /dev/null
+++ b/src/field/interpolators/TrilinearFieldInterpolator.h
@@ -0,0 +1,188 @@
+//======================================================================================================================
+//
+// 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 TrilinearFieldInterpolator.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::GhostLayerField with trilinear strategy
+ *
+ * \ingroup field
+ *
+ * This interpolator uses trilinear interpolation to obtain the interpolated value from the interpolation position.
+ * If at least one of the cells, that are required for trilinear interpolation, i.e. the 8 closest cells,
+ * is not contained in the evaluation mask, the nearest-neighbor interpolation (see NearestNeighborFieldInterpolator.h)
+ * will be used instead.
+ * Never construct this interpolator directly, but use the functionality from FieldInterpolatorCreators.h instead.
+ */
+template< typename Field_T, typename FlagField_T >
+class TrilinearFieldInterpolator
+{
+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 TrilinearFieldInterpolator<Field_T,FlagField_T> OwnType;
+
+ TrilinearFieldInterpolator( const shared_ptr<StructuredBlockStorage> & blockStorage, const IBlock & block,
+ const BaseField_T & baseField, const FlagField_T & flagField,
+ const flag_t & evaluationMask )
+ : blockStorage_( blockStorage ), block_( block ), baseField_( baseField ), flagField_( flagField ), evaluationMask_( evaluationMask ),
+ nearestNeighborInterpolator_( blockStorage, block, baseField, flagField, evaluationMask )
+ {
+ WALBERLA_ASSERT(baseField.nrOfGhostLayers() > uint_t(0), "field for trilinear interpolation needs at least one ghost layer");
+ }
+
+
+ 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() << " !");
+
+ 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_ ) );
+
+ Cell containingCell = blockStorage_->getBlockLocalCell( block_, x, y, z );
+ Vector3<real_t> containingCellCenter = blockStorage_->getBlockLocalCellCenter( block_, containingCell );
+
+ const cell_idx_t xNeighbor1 = cell_idx_c( floor( x - containingCellCenter[0] ) );
+ const cell_idx_t xNeighbor2 = xNeighbor1 + cell_idx_t(1);
+ const cell_idx_t yNeighbor1 = cell_idx_c( floor( y - containingCellCenter[1] ) );
+ const cell_idx_t yNeighbor2 = yNeighbor1 + cell_idx_t(1);
+ const cell_idx_t zNeighbor1 = cell_idx_c( floor( z - containingCellCenter[2] ) );
+ const cell_idx_t zNeighbor2 = zNeighbor1 + cell_idx_t(1);
+
+ // define the 8 nearest cells required for the trilienar interpolation
+ // the cell 'ccc' is the one with the smallest x-, y-, and z-indices
+ Cell ccc( containingCell.x() + xNeighbor1, containingCell.y() + yNeighbor1, containingCell.z() + zNeighbor1 );
+ Cell hcc( containingCell.x() + xNeighbor2, containingCell.y() + yNeighbor1, containingCell.z() + zNeighbor1 );
+ Cell chc( containingCell.x() + xNeighbor1, containingCell.y() + yNeighbor2, containingCell.z() + zNeighbor1 );
+ Cell hhc( containingCell.x() + xNeighbor2, containingCell.y() + yNeighbor2, containingCell.z() + zNeighbor1 );
+ Cell cch( containingCell.x() + xNeighbor1, containingCell.y() + yNeighbor1, containingCell.z() + zNeighbor2 );
+ Cell hch( containingCell.x() + xNeighbor2, containingCell.y() + yNeighbor1, containingCell.z() + zNeighbor2 );
+ Cell chh( containingCell.x() + xNeighbor1, containingCell.y() + yNeighbor2, containingCell.z() + zNeighbor2 );
+ Cell hhh( containingCell.x() + xNeighbor2, containingCell.y() + yNeighbor2, containingCell.z() + zNeighbor2 );
+
+ // check if all neighboring cells are domain cells
+ if( flagField_.isPartOfMaskSet( ccc, evaluationMask_ ) &&
+ flagField_.isPartOfMaskSet( hcc, evaluationMask_ ) &&
+ flagField_.isPartOfMaskSet( chc, evaluationMask_ ) &&
+ flagField_.isPartOfMaskSet( hhc, evaluationMask_ ) &&
+ flagField_.isPartOfMaskSet( cch, evaluationMask_ ) &&
+ flagField_.isPartOfMaskSet( hch, evaluationMask_ ) &&
+ flagField_.isPartOfMaskSet( chh, evaluationMask_ ) &&
+ flagField_.isPartOfMaskSet( hhh, evaluationMask_ ) )
+ {
+ // trilinear interpolation can be applied
+
+ const real_t inv_totalVolume = real_t(1) / ( dx * dy * dz );
+ Vector3<real_t> cccCellCenter = blockStorage_->getBlockLocalCellCenter( block_, ccc );
+
+ // weighting = volume of opposing volume element / total volume
+ real_t weighting(0.0);
+ // x: (c)---dxc--(X)--------------dxh-------------(h)
+ const real_t dxc = x - cccCellCenter[0];
+ const real_t dxh = cccCellCenter[0] + dx - x;
+ const real_t dyc = y - cccCellCenter[1];
+ const real_t dyh = cccCellCenter[1] + dy - y;
+ const real_t dzc = z - cccCellCenter[2];
+ const real_t dzh = cccCellCenter[2] + dz - z;
+
+
+ weighting = dxh * dyh * dzh * inv_totalVolume;
+ addWeightedCellValue( interpolationResultBegin, ccc, weighting );
+
+ weighting = dxc * dyh * dzh * inv_totalVolume;
+ addWeightedCellValue( interpolationResultBegin, hcc, weighting );
+
+ weighting = dxh * dyc * dzh * inv_totalVolume;
+ addWeightedCellValue( interpolationResultBegin, chc, weighting );
+
+ weighting = dxc * dyc * dzh * inv_totalVolume;
+ addWeightedCellValue( interpolationResultBegin, hhc, weighting );
+
+ weighting = dxh * dyh * dzc * inv_totalVolume;
+ addWeightedCellValue( interpolationResultBegin, cch, weighting );
+
+ weighting = dxc * dyh * dzc * inv_totalVolume;
+ addWeightedCellValue( interpolationResultBegin, hch, weighting );
+
+ weighting = dxh * dyc * dzc * inv_totalVolume;
+ addWeightedCellValue( interpolationResultBegin, chh, weighting );
+
+ weighting = dxc * dyc * dzc * inv_totalVolume;
+ addWeightedCellValue( interpolationResultBegin, hhh, weighting );
+
+ }
+ else
+ {
+ // revert to nearest neighbor interpolation
+ nearestNeighborInterpolator_.get( x, y, z, interpolationResultBegin );
+ }
+ }
+
+private:
+
+ template< typename ForwardIterator_T >
+ void addWeightedCellValue( ForwardIterator_T interpolationResultBegin, const Cell & curCell, const real_t & weighting )
+ {
+ 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 += weighting * baseField_( curCell, f);
+ ++interpolationResultBegin;
+ }
+ }
+
+ shared_ptr<StructuredBlockStorage> blockStorage_;
+ const IBlock & block_;
+ const BaseField_T & baseField_;
+ const FlagField_T & flagField_;
+ flag_t evaluationMask_;
+
+ NearestNeighborFieldInterpolator<Field_T,FlagField_T> nearestNeighborInterpolator_;
+
+};
+
+
+} // namespace field
+} // namespace walberla
diff --git a/src/field/interpolators/all.h b/src/field/interpolators/all.h
index 0e7cfa80d337842b0320ea205f0dc109e63094a7..92055c83158120b94f671722f30d2e1d46b01ce4 100644
--- a/src/field/interpolators/all.h
+++ b/src/field/interpolators/all.h
@@ -24,3 +24,8 @@
#include "NearestNeighborInterpolator.h"
#include "TrilinearInterpolator.h"
+
+#include "FieldInterpolatorCreators.h"
+#include "KernelFieldInterpolator.h"
+#include "NearestNeighborFieldInterpolator.h"
+#include "TrilinearFieldInterpolator.h"
diff --git a/tests/field/CMakeLists.txt b/tests/field/CMakeLists.txt
index 4307644a0dcbbdd7d803f6a91d91b8e02e1c869f..2e490ea32b22d0266e8a7aac742be338758057be 100644
--- a/tests/field/CMakeLists.txt
+++ b/tests/field/CMakeLists.txt
@@ -25,6 +25,9 @@ waLBerla_execute_test( NAME FlagFieldTest )
waLBerla_compile_test( FILES interpolators/InterpolationTest.cpp)
waLBerla_execute_test( NAME InterpolationTest )
+waLBerla_compile_test( FILES interpolators/FieldInterpolationTest.cpp)
+waLBerla_execute_test( NAME FieldInterpolationTest )
+
waLBerla_compile_test( FILES adaptors/AdaptorTest.cpp DEPENDS blockforest gui lbm )
waLBerla_execute_test( NAME AdaptorTest )
diff --git a/tests/field/interpolators/FieldInterpolationTest.cpp b/tests/field/interpolators/FieldInterpolationTest.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..32adf2898cf57a2ccbc988d01d0f9de02c8ebfd9
--- /dev/null
+++ b/tests/field/interpolators/FieldInterpolationTest.cpp
@@ -0,0 +1,430 @@
+//======================================================================================================================
+//
+// 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 FieldInterpolationTest.cpp
+//! \ingroup field
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//
+//======================================================================================================================
+
+#include "blockforest/all.h"
+
+#include "core/debug/TestSubsystem.h"
+#include "core/mpi/Environment.h"
+#include "core/math/all.h"
+
+#include "field/AddToStorage.h"
+#include "field/FlagField.h"
+#include "field/GhostLayerField.h"
+#include "field/interpolators/all.h"
+
+#include <vector>
+
+using namespace walberla;
+
+namespace field_interpolation_tests {
+
+const uint_t FieldGhostLayers( 1 );
+
+typedef walberla::uint8_t flag_t;
+typedef FlagField< flag_t > FlagField_T;
+
+typedef GhostLayerField< real_t, 1> ScalarField_T;
+typedef GhostLayerField< Vector3<real_t>, 1> Vec3Field_T;
+typedef GhostLayerField< real_t, 3> MultiComponentField_T;
+
+
+const FlagUID Domain_Flag ( "domain" );
+const FlagUID Boundary_Flag ( "boundary" );
+
+void initFlagField( FlagField_T * field, IBlock * const /*block*/ )
+{
+ auto domainFlag = field->getOrRegisterFlag( Domain_Flag );
+ WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ( field, field->addFlag( x, y, z, domainFlag ); );
+}
+
+void initScalarField( const shared_ptr<StructuredBlockStorage> & blocks, const BlockDataID & scalarFieldID )
+{
+ for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt )
+ {
+ auto field = blockIt->getData<ScalarField_T>( scalarFieldID );
+ WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(field,
+ const Vector3<real_t> p = blocks->getBlockLocalCellCenter(*blockIt, Cell(x,y,z)) - Vector3<real_t>(real_t(0.5));
+ //field->get(x,y,z) = real_t(2);
+ field->get(x,y,z) = p[0];
+ );
+ }
+}
+
+void initVectorField( const shared_ptr<StructuredBlockStorage> & blocks, const BlockDataID & vectorFieldID )
+{
+ for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt )
+ {
+ auto field = blockIt->getData<Vec3Field_T>( vectorFieldID );
+ WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(field,
+ const Vector3<real_t> p = blocks->getBlockLocalCellCenter(*blockIt, Cell(x,y,z)) - Vector3<real_t>(real_t(0.5));
+ field->get(x,y,z) = Vector3<real_t>(p[0], real_t(2), p[1]);
+ );
+ }
+}
+
+void initMultiComponentField( const shared_ptr<StructuredBlockStorage> & blocks, const BlockDataID & multiComponentFieldID )
+{
+ for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt )
+ {
+ auto field = blockIt->getData<MultiComponentField_T>( multiComponentFieldID );
+ WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(field,
+ const Vector3<real_t> p = blocks->getBlockLocalCellCenter(*blockIt, Cell(x,y,z)) - Vector3<real_t>(real_t(0.5));
+ field->get(x,y,z,0) = p[0];
+ field->get(x,y,z,1) = real_t(2);
+ field->get(x,y,z,2) = p[1];
+ );
+ }
+}
+
+void setBoundaryFlags( const shared_ptr<StructuredBlockStorage> & blocks,
+ const BlockDataID & flagFieldID, const BlockDataID & scalarFieldID )
+{
+ for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt )
+ {
+ auto flagField = blockIt->getData<FlagField_T>( flagFieldID );
+ auto valueField = blockIt->getData<ScalarField_T>( scalarFieldID );
+ auto domainFlag = flagField->getOrRegisterFlag( Domain_Flag );
+ auto boundaryFlag = flagField->getOrRegisterFlag( Boundary_Flag );
+ WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(flagField,
+ if( x == 2)
+ {
+ flagField->removeFlag(x,y,z,domainFlag);
+ flagField->addFlag(x,y,z,boundaryFlag);
+ valueField->get(x,y,z) = std::numeric_limits<real_t>::max();
+ }
+ );
+ }
+}
+
+void testNearestNeighborFieldInterpolator( const shared_ptr<StructuredBlockStorage> & blocks, const BlockDataID & flagFieldID,
+ const BlockDataID & scalarFieldID, const BlockDataID & vectorFieldID, const BlockDataID & multiComponentFieldID )
+{
+ // field interpolators
+ typedef field::NearestNeighborFieldInterpolator<ScalarField_T, FlagField_T> ScalarFieldInterpolator_T;
+ typedef field::NearestNeighborFieldInterpolator<Vec3Field_T, FlagField_T> Vec3FieldInterpolator_T;
+ typedef field::NearestNeighborFieldInterpolator<MultiComponentField_T, FlagField_T> MultiComponentFieldInterpolator_T;
+ BlockDataID scalarFieldInterpolatorID = field::addFieldInterpolator< ScalarFieldInterpolator_T, FlagField_T >( blocks, scalarFieldID, flagFieldID, Domain_Flag );
+ BlockDataID vectorFieldInterpolatorID = field::addFieldInterpolator< Vec3FieldInterpolator_T, FlagField_T >( blocks, vectorFieldID, flagFieldID, Domain_Flag );
+ BlockDataID multiComponentFieldInterpolatorID = field::addFieldInterpolator< MultiComponentFieldInterpolator_T, FlagField_T >( blocks, multiComponentFieldID, flagFieldID, Domain_Flag );
+
+ // check scalar interpolation
+ {
+ Vector3<real_t> interpolationPoint(real_t(1.9), real_t(2.1), real_t(2.1));
+ auto containingBlockID = blocks->getBlock(interpolationPoint);
+ if( containingBlockID != NULL )
+ {
+ real_t interpolatedValue(real_t(0));
+ auto interPtr = containingBlockID->getData<ScalarFieldInterpolator_T>(scalarFieldInterpolatorID);
+ interPtr->get(interpolationPoint, &interpolatedValue);
+ WALBERLA_CHECK_FLOAT_EQUAL(interpolatedValue, real_t(1), "NearestNeighborFieldInterpolator: Scalar interpolation failed");
+ }
+ }
+
+ // check vector interpolation
+ {
+ Vector3<real_t> interpolationPoint(real_t(5.4),real_t(2.1),real_t(3.2));
+ auto containingBlockID = blocks->getBlock( interpolationPoint );
+ if( containingBlockID != NULL ) {
+ Vector3<real_t> interpolatedValue(real_t(0));
+ auto interPtr = containingBlockID->getData<Vec3FieldInterpolator_T>(vectorFieldInterpolatorID);
+ interPtr->get(interpolationPoint, &interpolatedValue);
+ WALBERLA_CHECK_FLOAT_EQUAL(interpolatedValue[0], real_t(5), "NearestNeighborFieldInterpolator: Vec3[0] interpolation failed");
+ WALBERLA_CHECK_FLOAT_EQUAL(interpolatedValue[1], real_t(2), "NearestNeighborFieldInterpolator: Vec3[1] interpolation failed");
+ WALBERLA_CHECK_FLOAT_EQUAL(interpolatedValue[2], real_t(2), "NearestNeighborFieldInterpolator: Vec3[2] interpolation failed");
+ }
+ }
+
+ // check multi component interpolation
+ {
+ Vector3<real_t> interpolationPoint(real_t(4.4),real_t(2.1),real_t(3.2));
+ auto containingBlockID = blocks->getBlock( interpolationPoint );
+ if( containingBlockID != NULL ) {
+ std::vector<real_t> interpolatedValue(3, real_t(0));
+ auto interPtr = containingBlockID->getData<MultiComponentFieldInterpolator_T>(multiComponentFieldInterpolatorID);
+ interPtr->get(interpolationPoint, interpolatedValue.begin());
+ WALBERLA_CHECK_FLOAT_EQUAL(interpolatedValue[0], real_t(4), "NearestNeighborFieldInterpolator: Multi Component[0] interpolation failed");
+ WALBERLA_CHECK_FLOAT_EQUAL(interpolatedValue[1], real_t(2), "NearestNeighborFieldInterpolator: Multi Component[1] interpolation failed");
+ WALBERLA_CHECK_FLOAT_EQUAL(interpolatedValue[2], real_t(2), "NearestNeighborFieldInterpolator: Multi Component[2] interpolation failed");
+ }
+ }
+}
+
+void testTrilinearFieldInterpolator( const shared_ptr<StructuredBlockStorage> & blocks, const BlockDataID & flagFieldID,
+ const BlockDataID & scalarFieldID, const BlockDataID & vectorFieldID, const BlockDataID & multiComponentFieldID )
+{
+ // field interpolators
+ typedef field::TrilinearFieldInterpolator<ScalarField_T, FlagField_T> ScalarFieldInterpolator_T;
+ typedef field::TrilinearFieldInterpolator<Vec3Field_T, FlagField_T> Vec3FieldInterpolator_T;
+ typedef field::TrilinearFieldInterpolator<MultiComponentField_T, FlagField_T> MultiComponentFieldInterpolator_T;
+ BlockDataID scalarFieldInterpolatorID = field::addFieldInterpolator< ScalarFieldInterpolator_T, FlagField_T >( blocks, scalarFieldID, flagFieldID, Domain_Flag );
+ BlockDataID vectorFieldInterpolatorID = field::addFieldInterpolator< Vec3FieldInterpolator_T, FlagField_T >( blocks, vectorFieldID, flagFieldID, Domain_Flag );
+ BlockDataID multiComponentFieldInterpolatorID = field::addFieldInterpolator< MultiComponentFieldInterpolator_T, FlagField_T >( blocks, multiComponentFieldID, flagFieldID, Domain_Flag );
+
+ // check scalar interpolation
+ {
+ Vector3<real_t> interpolationPoint(real_t(1.9), real_t(2.1), real_t(2.1));
+ auto containingBlockID = blocks->getBlock(interpolationPoint);
+ if( containingBlockID != NULL )
+ {
+ real_t interpolatedValue(real_t(0));
+ auto interPtr = containingBlockID->getData<ScalarFieldInterpolator_T>(scalarFieldInterpolatorID);
+ interPtr->get(interpolationPoint, &interpolatedValue);
+ WALBERLA_CHECK_FLOAT_EQUAL(interpolatedValue, real_t(1.4), "TrilinearFieldInterpolator: Scalar interpolation failed");
+ }
+ }
+
+ // check vector interpolation
+ {
+ Vector3<real_t> interpolationPoint(real_t(5.4),real_t(2.1),real_t(3.2));
+ auto containingBlockID = blocks->getBlock( interpolationPoint );
+ if( containingBlockID != NULL ) {
+ Vector3<real_t> interpolatedValue(real_t(0));
+ auto interPtr = containingBlockID->getData<Vec3FieldInterpolator_T>(vectorFieldInterpolatorID);
+ interPtr->get(interpolationPoint, &interpolatedValue);
+ WALBERLA_CHECK_FLOAT_EQUAL(interpolatedValue[0], real_t(4.9), "TrilinearFieldInterpolator: Vec3[0] interpolation failed");
+ WALBERLA_CHECK_FLOAT_EQUAL(interpolatedValue[1], real_t(2.0), "TrilinearFieldInterpolator: Vec3[1] interpolation failed");
+ WALBERLA_CHECK_FLOAT_EQUAL(interpolatedValue[2], real_t(1.6), "TrilinearFieldInterpolator: Vec3[2] interpolation failed");
+ }
+ }
+
+ // check multi component interpolation
+ {
+ Vector3<real_t> interpolationPoint(real_t(4.4),real_t(2.1),real_t(3.2));
+ auto containingBlockID = blocks->getBlock( interpolationPoint );
+ if( containingBlockID != NULL ) {
+ std::vector<real_t> interpolatedValue(3, real_t(0));
+ auto interPtr = containingBlockID->getData<MultiComponentFieldInterpolator_T>(multiComponentFieldInterpolatorID);
+ interPtr->get(interpolationPoint, interpolatedValue.begin());
+ WALBERLA_CHECK_FLOAT_EQUAL(interpolatedValue[0], real_t(3.9), "TrilinearFieldInterpolator: Multi Component[0] interpolation failed");
+ WALBERLA_CHECK_FLOAT_EQUAL(interpolatedValue[1], real_t(2.0), "TrilinearFieldInterpolator: Multi Component[1] interpolation failed");
+ WALBERLA_CHECK_FLOAT_EQUAL(interpolatedValue[2], real_t(1.6), "TrilinearFieldInterpolator: Multi Component[2] interpolation failed");
+ }
+ }
+}
+
+void testKernelFieldInterpolator( const shared_ptr<StructuredBlockStorage> & blocks, const BlockDataID & flagFieldID,
+ const BlockDataID & scalarFieldID, const BlockDataID & vectorFieldID, const BlockDataID & multiComponentFieldID )
+{
+ // field interpolators
+ typedef field::KernelFieldInterpolator<ScalarField_T, FlagField_T> ScalarFieldInterpolator_T;
+ typedef field::KernelFieldInterpolator<Vec3Field_T, FlagField_T> Vec3FieldInterpolator_T;
+ typedef field::KernelFieldInterpolator<MultiComponentField_T, FlagField_T> MultiComponentFieldInterpolator_T;
+ BlockDataID scalarFieldInterpolatorID = field::addFieldInterpolator< ScalarFieldInterpolator_T, FlagField_T >( blocks, scalarFieldID, flagFieldID, Domain_Flag );
+ BlockDataID vectorFieldInterpolatorID = field::addFieldInterpolator< Vec3FieldInterpolator_T, FlagField_T >( blocks, vectorFieldID, flagFieldID, Domain_Flag );
+ BlockDataID multiComponentFieldInterpolatorID = field::addFieldInterpolator< MultiComponentFieldInterpolator_T, FlagField_T >( blocks, multiComponentFieldID, flagFieldID, Domain_Flag );
+
+ // check scalar interpolation
+ {
+ Vector3<real_t> interpolationPoint(real_t(1.9), real_t(2.1), real_t(2.1));
+ auto containingBlockID = blocks->getBlock(interpolationPoint);
+ if( containingBlockID != NULL )
+ {
+ real_t interpolatedValue(real_t(0));
+ auto interPtr = containingBlockID->getData<ScalarFieldInterpolator_T>(scalarFieldInterpolatorID);
+ interPtr->get(interpolationPoint, &interpolatedValue);
+ WALBERLA_CHECK_FLOAT_EQUAL(interpolatedValue, real_t(1.4), "KernelFieldInterpolator: Scalar interpolation failed");
+ }
+ }
+
+ // check vector interpolation
+ {
+ Vector3<real_t> interpolationPoint(real_t(5.4),real_t(2.1),real_t(3.2));
+ auto containingBlockID = blocks->getBlock( interpolationPoint );
+ if( containingBlockID != NULL ) {
+ Vector3<real_t> interpolatedValue(real_t(0));
+ auto interPtr = containingBlockID->getData<Vec3FieldInterpolator_T>(vectorFieldInterpolatorID);
+ interPtr->get(interpolationPoint, &interpolatedValue);
+ WALBERLA_CHECK_FLOAT_EQUAL(interpolatedValue[0], real_t(4.9), "KernelFieldInterpolator: Vec3[0] interpolation failed");
+ WALBERLA_CHECK_FLOAT_EQUAL(interpolatedValue[1], real_t(2.0), "KernelFieldInterpolator: Vec3[1] interpolation failed");
+ WALBERLA_CHECK_FLOAT_EQUAL(interpolatedValue[2], real_t(1.6), "KernelFieldInterpolator: Vec3[2] interpolation failed");
+ }
+ }
+
+ // check multi component interpolation
+ {
+ Vector3<real_t> interpolationPoint(real_t(4.4),real_t(2.1),real_t(3.2));
+ auto containingBlockID = blocks->getBlock( interpolationPoint );
+ if( containingBlockID != NULL ) {
+ std::vector<real_t> interpolatedValue(3, real_t(0));
+ auto interPtr = containingBlockID->getData<MultiComponentFieldInterpolator_T>(multiComponentFieldInterpolatorID);
+ interPtr->get(interpolationPoint, interpolatedValue.begin());
+ WALBERLA_CHECK_FLOAT_EQUAL(interpolatedValue[0], real_t(3.9), "KernelFieldInterpolator: Multi Component[0] interpolation failed");
+ WALBERLA_CHECK_FLOAT_EQUAL(interpolatedValue[1], real_t(2.0), "KernelFieldInterpolator: Multi Component[1] interpolation failed");
+ WALBERLA_CHECK_FLOAT_EQUAL(interpolatedValue[2], real_t(1.6), "KernelFieldInterpolator: Multi Component[2] interpolation failed");
+ }
+ }
+}
+
+
+void testNearestNeighborFieldInterpolatorAtBoundary( const shared_ptr<StructuredBlockStorage> & blocks,
+ const BlockDataID & flagFieldID, const BlockDataID & scalarFieldID ) {
+ // field interpolators
+ typedef field::NearestNeighborFieldInterpolator<ScalarField_T, FlagField_T> ScalarFieldInterpolator_T;
+ BlockDataID scalarFieldInterpolatorID = field::addFieldInterpolator<ScalarFieldInterpolator_T, FlagField_T>(blocks, scalarFieldID, flagFieldID, Domain_Flag);
+
+ // check scalar interpolation close to boundary
+ {
+ Vector3<real_t> interpolationPoint(real_t(1.9), real_t(2.1), real_t(2.1));
+ auto containingBlockID = blocks->getBlock(interpolationPoint);
+ if (containingBlockID != NULL) {
+ real_t interpolatedValue(real_t(0));
+ auto interPtr = containingBlockID->getData<ScalarFieldInterpolator_T>(scalarFieldInterpolatorID);
+ interPtr->get(interpolationPoint, &interpolatedValue);
+ WALBERLA_CHECK_FLOAT_EQUAL(interpolatedValue, real_t(1),
+ "NearestNeighborFieldInterpolator: Scalar interpolation near boundary failed");
+ }
+ }
+
+ // check scalar interpolation inside boundary
+ {
+ Vector3<real_t> interpolationPoint(real_t(2.7), real_t(2.1), real_t(1.1));
+ auto containingBlockID = blocks->getBlock(interpolationPoint);
+ if (containingBlockID != NULL) {
+ real_t interpolatedValue(real_t(0));
+ auto interPtr = containingBlockID->getData<ScalarFieldInterpolator_T>(scalarFieldInterpolatorID);
+ interPtr->get(interpolationPoint, &interpolatedValue);
+ WALBERLA_CHECK_FLOAT_EQUAL(interpolatedValue, real_t(3),
+ "NearestNeighborFieldInterpolator: Scalar interpolation inside boundary failed");
+ }
+ }
+}
+
+void testTrilinearFieldInterpolatorAtBoundary( const shared_ptr<StructuredBlockStorage> & blocks,
+ const BlockDataID & flagFieldID, const BlockDataID & scalarFieldID ) {
+ // field interpolators
+ typedef field::TrilinearFieldInterpolator<ScalarField_T, FlagField_T> ScalarFieldInterpolator_T;
+ BlockDataID scalarFieldInterpolatorID = field::addFieldInterpolator<ScalarFieldInterpolator_T, FlagField_T>(blocks, scalarFieldID, flagFieldID, Domain_Flag);
+
+ // check scalar interpolation close to boundary
+ {
+ Vector3<real_t> interpolationPoint(real_t(1.9), real_t(2.1), real_t(2.1));
+ auto containingBlockID = blocks->getBlock(interpolationPoint);
+ if (containingBlockID != NULL) {
+ real_t interpolatedValue(real_t(0));
+ auto interPtr = containingBlockID->getData<ScalarFieldInterpolator_T>(scalarFieldInterpolatorID);
+ interPtr->get(interpolationPoint, &interpolatedValue);
+ WALBERLA_CHECK_FLOAT_EQUAL(interpolatedValue, real_t(1),
+ "TrilinearFieldInterpolator: Scalar interpolation near boundary failed");
+ }
+ }
+
+ // check scalar interpolation inside boundary
+ {
+ Vector3<real_t> interpolationPoint(real_t(2.7), real_t(2.1), real_t(1.1));
+ auto containingBlockID = blocks->getBlock(interpolationPoint);
+ if (containingBlockID != NULL) {
+ real_t interpolatedValue(real_t(0));
+ auto interPtr = containingBlockID->getData<ScalarFieldInterpolator_T>(scalarFieldInterpolatorID);
+ interPtr->get(interpolationPoint, &interpolatedValue);
+ WALBERLA_CHECK_FLOAT_EQUAL(interpolatedValue, real_t(3),
+ "TrilinearFieldInterpolator: Scalar interpolation inside boundary failed");
+ }
+ }
+}
+
+void testKernelFieldInterpolatorAtBoundary( const shared_ptr<StructuredBlockStorage> & blocks,
+ const BlockDataID & flagFieldID, const BlockDataID & scalarFieldID ) {
+ // field interpolators
+ typedef field::KernelFieldInterpolator<ScalarField_T, FlagField_T> ScalarFieldInterpolator_T;
+ BlockDataID scalarFieldInterpolatorID = field::addFieldInterpolator<ScalarFieldInterpolator_T, FlagField_T>(blocks, scalarFieldID, flagFieldID, Domain_Flag);
+
+ // check scalar interpolation close to boundary
+ {
+ Vector3<real_t> interpolationPoint(real_t(1.9), real_t(2.1), real_t(2.1));
+ auto containingBlockID = blocks->getBlock(interpolationPoint);
+ if (containingBlockID != NULL) {
+ real_t interpolatedValue(real_t(0));
+ auto interPtr = containingBlockID->getData<ScalarFieldInterpolator_T>(scalarFieldInterpolatorID);
+ interPtr->get(interpolationPoint, &interpolatedValue);
+ // kernel interpolation can not extrapolate values from the available ones (see comments in KernelFieldInterpolator.h)
+ // it will thus yield a value between the available ones, which are 0 and 1
+ WALBERLA_CHECK(interpolatedValue < real_t(1),
+ "KernelFieldInterpolator: Scalar interpolation near boundary failed");
+ WALBERLA_CHECK(interpolatedValue > real_t(0),
+ "KernelFieldInterpolator: Scalar interpolation near boundary failed");
+ }
+ }
+
+ // check scalar interpolation inside boundary
+ {
+ Vector3<real_t> interpolationPoint(real_t(2.7), real_t(2.1), real_t(1.1));
+ auto containingBlockID = blocks->getBlock(interpolationPoint);
+ if (containingBlockID != NULL) {
+ real_t interpolatedValue(real_t(0));
+ auto interPtr = containingBlockID->getData<ScalarFieldInterpolator_T>(scalarFieldInterpolatorID);
+ interPtr->get(interpolationPoint, &interpolatedValue);
+ // values has to be between the available ones, i.e. 1 and 3
+ WALBERLA_CHECK(interpolatedValue > real_t(1),
+ "KernelFieldInterpolator: Scalar interpolation inside boundary failed");
+ WALBERLA_CHECK(interpolatedValue < real_t(3),
+ "KernelFieldInterpolator: Scalar interpolation inside boundary failed");
+ }
+ }
+}
+
+int main(int argc, char **argv) {
+
+ mpi::Environment mpiEnv(argc, argv);
+ debug::enterTestMode();
+
+ const uint_t numberOfBlocksInDirection = 2;
+ const uint_t numberOfCellsPerBlockInDirection = 4;
+ const real_t dx = real_t(1);
+
+ // block storage
+ auto blocks = blockforest::createUniformBlockGrid( numberOfBlocksInDirection, numberOfBlocksInDirection, numberOfBlocksInDirection,
+ numberOfCellsPerBlockInDirection, numberOfCellsPerBlockInDirection, numberOfCellsPerBlockInDirection,
+ dx, 0, false, false,
+ false, false, false,
+ false );
+ // flag field
+ BlockDataID flagFieldID = field::addFlagFieldToStorage<FlagField_T>( blocks, "flag field", FieldGhostLayers, false, initFlagField );
+
+ // data fields
+ BlockDataID scalarFieldID = field::addToStorage< ScalarField_T >( blocks, "scalar field", real_t(0), field::zyxf, FieldGhostLayers );
+ BlockDataID vectorFieldID = field::addToStorage< Vec3Field_T >( blocks, "vec3 field", Vector3<real_t>(real_t(0)), field::zyxf, FieldGhostLayers );
+ BlockDataID multiComponentFieldID = field::addToStorage< MultiComponentField_T >( blocks, "multi component field", real_t(0), field::zyxf, FieldGhostLayers );
+
+ initScalarField(blocks, scalarFieldID);
+ initVectorField(blocks, vectorFieldID );
+ initMultiComponentField(blocks, multiComponentFieldID );
+
+ // test all interpolators with domain flags everywhere, i.e. without special boundary treatment necessary
+ testNearestNeighborFieldInterpolator(blocks, flagFieldID, scalarFieldID, vectorFieldID, multiComponentFieldID);
+ testTrilinearFieldInterpolator(blocks, flagFieldID, scalarFieldID, vectorFieldID, multiComponentFieldID);
+ testKernelFieldInterpolator(blocks, flagFieldID, scalarFieldID, vectorFieldID, multiComponentFieldID);
+
+ // set some boundary flags in flag field and invalidate the corresponding scalar field values
+ setBoundaryFlags(blocks, flagFieldID, scalarFieldID);
+
+ // test all interpolators' behavior close to boundary cells
+ testNearestNeighborFieldInterpolatorAtBoundary(blocks, flagFieldID, scalarFieldID);
+ testTrilinearFieldInterpolatorAtBoundary(blocks, flagFieldID, scalarFieldID);
+ testKernelFieldInterpolatorAtBoundary(blocks, flagFieldID, scalarFieldID);
+
+ return 0;
+}
+
+} // namespace field_interpolation_tests
+
+int main( int argc, char **argv ){
+ field_interpolation_tests::main(argc, argv);
+}
\ No newline at end of file