src/OpenFOAM/interpolations/interpolateSplineXY/interpolateSplineXY.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
   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
  13     the Free Software Foundation, either version 3 of the License, or
  14     (at your 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, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "interpolateSplineXY.H"
  27 #include "primitiveFields.H"
  28
  29 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  30
  31 template<class Type>
  32 Foam::Field<Type> Foam::interpolateSplineXY
  33 (
  34     const scalarField& xNew,
  35     const scalarField& xOld,
  36     const Field<Type>& yOld
  37 )
  38 {
  39     Field<Type> yNew(xNew.size());
  40
  41     forAll(xNew, i)
  42     {
  43         yNew[i] = interpolateSmoothXY(xNew[i], xOld, yOld);
  44     }
  45
  46     return yNew;
  47 }
  48
  49
  50 template<class Type>
  51 Type Foam::interpolateSplineXY
  52 (
  53     const scalar x,
  54     const scalarField& xOld,
  55     const Field<Type>& yOld
  56 )
  57 {
  58     label n = xOld.size();
  59
  60     // early exit if out of bounds or only one value
  61     if (n == 1 || x < xOld[0])
  62     {
  63         return yOld[0];
  64     }
  65     if (x > xOld[n - 1])
  66     {
  67         return yOld[n - 1];
  68     }
  69
  70     // linear interpolation if only two values
  71     if (n == 2)
  72     {
  73         return (x - xOld[0])/(xOld[1] - xOld[0])*(yOld[1] - yOld[0]) + yOld[0];
  74     }
  75
  76     // find bounding knots
  77     label hi = 0;
  78     while (hi < n && xOld[hi] < x)
  79     {
  80         hi++;
  81     }
  82
  83     label lo = hi - 1;
  84
  85     const Type& y1 = yOld[lo];
  86     const Type& y2 = yOld[hi];
  87
  88     Type y0;
  89     if (lo == 0)
  90     {
  91         y0 = 2*y1 - y2;
  92     }
  93     else
  94     {
  95         y0 = yOld[lo - 1];
  96     }
  97
  98     Type y3;
  99     if (hi + 1 == n)
 100     {
 101         y3 = 2*y2 - y1;
 102     }
 103     else
 104     {
 105         y3 = yOld[hi + 1];
 106     }
 107
 108     // weighting
 109     scalar mu = (x - xOld[lo])/(xOld[hi] - xOld[lo]);
 110
 111     // interpolate
 112     return
 113         0.5
 114        *(
 115             2*y1
 116           + mu
 117            *(
 118                -y0 + y2
 119               + mu*((2*y0 - 5*y1 + 4*y2 - y3) + mu*(-y0 + 3*y1 - 3*y2 + y3))
 120             )
 121         );
 122 }
 123
 124
 125 // ************************************************************************* //