FreeFOAM The Cross-Platform CFD Toolkit
timeVaryingMappedFixedValueFvPatchField.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) 1991-2011 OpenCFD Ltd.
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 <OpenFOAM/Time.H>
29 #include <triSurface/triSurface.H>
30 #include <OpenFOAM/vector2D.H>
31 #include <OpenFOAM/OFstream.H>
33 #include <OpenFOAM/Random.H>
34 
35 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
36 
37 namespace Foam
38 {
39 
40 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
41 
42 template<class Type>
45 (
46  const fvPatch& p,
48 )
49 :
51  fieldTableName_(iF.name()),
52  setAverage_(false),
53  perturb_(0),
54  referenceCS_(NULL),
55  nearestVertex_(0),
56  nearestVertexWeight_(0),
57  sampleTimes_(0),
58  startSampleTime_(-1),
59  startSampledValues_(0),
60  startAverage_(pTraits<Type>::zero),
61  endSampleTime_(-1),
62  endSampledValues_(0),
63  endAverage_(pTraits<Type>::zero)
64 {}
65 
66 
67 template<class Type>
70 (
72  const fvPatch& p,
74  const fvPatchFieldMapper& mapper
75 )
76 :
77  fixedValueFvPatchField<Type>(ptf, p, iF, mapper),
78  fieldTableName_(ptf.fieldTableName_),
79  setAverage_(ptf.setAverage_),
80  perturb_(ptf.perturb_),
81  referenceCS_(NULL),
82  nearestVertex_(0),
83  nearestVertexWeight_(0),
84  sampleTimes_(0),
85  startSampleTime_(-1),
86  startSampledValues_(0),
87  startAverage_(pTraits<Type>::zero),
88  endSampleTime_(-1),
89  endSampledValues_(0),
90  endAverage_(pTraits<Type>::zero)
91 {}
92 
93 
94 template<class Type>
97 (
98  const fvPatch& p,
100  const dictionary& dict
101 )
102 :
104  fieldTableName_(iF.name()),
105  setAverage_(readBool(dict.lookup("setAverage"))),
106  perturb_(dict.lookupOrDefault("perturb", 1E-5)),
107  referenceCS_(NULL),
108  nearestVertex_(0),
109  nearestVertexWeight_(0),
110  sampleTimes_(0),
111  startSampleTime_(-1),
112  startSampledValues_(0),
113  startAverage_(pTraits<Type>::zero),
114  endSampleTime_(-1),
115  endSampledValues_(0),
116  endAverage_(pTraits<Type>::zero)
117 {
118  if (dict.found("fieldTableName"))
119  {
120  dict.lookup("fieldTableName") >> fieldTableName_;
121  }
122 
123  if (dict.found("value"))
124  {
125  fvPatchField<Type>::operator==(Field<Type>("value", dict, p.size()));
126  }
127  else
128  {
129  updateCoeffs();
130  }
131 }
132 
133 
134 template<class Type>
137 (
139 )
140 :
142  fieldTableName_(ptf.fieldTableName_),
143  setAverage_(ptf.setAverage_),
144  perturb_(ptf.perturb_),
145  referenceCS_(ptf.referenceCS_),
146  nearestVertex_(ptf.nearestVertex_),
147  nearestVertexWeight_(ptf.nearestVertexWeight_),
148  sampleTimes_(ptf.sampleTimes_),
149  startSampleTime_(ptf.startSampleTime_),
150  startSampledValues_(ptf.startSampledValues_),
151  startAverage_(ptf.startAverage_),
152  endSampleTime_(ptf.endSampleTime_),
153  endSampledValues_(ptf.endSampledValues_),
154  endAverage_(ptf.endAverage_)
155 {}
156 
157 
158 
159 template<class Type>
162 (
165 )
166 :
168  fieldTableName_(ptf.fieldTableName_),
169  setAverage_(ptf.setAverage_),
170  perturb_(ptf.perturb_),
171  referenceCS_(ptf.referenceCS_),
172  nearestVertex_(ptf.nearestVertex_),
173  nearestVertexWeight_(ptf.nearestVertexWeight_),
174  sampleTimes_(ptf.sampleTimes_),
175  startSampleTime_(ptf.startSampleTime_),
176  startSampledValues_(ptf.startSampledValues_),
177  startAverage_(ptf.startAverage_),
178  endSampleTime_(ptf.endSampleTime_),
179  endSampledValues_(ptf.endSampledValues_),
180  endAverage_(ptf.endAverage_)
181 {}
182 
183 
184 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
185 
186 template<class Type>
188 (
189  const fvPatchFieldMapper& m
190 )
191 {
193  if (startSampledValues_.size())
194  {
195  startSampledValues_.autoMap(m);
196  endSampledValues_.autoMap(m);
197  }
198 }
199 
200 
201 template<class Type>
203 (
204  const fvPatchField<Type>& ptf,
205  const labelList& addr
206 )
207 {
209 
211  refCast<const timeVaryingMappedFixedValueFvPatchField<Type> >(ptf);
212 
213  startSampledValues_.rmap(tiptf.startSampledValues_, addr);
214  endSampledValues_.rmap(tiptf.endSampledValues_, addr);
215 }
216 
217 
218 template<class Type>
220 {
221  // Read the sample points
222 
223  pointIOField samplePoints
224  (
225  IOobject
226  (
227  "points",
228  this->db().time().constant(),
229  "boundaryData"/this->patch().name(),
230  this->db(),
233  false
234  )
235  );
236 
237  const fileName samplePointsFile = samplePoints.filePath();
238 
239  if (debug)
240  {
241  Info<< "timeVaryingMappedFixedValueFvPatchField :"
242  << " Read " << samplePoints.size() << " sample points from "
243  << samplePointsFile << endl;
244  }
245 
246  // Determine coordinate system from samplePoints
247 
248  if (samplePoints.size() < 3)
249  {
251  (
252  "timeVaryingMappedFixedValueFvPatchField<Type>::readSamplePoints()"
253  ) << "Only " << samplePoints.size() << " points read from file "
254  << samplePoints.objectPath() << nl
255  << "Need at least three non-colinear samplePoints"
256  << " to be able to interpolate."
257  << "\n on patch " << this->patch().name()
258  << " of points " << samplePoints.name()
259  << " in file " << samplePoints.objectPath()
260  << exit(FatalError);
261  }
262 
263  const point& p0 = samplePoints[0];
264 
265  // Find furthest away point
266  vector e1;
267  label index1 = -1;
268  scalar maxDist = -GREAT;
269 
270  for (label i = 1; i < samplePoints.size(); i++)
271  {
272  const vector d = samplePoints[i] - p0;
273  scalar magD = mag(d);
274 
275  if (magD > maxDist)
276  {
277  e1 = d/magD;
278  index1 = i;
279  maxDist = magD;
280  }
281  }
282 
283  // Find point that is furthest away from line p0-p1
284  const point& p1 = samplePoints[index1];
285 
286  label index2 = -1;
287  maxDist = -GREAT;
288  for (label i = 1; i < samplePoints.size(); i++)
289  {
290  if (i != index1)
291  {
292  const point& p2 = samplePoints[i];
293  vector e2(p2 - p0);
294  e2 -= (e2&e1)*e1;
295  scalar magE2 = mag(e2);
296 
297  if (magE2 > maxDist)
298  {
299  index2 = i;
300  maxDist = magE2;
301  }
302  }
303  }
304 
305  if (index2 == -1)
306  {
308  (
309  "timeVaryingMappedFixedValueFvPatchField<Type>::readSamplePoints()"
310  ) << "Cannot find points that make valid normal." << nl
311  << "Have so far points " << p0 << " and " << p1
312  << "Need at least three sample points which are not in a line."
313  << "\n on patch " << this->patch().name()
314  << " of points " << samplePoints.name()
315  << " in file " << samplePoints.objectPath()
316  << exit(FatalError);
317  }
318 
319  vector n = e1^(samplePoints[index2]-p0);
320 
321  n /= mag(n);
322 
323  if (debug)
324  {
325  Info<< "timeVaryingMappedFixedValueFvPatchField :"
326  << " Used points " << p0 << ' ' << samplePoints[index1]
327  << ' ' << samplePoints[index2]
328  << " to define coordinate system with normal " << n << endl;
329  }
330 
331  referenceCS_.reset
332  (
333  new coordinateSystem
334  (
335  "reference",
336  p0, // origin
337  n, // normal
338  e1 // 0-axis
339  )
340  );
341 
342  tmp<vectorField> tlocalVertices
343  (
344  referenceCS().localPosition(samplePoints)
345  );
346  vectorField& localVertices = tlocalVertices();
347 
348  const boundBox bb(localVertices, true);
349  const point bbMid(bb.midpoint());
350 
351  if (debug)
352  {
353  Info<< "timeVaryingMappedFixedValueFvPatchField :"
354  << " Perturbing points with " << perturb_
355  << " fraction of a random position inside " << bb
356  << " to break any ties on regular meshes."
357  << nl << endl;
358  }
359 
360  Random rndGen(123456);
361  forAll(localVertices, i)
362  {
363  localVertices[i] +=
364  perturb_
365  *(rndGen.position(bb.min(), bb.max())-bbMid);
366  }
367 
368  // Determine triangulation
369  List<vector2D> localVertices2D(localVertices.size());
370  forAll(localVertices, i)
371  {
372  localVertices2D[i][0] = localVertices[i][0];
373  localVertices2D[i][1] = localVertices[i][1];
374  }
375 
376  triSurface s(triSurfaceTools::delaunay2D(localVertices2D));
377 
378  tmp<pointField> tlocalFaceCentres
379  (
380  referenceCS().localPosition
381  (
382  this->patch().patch().faceCentres()
383  )
384  );
385  const pointField& localFaceCentres = tlocalFaceCentres();
386 
387  if (debug)
388  {
389  Pout<< "readSamplePoints :"
390  <<" Dumping triangulated surface to triangulation.stl" << endl;
391  s.write(this->db().time().path()/"triangulation.stl");
392 
393  OFstream str(this->db().time().path()/"localFaceCentres.obj");
394  Pout<< "readSamplePoints :"
395  << " Dumping face centres to " << str.name() << endl;
396 
397  forAll(localFaceCentres, i)
398  {
399  const point& p = localFaceCentres[i];
400  str<< "v " << p.x() << ' ' << p.y() << ' ' << p.z() << nl;
401  }
402  }
403 
404  // Determine interpolation onto face centres.
406  (
407  s,
408  localFaceCentres, // points to interpolate to
409  nearestVertex_,
410  nearestVertexWeight_
411  );
412 
413 
414 
415  // Read the times for which data is available
416 
417  const fileName samplePointsDir = samplePointsFile.path();
418 
419  sampleTimes_ = Time::findTimes(samplePointsDir);
420 
421  if (debug)
422  {
423  Info<< "timeVaryingMappedFixedValueFvPatchField : In directory "
424  << samplePointsDir << " found times " << timeNames(sampleTimes_)
425  << endl;
426  }
427 }
428 
429 
430 template<class Type>
431 wordList timeVaryingMappedFixedValueFvPatchField<Type>::timeNames
432 (
433  const instantList& times
434 )
435 {
436  wordList names(times.size());
437 
438  forAll(times, i)
439  {
440  names[i] = times[i].name();
441  }
442  return names;
443 }
444 
445 
446 template<class Type>
447 void timeVaryingMappedFixedValueFvPatchField<Type>::findTime
448 (
449  const fileName& instance,
450  const fileName& local,
451  const scalar timeVal,
452  label& lo,
453  label& hi
454 ) const
455 {
456  lo = startSampleTime_;
457  hi = -1;
458 
459  for (label i = startSampleTime_+1; i < sampleTimes_.size(); i++)
460  {
461  if (sampleTimes_[i].value() > timeVal)
462  {
463  break;
464  }
465  else
466  {
467  lo = i;
468  }
469  }
470 
471  if (lo == -1)
472  {
473  FatalErrorIn("findTime")
474  << "Cannot find starting sampling values for current time "
475  << timeVal << nl
476  << "Have sampling values for times "
477  << timeNames(sampleTimes_) << nl
478  << "In directory "
479  << this->db().time().constant()/"boundaryData"/this->patch().name()
480  << "\n on patch " << this->patch().name()
481  << " of field " << fieldTableName_
482  << exit(FatalError);
483  }
484 
485  if (lo < sampleTimes_.size()-1)
486  {
487  hi = lo+1;
488  }
489 
490 
491  if (debug)
492  {
493  if (hi == -1)
494  {
495  Pout<< "findTime : Found time " << timeVal << " after"
496  << " index:" << lo << " time:" << sampleTimes_[lo].value()
497  << endl;
498  }
499  else
500  {
501  Pout<< "findTime : Found time " << timeVal << " inbetween"
502  << " index:" << lo << " time:" << sampleTimes_[lo].value()
503  << " and index:" << hi << " time:" << sampleTimes_[hi].value()
504  << endl;
505  }
506  }
507 }
508 
509 
510 template<class Type>
512 {
513  // Initialise
514  if (startSampleTime_ == -1 && endSampleTime_ == -1)
515  {
516  readSamplePoints();
517  }
518 
519  // Find current time in sampleTimes
520  label lo = -1;
521  label hi = -1;
522 
523  findTime
524  (
525  this->db().time().constant(),
526  "boundaryData"/this->patch().name(),
527  this->db().time().value(),
528  lo,
529  hi
530  );
531 
532  // Update sampled data fields.
533 
534  if (lo != startSampleTime_)
535  {
536  startSampleTime_ = lo;
537 
538  if (startSampleTime_ == endSampleTime_)
539  {
540  // No need to reread since are end values
541  if (debug)
542  {
543  Pout<< "checkTable : Setting startValues to (already read) "
544  << "boundaryData"
545  /this->patch().name()
546  /sampleTimes_[startSampleTime_].name()
547  << endl;
548  }
549  startSampledValues_ = endSampledValues_;
550  startAverage_ = endAverage_;
551  }
552  else
553  {
554  if (debug)
555  {
556  Pout<< "checkTable : Reading startValues from "
557  << "boundaryData"
558  /this->patch().name()
559  /sampleTimes_[lo].name()
560  << endl;
561  }
562 
563 
564  // Reread values and interpolate
566  (
567  IOobject
568  (
569  fieldTableName_,
570  this->db().time().constant(),
571  "boundaryData"
572  /this->patch().name()
573  /sampleTimes_[startSampleTime_].name(),
574  this->db(),
577  false
578  )
579  );
580 
581  startAverage_ = vals.average();
582  startSampledValues_ = interpolate(vals);
583  }
584  }
585 
586  if (hi != endSampleTime_)
587  {
588  endSampleTime_ = hi;
589 
590  if (endSampleTime_ == -1)
591  {
592  // endTime no longer valid. Might as well clear endValues.
593  if (debug)
594  {
595  Pout<< "checkTable : Clearing endValues" << endl;
596  }
597  endSampledValues_.clear();
598  }
599  else
600  {
601  if (debug)
602  {
603  Pout<< "checkTable : Reading endValues from "
604  << "boundaryData"
605  /this->patch().name()
606  /sampleTimes_[endSampleTime_].name()
607  << endl;
608  }
609  // Reread values and interpolate
611  (
612  IOobject
613  (
614  fieldTableName_,
615  this->db().time().constant(),
616  "boundaryData"
617  /this->patch().name()
618  /sampleTimes_[endSampleTime_].name(),
619  this->db(),
622  false
623  )
624  );
625  endAverage_ = vals.average();
626  endSampledValues_ = interpolate(vals);
627  }
628  }
629 }
630 
631 
632 template<class Type>
634 (
635  const Field<Type>& sourceFld
636 ) const
637 {
638  tmp<Field<Type> > tfld(new Field<Type>(nearestVertex_.size()));
639  Field<Type>& fld = tfld();
640 
641  forAll(fld, i)
642  {
643  const FixedList<label, 3>& verts = nearestVertex_[i];
644  const FixedList<scalar, 3>& w = nearestVertexWeight_[i];
645 
646  if (verts[2] == -1)
647  {
648  if (verts[1] == -1)
649  {
650  // Use vertex0 only
651  fld[i] = sourceFld[verts[0]];
652  }
653  else
654  {
655  // Use vertex 0,1
656  fld[i] =
657  w[0]*sourceFld[verts[0]]
658  + w[1]*sourceFld[verts[1]];
659  }
660  }
661  else
662  {
663  fld[i] =
664  w[0]*sourceFld[verts[0]]
665  + w[1]*sourceFld[verts[1]]
666  + w[2]*sourceFld[verts[2]];
667  }
668  }
669  return tfld;
670 }
671 
672 
673 template<class Type>
675 {
676  if (this->updated())
677  {
678  return;
679  }
680 
681  checkTable();
682 
683  // Interpolate between the sampled data
684 
685  Type wantedAverage;
686 
687  if (endSampleTime_ == -1)
688  {
689  // only start value
690  if (debug)
691  {
692  Pout<< "updateCoeffs : Sampled, non-interpolated values"
693  << " from start time:"
694  << sampleTimes_[startSampleTime_].name() << nl;
695  }
696 
697  this->operator==(startSampledValues_);
698  wantedAverage = startAverage_;
699  }
700  else
701  {
702  scalar start = sampleTimes_[startSampleTime_].value();
703  scalar end = sampleTimes_[endSampleTime_].value();
704 
705  scalar s = (this->db().time().value()-start)/(end-start);
706 
707  if (debug)
708  {
709  Pout<< "updateCoeffs : Sampled, interpolated values"
710  << " between start time:"
711  << sampleTimes_[startSampleTime_].name()
712  << " and end time:" << sampleTimes_[endSampleTime_].name()
713  << " with weight:" << s << endl;
714  }
715 
716  this->operator==((1-s)*startSampledValues_ + s*endSampledValues_);
717  wantedAverage = (1-s)*startAverage_ + s*endAverage_;
718  }
719 
720  // Enforce average. Either by scaling (if scaling factor > 0.5) or by
721  // offsetting.
722  if (setAverage_)
723  {
724  const Field<Type>& fld = *this;
725 
726  Type averagePsi =
727  gSum(this->patch().magSf()*fld)
728  /gSum(this->patch().magSf());
729 
730  if (debug)
731  {
732  Pout<< "updateCoeffs :"
733  << " actual average:" << averagePsi
734  << " wanted average:" << wantedAverage
735  << endl;
736  }
737 
738  if (mag(averagePsi) < VSMALL)
739  {
740  // Field too small to scale. Offset instead.
741  const Type offset = wantedAverage - averagePsi;
742  if (debug)
743  {
744  Pout<< "updateCoeffs :"
745  << " offsetting with:" << offset << endl;
746  }
747  this->operator==(fld+offset);
748  }
749  else
750  {
751  const scalar scale = mag(wantedAverage)/mag(averagePsi);
752 
753  if (debug)
754  {
755  Pout<< "updateCoeffs :"
756  << " scaling with:" << scale << endl;
757  }
758  this->operator==(scale*fld);
759  }
760  }
761 
762  if (debug)
763  {
764  Pout<< "updateCoeffs : set fixedValue to min:" << gMin(*this)
765  << " max:" << gMax(*this) << endl;
766  }
767 
769 }
770 
771 
772 template<class Type>
774 {
776  os.writeKeyword("setAverage") << setAverage_ << token::END_STATEMENT << nl;
777  os.writeKeyword("peturb") << perturb_ << token::END_STATEMENT << nl;
778 
779  if (fieldTableName_ != this->dimensionedInternalField().name())
780  {
781  os.writeKeyword("fieldTableName") << fieldTableName_
782  << token::END_STATEMENT << nl;
783  }
784 
785  this->writeEntry("value", os);
786 }
787 
788 
789 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
790 
791 } // End namespace Foam
792 
793 // ************************ vim: set sw=4 sts=4 et: ************************ //