src/finiteArea/interpolation/edgeInterpolation/edgeInterpolation.C

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