src/dynamicMesh/motionSmoother/motionSmoother.C

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