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-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 "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 PackedList<1>& 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 PackedList<1>& 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     PackedList<1>& 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 Foam::pointVectorField& Foam::motionSmoother::displacement()
 506 {
 507     return displacement_;
 508 }
 509
 510
 511 const Foam::pointVectorField& Foam::motionSmoother::displacement() const
 512 {
 513     return displacement_;
 514 }
 515
 516
 517 const Foam::pointScalarField& Foam::motionSmoother::scale() const
 518 {
 519     return scale_;
 520 }
 521
 522
 523 const Foam::pointField& Foam::motionSmoother::oldPoints() const
 524 {
 525     return oldPoints_;
 526 }
 527
 528
 529 void Foam::motionSmoother::correct()
 530 {
 531     oldPoints_ = mesh_.points();
 532
 533     scale_ = 1.0;
 534
 535     // No need to update twoDmotion corrector since only holds edge labels
 536     // which will remain the same as before. So unless the mesh was distorted
 537     // severely outside of motionSmoother there will be no need.
 538 }
 539
 540
 541 void Foam::motionSmoother::setDisplacement(pointField& patchDisp)
 542 {
 543     // See comment in .H file about shared points.
 544     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 545
 546     forAll(patches, patchI)
 547     {
 548         const polyPatch& pp = patches[patchI];
 549
 550         if (pp.coupled())
 551         {
 552             const labelList& meshPoints = pp.meshPoints();
 553
 554             forAll(meshPoints, i)
 555             {
 556                 displacement_[meshPoints[i]] = vector::zero;
 557             }
 558         }
 559     }
 560
 561     const labelList& ppMeshPoints = pp_.meshPoints();
 562
 563     // Set internal point data from displacement on combined patch points.
 564     forAll(ppMeshPoints, patchPointI)
 565     {
 566         displacement_[ppMeshPoints[patchPointI]] = patchDisp[patchPointI];
 567     }
 568
 569     // Adapt the fixedValue bc's (i.e. copy internal point data to
 570     // boundaryField for all affected patches)
 571     forAll(adaptPatchIDs_, i)
 572     {
 573         label patchI = adaptPatchIDs_[i];
 574
 575         displacement_.boundaryField()[patchI] ==
 576             displacement_.boundaryField()[patchI].patchInternalField();
 577     }
 578
 579     // Make consistent with non-adapted bc's by evaluating those now and
 580     // resetting the displacement from the values.
 581     // Note that we're just doing a correctBoundaryConditions with
 582     // fixedValue bc's first.
 583     labelHashSet adaptPatchSet(adaptPatchIDs_);
 584
 585     const lduSchedule& patchSchedule = mesh_.globalData().patchSchedule();
 586
 587     forAll(patchSchedule, patchEvalI)
 588     {
 589         label patchI = patchSchedule[patchEvalI].patch;
 590
 591         if (!adaptPatchSet.found(patchI))
 592         {
 593             if (patchSchedule[patchEvalI].init)
 594             {
 595                 displacement_.boundaryField()[patchI]
 596                     .initEvaluate(Pstream::scheduled);
 597             }
 598             else
 599             {
 600                 displacement_.boundaryField()[patchI]
 601                     .evaluate(Pstream::scheduled);
 602             }
 603         }
 604     }
 605
 606     // Multi-patch constraints
 607     applyCornerConstraints(displacement_);
 608
 609     // Correct for problems introduced by corner constraints
 610     syncTools::syncPointList
 611     (
 612         mesh_,
 613         displacement_,
 614         maxMagEqOp(),   // combine op
 615         vector::zero,   // null value
 616         false           // no separation
 617     );
 618
 619     // Adapt the fixedValue bc's (i.e. copy internal point data to
 620     // boundaryField for all affected patches) to take the changes caused
 621     // by multi-corner constraints into account.
 622     forAll(adaptPatchIDs_, i)
 623     {
 624         label patchI = adaptPatchIDs_[i];
 625
 626         displacement_.boundaryField()[patchI] ==
 627             displacement_.boundaryField()[patchI].patchInternalField();
 628     }
 629
 630     if (debug)
 631     {
 632         OFstream str(mesh_.db().path()/"changedPoints.obj");
 633         label nVerts = 0;
 634         forAll(ppMeshPoints, patchPointI)
 635         {
 636             const vector& newDisp = displacement_[ppMeshPoints[patchPointI]];
 637
 638             if (mag(newDisp-patchDisp[patchPointI]) > SMALL)
 639             {
 640                 const point& pt = mesh_.points()[ppMeshPoints[patchPointI]];
 641
 642                 meshTools::writeOBJ(str, pt);
 643                 nVerts++;
 644                 //Pout<< "Point:" << pt
 645                 //    << " oldDisp:" << patchDisp[patchPointI]
 646                 //    << " newDisp:" << newDisp << endl;
 647             }
 648         }
 649         Pout<< "Written " << nVerts << " points that are changed to file "
 650             << str.name() << endl;
 651     }
 652
 653     // Now reset input displacement
 654     forAll(ppMeshPoints, patchPointI)
 655     {
 656         patchDisp[patchPointI] = displacement_[ppMeshPoints[patchPointI]];
 657     }
 658 }
 659
 660
 661 // correctBoundaryConditions with fixedValue bc's first.
 662 void Foam::motionSmoother::correctBoundaryConditions
 663 (
 664     pointVectorField& displacement
 665 ) const
 666 {
 667     labelHashSet adaptPatchSet(adaptPatchIDs_);
 668
 669     const lduSchedule& patchSchedule = mesh_.globalData().patchSchedule();
 670
 671     // 1. evaluate on adaptPatches
 672     forAll(patchSchedule, patchEvalI)
 673     {
 674         label patchI = patchSchedule[patchEvalI].patch;
 675
 676         if (adaptPatchSet.found(patchI))
 677         {
 678             if (patchSchedule[patchEvalI].init)
 679             {
 680                 displacement.boundaryField()[patchI]
 681                     .initEvaluate(Pstream::blocking);
 682             }
 683             else
 684             {
 685                 displacement.boundaryField()[patchI]
 686                     .evaluate(Pstream::blocking);
 687             }
 688         }
 689     }
 690
 691
 692     // 2. evaluate on non-AdaptPatches
 693     forAll(patchSchedule, patchEvalI)
 694     {
 695         label patchI = patchSchedule[patchEvalI].patch;
 696
 697         if (!adaptPatchSet.found(patchI))
 698         {
 699             if (patchSchedule[patchEvalI].init)
 700             {
 701                 displacement.boundaryField()[patchI]
 702                     .initEvaluate(Pstream::blocking);
 703             }
 704             else
 705             {
 706                 displacement.boundaryField()[patchI]
 707                     .evaluate(Pstream::blocking);
 708             }
 709         }
 710     }
 711
 712     // Multi-patch constraints
 713     applyCornerConstraints(displacement);
 714
 715     // Correct for problems introduced by corner constraints
 716     syncTools::syncPointList
 717     (
 718         mesh_,
 719         displacement,
 720         maxMagEqOp(),           // combine op
 721         vector::zero,           // null value
 722         false                   // no separation
 723     );
 724 }
 725
 726
 727 Foam::tmp<Foam::scalarField> Foam::motionSmoother::movePoints
 728 (
 729     pointField& newPoints
 730 )
 731 {
 732     // Correct for 2-D motion
 733     if (twoDCorrector_.required())
 734     {
 735         Info<< "Correct-ing 2-D mesh motion";
 736
 737         if (mesh_.globalData().parallel())
 738         {
 739             WarningIn("motionSmoother::movePoints(pointField& newPoints)")
 740                 << "2D mesh-motion probably not correct in parallel" << endl;
 741         }
 742
 743         // We do not want to move 3D planes so project all points onto those
 744         const pointField& oldPoints = mesh_.points();
 745         const edgeList& edges = mesh_.edges();
 746
 747         const labelList& neIndices = twoDCorrector().normalEdgeIndices();
 748         const vector& pn = twoDCorrector().planeNormal();
 749
 750         forAll(neIndices, i)
 751         {
 752             const edge& e = edges[neIndices[i]];
 753
 754             point& pStart = newPoints[e.start()];
 755
 756             pStart += pn*(pn & (oldPoints[e.start()] - pStart));
 757
 758             point& pEnd = newPoints[e.end()];
 759
 760             pEnd += pn*(pn & (oldPoints[e.end()] - pEnd));
 761         }
 762
 763         // Correct tangentially
 764         twoDCorrector_.correctPoints(newPoints);
 765         Info<< " ...done" << endl;
 766     }
 767
 768     if (debug)
 769     {
 770         Pout<< "motionSmoother::movePoints : testing sync of newPoints."
 771             << endl;
 772         testSyncField
 773         (
 774             newPoints,
 775             minEqOp<point>(),           // combine op
 776             vector(GREAT,GREAT,GREAT),  // null
 777             true                        // separation
 778         );
 779     }
 780
 781     tmp<scalarField> tsweptVol = mesh_.movePoints(newPoints);
 782
 783     //!!! Workaround for movePoints bug
 784     const_cast<polyBoundaryMesh&>(mesh_.boundaryMesh()).movePoints(newPoints);
 785
 786     pp_.movePoints(mesh_.points());
 787
 788     return tsweptVol;
 789 }
 790
 791
 792 Foam::scalar Foam::motionSmoother::setErrorReduction
 793 (
 794     const scalar errorReduction
 795 )
 796 {
 797     scalar oldErrorReduction = readScalar(paramDict_.lookup("errorReduction"));
 798
 799     paramDict_.remove("errorReduction");
 800     paramDict_.add("errorReduction", errorReduction);
 801
 802     return oldErrorReduction;
 803 }
 804
 805
 806 bool Foam::motionSmoother::scaleMesh
 807 (
 808     labelList& checkFaces,
 809     const bool smoothMesh,
 810     const label nAllowableErrors
 811 )
 812 {
 813     List<labelPair> emptyBaffles;
 814     return scaleMesh
 815     (
 816         checkFaces,
 817         emptyBaffles,
 818         smoothMesh,
 819         nAllowableErrors
 820     );
 821 }
 822
 823
 824 bool Foam::motionSmoother::scaleMesh
 825 (
 826     labelList& checkFaces,
 827     const List<labelPair>& baffles,
 828     const bool smoothMesh,
 829     const label nAllowableErrors
 830 )
 831 {
 832     if (!smoothMesh && adaptPatchIDs_.size() == 0)
 833     {
 834         FatalErrorIn("motionSmoother::scaleMesh(const bool")
 835             << "You specified both no movement on the internal mesh points"
 836             << " (smoothMesh = false)" << nl
 837             << "and no movement on the patch (adaptPatchIDs is empty)" << nl
 838             << "Hence nothing to adapt."
 839             << exit(FatalError);
 840     }
 841
 842     if (debug)
 843     {
 844         // Had a problem with patches moved non-synced. Check transformations.
 845         const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 846
 847         Pout<< "Entering scaleMesh : coupled patches:" << endl;
 848         forAll(patches, patchI)
 849         {
 850             if (patches[patchI].coupled())
 851             {
 852                 const coupledPolyPatch& pp =
 853                     refCast<const coupledPolyPatch>(patches[patchI]);
 854
 855                 Pout<< '\t' << patchI << '\t' << pp.name()
 856                     << " parallel:" << pp.parallel()
 857                     << " separated:" << pp.separated()
 858                     << " forwardT:" << pp.forwardT().size()
 859                     << endl;
 860             }
 861         }
 862     }
 863
 864     const scalar errorReduction =
 865         readScalar(paramDict_.lookup("errorReduction"));
 866     const label nSmoothScale =
 867         readLabel(paramDict_.lookup("nSmoothScale"));
 868
 869
 870     // Note: displacement_ should already be synced already from setDisplacement
 871     // but just to make sure.
 872     syncTools::syncPointList
 873     (
 874         mesh_,
 875         displacement_,
 876         maxMagEqOp(),
 877         vector::zero,   // null value
 878         false           // no separation
 879     );
 880
 881     // Set newPoints as old + scale*displacement
 882     pointField newPoints;
 883     {
 884         // Create overall displacement with same b.c.s as displacement_
 885         pointVectorField totalDisplacement
 886         (
 887             IOobject
 888             (
 889                 "totalDisplacement",
 890                 mesh_.time().timeName(),
 891                 mesh_,
 892                 IOobject::NO_READ,
 893                 IOobject::NO_WRITE,
 894                 false
 895             ),
 896             scale_*displacement_,
 897             displacement_.boundaryField().types()
 898         );
 899         correctBoundaryConditions(totalDisplacement);
 900
 901         if (debug)
 902         {
 903             Pout<< "scaleMesh : testing sync of totalDisplacement" << endl;
 904             testSyncField
 905             (
 906                 totalDisplacement,
 907                 maxMagEqOp(),
 908                 vector::zero,   // null value
 909                 false           // separation
 910             );
 911         }
 912
 913         newPoints = oldPoints_ + totalDisplacement.internalField();
 914     }
 915
 916     Info<< "Moving mesh using diplacement scaling :"
 917         << " min:" << gMin(scale_.internalField())
 918         << "  max:" << gMax(scale_.internalField())
 919         << endl;
 920
 921
 922     // Move
 923     movePoints(newPoints);
 924
 925     // Check. Returns parallel number of incorrect faces.
 926     faceSet wrongFaces(mesh_, "wrongFaces", mesh_.nFaces()/100+100);
 927     checkMesh(false, mesh_, paramDict_, checkFaces, baffles, wrongFaces);
 928
 929     if (returnReduce(wrongFaces.size(), sumOp<label>()) <= nAllowableErrors)
 930     {
 931         return true;
 932     }
 933     else
 934     {
 935         // Sync across coupled faces by extending the set.
 936         wrongFaces.sync(mesh_);
 937
 938         // Special case:
 939         // if errorReduction is set to zero, extend wrongFaces
 940         // to face-Cell-faces to ensure quick return to previously valid mesh
 941
 942         if (mag(errorReduction) < SMALL)
 943         {
 944             labelHashSet newWrongFaces(wrongFaces);
 945             forAllConstIter(labelHashSet, wrongFaces, iter)
 946             {
 947                 label own = mesh_.faceOwner()[iter.key()];
 948                 const cell& ownFaces = mesh_.cells()[own];
 949
 950                 forAll(ownFaces, cfI)
 951                 {
 952                     newWrongFaces.insert(ownFaces[cfI]);
 953                 }
 954
 955                 if (iter.key() < mesh_.nInternalFaces())
 956                 {
 957                     label nei = mesh_.faceNeighbour()[iter.key()];
 958                     const cell& neiFaces = mesh_.cells()[nei];
 959
 960                     forAll(neiFaces, cfI)
 961                     {
 962                         newWrongFaces.insert(neiFaces[cfI]);
 963                     }
 964                 }
 965             }
 966             wrongFaces.transfer(newWrongFaces);
 967             wrongFaces.sync(mesh_);
 968         }
 969
 970
 971         // Find out points used by wrong faces and scale displacement.
 972         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 973
 974         pointSet usedPoints(mesh_, "usedPoints", getPoints(wrongFaces));
 975         usedPoints.sync(mesh_);
 976
 977
 978
 979         // Grow a few layers to determine
 980         // - points to be smoothed
 981         // - faces to be checked in next iteration
 982         PackedList<1> isAffectedPoint(mesh_.nPoints(), 0);
 983         getAffectedFacesAndPoints
 984         (
 985             nSmoothScale,       // smoothing iterations
 986             wrongFaces,         // error faces
 987             checkFaces,
 988             isAffectedPoint
 989         );
 990
 991         if (debug)
 992         {
 993             Pout<< "Faces in error:" << wrongFaces.size()
 994                 << "  with points:" << usedPoints.size()
 995                 << endl;
 996         }
 997
 998         if (adaptPatchIDs_.size() != 0)
 999         {
1000             // Scale conflicting patch points
1001             scaleField(pp_.meshPoints(), usedPoints, errorReduction, scale_);
1002         }
1003         if (smoothMesh)
1004         {
1005             // Scale conflicting internal points
1006             scaleField(usedPoints, errorReduction, scale_);
1007         }
1008
1009         for (label i = 0; i < nSmoothScale; i++)
1010         {
1011             if (adaptPatchIDs_.size() != 0)
1012             {
1013                 // Smooth patch values
1014                 pointScalarField oldScale(scale_);
1015                 minSmooth
1016                 (
1017                     isAffectedPoint,
1018                     pp_.meshPoints(),
1019                     oldScale,
1020                     scale_
1021                 );
1022                 checkFld(scale_);
1023             }
1024             if (smoothMesh)
1025             {
1026                 // Smooth internal values
1027                 pointScalarField oldScale(scale_);
1028                 minSmooth(isAffectedPoint, oldScale, scale_);
1029                 checkFld(scale_);
1030             }
1031         }
1032
1033         syncTools::syncPointList
1034         (
1035             mesh_,
1036             scale_,
1037             maxEqOp<scalar>(),
1038             -GREAT,             // null value
1039             false               // no separation
1040         );
1041
1042
1043         if (debug)
1044         {
1045             Pout<< "scale_ after smoothing :"
1046                 << " min:" << Foam::gMin(scale_)
1047                 << " max:" << Foam::gMax(scale_)
1048                 << endl;
1049         }
1050
1051         return false;
1052     }
1053 }
1054
1055
1056 void Foam::motionSmoother::updateMesh()
1057 {
1058     const pointBoundaryMesh& patches = pMesh_.boundary();
1059
1060     // Check whether displacement has fixed value b.c. on adaptPatchID
1061     forAll(adaptPatchIDs_, i)
1062     {
1063         label patchI = adaptPatchIDs_[i];
1064
1065         if
1066         (
1067            !isA<fixedValuePointPatchVectorField>
1068             (
1069                 displacement_.boundaryField()[patchI]
1070             )
1071         )
1072         {
1073             FatalErrorIn
1074             (
1075                 "motionSmoother::motionSmoother"
1076             )   << "Patch " << patches[patchI].name()
1077                 << " has wrong boundary condition "
1078                 << displacement_.boundaryField()[patchI].type()
1079                 << " on field " << displacement_.name() << nl
1080                 << "Only type allowed is "
1081                 << fixedValuePointPatchVectorField::typeName
1082                 << exit(FatalError);
1083         }
1084     }
1085
1086
1087     // Determine internal points. Note that for twoD there are no internal
1088     // points so we use the points of adaptPatchIDs instead
1089
1090     twoDCorrector_.updateMesh();
1091
1092     const labelList& meshPoints = pp_.meshPoints();
1093
1094     forAll(meshPoints, i)
1095     {
1096         isInternalPoint_.set(meshPoints[i], 0);
1097     }
1098
1099     // Calculate master edge addressing
1100     isMasterEdge_ = syncTools::getMasterEdges(mesh_);
1101
1102     makePatchPatchAddressing();
1103 }
1104
1105
1106 // Specialisation of applyCornerConstraints for scalars because
1107 // no constraint need be applied
1108 template<>
1109 void Foam::motionSmoother::applyCornerConstraints<Foam::scalar>
1110 (
1111     GeometricField<scalar, pointPatchField, pointMesh>& pf
1112 ) const
1113 {}
1114
1115
1116 // ************************************************************************* //