interpolateSplineXY.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | www.openfoam.com
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8  Copyright (C) 2011 OpenFOAM Foundation
9 -------------------------------------------------------------------------------
10 License
11  This file is part of OpenFOAM.
12 
13  OpenFOAM is free software: you can redistribute it and/or modify it
14  under the terms of the GNU General Public License as published by
15  the Free Software Foundation, either version 3 of the License, or
16  (at your option) any later version.
17 
18  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
19  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
20  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
21  for more details.
22 
23  You should have received a copy of the GNU General Public License
24  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
25 
26 \*---------------------------------------------------------------------------*/
27 
28 #include "interpolateSplineXY.H"
29 #include "primitiveFields.H"
30 
31 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
32 
33 template<class Type>
35 (
36  const scalarField& xNew,
37  const scalarField& xOld,
38  const Field<Type>& yOld
39 )
40 {
41  Field<Type> yNew(xNew.size());
42 
43  forAll(xNew, i)
44  {
45  yNew[i] = interpolateSplineXY(xNew[i], xOld, yOld);
46  }
47 
48  return yNew;
49 }
50 
51 
52 template<class Type>
54 (
55  const scalar x,
56  const scalarField& xOld,
57  const Field<Type>& yOld
58 )
59 {
60  label n = xOld.size();
61 
62  // early exit if out of bounds or only one value
63  if (n == 1 || x < xOld[0])
64  {
65  return yOld[0];
66  }
67  if (x > xOld[n - 1])
68  {
69  return yOld[n - 1];
70  }
71 
72  // linear interpolation if only two values
73  if (n == 2)
74  {
75  return (x - xOld[0])/(xOld[1] - xOld[0])*(yOld[1] - yOld[0]) + yOld[0];
76  }
77 
78  // find bounding knots
79  label hi = 0;
80  while (hi < n && xOld[hi] < x)
81  {
82  hi++;
83  }
84 
85  label lo = hi - 1;
86 
87  const Type& y1 = yOld[lo];
88  const Type& y2 = yOld[hi];
89 
90  Type y0;
91  if (lo == 0)
92  {
93  y0 = 2*y1 - y2;
94  }
95  else
96  {
97  y0 = yOld[lo - 1];
98  }
99 
100  Type y3;
101  if (hi + 1 == n)
102  {
103  y3 = 2*y2 - y1;
104  }
105  else
106  {
107  y3 = yOld[hi + 1];
108  }
109 
110  // weighting
111  scalar mu = (x - xOld[lo])/(xOld[hi] - xOld[lo]);
112 
113  // interpolate
114  return
115  0.5
116  *(
117  2*y1
118  + mu
119  *(
120  -y0 + y2
121  + mu*((2*y0 - 5*y1 + 4*y2 - y3) + mu*(-y0 + 3*y1 - 3*y2 + y3))
122  )
123  );
124 }
125 
126 
127 // ************************************************************************* //
Foam::scalarField
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Definition: primitiveFieldsFwd.H:52
Foam::y1
dimensionedScalar y1(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:282
Foam::constant::physicoChemical::mu
const dimensionedScalar mu
Atomic mass unit.
Definition: createFieldRefs.H:4
primitiveFields.H
Specialisations of Field<T> for scalar, vector and tensor.
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:296
n
label n
Definition: TABSMDCalcMethod2.H:31
Foam::Field
Generic templated field type.
Definition: Field.H:63
Foam::y0
dimensionedScalar y0(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:281
x
x
Definition: LISASMDCalcMethod2.H:52
interpolateSplineXY.H
Interpolates y values from one curve to another with a different x distribution.
Foam::interpolateSplineXY
Field< Type > interpolateSplineXY(const scalarField &xNew, const scalarField &xOld, const Field< Type > &yOld)