src/dynamicMesh/motionSmoother/motionSmoother.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2009 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 "motionSmoother.H"
  28 #include "twoDPointCorrector.H"
  29 #include "faceSet.H"
  30 #include "pointSet.H"
  31 #include "fixedValuePointPatchFields.H"
  32 #include "pointConstraint.H"
  33 #include "syncTools.H"
  34 #include "meshTools.H"
  35 #include "OFstream.H"
  36
  37 namespace Foam
  38 {
  39
  40 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  41
  42 defineTypeNameAndDebug(motionSmoother, 0);
  43
  44 }
  45
  46
  47 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  48
  49 // From pointPatchInterpolation
  50 void Foam::motionSmoother::makePatchPatchAddressing()
  51 {
  52     if (debug)
  53     {
  54         Pout<< "motionSmoother::makePatchPatchAddressing() : "
  55             << "constructing boundary addressing"
  56             << endl;
  57     }
  58
  59     const polyBoundaryMesh& bm = mesh_.boundaryMesh();
  60     const pointBoundaryMesh& pbm = pMesh_.boundary();
  61
  62     // first count the total number of patch-patch points
  63
  64     label nPatchPatchPoints = 0;
  65
  66     forAll(bm, patchi)
  67     {
  68         if(!isA<emptyPolyPatch>(bm[patchi]))
  69         {
  70             nPatchPatchPoints += bm[patchi].boundaryPoints().size();
  71         }
  72     }
  73
  74
  75     // Go through all patches and mark up the external edge points
  76     Map<label> patchPatchPointSet(2*nPatchPatchPoints);
  77
  78     labelList patchPatchPoints(nPatchPatchPoints);
  79
  80     List<pointConstraint> patchPatchPointConstraints(nPatchPatchPoints);
  81
  82     label pppi = 0;
  83
  84     forAll(bm, patchi)
  85     {
  86         if(!isA<emptyPolyPatch>(bm[patchi]))
  87         {
  88             const labelList& bp = bm[patchi].boundaryPoints();
  89             const labelList& meshPoints = bm[patchi].meshPoints();
  90
  91             forAll(bp, pointi)
  92             {
  93                 label ppp = meshPoints[bp[pointi]];
  94
  95                 Map<label>::iterator iter = patchPatchPointSet.find(ppp);
  96
  97                 if (iter == patchPatchPointSet.end())
  98                 {
  99                     patchPatchPointSet.insert(ppp, pppi);
 100                     patchPatchPoints[pppi] = ppp;
 101                     pbm[patchi].applyConstraint
 102                     (
 103                         bp[pointi],
 104                         patchPatchPointConstraints[pppi]
 105                     );
 106                     pppi++;
 107                 }
 108                 else
 109                 {
 110                     pbm[patchi].applyConstraint
 111                     (
 112                         bp[pointi],
 113                         patchPatchPointConstraints[iter()]
 114                     );
 115                 }
 116             }
 117         }
 118     }
 119
 120     nPatchPatchPoints = pppi;
 121     patchPatchPoints.setSize(nPatchPatchPoints);
 122     patchPatchPointConstraints.setSize(nPatchPatchPoints);
 123
 124     patchPatchPointConstraintPoints_.setSize(nPatchPatchPoints);
 125     patchPatchPointConstraintTensors_.setSize(nPatchPatchPoints);
 126
 127     label nConstraints = 0;
 128
 129     forAll(patchPatchPointConstraints, i)
 130     {
 131         if (patchPatchPointConstraints[i].first() != 0)
 132         {
 133             patchPatchPointConstraintPoints_[nConstraints] =
 134                 patchPatchPoints[i];
 135
 136             patchPatchPointConstraintTensors_[nConstraints] =
 137                 patchPatchPointConstraints[i].constraintTransformation();
 138
 139             nConstraints++;
 140         }
 141     }
 142
 143     patchPatchPointConstraintPoints_.setSize(nConstraints);
 144     patchPatchPointConstraintTensors_.setSize(nConstraints);
 145
 146
 147     if (debug)
 148     {
 149         OFstream str(mesh_.db().path()/"constraintPoints.obj");
 150
 151         Pout<< "Dumping " << patchPatchPointConstraintPoints_.size()
 152             << " constraintPoints to " << str.name() << endl;
 153         forAll(patchPatchPointConstraintPoints_, i)
 154         {
 155             label pointI = patchPatchPointConstraintPoints_[i];
 156
 157             meshTools::writeOBJ(str, mesh_.points()[pointI]);
 158         }
 159
 160         Pout<< "motionSmoother::makePatchPatchAddressing() : "
 161             << "finished constructing boundary addressing"
 162             << endl;
 163     }
 164 }
 165
 166
 167 void Foam::motionSmoother::checkFld(const pointScalarField& fld)
 168 {
 169     forAll(fld, pointI)
 170     {
 171         const scalar val = fld[pointI];
 172
 173         if ((val > -GREAT) && (val < GREAT))
 174         {}
 175         else
 176         {
 177             FatalErrorIn("motionSmoother::checkFld")
 178                 << "Problem : point:" << pointI << " value:" << val
 179                 << abort(FatalError);
 180         }
 181     }
 182 }
 183
 184
 185 Foam::labelHashSet Foam::motionSmoother::getPoints
 186 (
 187     const labelHashSet& faceLabels
 188 ) const
 189 {
 190     labelHashSet usedPoints(mesh_.nPoints()/100);
 191
 192     forAllConstIter(labelHashSet, faceLabels, iter)
 193     {
 194         const face& f = mesh_.faces()[iter.key()];
 195
 196         forAll(f, fp)
 197         {
 198             usedPoints.insert(f[fp]);
 199         }
 200     }
 201
 202     return usedPoints;
 203 }
 204
 205
 206 // Smooth on selected points (usually patch points)
 207 void Foam::motionSmoother::minSmooth
 208 (
 209     const PackedBoolList& isAffectedPoint,
 210     const labelList& meshPoints,
 211     const pointScalarField& fld,
 212     pointScalarField& newFld
 213 ) const
 214 {
 215     tmp<pointScalarField> tavgFld = avg
 216     (
 217         fld,
 218         scalarField(mesh_.nEdges(), 1.0),   // uniform weighting
 219         false                               // fld is not position.
 220     );
 221     const pointScalarField& avgFld = tavgFld();
 222
 223     forAll(meshPoints, i)
 224     {
 225         label pointI = meshPoints[i];
 226         if (isAffectedPoint.get(pointI) == 1)
 227         {
 228             newFld[pointI] = min
 229             (
 230                 fld[pointI],
 231                 0.5*fld[pointI] + 0.5*avgFld[pointI]
 232             );
 233         }
 234     }
 235
 236     newFld.correctBoundaryConditions();
 237     applyCornerConstraints(newFld);
 238 }
 239
 240
 241 // Smooth on all internal points
 242 void Foam::motionSmoother::minSmooth
 243 (
 244     const PackedBoolList& isAffectedPoint,
 245     const pointScalarField& fld,
 246     pointScalarField& newFld
 247 ) const
 248 {
 249     tmp<pointScalarField> tavgFld = avg
 250     (
 251         fld,
 252         scalarField(mesh_.nEdges(), 1.0),   // uniform weighting
 253         false                               // fld is not position.
 254     );
 255     const pointScalarField& avgFld = tavgFld();
 256
 257     forAll(fld, pointI)
 258     {
 259         if (isAffectedPoint.get(pointI) == 1 && isInternalPoint(pointI))
 260         {
 261             newFld[pointI] = min
 262             (
 263                 fld[pointI],
 264                 0.5*fld[pointI] + 0.5*avgFld[pointI]
 265             );
 266         }
 267     }
 268
 269     newFld.correctBoundaryConditions();
 270     applyCornerConstraints(newFld);
 271 }
 272
 273
 274 // Scale on selected points
 275 void Foam::motionSmoother::scaleField
 276 (
 277     const labelHashSet& pointLabels,
 278     const scalar scale,
 279     pointScalarField& fld
 280 ) const
 281 {
 282     forAllConstIter(labelHashSet, pointLabels, iter)
 283     {
 284         if (isInternalPoint(iter.key()))
 285         {
 286             fld[iter.key()] *= scale;
 287         }
 288     }
 289     fld.correctBoundaryConditions();
 290     applyCornerConstraints(fld);
 291 }
 292
 293
 294 // Scale on selected points (usually patch points)
 295 void Foam::motionSmoother::scaleField
 296 (
 297     const labelList& meshPoints,
 298     const labelHashSet& pointLabels,
 299     const scalar scale,
 300     pointScalarField& fld
 301 ) const
 302 {
 303     forAll(meshPoints, i)
 304     {
 305         label pointI = meshPoints[i];
 306
 307         if (pointLabels.found(pointI))
 308         {
 309             fld[pointI] *= scale;
 310         }
 311     }
 312 }
 313
 314
 315 bool Foam::motionSmoother::isInternalPoint(const label pointI) const
 316 {
 317     return isInternalPoint_.get(pointI) == 1;
 318 }
 319
 320
 321 void Foam::motionSmoother::getAffectedFacesAndPoints
 322 (
 323     const label nPointIter,
 324     const faceSet& wrongFaces,
 325
 326     labelList& affectedFaces,
 327     PackedBoolList& isAffectedPoint
 328 ) const
 329 {
 330     isAffectedPoint.setSize(mesh_.nPoints());
 331     isAffectedPoint = 0;
 332
 333     faceSet nbrFaces(mesh_, "checkFaces", wrongFaces);
 334
 335     // Find possible points influenced by nPointIter iterations of
 336     // scaling and smoothing by doing pointCellpoint walk.
 337     // Also update faces-to-be-checked to extend one layer beyond the points
 338     // that will get updated.
 339
 340     for (label i = 0; i < nPointIter; i++)
 341     {
 342         pointSet nbrPoints(mesh_, "grownPoints", getPoints(nbrFaces.toc()));
 343
 344         forAllConstIter(pointSet, nbrPoints, iter)
 345         {
 346             const labelList& pCells = mesh_.pointCells(iter.key());
 347
 348             forAll(pCells, pCellI)
 349             {
 350                 const cell& cFaces = mesh_.cells()[pCells[pCellI]];
 351
 352                 forAll(cFaces, cFaceI)
 353                 {
 354                     nbrFaces.insert(cFaces[cFaceI]);
 355                 }
 356             }
 357         }
 358         nbrFaces.sync(mesh_);
 359
 360         if (i == nPointIter - 2)
 361         {
 362             forAllConstIter(faceSet, nbrFaces, iter)
 363             {
 364                 const face& f = mesh_.faces()[iter.key()];
 365                 forAll(f, fp)
 366                 {
 367                     isAffectedPoint.set(f[fp], 1);
 368                 }
 369             }
 370         }
 371     }
 372
 373     affectedFaces = nbrFaces.toc();
 374 }
 375
 376
 377 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 378
 379 Foam::motionSmoother::motionSmoother
 380 (
 381     polyMesh& mesh,
 382     pointMesh& pMesh,
 383     indirectPrimitivePatch& pp,
 384     const labelList& adaptPatchIDs,
 385     const dictionary& paramDict
 386 )
 387 :
 388     mesh_(mesh),
 389     pMesh_(pMesh),
 390     pp_(pp),
 391     adaptPatchIDs_(adaptPatchIDs),
 392     paramDict_(paramDict),
 393     displacement_
 394     (
 395         IOobject
 396         (
 397             "displacement",
 398             mesh_.time().timeName(),
 399             mesh_,
 400             IOobject::MUST_READ,
 401             IOobject::AUTO_WRITE
 402         ),
 403         pMesh_
 404     ),
 405     scale_
 406     (
 407         IOobject
 408         (
 409             "scale",
 410             mesh_.time().timeName(),
 411             mesh_,
 412             IOobject::NO_READ,
 413             IOobject::AUTO_WRITE
 414         ),
 415         pMesh_,
 416         dimensionedScalar("scale", dimless, 1.0)
 417     ),
 418     oldPoints_(mesh_.points()),
 419     isInternalPoint_(mesh_.nPoints(), 1),
 420     twoDCorrector_(mesh_)
 421 {
 422     updateMesh();
 423 }
 424
 425
 426 Foam::motionSmoother::motionSmoother
 427 (
 428     polyMesh& mesh,
 429     indirectPrimitivePatch& pp,
 430     const labelList& adaptPatchIDs,
 431     const pointVectorField& displacement,
 432     const dictionary& paramDict
 433 )
 434 :
 435     mesh_(mesh),
 436     pMesh_(const_cast<pointMesh&>(displacement.mesh())),
 437     pp_(pp),
 438     adaptPatchIDs_(adaptPatchIDs),
 439     paramDict_(paramDict),
 440     displacement_
 441     (
 442         IOobject
 443         (
 444             "displacement",
 445             mesh_.time().timeName(),
 446             mesh_,
 447             IOobject::NO_READ,
 448             IOobject::AUTO_WRITE
 449         ),
 450         displacement
 451     ),
 452     scale_
 453     (
 454         IOobject
 455         (
 456             "scale",
 457             mesh_.time().timeName(),
 458             mesh_,
 459             IOobject::NO_READ,
 460             IOobject::AUTO_WRITE
 461         ),
 462         pMesh_,
 463         dimensionedScalar("scale", dimless, 1.0)
 464     ),
 465     oldPoints_(mesh_.points()),
 466     isInternalPoint_(mesh_.nPoints(), 1),
 467     twoDCorrector_(mesh_)
 468 {
 469     updateMesh();
 470 }
 471
 472
 473 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 474
 475 Foam::motionSmoother::~motionSmoother()
 476 {}
 477
 478
 479 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 480
 481 const Foam::polyMesh& Foam::motionSmoother::mesh() const
 482 {
 483     return mesh_;
 484 }
 485
 486
 487 const Foam::pointMesh& Foam::motionSmoother::pMesh() const
 488 {
 489     return pMesh_;
 490 }
 491
 492
 493 const Foam::indirectPrimitivePatch& Foam::motionSmoother::patch() const
 494 {
 495     return pp_;
 496 }
 497
 498
 499 const Foam::labelList& Foam::motionSmoother::adaptPatchIDs() const
 500 {
 501     return adaptPatchIDs_;
 502 }
 503
 504
 505 const Foam::dictionary& Foam::motionSmoother::paramDict() const
 506 {
 507     return paramDict_;
 508 }
 509
 510
 511 Foam::pointVectorField& Foam::motionSmoother::displacement()
 512 {
 513     return displacement_;
 514 }
 515
 516
 517 const Foam::pointVectorField& Foam::motionSmoother::displacement() const
 518 {
 519     return displacement_;
 520 }
 521
 522
 523 const Foam::pointScalarField& Foam::motionSmoother::scale() const
 524 {
 525     return scale_;
 526 }
 527
 528
 529 const Foam::pointField& Foam::motionSmoother::oldPoints() const
 530 {
 531     return oldPoints_;
 532 }
 533
 534
 535 void Foam::motionSmoother::correct()
 536 {
 537     oldPoints_ = mesh_.points();
 538
 539     scale_ = 1.0;
 540
 541     // No need to update twoDmotion corrector since only holds edge labels
 542     // which will remain the same as before. So unless the mesh was distorted
 543     // severely outside of motionSmoother there will be no need.
 544 }
 545
 546
 547 void Foam::motionSmoother::setDisplacement(pointField& patchDisp)
 548 {
 549     // See comment in .H file about shared points.
 550     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 551
 552     forAll(patches, patchI)
 553     {
 554         const polyPatch& pp = patches[patchI];
 555
 556         if (pp.coupled())
 557         {
 558             const labelList& meshPoints = pp.meshPoints();
 559
 560             forAll(meshPoints, i)
 561             {
 562                 displacement_[meshPoints[i]] = vector::zero;
 563             }
 564         }
 565     }
 566
 567     const labelList& ppMeshPoints = pp_.meshPoints();
 568
 569     // Set internal point data from displacement on combined patch points.
 570     forAll(ppMeshPoints, patchPointI)
 571     {
 572         displacement_[ppMeshPoints[patchPointI]] = patchDisp[patchPointI];
 573     }
 574
 575     // Adapt the fixedValue bc's (i.e. copy internal point data to
 576     // boundaryField for all affected patches)
 577     forAll(adaptPatchIDs_, i)
 578     {
 579         label patchI = adaptPatchIDs_[i];
 580
 581         displacement_.boundaryField()[patchI] ==
 582             displacement_.boundaryField()[patchI].patchInternalField();
 583     }
 584
 585     // Make consistent with non-adapted bc's by evaluating those now and
 586     // resetting the displacement from the values.
 587     // Note that we're just doing a correctBoundaryConditions with
 588     // fixedValue bc's first.
 589     labelHashSet adaptPatchSet(adaptPatchIDs_);
 590
 591     const lduSchedule& patchSchedule = mesh_.globalData().patchSchedule();
 592
 593     forAll(patchSchedule, patchEvalI)
 594     {
 595         label patchI = patchSchedule[patchEvalI].patch;
 596
 597         if (!adaptPatchSet.found(patchI))
 598         {
 599             if (patchSchedule[patchEvalI].init)
 600             {
 601                 displacement_.boundaryField()[patchI]
 602                     .initEvaluate(Pstream::scheduled);
 603             }
 604             else
 605             {
 606                 displacement_.boundaryField()[patchI]
 607                     .evaluate(Pstream::scheduled);
 608             }
 609         }
 610     }
 611
 612     // Multi-patch constraints
 613     applyCornerConstraints(displacement_);
 614
 615     // Correct for problems introduced by corner constraints
 616     syncTools::syncPointList
 617     (
 618         mesh_,
 619         displacement_,
 620         maxMagEqOp(),   // combine op
 621         vector::zero,   // null value
 622         false           // no separation
 623     );
 624
 625     // Adapt the fixedValue bc's (i.e. copy internal point data to
 626     // boundaryField for all affected patches) to take the changes caused
 627     // by multi-corner constraints into account.
 628     forAll(adaptPatchIDs_, i)
 629     {
 630         label patchI = adaptPatchIDs_[i];
 631
 632         displacement_.boundaryField()[patchI] ==
 633             displacement_.boundaryField()[patchI].patchInternalField();
 634     }
 635
 636     if (debug)
 637     {
 638         OFstream str(mesh_.db().path()/"changedPoints.obj");
 639         label nVerts = 0;
 640         forAll(ppMeshPoints, patchPointI)
 641         {
 642             const vector& newDisp = displacement_[ppMeshPoints[patchPointI]];
 643
 644             if (mag(newDisp-patchDisp[patchPointI]) > SMALL)
 645             {
 646                 const point& pt = mesh_.points()[ppMeshPoints[patchPointI]];
 647
 648                 meshTools::writeOBJ(str, pt);
 649                 nVerts++;
 650                 //Pout<< "Point:" << pt
 651                 //    << " oldDisp:" << patchDisp[patchPointI]
 652                 //    << " newDisp:" << newDisp << endl;
 653             }
 654         }
 655         Pout<< "Written " << nVerts << " points that are changed to file "
 656             << str.name() << endl;
 657     }
 658
 659     // Now reset input displacement
 660     forAll(ppMeshPoints, patchPointI)
 661     {
 662         patchDisp[patchPointI] = displacement_[ppMeshPoints[patchPointI]];
 663     }
 664 }
 665
 666
 667 // correctBoundaryConditions with fixedValue bc's first.
 668 void Foam::motionSmoother::correctBoundaryConditions
 669 (
 670     pointVectorField& displacement
 671 ) const
 672 {
 673     labelHashSet adaptPatchSet(adaptPatchIDs_);
 674
 675     const lduSchedule& patchSchedule = mesh_.globalData().patchSchedule();
 676
 677     // 1. evaluate on adaptPatches
 678     forAll(patchSchedule, patchEvalI)
 679     {
 680         label patchI = patchSchedule[patchEvalI].patch;
 681
 682         if (adaptPatchSet.found(patchI))
 683         {
 684             if (patchSchedule[patchEvalI].init)
 685             {
 686                 displacement.boundaryField()[patchI]
 687                     .initEvaluate(Pstream::blocking);
 688             }
 689             else
 690             {
 691                 displacement.boundaryField()[patchI]
 692                     .evaluate(Pstream::blocking);
 693             }
 694         }
 695     }
 696
 697
 698     // 2. evaluate on non-AdaptPatches
 699     forAll(patchSchedule, patchEvalI)
 700     {
 701         label patchI = patchSchedule[patchEvalI].patch;
 702
 703         if (!adaptPatchSet.found(patchI))
 704         {
 705             if (patchSchedule[patchEvalI].init)
 706             {
 707                 displacement.boundaryField()[patchI]
 708                     .initEvaluate(Pstream::blocking);
 709             }
 710             else
 711             {
 712                 displacement.boundaryField()[patchI]
 713                     .evaluate(Pstream::blocking);
 714             }
 715         }
 716     }
 717
 718     // Multi-patch constraints
 719     applyCornerConstraints(displacement);
 720
 721     // Correct for problems introduced by corner constraints
 722     syncTools::syncPointList
 723     (
 724         mesh_,
 725         displacement,
 726         maxMagEqOp(),           // combine op
 727         vector::zero,           // null value
 728         false                   // no separation
 729     );
 730 }
 731
 732
 733 Foam::tmp<Foam::scalarField> Foam::motionSmoother::movePoints
 734 (
 735     pointField& newPoints
 736 )
 737 {
 738     // Correct for 2-D motion
 739     if (twoDCorrector_.required())
 740     {
 741         Info<< "Correct-ing 2-D mesh motion";
 742
 743         if (mesh_.globalData().parallel())
 744         {
 745             WarningIn("motionSmoother::movePoints(pointField& newPoints)")
 746                 << "2D mesh-motion probably not correct in parallel" << endl;
 747         }
 748
 749         // We do not want to move 3D planes so project all points onto those
 750         const pointField& oldPoints = mesh_.points();
 751         const edgeList& edges = mesh_.edges();
 752
 753         const labelList& neIndices = twoDCorrector().normalEdgeIndices();
 754         const vector& pn = twoDCorrector().planeNormal();
 755
 756         forAll(neIndices, i)
 757         {
 758             const edge& e = edges[neIndices[i]];
 759
 760             point& pStart = newPoints[e.start()];
 761
 762             pStart += pn*(pn & (oldPoints[e.start()] - pStart));
 763
 764             point& pEnd = newPoints[e.end()];
 765
 766             pEnd += pn*(pn & (oldPoints[e.end()] - pEnd));
 767         }
 768
 769         // Correct tangentially
 770         twoDCorrector_.correctPoints(newPoints);
 771         Info<< " ...done" << endl;
 772     }
 773
 774     if (debug)
 775     {
 776         Pout<< "motionSmoother::movePoints : testing sync of newPoints."
 777             << endl;
 778         testSyncField
 779         (
 780             newPoints,
 781             minEqOp<point>(),           // combine op
 782             vector(GREAT,GREAT,GREAT),  // null
 783             true                        // separation
 784         );
 785     }
 786
 787     tmp<scalarField> tsweptVol = mesh_.movePoints(newPoints);
 788
 789     pp_.movePoints(mesh_.points());
 790
 791     return tsweptVol;
 792 }
 793
 794
 795 Foam::scalar Foam::motionSmoother::setErrorReduction
 796 (
 797     const scalar errorReduction
 798 )
 799 {
 800     scalar oldErrorReduction = readScalar(paramDict_.lookup("errorReduction"));
 801
 802     paramDict_.remove("errorReduction");
 803     paramDict_.add("errorReduction", errorReduction);
 804
 805     return oldErrorReduction;
 806 }
 807
 808
 809 bool Foam::motionSmoother::scaleMesh
 810 (
 811     labelList& checkFaces,
 812     const bool smoothMesh,
 813     const label nAllowableErrors
 814 )
 815 {
 816     List<labelPair> emptyBaffles;
 817     return scaleMesh
 818     (
 819         checkFaces,
 820         emptyBaffles,
 821         smoothMesh,
 822         nAllowableErrors
 823     );
 824 }
 825
 826
 827 bool Foam::motionSmoother::scaleMesh
 828 (
 829     labelList& checkFaces,
 830     const List<labelPair>& baffles,
 831     const bool smoothMesh,
 832     const label nAllowableErrors
 833 )
 834 {
 835     return scaleMesh
 836     (
 837         checkFaces,
 838         baffles,
 839         paramDict_,
 840         paramDict_,
 841         smoothMesh,
 842         nAllowableErrors
 843     );
 844 }
 845
 846
 847 bool Foam::motionSmoother::scaleMesh
 848 (
 849     labelList& checkFaces,
 850     const List<labelPair>& baffles,
 851     const dictionary& paramDict,
 852     const dictionary& meshQualityDict,
 853     const bool smoothMesh,
 854     const label nAllowableErrors
 855 )
 856 {
 857     if (!smoothMesh && adaptPatchIDs_.empty())
 858     {
 859         FatalErrorIn("motionSmoother::scaleMesh(const bool")
 860             << "You specified both no movement on the internal mesh points"
 861             << " (smoothMesh = false)" << nl
 862             << "and no movement on the patch (adaptPatchIDs is empty)" << nl
 863             << "Hence nothing to adapt."
 864             << exit(FatalError);
 865     }
 866
 867     if (debug)
 868     {
 869         // Had a problem with patches moved non-synced. Check transformations.
 870         const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 871
 872         Pout<< "Entering scaleMesh : coupled patches:" << endl;
 873         forAll(patches, patchI)
 874         {
 875             if (patches[patchI].coupled())
 876             {
 877                 const coupledPolyPatch& pp =
 878                     refCast<const coupledPolyPatch>(patches[patchI]);
 879
 880                 Pout<< '\t' << patchI << '\t' << pp.name()
 881                     << " parallel:" << pp.parallel()
 882                     << " separated:" << pp.separated()
 883                     << " forwardT:" << pp.forwardT().size()
 884                     << endl;
 885             }
 886         }
 887     }
 888
 889     const scalar errorReduction =
 890         readScalar(paramDict.lookup("errorReduction"));
 891     const label nSmoothScale =
 892         readLabel(paramDict.lookup("nSmoothScale"));
 893
 894
 895     // Note: displacement_ should already be synced already from setDisplacement
 896     // but just to make sure.
 897     syncTools::syncPointList
 898     (
 899         mesh_,
 900         displacement_,
 901         maxMagEqOp(),
 902         vector::zero,   // null value
 903         false           // no separation
 904     );
 905
 906     // Set newPoints as old + scale*displacement
 907     pointField newPoints;
 908     {
 909         // Create overall displacement with same b.c.s as displacement_
 910         pointVectorField totalDisplacement
 911         (
 912             IOobject
 913             (
 914                 "totalDisplacement",
 915                 mesh_.time().timeName(),
 916                 mesh_,
 917                 IOobject::NO_READ,
 918                 IOobject::NO_WRITE,
 919                 false
 920             ),
 921             scale_*displacement_,
 922             displacement_.boundaryField().types()
 923         );
 924         correctBoundaryConditions(totalDisplacement);
 925
 926         if (debug)
 927         {
 928             Pout<< "scaleMesh : testing sync of totalDisplacement" << endl;
 929             testSyncField
 930             (
 931                 totalDisplacement,
 932                 maxMagEqOp(),
 933                 vector::zero,   // null value
 934                 false           // separation
 935             );
 936         }
 937
 938         newPoints = oldPoints_ + totalDisplacement.internalField();
 939     }
 940
 941     Info<< "Moving mesh using diplacement scaling :"
 942         << " min:" << gMin(scale_.internalField())
 943         << "  max:" << gMax(scale_.internalField())
 944         << endl;
 945
 946
 947     // Move
 948     movePoints(newPoints);
 949
 950     // Check. Returns parallel number of incorrect faces.
 951     faceSet wrongFaces(mesh_, "wrongFaces", mesh_.nFaces()/100+100);
 952     checkMesh(false, mesh_, meshQualityDict, checkFaces, baffles, wrongFaces);
 953
 954     if (returnReduce(wrongFaces.size(), sumOp<label>()) <= nAllowableErrors)
 955     {
 956         return true;
 957     }
 958     else
 959     {
 960         // Sync across coupled faces by extending the set.
 961         wrongFaces.sync(mesh_);
 962
 963         // Special case:
 964         // if errorReduction is set to zero, extend wrongFaces
 965         // to face-Cell-faces to ensure quick return to previously valid mesh
 966
 967         if (mag(errorReduction) < SMALL)
 968         {
 969             labelHashSet newWrongFaces(wrongFaces);
 970             forAllConstIter(labelHashSet, wrongFaces, iter)
 971             {
 972                 label own = mesh_.faceOwner()[iter.key()];
 973                 const cell& ownFaces = mesh_.cells()[own];
 974
 975                 forAll(ownFaces, cfI)
 976                 {
 977                     newWrongFaces.insert(ownFaces[cfI]);
 978                 }
 979
 980                 if (iter.key() < mesh_.nInternalFaces())
 981                 {
 982                     label nei = mesh_.faceNeighbour()[iter.key()];
 983                     const cell& neiFaces = mesh_.cells()[nei];
 984
 985                     forAll(neiFaces, cfI)
 986                     {
 987                         newWrongFaces.insert(neiFaces[cfI]);
 988                     }
 989                 }
 990             }
 991             wrongFaces.transfer(newWrongFaces);
 992             wrongFaces.sync(mesh_);
 993         }
 994
 995
 996         // Find out points used by wrong faces and scale displacement.
 997         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 998
 999         pointSet usedPoints(mesh_, "usedPoints", getPoints(wrongFaces));
1000         usedPoints.sync(mesh_);
1001
1002
1003
1004         // Grow a few layers to determine
1005         // - points to be smoothed
1006         // - faces to be checked in next iteration
1007         PackedBoolList isAffectedPoint(mesh_.nPoints());
1008         getAffectedFacesAndPoints
1009         (
1010             nSmoothScale,       // smoothing iterations
1011             wrongFaces,         // error faces
1012             checkFaces,
1013             isAffectedPoint
1014         );
1015
1016         if (debug)
1017         {
1018             Pout<< "Faces in error:" << wrongFaces.size()
1019                 << "  with points:" << usedPoints.size()
1020                 << endl;
1021         }
1022
1023         if (adaptPatchIDs_.size())
1024         {
1025             // Scale conflicting patch points
1026             scaleField(pp_.meshPoints(), usedPoints, errorReduction, scale_);
1027         }
1028         if (smoothMesh)
1029         {
1030             // Scale conflicting internal points
1031             scaleField(usedPoints, errorReduction, scale_);
1032         }
1033
1034         for (label i = 0; i < nSmoothScale; i++)
1035         {
1036             if (adaptPatchIDs_.size())
1037             {
1038                 // Smooth patch values
1039                 pointScalarField oldScale(scale_);
1040                 minSmooth
1041                 (
1042                     isAffectedPoint,
1043                     pp_.meshPoints(),
1044                     oldScale,
1045                     scale_
1046                 );
1047                 checkFld(scale_);
1048             }
1049             if (smoothMesh)
1050             {
1051                 // Smooth internal values
1052                 pointScalarField oldScale(scale_);
1053                 minSmooth(isAffectedPoint, oldScale, scale_);
1054                 checkFld(scale_);
1055             }
1056         }
1057
1058         syncTools::syncPointList
1059         (
1060             mesh_,
1061             scale_,
1062             maxEqOp<scalar>(),
1063             -GREAT,             // null value
1064             false               // no separation
1065         );
1066
1067
1068         if (debug)
1069         {
1070             Pout<< "scale_ after smoothing :"
1071                 << " min:" << Foam::gMin(scale_)
1072                 << " max:" << Foam::gMax(scale_)
1073                 << endl;
1074         }
1075
1076         return false;
1077     }
1078 }
1079
1080
1081 void Foam::motionSmoother::updateMesh()
1082 {
1083     const pointBoundaryMesh& patches = pMesh_.boundary();
1084
1085     // Check whether displacement has fixed value b.c. on adaptPatchID
1086     forAll(adaptPatchIDs_, i)
1087     {
1088         label patchI = adaptPatchIDs_[i];
1089
1090         if
1091         (
1092            !isA<fixedValuePointPatchVectorField>
1093             (
1094                 displacement_.boundaryField()[patchI]
1095             )
1096         )
1097         {
1098             FatalErrorIn
1099             (
1100                 "motionSmoother::motionSmoother"
1101             )   << "Patch " << patches[patchI].name()
1102                 << " has wrong boundary condition "
1103                 << displacement_.boundaryField()[patchI].type()
1104                 << " on field " << displacement_.name() << nl
1105                 << "Only type allowed is "
1106                 << fixedValuePointPatchVectorField::typeName
1107                 << exit(FatalError);
1108         }
1109     }
1110
1111
1112     // Determine internal points. Note that for twoD there are no internal
1113     // points so we use the points of adaptPatchIDs instead
1114
1115     twoDCorrector_.updateMesh();
1116
1117     const labelList& meshPoints = pp_.meshPoints();
1118
1119     forAll(meshPoints, i)
1120     {
1121         isInternalPoint_.set(meshPoints[i], 0);
1122     }
1123
1124     // Calculate master edge addressing
1125     isMasterEdge_ = syncTools::getMasterEdges(mesh_);
1126
1127     makePatchPatchAddressing();
1128 }
1129
1130
1131 // Specialisation of applyCornerConstraints for scalars because
1132 // no constraint need be applied
1133 template<>
1134 void Foam::motionSmoother::applyCornerConstraints<Foam::scalar>
1135 (
1136     GeometricField<scalar, pointPatchField, pointMesh>& pf
1137 ) const
1138 {}
1139
1140
1141 // ************************************************************************* //