linearInterpolationWeights.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) 2012-2019 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 
28 #include "ListOps.H"
29 #include "Pair.H"
30 
31 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
32 
33 namespace Foam
34 {
35 
36 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
37 
38 defineTypeNameAndDebug(linearInterpolationWeights, 0);
40 (
41  interpolationWeights,
42  linearInterpolationWeights,
43  word
44 );
45 
46 
47 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
48 
50 (
51  const scalarField& samples
52 )
53 :
54  interpolationWeights(samples),
55  index_(-1)
56 {}
57 
58 
59 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
60 
62 (
63  const scalar t,
64  labelList& indices,
65  scalarField& weights
66 ) const
67 {
68  // Check if current index is still valid
69  bool changed = false;
70  if
71  (
72  (
73  index_ == -1
74  && t <= samples_.first()
75  )
76  || (
77  index_ >= 0
78  && index_ < samples_.size() - 1
79  && t >= samples_[index_]
80  && t <= samples_[index_ + 1]
81  )
82  || (
83  index_ == samples_.size() - 1
84  && t >= samples_.last()
85  )
86  )
87  {
88  // The index is still in the correct slot
89  }
90  else
91  {
92  // The index is no longer in the correct slot, so search for a new one
93  index_ = findLower(samples_, t);
94  changed = true;
95  }
96 
97  // Calculate the number of indices and weights and resize the result
98  const label n = index_ == -1 || index_ == samples_.size() - 1 ? 1 : 2;
99  indices.resize(n);
100  weights.resize(n);
101 
102  // Compute the value
103  if (index_ == -1)
104  {
105  // Use the first value
106  indices[0] = 0;
107  weights[0] = 1;
108  }
109  else if (index_ == samples_.size() - 1)
110  {
111  // Use the last value
112  indices[0] = samples_.size() - 1;
113  weights[0] = 1;
114  }
115  else
116  {
117  // Interpolate within the interval
118  const scalar f =
119  (t - samples_[index_])/(samples_[index_ + 1] - samples_[index_]);
120  indices[0] = index_;
121  weights[0] = 1 - f;
122  indices[1] = index_ + 1;
123  weights[1] = f;
124  }
125 
126  return changed;
127 }
128 
129 
131 (
132  scalar t1,
133  scalar t2,
134  labelList& indices,
135  scalarField& weights
136 ) const
137 {
138  // If the arguments are in descending order, then swap them and set the
139  // weights' sign negative
140  label sign = +1;
141  if (t1 > t2)
142  {
143  Swap(t1, t2);
144  sign = -1;
145  }
146 
147  //- Search for lower indices
148  // Note: currently there is no caching of this search like in valueWeights
149  const label i1 = findLower(samples_, t1);
150  const label i2 = findLower(samples_, t2);
151  const label iClip1 = min(max(i1, 0), samples_.size() - 2);
152  const label iClip2 = min(max(i2, 0), samples_.size() - 2);
153  const label n = max(i2 - i1 + (i1 == iClip1) + (i2 == iClip2), 1);
154 
155  // Check if anything changed
156  bool changed = false;
157  if (indices.size() != n)
158  {
159  changed = true;
160  }
161  else
162  {
163  forAll(indices, indexi)
164  {
165  if (indices[indexi] == indexi + max(i1, 0))
166  {
167  changed = true;
168  break;
169  }
170  }
171  }
172 
173  // Resize the result arrays
174  indices.resize(n);
175  indices = -1;
176  weights.resize(n);
177  weights = 0;
178 
179  // Add out of bounds interval below the table
180  if (i1 == -1)
181  {
182  indices[0] = 0;
183  weights[0] += sign*(samples_[0] - t1);
184  }
185  if (i2 == -1)
186  {
187  indices[0] = 0;
188  weights[0] -= sign*(samples_[0] - t2);
189  }
190 
191  // Add partial interval from t1 to i1 + 1
192  if (i1 == iClip1)
193  {
194  const scalar f = (t1 - samples_[i1])/(samples_[i1 + 1] - samples_[i1]);
195  const scalar d = samples_[i1 + 1] - t1;
196  indices[0] = i1;
197  weights[0] += sign*(1 - f)*d/2;
198  indices[1] = i1 + 1;
199  weights[1] += sign*(1 + f)*d/2;
200  }
201 
202  // Sum whole intervals from i1 + 1 to i2
203  if (i1 != i2) for (label i = i1 + 1; i <= iClip2; i ++)
204  {
205  const scalar d = samples_[i + 1] - samples_[i];
206  indices[i - iClip1] = i;
207  weights[i - iClip1] += sign*d/2;
208  indices[i - iClip1 + 1] = i + 1;
209  weights[i - iClip1 + 1] += sign*d/2;
210  }
211 
212  // Subtract partial interval from t2 to i2 + 1
213  if (i2 == iClip2)
214  {
215  const scalar f = (t2 - samples_[i2])/(samples_[i2 + 1] - samples_[i2]);
216  const scalar d = samples_[i2 + 1] - t2;
217  indices[n - 2] = i2;
218  weights[n - 2] -= sign*(1 - f)*d/2;
219  indices[n - 1] = i2 + 1;
220  weights[n - 1] -= sign*(1 + f)*d/2;
221  }
222 
223  // Add out of bounds interval above the table
224  if (i1 == samples_.size() - 1)
225  {
226  indices[n - 1] = samples_.size() - 1;
227  weights[n - 1] -= sign*(t1 - samples_.last());
228  }
229  if (i2 == samples_.size() - 1)
230  {
231  indices[n - 1] = samples_.size() - 1;
232  weights[n - 1] += sign*(t2 - samples_.last());
233  }
234 
235  return changed;
236 }
237 
238 
239 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
240 
241 } // End namespace Foam
242 
243 // ************************************************************************* //
dimensionedScalar sign(const dimensionedScalar &ds)
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
void resize(const label)
Alias for setSize(const label)
Definition: ListI.H:138
Abstract base class for interpolating in 1D.
Macros for easy insertion into run-time selection tables.
Various functions to operate on Lists.
virtual bool valueWeights(const scalar t, labelList &indices, scalarField &weights) const
Calculate weights and indices to calculate t from samples.
void Swap(T &a, T &b)
Definition: Swap.H:43
linearInterpolationWeights(const scalarField &samples)
Construct from components.
layerAndWeight min(const layerAndWeight &a, const layerAndWeight &b)
addToRunTimeSelectionTable(ensightPart, ensightPartCells, istream)
defineTypeNameAndDebug(combustionModel, 0)
labelList f(nPoints)
virtual bool integrationWeights(scalar t1, scalar t2, labelList &indices, scalarField &weights) const
Calculate weights and indices to calculate integrand of t1..t2.
label findLower(const ListType &, typename ListType::const_reference, const label stary, const BinaryOp &bop)
Find last element < given value in sorted list and return index,.
label n
Namespace for OpenFOAM.