src/finiteArea/interpolation/edgeInterpolation/edgeInterpolation.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright held by original author
   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 Description
  26     Face to edge interpolation scheme. Included in faMesh.
  27
  28 \*---------------------------------------------------------------------------*/
  29
  30 #include "faMesh.H"
  31 #include "areaFields.H"
  32 #include "edgeFields.H"
  33 #include "demandDrivenData.H"
  34 #include "faPatchFields.H"
  35
  36 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  37
  38 namespace Foam
  39 {
  40
  41 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  42
  43 defineTypeNameAndDebug(edgeInterpolation, 0);
  44
  45
  46 // * * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * //
  47
  48 void edgeInterpolation::clearOut()
  49 {
  50     deleteDemandDrivenData(lPN_);
  51     deleteDemandDrivenData(weightingFactors_);
  52     deleteDemandDrivenData(differenceFactors_);
  53     deleteDemandDrivenData(correctionVectors_);
  54     deleteDemandDrivenData(skewCorrectionVectors_);
  55 //     deleteDemandDrivenData(leastSquarePvectors_);
  56 //     deleteDemandDrivenData(leastSquareNvectors_);
  57 }
  58
  59
  60 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
  61
  62 edgeInterpolation::edgeInterpolation(const faMesh& fam)
  63 :
  64     faSchemes(fam()),
  65     faSolution(fam()),
  66     faMesh_(fam),
  67     lPN_(NULL),
  68     weightingFactors_(NULL),
  69     differenceFactors_(NULL),
  70     orthogonal_(false),
  71     correctionVectors_(NULL),
  72     skew_(true),
  73     skewCorrectionVectors_(NULL)
  74 //     leastSquarePvectors_(NULL),
  75 //     leastSquareNvectors_(NULL)
  76 {}
  77
  78
  79 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
  80
  81 edgeInterpolation::~edgeInterpolation()
  82 {
  83     clearOut();
  84 }
  85
  86
  87 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  88
  89 const edgeScalarField& edgeInterpolation::lPN() const
  90 {
  91     if (!lPN_)
  92     {
  93         makeLPN();
  94     }
  95
  96     return (*lPN_);
  97 }
  98
  99
 100 const edgeScalarField& edgeInterpolation::weights() const
 101 {
 102     if (!weightingFactors_)
 103     {
 104         makeWeights();
 105     }
 106
 107     return (*weightingFactors_);
 108 }
 109
 110
 111 const edgeScalarField& edgeInterpolation::deltaCoeffs() const
 112 {
 113     if (!differenceFactors_)
 114     {
 115         makeDeltaCoeffs();
 116     }
 117
 118     return (*differenceFactors_);
 119 }
 120
 121
 122 bool edgeInterpolation::orthogonal() const
 123 {
 124     if (orthogonal_ == false && !correctionVectors_)
 125     {
 126         makeCorrectionVectors();
 127     }
 128
 129     return orthogonal_;
 130 }
 131
 132
 133 const edgeVectorField& edgeInterpolation::correctionVectors() const
 134 {
 135     if (orthogonal())
 136     {
 137         FatalErrorIn("edgeInterpolation::correctionVectors()")
 138             << "cannot return correctionVectors; mesh is orthogonal"
 139             << abort(FatalError);
 140     }
 141
 142     return (*correctionVectors_);
 143 }
 144
 145
 146 bool edgeInterpolation::skew() const
 147 {
 148     if (skew_ == true && !skewCorrectionVectors_)
 149     {
 150         makeSkewCorrectionVectors();
 151     }
 152
 153     return skew_;
 154 }
 155
 156
 157 const edgeVectorField& edgeInterpolation::skewCorrectionVectors() const
 158 {
 159     if (!skew())
 160     {
 161         FatalErrorIn("edgeInterpolation::skewCorrectionVectors()")
 162             << "cannot return skewCorrectionVectors; mesh is now skewed"
 163             << abort(FatalError);
 164     }
 165
 166     return (*skewCorrectionVectors_);
 167 }
 168
 169
 170 // const edgeVectorField& edgeInterpolation::leastSquarePvectors() const
 171 // {
 172 //     if (!leastSquarePvectors_)
 173 //     {
 174 //         makeLeastSquareVectors();
 175 //     }
 176
 177 //     return (*leastSquarePvectors_);
 178 // }
 179
 180
 181 // const edgeVectorField& edgeInterpolation::leastSquareNvectors() const
 182 // {
 183 //     if (!leastSquareNvectors_)
 184 //     {
 185 //         makeLeastSquareVectors();
 186 //     }
 187
 188 //     return (*leastSquareNvectors_);
 189 // }
 190
 191
 192 // Do what is neccessary if the mesh has moved
 193 bool edgeInterpolation::movePoints()
 194 {
 195     deleteDemandDrivenData(lPN_);
 196     deleteDemandDrivenData(weightingFactors_);
 197     deleteDemandDrivenData(differenceFactors_);
 198
 199     orthogonal_ = false;
 200     deleteDemandDrivenData(correctionVectors_);
 201
 202     skew_ = true;
 203     deleteDemandDrivenData(skewCorrectionVectors_);
 204
 205 //     deleteDemandDrivenData(leastSquarePvectors_);
 206 //     deleteDemandDrivenData(leastSquareNvectors_);
 207
 208     return true;
 209 }
 210
 211
 212 void edgeInterpolation::makeLPN() const
 213 {
 214     if (debug)
 215     {
 216         Info<< "edgeInterpolation::makeLPN() : "
 217             << "Constructing geodesic distance between points P and N"
 218             << endl;
 219     }
 220
 221
 222     lPN_ = new edgeScalarField
 223     (
 224         IOobject
 225         (
 226             "lPN",
 227             faSolution::time().constant(),
 228             faSolution::db(),
 229             IOobject::NO_READ,
 230             IOobject::NO_WRITE,
 231             false
 232         ),
 233         mesh(),
 234         dimLength
 235     );
 236     edgeScalarField& lPN = *lPN_;
 237
 238
 239     // Set local references to mesh data
 240     const edgeVectorField& edgeCentres = mesh().edgeCentres();
 241     const areaVectorField& faceCentres = mesh().areaCentres();
 242     const unallocLabelList& owner = mesh().owner();
 243     const unallocLabelList& neighbour = mesh().neighbour();
 244
 245     forAll(owner, edgeI)
 246     {
 247         vector curSkewCorrVec = vector::zero;
 248
 249         if (skew())
 250         {
 251             curSkewCorrVec = skewCorrectionVectors()[edgeI];
 252         }
 253
 254         scalar lPE =
 255             mag
 256             (
 257                 edgeCentres[edgeI]
 258               - curSkewCorrVec
 259               - faceCentres[owner[edgeI]]
 260             );
 261
 262         scalar lEN =
 263             mag
 264             (
 265                 faceCentres[neighbour[edgeI]]
 266               - edgeCentres[edgeI]
 267               + curSkewCorrVec
 268             );
 269
 270         lPN.internalField()[edgeI] = (lPE + lEN);
 271     }
 272
 273
 274     forAll(lPN.boundaryField(), patchI)
 275     {
 276         mesh().boundary()[patchI].makeDeltaCoeffs
 277         (
 278             lPN.boundaryField()[patchI]
 279         );
 280
 281         lPN.boundaryField()[patchI] = 1.0/lPN.boundaryField()[patchI];
 282     }
 283
 284
 285     if (debug)
 286     {
 287         Info<< "edgeInterpolation::makeLPN() : "
 288             << "Finished constructing geodesic distance PN"
 289             << endl;
 290     }
 291 }
 292
 293
 294 void edgeInterpolation::makeWeights() const
 295 {
 296     if (debug)
 297     {
 298         Info<< "edgeInterpolation::makeWeights() : "
 299             << "Constructing weighting factors for edge interpolation"
 300             << endl;
 301     }
 302
 303
 304     weightingFactors_ = new edgeScalarField
 305     (
 306         IOobject
 307         (
 308             "weightingFactors",
 309             mesh()().pointsInstance(),
 310             mesh()(),
 311             IOobject::NO_READ,
 312             IOobject::NO_WRITE,
 313             false
 314         ),
 315         mesh(),
 316         dimless
 317     );
 318     edgeScalarField& weightingFactors = *weightingFactors_;
 319
 320
 321     // Set local references to mesh data
 322     const edgeVectorField& edgeCentres = mesh().edgeCentres();
 323     const areaVectorField& faceCentres = mesh().areaCentres();
 324     const unallocLabelList& owner = mesh().owner();
 325     const unallocLabelList& neighbour = mesh().neighbour();
 326
 327     forAll(owner, edgeI)
 328     {
 329         vector curSkewCorrVec = vector::zero;
 330
 331         if (skew())
 332         {
 333             curSkewCorrVec = skewCorrectionVectors()[edgeI];
 334         }
 335
 336         scalar lPE =
 337             mag
 338             (
 339                 edgeCentres[edgeI]
 340               - curSkewCorrVec
 341               - faceCentres[owner[edgeI]]
 342             );
 343
 344         scalar lEN =
 345             mag
 346             (
 347                 faceCentres[neighbour[edgeI]]
 348               - edgeCentres[edgeI]
 349               + curSkewCorrVec
 350             );
 351
 352         weightingFactors.internalField()[edgeI] =
 353             lEN
 354             /(
 355                 lPE
 356 #               ifdef BAD_MESH_STABILISATION
 357               + VSMALL
 358 #               endif
 359               + lEN
 360             );
 361     }
 362
 363     forAll(mesh().boundary(), patchI)
 364     {
 365         mesh().boundary()[patchI].makeWeights
 366         (
 367             weightingFactors.boundaryField()[patchI]
 368         );
 369     }
 370
 371     if (debug)
 372     {
 373         Info<< "edgeInterpolation::makeWeights() : "
 374             << "Finished constructing weighting factors for face interpolation"
 375             << endl;
 376     }
 377 }
 378
 379
 380 void edgeInterpolation::makeDeltaCoeffs() const
 381 {
 382     if (debug)
 383     {
 384         Info<< "edgeInterpolation::makeDeltaCoeffs() : "
 385             << "Constructing differencing factors array for edge gradient"
 386             << endl;
 387     }
 388
 389     // Force the construction of the weighting factors
 390     // needed to make sure deltaCoeffs are calculated for parallel runs.
 391     weights();
 392
 393     differenceFactors_ = new edgeScalarField
 394     (
 395         IOobject
 396         (
 397             "differenceFactors_",
 398             mesh()().pointsInstance(),
 399             mesh()(),
 400             IOobject::NO_READ,
 401             IOobject::NO_WRITE,
 402             false
 403         ),
 404         mesh(),
 405         dimless/dimLength
 406     );
 407     edgeScalarField& DeltaCoeffs = *differenceFactors_;
 408     scalarField& dc = DeltaCoeffs.internalField();
 409
 410
 411     // Set local references to mesh data
 412     const edgeVectorField& edgeCentres = mesh().edgeCentres();
 413     const areaVectorField& faceCentres = mesh().areaCentres();
 414     const unallocLabelList& owner = mesh().owner();
 415     const unallocLabelList& neighbour = mesh().neighbour();
 416     const edgeVectorField& lengths = mesh().Le();
 417
 418     const edgeList& edges = mesh().edges();
 419     const pointField& points = mesh().points();
 420
 421
 422     forAll(owner, edgeI)
 423     {
 424         // Edge normal - area normal
 425         vector edgeNormal = lengths[edgeI]^edges[edgeI].vec(points);
 426         edgeNormal /= mag(edgeNormal);
 427
 428
 429         // Unit delta vector
 430         vector unitDelta =
 431             faceCentres[neighbour[edgeI]]
 432           - faceCentres[owner[edgeI]];
 433
 434         unitDelta -=
 435             edgeNormal*(edgeNormal&unitDelta);
 436
 437         unitDelta /= mag(unitDelta);
 438
 439
 440         // Calc PN arc length
 441         vector curSkewCorrVec = vector::zero;
 442
 443         if (skew())
 444         {
 445             curSkewCorrVec = skewCorrectionVectors()[edgeI];
 446         }
 447
 448         scalar lPE =
 449             mag
 450             (
 451                 edgeCentres[edgeI]
 452               - curSkewCorrVec
 453               - faceCentres[owner[edgeI]]
 454             );
 455
 456         scalar lEN =
 457             mag
 458             (
 459                 faceCentres[neighbour[edgeI]]
 460               - edgeCentres[edgeI]
 461               + curSkewCorrVec
 462             );
 463
 464         scalar lPN = lPE + lEN;
 465
 466
 467         // Edge normal - area tangent
 468         edgeNormal = lengths[edgeI]/mag(lengths[edgeI]);
 469
 470         // Delta coefficient as per definition
 471 //         dc[edgeI] = 1.0/(lPN*(unitDelta & edgeNormal));
 472
 473         // Stabilised form for bad meshes.  HJ, 23/Jul/2009
 474         dc[edgeI] = 1.0/max((lPN*(unitDelta & edgeNormal)), 0.05*lPN);
 475     }
 476
 477
 478     forAll(DeltaCoeffs.boundaryField(), patchI)
 479     {
 480         mesh().boundary()[patchI].makeDeltaCoeffs
 481         (
 482             DeltaCoeffs.boundaryField()[patchI]
 483         );
 484     }
 485 }
 486
 487
 488 void edgeInterpolation::makeCorrectionVectors() const
 489 {
 490     if (debug)
 491     {
 492         Info<< "edgeInterpolation::makeCorrectionVectors() : "
 493             << "Constructing non-orthogonal correction vectors"
 494             << endl;
 495     }
 496
 497     correctionVectors_ = new edgeVectorField
 498     (
 499         IOobject
 500         (
 501             "correctionVectors",
 502             mesh()().pointsInstance(),
 503             mesh()(),
 504             IOobject::NO_READ,
 505             IOobject::NO_WRITE,
 506             false
 507         ),
 508         mesh(),
 509         dimless
 510     );
 511     edgeVectorField& CorrVecs = *correctionVectors_;
 512
 513     // Set local references to mesh data
 514     const areaVectorField& faceCentres = mesh().areaCentres();
 515
 516     const unallocLabelList& owner = mesh().owner();
 517     const unallocLabelList& neighbour = mesh().neighbour();
 518
 519     const edgeVectorField& lengths = mesh().Le();
 520     const edgeScalarField& magLengths = mesh().magLe();
 521
 522     const edgeList& edges = mesh().edges();
 523     const pointField& points = mesh().points();
 524
 525     scalarField deltaCoeffs(owner.size());
 526
 527     forAll(owner, edgeI)
 528     {
 529         // Edge normal - area normal
 530         vector edgeNormal = lengths[edgeI] ^ edges[edgeI].vec(points);
 531         edgeNormal /= mag(edgeNormal);
 532
 533         // Unit delta vector
 534         vector unitDelta =
 535             faceCentres[neighbour[edgeI]]
 536           - faceCentres[owner[edgeI]];
 537
 538         unitDelta -= edgeNormal*(edgeNormal & unitDelta);
 539
 540         unitDelta /= mag(unitDelta);
 541
 542         // Edge normal - area tangent
 543         edgeNormal = lengths[edgeI]/magLengths[edgeI];
 544
 545         // Delta coeffs
 546         deltaCoeffs[edgeI] = 1.0/(unitDelta & edgeNormal);
 547
 548         // Edge correction vector
 549         CorrVecs.internalField()[edgeI] =
 550             edgeNormal
 551           - deltaCoeffs[edgeI]*unitDelta;
 552     }
 553
 554     forAll(CorrVecs.boundaryField(), patchI)
 555     {
 556         faePatchVectorField& patchCorrVecs = CorrVecs.boundaryField()[patchI];
 557
 558         patchCorrVecs = vector::zero;
 559     }
 560
 561     scalar NonOrthogCoeff = 0.0;
 562
 563     if (owner.size() > 0)
 564     {
 565         scalarField sinAlpha = deltaCoeffs*mag(CorrVecs.internalField());
 566
 567         forAll(sinAlpha, edgeI)
 568         {
 569             sinAlpha[edgeI] = max(-1, min(sinAlpha[edgeI], 1));
 570         }
 571
 572         NonOrthogCoeff = max(Foam::asin(sinAlpha)*180.0/M_PI);
 573     }
 574
 575     if (debug)
 576     {
 577         Info<< "edgeInterpolation::makeCorrectionVectors() : "
 578             << "non-orthogonality coefficient = " << NonOrthogCoeff << " deg."
 579             << endl;
 580     }
 581
 582     if (NonOrthogCoeff < 0.1)
 583     {
 584         orthogonal_ = true;
 585         deleteDemandDrivenData(correctionVectors_);
 586     }
 587     else
 588     {
 589         orthogonal_ = false;
 590     }
 591
 592     if (debug)
 593     {
 594         Info<< "edgeInterpolation::makeCorrectionVectors() : "
 595             << "Finished constructing non-orthogonal correction vectors"
 596             << endl;
 597     }
 598 }
 599
 600
 601 void edgeInterpolation::makeSkewCorrectionVectors() const
 602 {
 603     if (debug)
 604     {
 605         Info<< "edgeInterpolation::makeSkewCorrectionVectors() : "
 606             << "Constructing skew correction vectors"
 607             << endl;
 608     }
 609
 610     skewCorrectionVectors_ = new edgeVectorField
 611     (
 612         IOobject
 613         (
 614             "skewCorrectionVectors",
 615             mesh()().pointsInstance(),
 616             mesh()(),
 617             IOobject::NO_READ,
 618             IOobject::NO_WRITE,
 619             false
 620         ),
 621         mesh(),
 622         dimless
 623     );
 624     edgeVectorField& SkewCorrVecs = *skewCorrectionVectors_;
 625
 626     // Set local references to mesh data
 627     const areaVectorField& C = mesh().areaCentres();
 628     const edgeVectorField& Ce = mesh().edgeCentres();
 629
 630     const unallocLabelList& owner = mesh().owner();
 631     const unallocLabelList& neighbour = mesh().neighbour();
 632
 633     const pointField& points = mesh().points();
 634     const edgeList& edges = mesh().edges();
 635
 636
 637     forAll(neighbour, edgeI)
 638     {
 639         vector P = C[owner[edgeI]];
 640         vector N = C[neighbour[edgeI]];
 641         vector S = points[edges[edgeI].start()];
 642         vector e = edges[edgeI].vec(points);
 643
 644         scalar alpha = - ( ( (N - P)^(S - P) )&( (N - P)^e ) )/
 645             ( ( (N - P)^e )&( (N - P)^e ) );
 646
 647         vector E = S + alpha*e;
 648
 649         SkewCorrVecs[edgeI] = Ce[edgeI] - E;
 650     }
 651
 652
 653     forAll(SkewCorrVecs.boundaryField(), patchI)
 654     {
 655         faePatchVectorField& patchSkewCorrVecs =
 656             SkewCorrVecs.boundaryField()[patchI];
 657
 658         if (patchSkewCorrVecs.coupled())
 659         {
 660             const unallocLabelList& edgeFaces =
 661                 mesh().boundary()[patchI].edgeFaces();
 662
 663             const edgeList::subList patchEdges =
 664                 mesh().boundary()[patchI].patchSlice(edges);
 665
 666             vectorField ngbC =
 667                 C.boundaryField()[patchI].patchNeighbourField();
 668
 669             forAll (patchSkewCorrVecs, edgeI)
 670             {
 671                 vector P = C[edgeFaces[edgeI]];
 672                 vector N = ngbC[edgeI];
 673                 vector S = points[patchEdges[edgeI].start()];
 674                 vector e = patchEdges[edgeI].vec(points);
 675
 676                 scalar alpha = - ( ( (N - P)^(S - P) )&( (N - P)^e ) )/
 677                     ( ( (N - P)^e )&( (N - P)^e ) );
 678
 679                 vector E = S + alpha*e;
 680
 681                 patchSkewCorrVecs[edgeI] =
 682                     Ce.boundaryField()[patchI][edgeI] - E;
 683             }
 684         }
 685         else
 686         {
 687             patchSkewCorrVecs = vector::zero;
 688         }
 689     }
 690
 691
 692     scalar skewCoeff = 0.0;
 693
 694     // Calculating PN arc length
 695     scalarField lPN(owner.size());
 696
 697     forAll (owner, edgeI)
 698     {
 699         lPN[edgeI] =
 700             mag
 701             (
 702                 Ce[edgeI]
 703               - SkewCorrVecs[edgeI]
 704               - C[owner[edgeI]]
 705             )
 706           + mag
 707             (
 708                 C[neighbour[edgeI]]
 709               - Ce[edgeI]
 710               + SkewCorrVecs[edgeI]
 711             );
 712     }
 713
 714     if (lPN.size() > 0)
 715     {
 716         skewCoeff = max(mag(SkewCorrVecs.internalField())/mag(lPN));
 717     }
 718
 719     if (debug)
 720     {
 721         Info<< "edgeInterpolation::makeSkewCorrectionVectors() : "
 722             << "skew coefficient = " << skewCoeff << endl;
 723     }
 724
 725     if (skewCoeff < 0.1)
 726     {
 727         skew_ = false;
 728         deleteDemandDrivenData(skewCorrectionVectors_);
 729     }
 730     else
 731     {
 732         skew_ = true;
 733     }
 734
 735     if (debug)
 736     {
 737         Info<< "edgeInterpolation::makeSkewCorrectionVectors() : "
 738             << "Finished constructing skew correction vectors"
 739             << endl;
 740     }
 741 }
 742
 743
 744 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 745
 746 } // End namespace Foam
 747
 748 // ************************************************************************* //