src/finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValueFvPatchField.C

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