src/OpenFOAM/interpolations/RBFInterpolation/RBFInterpolationTemplates.C

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