timeVaryingMappedFixedValuePointPatchField.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-2014 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 
27 #include "Time.H"
28 #include "AverageIOField.H"
29 
30 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
31 
32 template<class Type>
33 Foam::
36 (
37  const pointPatch& p,
39 )
40 :
42  fieldTableName_(iF.name()),
43  setAverage_(false),
44  perturb_(0),
45  mapperPtr_(NULL),
46  sampleTimes_(0),
47  startSampleTime_(-1),
48  startSampledValues_(0),
49  startAverage_(pTraits<Type>::zero),
50  endSampleTime_(-1),
51  endSampledValues_(0),
52  endAverage_(pTraits<Type>::zero),
53  offset_()
54 {}
55 
56 
57 template<class Type>
58 Foam::
61 (
63  const pointPatch& p,
65  const pointPatchFieldMapper& mapper
66 )
67 :
68  fixedValuePointPatchField<Type>(ptf, p, iF, mapper),
69  fieldTableName_(ptf.fieldTableName_),
70  setAverage_(ptf.setAverage_),
71  perturb_(ptf.perturb_),
72  mapMethod_(ptf.mapMethod_),
73  mapperPtr_(NULL),
74  sampleTimes_(0),
75  startSampleTime_(-1),
76  startSampledValues_(0),
77  startAverage_(pTraits<Type>::zero),
78  endSampleTime_(-1),
79  endSampledValues_(0),
80  endAverage_(pTraits<Type>::zero),
81  offset_
82  (
83  ptf.offset_.valid()
84  ? ptf.offset_().clone().ptr()
85  : NULL
86  )
87 {}
88 
89 
90 template<class Type>
91 Foam::
94 (
95  const pointPatch& p,
97  const dictionary& dict
98 )
99 :
101  fieldTableName_(iF.name()),
102  setAverage_(readBool(dict.lookup("setAverage"))),
103  perturb_(dict.lookupOrDefault("perturb", 1e-5)),
104  mapMethod_
105  (
106  dict.lookupOrDefault<word>
107  (
108  "mapMethod",
109  "planarInterpolation"
110  )
111  ),
112  mapperPtr_(NULL),
113  sampleTimes_(0),
114  startSampleTime_(-1),
115  startSampledValues_(0),
116  startAverage_(pTraits<Type>::zero),
117  endSampleTime_(-1),
118  endSampledValues_(0),
119  endAverage_(pTraits<Type>::zero),
120  offset_()
121 {
122  if (dict.found("offset"))
123  {
124  offset_ = DataEntry<Type>::New("offset", dict);
125  }
126 
127  dict.readIfPresent("fieldTableName", fieldTableName_);
128 
129  if (dict.found("value"))
130  {
132  (
133  Field<Type>("value", dict, p.size())
134  );
135  }
136  else
137  {
138  // Note: use evaluate to do updateCoeffs followed by a reset
139  // of the pointPatchField::updated_ flag. This is
140  // so if first use is in the next time step it retriggers
141  // a new update.
142  pointPatchField<Type>::evaluate(Pstream::blocking);
143  }
144 }
145 
146 
147 template<class Type>
148 Foam::
151 (
153 )
154 :
156  fieldTableName_(ptf.fieldTableName_),
157  setAverage_(ptf.setAverage_),
158  perturb_(ptf.perturb_),
159  mapMethod_(ptf.mapMethod_),
160  mapperPtr_(ptf.mapperPtr_),
161  sampleTimes_(ptf.sampleTimes_),
162  startSampleTime_(ptf.startSampleTime_),
163  startSampledValues_(ptf.startSampledValues_),
164  startAverage_(ptf.startAverage_),
165  endSampleTime_(ptf.endSampleTime_),
166  endSampledValues_(ptf.endSampledValues_),
167  endAverage_(ptf.endAverage_),
168  offset_
169  (
170  ptf.offset_.valid()
171  ? ptf.offset_().clone().ptr()
172  : NULL
173  )
174 {}
175 
176 
177 template<class Type>
178 Foam::
181 (
184 )
185 :
187  fieldTableName_(ptf.fieldTableName_),
188  setAverage_(ptf.setAverage_),
189  perturb_(ptf.perturb_),
190  mapMethod_(ptf.mapMethod_),
191  mapperPtr_(ptf.mapperPtr_),
192  sampleTimes_(ptf.sampleTimes_),
193  startSampleTime_(ptf.startSampleTime_),
194  startSampledValues_(ptf.startSampledValues_),
195  startAverage_(ptf.startAverage_),
196  endSampleTime_(ptf.endSampleTime_),
197  endSampledValues_(ptf.endSampledValues_),
198  endAverage_(ptf.endAverage_),
199  offset_
200  (
201  ptf.offset_.valid()
202  ? ptf.offset_().clone().ptr()
203  : NULL
204  )
205 {}
206 
207 
208 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
209 
210 template<class Type>
212 (
213  const pointPatchFieldMapper& m
214 )
215 {
217  if (startSampledValues_.size())
218  {
219  startSampledValues_.autoMap(m);
220  endSampledValues_.autoMap(m);
221  }
222  // Clear interpolator
223  mapperPtr_.clear();
224  startSampleTime_ = -1;
225  endSampleTime_ = -1;
226 }
227 
228 
229 template<class Type>
231 (
232  const pointPatchField<Type>& ptf,
233  const labelList& addr
234 )
235 {
237 
239  refCast<const timeVaryingMappedFixedValuePointPatchField<Type> >(ptf);
240 
241  startSampledValues_.rmap(tiptf.startSampledValues_, addr);
242  endSampledValues_.rmap(tiptf.endSampledValues_, addr);
243 
244  // Clear interpolator
245  mapperPtr_.clear();
246  startSampleTime_ = -1;
247  endSampleTime_ = -1;
248 }
249 
250 
251 template<class Type>
253 {
254  // Initialise
255  if (startSampleTime_ == -1 && endSampleTime_ == -1)
256  {
257  const polyMesh& pMesh = this->patch().boundaryMesh().mesh()();
258 
259  // Read the initial point position
260  pointField meshPts;
261 
262  if (pMesh.pointsInstance() == pMesh.facesInstance())
263  {
264  meshPts = pointField(pMesh.points(), this->patch().meshPoints());
265  }
266  else
267  {
268  // Load points from facesInstance
269  if (debug)
270  {
271  Info<< "Reloading points0 from " << pMesh.facesInstance()
272  << endl;
273  }
274 
275  pointIOField points0
276  (
277  IOobject
278  (
279  "points",
280  pMesh.facesInstance(),
281  polyMesh::meshSubDir,
282  pMesh,
283  IOobject::MUST_READ,
284  IOobject::NO_WRITE,
285  false
286  )
287  );
288  meshPts = pointField(points0, this->patch().meshPoints());
289  }
290 
291  pointIOField samplePoints
292  (
293  IOobject
294  (
295  "points",
296  this->db().time().constant(),
297  "boundaryData"/this->patch().name(),
298  this->db(),
299  IOobject::MUST_READ,
300  IOobject::AUTO_WRITE,
301  false
302  )
303  );
304 
305  // tbd: run-time selection
306  bool nearestOnly =
307  (
308  !mapMethod_.empty()
309  && mapMethod_ != "planarInterpolation"
310  );
311 
312  // Allocate the interpolator
313  mapperPtr_.reset
314  (
316  (
317  samplePoints,
318  meshPts,
319  perturb_,
320  nearestOnly
321  )
322  );
323 
324  // Read the times for which data is available
325 
326  const fileName samplePointsFile = samplePoints.filePath();
327  const fileName samplePointsDir = samplePointsFile.path();
328  sampleTimes_ = Time::findTimes(samplePointsDir);
329 
330  if (debug)
331  {
332  Info<< "timeVaryingMappedFixedValuePointPatchField : In directory "
333  << samplePointsDir << " found times "
334  << pointToPointPlanarInterpolation::timeNames(sampleTimes_)
335  << endl;
336  }
337  }
338 
339  // Find current time in sampleTimes
340  label lo = -1;
341  label hi = -1;
342 
343  bool foundTime = mapperPtr_().findTime
344  (
345  sampleTimes_,
346  startSampleTime_,
347  this->db().time().value(),
348  lo,
349  hi
350  );
351 
352  if (!foundTime)
353  {
355  (
356  "timeVaryingMappedFixedValuePointPatchField<Type>::checkTable"
357  ) << "Cannot find starting sampling values for current time "
358  << this->db().time().value() << nl
359  << "Have sampling values for times "
360  << pointToPointPlanarInterpolation::timeNames(sampleTimes_) << nl
361  << "In directory "
362  << this->db().time().constant()/"boundaryData"/this->patch().name()
363  << "\n on patch " << this->patch().name()
364  << " of field " << fieldTableName_
365  << exit(FatalError);
366  }
367 
368 
369  // Update sampled data fields.
370 
371  if (lo != startSampleTime_)
372  {
373  startSampleTime_ = lo;
374 
375  if (startSampleTime_ == endSampleTime_)
376  {
377  // No need to reread since are end values
378  if (debug)
379  {
380  Pout<< "checkTable : Setting startValues to (already read) "
381  << "boundaryData"
382  /this->patch().name()
383  /sampleTimes_[startSampleTime_].name()
384  << endl;
385  }
386  startSampledValues_ = endSampledValues_;
387  startAverage_ = endAverage_;
388  }
389  else
390  {
391  if (debug)
392  {
393  Pout<< "checkTable : Reading startValues from "
394  << "boundaryData"
395  /this->patch().name()
396  /sampleTimes_[lo].name()
397  << endl;
398  }
399 
400  // Reread values and interpolate
402  (
403  IOobject
404  (
405  fieldTableName_,
406  this->db().time().constant(),
407  "boundaryData"
408  /this->patch().name()
409  /sampleTimes_[startSampleTime_].name(),
410  this->db(),
411  IOobject::MUST_READ,
412  IOobject::AUTO_WRITE,
413  false
414  )
415  );
416 
417  if (vals.size() != mapperPtr_().sourceSize())
418  {
420  (
421  "timeVaryingMappedFixedValuePointPatchField<Type>::"
422  "checkTable()"
423  ) << "Number of values (" << vals.size()
424  << ") differs from the number of points ("
425  << mapperPtr_().sourceSize()
426  << ") in file " << vals.objectPath() << exit(FatalError);
427  }
428 
429  startAverage_ = vals.average();
430  startSampledValues_ = mapperPtr_().interpolate(vals);
431  }
432  }
433 
434  if (hi != endSampleTime_)
435  {
436  endSampleTime_ = hi;
437 
438  if (endSampleTime_ == -1)
439  {
440  // endTime no longer valid. Might as well clear endValues.
441  if (debug)
442  {
443  Pout<< "checkTable : Clearing endValues" << endl;
444  }
445  endSampledValues_.clear();
446  }
447  else
448  {
449  if (debug)
450  {
451  Pout<< "checkTable : Reading endValues from "
452  << "boundaryData"
453  /this->patch().name()
454  /sampleTimes_[endSampleTime_].name()
455  << endl;
456  }
457  // Reread values and interpolate
459  (
460  IOobject
461  (
462  fieldTableName_,
463  this->db().time().constant(),
464  "boundaryData"
465  /this->patch().name()
466  /sampleTimes_[endSampleTime_].name(),
467  this->db(),
468  IOobject::MUST_READ,
469  IOobject::AUTO_WRITE,
470  false
471  )
472  );
473 
474  if (vals.size() != mapperPtr_().sourceSize())
475  {
477  (
478  "timeVaryingMappedFixedValuePointPatchField<Type>::"
479  "checkTable()"
480  ) << "Number of values (" << vals.size()
481  << ") differs from the number of points ("
482  << mapperPtr_().sourceSize()
483  << ") in file " << vals.objectPath() << exit(FatalError);
484  }
485 
486  endAverage_ = vals.average();
487  endSampledValues_ = mapperPtr_().interpolate(vals);
488  }
489  }
490 }
491 
492 
493 template<class Type>
495 {
496  if (this->updated())
497  {
498  return;
499  }
500 
501  checkTable();
502 
503  // Interpolate between the sampled data
504 
505  Type wantedAverage;
506 
507  if (endSampleTime_ == -1)
508  {
509  // only start value
510  if (debug)
511  {
512  Pout<< "updateCoeffs : Sampled, non-interpolated values"
513  << " from start time:"
514  << sampleTimes_[startSampleTime_].name() << nl;
515  }
516 
517  this->operator==(startSampledValues_);
518  wantedAverage = startAverage_;
519  }
520  else
521  {
522  scalar start = sampleTimes_[startSampleTime_].value();
523  scalar end = sampleTimes_[endSampleTime_].value();
524 
525  scalar s = (this->db().time().value()-start)/(end-start);
526 
527  if (debug)
528  {
529  Pout<< "updateCoeffs : Sampled, interpolated values"
530  << " between start time:"
531  << sampleTimes_[startSampleTime_].name()
532  << " and end time:" << sampleTimes_[endSampleTime_].name()
533  << " with weight:" << s << endl;
534  }
535 
536  this->operator==((1-s)*startSampledValues_ + s*endSampledValues_);
537  wantedAverage = (1-s)*startAverage_ + s*endAverage_;
538  }
539 
540  // Enforce average. Either by scaling (if scaling factor > 0.5) or by
541  // offsetting.
542  if (setAverage_)
543  {
544  const Field<Type>& fld = *this;
545 
546  Type averagePsi = gAverage(fld);
547 
548  if (debug)
549  {
550  Pout<< "updateCoeffs :"
551  << " actual average:" << averagePsi
552  << " wanted average:" << wantedAverage
553  << endl;
554  }
555 
556  if (mag(averagePsi) < VSMALL)
557  {
558  // Field too small to scale. Offset instead.
559  const Type offset = wantedAverage - averagePsi;
560  if (debug)
561  {
562  Pout<< "updateCoeffs :"
563  << " offsetting with:" << offset << endl;
564  }
565  this->operator==(fld+offset);
566  }
567  else
568  {
569  const scalar scale = mag(wantedAverage)/mag(averagePsi);
570 
571  if (debug)
572  {
573  Pout<< "updateCoeffs :"
574  << " scaling with:" << scale << endl;
575  }
576  this->operator==(scale*fld);
577  }
578  }
579 
580  // apply offset to mapped values
581  if (offset_.valid())
582  {
583  const scalar t = this->db().time().timeOutputValue();
584  this->operator==(*this + offset_->value(t));
585  }
586 
587  if (debug)
588  {
589  Pout<< "updateCoeffs : set fixedValue to min:" << gMin(*this)
590  << " max:" << gMax(*this)
591  << " avg:" << gAverage(*this) << endl;
592  }
593 
595 }
596 
597 
598 template<class Type>
600 (
601  Ostream& os
602 ) const
603 {
605  os.writeKeyword("setAverage") << setAverage_ << token::END_STATEMENT << nl;
606  if (perturb_ != 1e-5)
607  {
608  os.writeKeyword("perturb") << perturb_ << token::END_STATEMENT << nl;
609  }
610 
611  if (fieldTableName_ != this->dimensionedInternalField().name())
612  {
613  os.writeKeyword("fieldTableName") << fieldTableName_
614  << token::END_STATEMENT << nl;
615  }
616 
617  if
618  (
619  (
620  !mapMethod_.empty()
621  && mapMethod_ != "planarInterpolation"
622  )
623  )
624  {
625  os.writeKeyword("mapMethod") << mapMethod_
626  << token::END_STATEMENT << nl;
627  }
628 
629  if (offset_.valid())
630  {
631  offset_->writeData(os);
632  }
633 }
634 
635 
636 // ************************************************************************* //
Interpolates between two sets of unstructured points using 2D Delaunay triangulation. Used in e.g. timeVaryingMapped bcs.
bool found(const word &, bool recursive=false, bool patternMatch=true) const
Search dictionary for given keyword.
Definition: dictionary.C:306
Basic pointPatch represents a set of points from the mesh.
Definition: pointPatch.H:56
bool readBool(Istream &)
Definition: boolIO.C:60
word name() const
Return file name (part beyond last /)
Definition: fileName.C:206
tmp< fvMatrix< Type > > operator==(const fvMatrix< Type > &, const fvMatrix< Type > &)
Field with dimensions and associated with geometry type GeoMesh which is used to size the field and a...
fileName filePath() const
Return complete path + object name if the file exists.
Definition: IOobject.C:323
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 ))
virtual void updateCoeffs()
Update the coefficients associated with the patch field.
dimensioned< scalar > mag(const dimensioned< Type > &)
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
virtual void autoMap(const pointPatchFieldMapper &)
Map (and resize as needed) from self given a mapping object.
A primitive field + average with IO.
virtual autoPtr< pointPatchField< Type > > clone() const
Construct and return a clone.
A class for handling words, derived from string.
Definition: word.H:59
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.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
messageStream Info
fileName path() const
Return directory path name (part before last /)
Definition: fileName.C:287
const fileName & pointsInstance() const
Return the current instance directory for points.
Definition: polyMesh.C:818
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:137
void checkTable()
Find boundary data inbetween current time and interpolate.
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
autoPtr< BasicCompressibleTurbulenceModel > New(const volScalarField &rho, const volVectorField &U, const surfaceScalarField &phi, const typename BasicCompressibleTurbulenceModel::transportModel &transport, const word &propertiesName)
runTime write()
const Type & average() const
dictionary dict
Abstract base class for point-mesh patch fields.
static const char nl
Definition: Ostream.H:260
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:91
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:1035
volScalarField & p
Definition: createFields.H:51
A time-varying form of a mapped fixed value boundary condition.
const dimensionedScalar e
Elementary charge.
Definition: doubleFloat.H:78
virtual void rmap(const pointPatchField< Type > &, const labelList &)
Reverse map the given PointPatchField onto.
fileName objectPath() const
Return complete path + object name.
Definition: IOobject.H:363
Pre-declare SubField and related Field type.
Definition: Field.H:57
const word & name() const
Return name.
Definition: IOobject.H:260
Type gMin(const FieldField< Field, Type > &f)
virtual const fileName & name() const
Return the name of the stream.
Definition: OSstream.H:88
Ostream & writeKeyword(const keyType &)
Write the keyword followed by an appropriate indentation.
Definition: Ostream.C:59
Foam::pointPatchFieldMapper.
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){const word &name=lagrangianScalarNames[i];IOField< scalar > fld(IOobject( name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE ))
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:421
#define FatalErrorIn(functionName)
Report an error message using Foam::FatalError.
Definition: error.H:314
rDeltaT dimensionedInternalField()
const fileName & facesInstance() const
Return the current instance directory for faces.
Definition: polyMesh.C:824
timeVaryingMappedFixedValuePointPatchField(const pointPatch &, const DimensionedField< Type, pointMesh > &)
Construct from patch and internal field.
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
Traits class for primitives.
Definition: pTraits.H:50
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
Definition: dictionary.C:452
error FatalError
Type gAverage(const FieldField< Field, Type > &f)
Constant dispersed-phase particle diameter model.
T lookupOrDefault(const word &, const T &, bool recursive=false, bool patternMatch=true) const
Find and return a T,.
const polyMesh & mesh() const
Return the mesh reference.
A class for handling file names.
Definition: fileName.H:69
bool readIfPresent(const word &, T &, bool recursive=false, bool patternMatch=true) const
Find an entry if present, and assign to T.
A FixedValue boundary condition for pointField.
An Ostream is an abstract base class for all output systems (streams, files, token lists...
Definition: Ostream.H:53
Type gMax(const FieldField< Field, Type > &f)
label size() const
Return size.
prefixOSstream Pout(cout,"Pout")
Definition: IOstreams.H:53