interpolateSplineXY.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2018 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>
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] = interpolateSplineXY(xNew[i], xOld, yOld);
44  }
45 
46  return yNew;
47 }
48 
49 
50 template<class Type>
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 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
Interpolates y values from one curve to another with a different x distribution.
dimensionedScalar y0(const dimensionedScalar &ds)
Field< Type > interpolateSplineXY(const scalarField &xNew, const scalarField &xOld, const Field< Type > &yOld)
Pre-declare SubField and related Field type.
Definition: Field.H:56
dimensionedScalar y1(const dimensionedScalar &ds)
volScalarField scalarField(fieldObject, mesh)
Specialisations of Field<T> for scalar, vector and tensor.
const dimensionedScalar & mu
Atomic mass unit.
label n