linearInterpolationWeights.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 | Copyright (C) 2012 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 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
48 
49 Foam::Pair<Foam::scalar> linearInterpolationWeights::integrationWeights
50 (
51  const label i,
52  const scalar t
53 ) const
54 {
55  // t is in range samples_[i] .. samples_[i+1]
56 
57  scalar s = (t-samples_[i])/(samples_[i+1]-samples_[i]);
58 
59  if (s < -SMALL || s > 1+SMALL)
60  {
61  FatalErrorIn("linearInterpolationWeights::integrationWeights(..)")
62  << "Value " << t << " outside range " << samples_[i]
63  << " .. " << samples_[i+1]
64  << exit(FatalError);
65  }
66 
67  scalar d = samples_[i+1]-t;
68 
69  return Pair<scalar>(d*0.5*(1-s), d*0.5*(1+s));
70 }
71 
72 
73 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
74 
76 (
77  const scalarField& samples
78 )
79 :
80  interpolationWeights(samples),
81  index_(-1)
82 {}
83 
84 
85 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
86 
88 (
89  const scalar t,
90  labelList& indices,
91  scalarField& weights
92 ) const
93 {
94  bool indexChanged = false;
95 
96  // Check if current timeIndex is still valid
97  if
98  (
99  index_ >= 0
100  && index_ < samples_.size()
101  && (
102  samples_[index_] <= t
103  && (index_ == samples_.size()-1 || t <= samples_[index_+1])
104  )
105  )
106  {
107  // index_ still at correct slot
108  }
109  else
110  {
111  // search for correct index
112  index_ = findLower(samples_, t);
113  indexChanged = true;
114  }
115 
116 
117  if (index_ == -1)
118  {
119  // Use first element only
120  indices.setSize(1);
121  weights.setSize(1);
122  indices[0] = 0;
123  weights[0] = 1.0;
124  }
125  else if (index_ == samples_.size()-1)
126  {
127  // Use last element only
128  indices.setSize(1);
129  weights.setSize(1);
130  indices[0] = samples_.size()-1;
131  weights[0] = 1.0;
132  }
133  else
134  {
135  // Interpolate
136  indices.setSize(2);
137  weights.setSize(2);
138 
139  indices[0] = index_;
140  indices[1] = index_+1;
141 
142  scalar t0 = samples_[indices[0]];
143  scalar t1 = samples_[indices[1]];
144  scalar deltaT = t1-t0;
145 
146  weights[0] = (t1-t)/deltaT;
147  weights[1] = 1.0-weights[0];
148  }
149 
150  return indexChanged;
151 }
152 
153 
154 bool linearInterpolationWeights::integrationWeights
155 (
156  const scalar t1,
157  const scalar t2,
158  labelList& indices,
159  scalarField& weights
160 ) const
161 {
162  if (t2 < t1-VSMALL)
163  {
164  FatalErrorIn("linearInterpolationWeights::integrationWeights(..)")
165  << "Integration should be in positive direction."
166  << " t1:" << t1 << " t2:" << t2
167  << exit(FatalError);
168  }
169 
170  // Currently no fancy logic on cached index like in value
171 
172  //- Find lower or equal index
173  label i1 = findLower(samples_, t1, 0, lessEqOp<scalar>());
174  //- Find lower index
175  label i2 = findLower(samples_, t2);
176 
177  // For now just fail if any outside table
178  if (i1 == -1 || i2 == samples_.size()-1)
179  {
180  FatalErrorIn("linearInterpolationWeights::integrationWeights(..)")
181  << "Integrating outside table " << samples_[0] << ".."
182  << samples_.last() << " not implemented."
183  << " t1:" << t1 << " t2:" << t2 << exit(FatalError);
184  }
185 
186  label nIndices = i2-i1+2;
187 
188 
189  // Determine if indices already correct
190  bool anyChanged = false;
191 
192  if (nIndices != indices.size())
193  {
194  anyChanged = true;
195  }
196  else
197  {
198  // Closer check
199 
200  label index = i1;
201  forAll(indices, i)
202  {
203  if (indices[i] != index)
204  {
205  anyChanged = true;
206  break;
207  }
208  index++;
209  }
210  }
211 
212  indices.setSize(nIndices);
213  weights.setSize(nIndices);
214  weights = 0.0;
215 
216  // Sum from i1+1 to i2+1
217  for (label i = i1+1; i <= i2; i++)
218  {
219  scalar d = samples_[i+1]-samples_[i];
220  indices[i-i1] = i;
221  weights[i-i1] += 0.5*d;
222  indices[i+1-i1] = i+1;
223  weights[i+1-i1] += 0.5*d;
224  }
225 
226  // Add from i1 to t1
227  {
228  Pair<scalar> i1Tot1 = integrationWeights(i1, t1);
229  indices[0] = i1;
230  weights[0] += i1Tot1.first();
231  indices[1] = i1+1;
232  weights[1] += i1Tot1.second();
233  }
234 
235  // Subtract from t2 to i2+1
236  {
237  Pair<scalar> wghts = integrationWeights(i2, t2);
238  indices[i2-i1] = i2;
239  weights[i2-i1] += -wghts.first();
240  indices[i2-i1+1] = i2+1;
241  weights[i2-i1+1] += -wghts.second();
242  }
243 
244  return anyChanged;
245 }
246 
247 
248 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
249 
250 } // End namespace Foam
251 
252 // ************************************************************************* //
virtual bool valueWeights(const scalar t, labelList &indices, scalarField &weights) const
Calculate weights and indices to calculate t from samples.
An ordered pair of two objects of type <T> with first() and second() elements.
Definition: contiguous.H:49
gmvFile<< "tracers "<< particles.size()<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){gmvFile<< iter().position().x()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){gmvFile<< iter().position().y()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){gmvFile<< iter().position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject( name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE ))
T & last()
Return the last element of the list.
Definition: UListI.H:131
const Type & second() const
Return second.
Definition: Pair.H:99
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:76
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
Various functions to operate on Lists.
Namespace for OpenFOAM.
void setSize(const label)
Reset size of List.
Definition: List.C:318
linearInterpolationWeights(const scalarField &samples)
Construct from components.
#define forAll(list, i)
Definition: UList.H:421
Macros for easy insertion into run-time selection tables.
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,.
#define FatalErrorIn(functionName)
Report an error message using Foam::FatalError.
Definition: error.H:314
error FatalError
const Type & first() const
Return first.
Definition: Pair.H:87
Abstract base class for interpolating in 1D.
addToRunTimeSelectionTable(ensightPart, ensightPartCells, istream)
defineTypeNameAndDebug(combustionModel, 0)