src/foam/interpolations/RBFInterpolation/RBFInterpolationTemplates.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend is free software: you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by the
  13     Free Software Foundation, either version 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 Description
  25     RBF interpolation templates
  26
  27 Author
  28     Frank Bos, TU Delft.  All rights reserved.
  29     Dubravko Matijasevic, FSB Zagreb.
  30
  31 \*---------------------------------------------------------------------------*/
  32
  33 #include "RBFInterpolation.H"
  34
  35 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  36
  37 template<class Type>
  38 Foam::tmp<Foam::Field<Type> > Foam::RBFInterpolation::interpolate
  39 (
  40     const Field<Type>& ctrlField
  41 ) const
  42 {
  43     // HJ and FB (05 Jan 2009)
  44     // Collect the values from ALL control points to all CPUs
  45     // Then, each CPU will do interpolation only on local dataPoints_
  46
  47     if (ctrlField.size() != controlPoints_.size())
  48     {
  49         FatalErrorIn
  50         (
  51             "tmp<Field<Type> > RBFInterpolation::interpolate\n"
  52             "(\n"
  53             "    const Field<Type>& ctrlField\n"
  54             ") const"
  55         )   << "Incorrect size of source field.  Size = " << ctrlField.size()
  56             << " nControlPoints = " << controlPoints_.size()
  57             << abort(FatalError);
  58     }
  59
  60     tmp<Field<Type> > tresult
  61     (
  62         new Field<Type>(dataPoints_.size(), pTraits<Type>::zero)
  63     );
  64
  65     Field<Type>& result = tresult();
  66
  67     // FB 21-12-2008
  68     // 1) Calculate alpha and beta coefficients using the Inverse
  69     // 2) Calculate displacements of internal nodes using RBF values,
  70     //    alpha's and beta's
  71     // 3) Return displacements using tresult()
  72
  73     const label nControlPoints = controlPoints_.size();
  74     const scalarSquareMatrix& mat = this->B();
  75
  76     // Determine interpolation coefficients
  77     Field<Type> alpha(nControlPoints, pTraits<Type>::zero);
  78     Field<Type> beta(4, pTraits<Type>::zero);
  79
  80     for (label row = 0; row < nControlPoints; row++)
  81     {
  82         for (label col = 0; col < nControlPoints; col++)
  83         {
  84             alpha[row] += mat[row][col]*ctrlField[col];
  85         }
  86     }
  87
  88     if (polynomials_)
  89     {
  90         for
  91         (
  92             label row = nControlPoints;
  93             row < nControlPoints + 4;
  94             row++
  95         )
  96         {
  97             for (label col = 0; col < nControlPoints; col++)
  98             {
  99                 beta[row - nControlPoints] += mat[row][col]*ctrlField[col];
 100             }
 101         }
 102     }
 103
 104     // Evaluation
 105     scalar t;
 106
 107     // Algorithmic improvement, Matteo Lombardi.  21/Mar/2011
 108
 109     forAll (dataPoints_, flPoint)
 110     {
 111         // Cut-off function to justify neglecting outer boundary points
 112         t = (mag(dataPoints_[flPoint] - focalPoint_) - innerRadius_)/
 113             (outerRadius_ - innerRadius_);
 114
 115         if (t >= 1)
 116         {
 117             // Increment is zero: w = 0
 118             result[flPoint] = 0*result[flPoint];
 119         }
 120         else
 121         {
 122             // Full calculation of weights
 123             scalarField weights =
 124                 RBF_->weights(controlPoints_, dataPoints_[flPoint]);
 125
 126             forAll (controlPoints_, i)
 127             {
 128                 result[flPoint] += weights[i]*alpha[i];
 129             }
 130
 131             if (polynomials_)
 132             {
 133                 result[flPoint] +=
 134                     beta[0]
 135                   + beta[1]*dataPoints_[flPoint].x()
 136                   + beta[2]*dataPoints_[flPoint].y()
 137                   + beta[3]*dataPoints_[flPoint].z();
 138             }
 139
 140             scalar w;
 141
 142             if (t <= 0)
 143             {
 144                 w = 1.0;
 145             }
 146             else
 147             {
 148                 w = 1 - sqr(t)*(3 - 2*t);
 149             }
 150
 151             result[flPoint] = w*result[flPoint];
 152         }
 153     }
 154
 155     return tresult;
 156 }
 157
 158
 159 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 160