applications/utilities/mesh/manipulation/createPatch/createPatch.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright held by original author
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of OpenFOAM.
  10
  11     OpenFOAM is free software; you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by the
  13     Free Software Foundation; either version 2 of the License, or (at your
  14     option) any later version.
  15
  16     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with OpenFOAM; if not, write to the Free Software Foundation,
  23     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
  24
  25 Description
  26     Utility to create patches out of selected boundary faces. Faces come either
  27     from existing patches or from a faceSet.
  28
  29     More specifically it:
  30     - creates new patches (from selected boundary faces). Synchronise faces
  31       on coupled patches.
  32     - synchronises points on coupled boundaries
  33     - remove patches with 0 faces in them
  34
  35 \*---------------------------------------------------------------------------*/
  36
  37 #include "cyclicPolyPatch.H"
  38 #include "syncTools.H"
  39 #include "argList.H"
  40 #include "polyMesh.H"
  41 #include "Time.H"
  42 #include "SortableList.H"
  43 #include "OFstream.H"
  44 #include "meshTools.H"
  45 #include "faceSet.H"
  46 #include "IOPtrList.H"
  47 #include "mapPolyMesh.H"
  48 #include "directTopoChange.H"
  49 #include "polyModifyFace.H"
  50
  51 using namespace Foam;
  52
  53 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  54
  55 namespace Foam
  56 {
  57     defineTemplateTypeNameAndDebug(IOPtrList<dictionary>, 0);
  58 }
  59
  60 // Combine operator to synchronise points. We choose point nearest to origin so
  61 // we can use e.g. great,great,great as null value.
  62 class nearestEqOp
  63 {
  64
  65 public:
  66
  67     void operator()(vector& x, const vector& y) const
  68     {
  69         if (magSqr(y) < magSqr(x))
  70         {
  71             x = y;
  72         }
  73     }
  74 };
  75
  76
  77 void changePatchID
  78 (
  79     const polyMesh& mesh,
  80     const label faceID,
  81     const label patchID,
  82     directTopoChange& meshMod
  83 )
  84 {
  85     const label zoneID = mesh.faceZones().whichZone(faceID);
  86
  87     bool zoneFlip = false;
  88
  89     if (zoneID >= 0)
  90     {
  91         const faceZone& fZone = mesh.faceZones()[zoneID];
  92
  93         zoneFlip = fZone.flipMap()[fZone.whichFace(faceID)];
  94     }
  95
  96     meshMod.setAction
  97     (
  98         polyModifyFace
  99         (
 100             mesh.faces()[faceID],               // face
 101             faceID,                             // face ID
 102             mesh.faceOwner()[faceID],           // owner
 103             -1,                                 // neighbour
 104             false,                              // flip flux
 105             patchID,                            // patch ID
 106             false,                              // remove from zone
 107             zoneID,                             // zone ID
 108             zoneFlip                            // zone flip
 109         )
 110     );
 111 }
 112
 113
 114 // Filter out the empty patches.
 115 void filterPatches(polyMesh& mesh)
 116 {
 117     const polyBoundaryMesh& patches = mesh.boundaryMesh();
 118
 119     // Patches to keep
 120     DynamicList<polyPatch*> allPatches(patches.size());
 121
 122     label nOldPatches = returnReduce(patches.size(), sumOp<label>());
 123
 124     // Copy old patches.
 125     forAll(patches, patchI)
 126     {
 127         const polyPatch& pp = patches[patchI];
 128
 129         // Note: reduce possible since non-proc patches guaranteed in same order
 130         if (!isA<processorPolyPatch>(pp))
 131         {
 132             if (returnReduce(pp.size(), sumOp<label>()) > 0)
 133             {
 134                 allPatches.append
 135                 (
 136                     pp.clone
 137                     (
 138                         patches,
 139                         allPatches.size(),
 140                         pp.size(),
 141                         pp.start()
 142                     ).ptr()
 143                 );
 144             }
 145             else
 146             {
 147                 Info<< "Removing empty patch " << pp.name() << " at position "
 148                     << patchI << endl;
 149             }
 150         }
 151     }
 152     // Copy non-empty processor patches
 153     forAll(patches, patchI)
 154     {
 155         const polyPatch& pp = patches[patchI];
 156
 157         if (isA<processorPolyPatch>(pp))
 158         {
 159             if (pp.size())
 160             {
 161                 allPatches.append
 162                 (
 163                     pp.clone
 164                     (
 165                         patches,
 166                         allPatches.size(),
 167                         pp.size(),
 168                         pp.start()
 169                     ).ptr()
 170                 );
 171             }
 172             else
 173             {
 174                 Info<< "Removing empty processor patch " << pp.name()
 175                     << " at position " << patchI << endl;
 176             }
 177         }
 178     }
 179
 180     label nAllPatches = returnReduce(allPatches.size(), sumOp<label>());
 181     if (nAllPatches != nOldPatches)
 182     {
 183         Info<< "Removing patches." << endl;
 184         allPatches.shrink();
 185         mesh.removeBoundary();
 186         mesh.addPatches(allPatches);
 187     }
 188     else
 189     {
 190         Info<< "No patches removed." << endl;
 191     }
 192 }
 193
 194
 195 // Dump for all patches the current match
 196 void dumpCyclicMatch(const fileName& prefix, const polyMesh& mesh)
 197 {
 198     const polyBoundaryMesh& patches = mesh.boundaryMesh();
 199
 200     forAll(patches, patchI)
 201     {
 202         if (isA<cyclicPolyPatch>(patches[patchI]))
 203         {
 204             const cyclicPolyPatch& cycPatch =
 205                 refCast<const cyclicPolyPatch>(patches[patchI]);
 206
 207             label halfSize = cycPatch.size()/2;
 208
 209             // Dump halves
 210             {
 211                 OFstream str(prefix+cycPatch.name()+"_half0.obj");
 212                 Pout<< "Dumping " << cycPatch.name()
 213                     << " half0 faces to " << str.name() << endl;
 214                 meshTools::writeOBJ
 215                 (
 216                     str,
 217                     static_cast<faceList>
 218                     (
 219                         SubList<face>
 220                         (
 221                             cycPatch,
 222                             halfSize
 223                         )
 224                     ),
 225                     cycPatch.points()
 226                 );
 227             }
 228             {
 229                 OFstream str(prefix+cycPatch.name()+"_half1.obj");
 230                 Pout<< "Dumping " << cycPatch.name()
 231                     << " half1 faces to " << str.name() << endl;
 232                 meshTools::writeOBJ
 233                 (
 234                     str,
 235                     static_cast<faceList>
 236                     (
 237                         SubList<face>
 238                         (
 239                             cycPatch,
 240                             halfSize,
 241                             halfSize
 242                         )
 243                     ),
 244                     cycPatch.points()
 245                 );
 246             }
 247
 248
 249             // Lines between corresponding face centres
 250             OFstream str(prefix+cycPatch.name()+"_match.obj");
 251             label vertI = 0;
 252
 253             Pout<< "Dumping cyclic match as lines between face centres to "
 254                 << str.name() << endl;
 255
 256             for (label faceI = 0; faceI < halfSize; faceI++)
 257             {
 258                 const point& fc0 = mesh.faceCentres()[cycPatch.start()+faceI];
 259                 meshTools::writeOBJ(str, fc0);
 260                 vertI++;
 261
 262                 label nbrFaceI = halfSize + faceI;
 263                 const point& fc1 =
 264                     mesh.faceCentres()[cycPatch.start()+nbrFaceI];
 265                 meshTools::writeOBJ(str, fc1);
 266                 vertI++;
 267
 268                 str<< "l " << vertI-1 << ' ' << vertI << nl;
 269             }
 270         }
 271     }
 272 }
 273
 274
 275 void separateList
 276 (
 277     const vectorField& separation,
 278     UList<vector>& field
 279 )
 280 {
 281     if (separation.size() == 1)
 282     {
 283         // Single value for all.
 284
 285         forAll(field, i)
 286         {
 287             field[i] += separation[0];
 288         }
 289     }
 290     else if (separation.size() == field.size())
 291     {
 292         forAll(field, i)
 293         {
 294             field[i] += separation[i];
 295         }
 296     }
 297     else
 298     {
 299         FatalErrorIn
 300         (
 301             "separateList(const vectorField&, UList<vector>&)"
 302         )   << "Sizes of field and transformation not equal. field:"
 303             << field.size() << " transformation:" << separation.size()
 304             << abort(FatalError);
 305     }
 306 }
 307
 308
 309 // Synchronise points on both sides of coupled boundaries.
 310 template <class CombineOp>
 311 void syncPoints
 312 (
 313     const polyMesh& mesh,
 314     pointField& points,
 315     const CombineOp& cop,
 316     const point& nullValue
 317 )
 318 {
 319     if (points.size() != mesh.nPoints())
 320     {
 321         FatalErrorIn
 322         (
 323             "syncPoints"
 324             "(const polyMesh&, pointField&, const CombineOp&, const point&)"
 325         )   << "Number of values " << points.size()
 326             << " is not equal to the number of points in the mesh "
 327             << mesh.nPoints() << abort(FatalError);
 328     }
 329
 330     const polyBoundaryMesh& patches = mesh.boundaryMesh();
 331
 332     // Is there any coupled patch with transformation?
 333     bool hasTransformation = false;
 334
 335     if (Pstream::parRun())
 336     {
 337         // Send
 338
 339         forAll(patches, patchI)
 340         {
 341             const polyPatch& pp = patches[patchI];
 342
 343             if
 344             (
 345                 isA<processorPolyPatch>(pp)
 346              && pp.nPoints() > 0
 347              && refCast<const processorPolyPatch>(pp).owner()
 348             )
 349             {
 350                 const processorPolyPatch& procPatch =
 351                     refCast<const processorPolyPatch>(pp);
 352
 353                 // Get data per patchPoint in neighbouring point numbers.
 354                 pointField patchInfo(procPatch.nPoints(), nullValue);
 355
 356                 const labelList& meshPts = procPatch.meshPoints();
 357                 const labelList& nbrPts = procPatch.neighbPoints();
 358
 359                 forAll(nbrPts, pointI)
 360                 {
 361                     label nbrPointI = nbrPts[pointI];
 362                     if (nbrPointI >= 0 && nbrPointI < patchInfo.size())
 363                     {
 364                         patchInfo[nbrPointI] = points[meshPts[pointI]];
 365                     }
 366                 }
 367
 368                 OPstream toNbr(Pstream::blocking, procPatch.neighbProcNo());
 369                 toNbr << patchInfo;
 370             }
 371         }
 372
 373
 374         // Receive and set.
 375
 376         forAll(patches, patchI)
 377         {
 378             const polyPatch& pp = patches[patchI];
 379
 380             if
 381             (
 382                 isA<processorPolyPatch>(pp)
 383              && pp.nPoints() > 0
 384              && !refCast<const processorPolyPatch>(pp).owner()
 385             )
 386             {
 387                 const processorPolyPatch& procPatch =
 388                     refCast<const processorPolyPatch>(pp);
 389
 390                 pointField nbrPatchInfo(procPatch.nPoints());
 391                 {
 392                     // We do not know the number of points on the other side
 393                     // so cannot use Pstream::read.
 394                     IPstream fromNbr
 395                     (
 396                         Pstream::blocking,
 397                         procPatch.neighbProcNo()
 398                     );
 399                     fromNbr >> nbrPatchInfo;
 400                 }
 401                 // Null any value which is not on neighbouring processor
 402                 nbrPatchInfo.setSize(procPatch.nPoints(), nullValue);
 403
 404                 if (!procPatch.parallel())
 405                 {
 406                     hasTransformation = true;
 407                     transformList(procPatch.forwardT(), nbrPatchInfo);
 408                 }
 409                 else if (procPatch.separated())
 410                 {
 411                     hasTransformation = true;
 412                     separateList(-procPatch.separation(), nbrPatchInfo);
 413                 }
 414
 415                 const labelList& meshPts = procPatch.meshPoints();
 416
 417                 forAll(meshPts, pointI)
 418                 {
 419                     label meshPointI = meshPts[pointI];
 420                     points[meshPointI] = nbrPatchInfo[pointI];
 421                 }
 422             }
 423         }
 424     }
 425
 426     // Do the cyclics.
 427     forAll(patches, patchI)
 428     {
 429         const polyPatch& pp = patches[patchI];
 430
 431         if (isA<cyclicPolyPatch>(pp))
 432         {
 433             const cyclicPolyPatch& cycPatch =
 434                 refCast<const cyclicPolyPatch>(pp);
 435
 436             const edgeList& coupledPoints = cycPatch.coupledPoints();
 437             const labelList& meshPts = cycPatch.meshPoints();
 438
 439             pointField half0Values(coupledPoints.size());
 440
 441             forAll(coupledPoints, i)
 442             {
 443                 const edge& e = coupledPoints[i];
 444                 label point0 = meshPts[e[0]];
 445                 half0Values[i] = points[point0];
 446             }
 447
 448             if (!cycPatch.parallel())
 449             {
 450                 hasTransformation = true;
 451                 transformList(cycPatch.reverseT(), half0Values);
 452             }
 453             else if (cycPatch.separated())
 454             {
 455                 hasTransformation = true;
 456                 const vectorField& v = cycPatch.coupledPolyPatch::separation();
 457                 separateList(v, half0Values);
 458             }
 459
 460             forAll(coupledPoints, i)
 461             {
 462                 const edge& e = coupledPoints[i];
 463                 label point1 = meshPts[e[1]];
 464                 points[point1] = half0Values[i];
 465             }
 466         }
 467     }
 468
 469     //- Note: hasTransformation is only used for warning messages so
 470     //  reduction not strictly nessecary.
 471     //reduce(hasTransformation, orOp<bool>());
 472
 473     // Synchronize multiple shared points.
 474     const globalMeshData& pd = mesh.globalData();
 475
 476     if (pd.nGlobalPoints() > 0)
 477     {
 478         if (hasTransformation)
 479         {
 480             WarningIn
 481             (
 482                 "syncPoints"
 483                 "(const polyMesh&, pointField&, const CombineOp&, const point&)"
 484             )   << "There are decomposed cyclics in this mesh with"
 485                 << " transformations." << endl
 486                 << "This is not supported. The result will be incorrect"
 487                 << endl;
 488         }
 489
 490
 491         // Values on shared points.
 492         pointField sharedPts(pd.nGlobalPoints(), nullValue);
 493
 494         forAll(pd.sharedPointLabels(), i)
 495         {
 496             label meshPointI = pd.sharedPointLabels()[i];
 497             // Fill my entries in the shared points
 498             sharedPts[pd.sharedPointAddr()[i]] = points[meshPointI];
 499         }
 500
 501         // Combine on master.
 502         Pstream::listCombineGather(sharedPts, cop);
 503         Pstream::listCombineScatter(sharedPts);
 504
 505         // Now we will all have the same information. Merge it back with
 506         // my local information.
 507         forAll(pd.sharedPointLabels(), i)
 508         {
 509             label meshPointI = pd.sharedPointLabels()[i];
 510             points[meshPointI] = sharedPts[pd.sharedPointAddr()[i]];
 511         }
 512     }
 513 }
 514
 515
 516 // Main program:
 517
 518 int main(int argc, char *argv[])
 519 {
 520 #   include "addRegionOption.H"
 521     argList::validOptions.insert("overwrite", "");
 522
 523 #   include "setRootCase.H"
 524 #   include "createTime.H"
 525     runTime.functionObjects().off();
 526
 527     Foam::word meshRegionName = polyMesh::defaultRegion;
 528     args.optionReadIfPresent("region", meshRegionName);
 529
 530     const bool overwrite = args.optionFound("overwrite");
 531
 532     Info<< "Reading createPatchDict." << nl << endl;
 533
 534     IOdictionary dict
 535     (
 536         IOobject
 537         (
 538             "createPatchDict",
 539             runTime.system(),
 540             (
 541                 meshRegionName != polyMesh::defaultRegion
 542               ? meshRegionName
 543               : word::null
 544             ),
 545             runTime,
 546             IOobject::MUST_READ,
 547             IOobject::NO_WRITE,
 548             false
 549         )
 550     );
 551
 552
 553     // Whether to synchronise points
 554     const Switch pointSync(dict.lookup("pointSync"));
 555
 556
 557     // Change tolerance in controlDict instead.  HJ, 22/Oct/2008
 558
 559     // Set the matching tolerance so we can read illegal meshes
 560 //     scalar tol = readScalar(dict.lookup("matchTolerance"));
 561 //     Info<< "Using relative tolerance " << tol
 562 //         << " to match up faces and points" << nl << endl;
 563 //      polyPatch::matchTol_ = tol;
 564
 565 #   include "createNamedPolyMesh.H"
 566
 567     const word oldInstance = mesh.pointsInstance();
 568
 569     const polyBoundaryMesh& patches = mesh.boundaryMesh();
 570
 571     // If running parallel check same patches everywhere
 572     patches.checkParallelSync(true);
 573
 574
 575     dumpCyclicMatch("initial_", mesh);
 576
 577     // Read patch construct info from dictionary
 578     PtrList<dictionary> patchSources(dict.lookup("patchInfo"));
 579
 580
 581
 582     // 1. Add all new patches
 583     // ~~~~~~~~~~~~~~~~~~~~~~
 584
 585     if (patchSources.size())
 586     {
 587         // Old and new patches.
 588         DynamicList<polyPatch*> allPatches(patches.size()+patchSources.size());
 589
 590         label startFaceI = mesh.nInternalFaces();
 591
 592         // Copy old patches.
 593         forAll(patches, patchI)
 594         {
 595             const polyPatch& pp = patches[patchI];
 596
 597             if (!isA<processorPolyPatch>(pp))
 598             {
 599                 allPatches.append
 600                 (
 601                     pp.clone
 602                     (
 603                         patches,
 604                         patchI,
 605                         pp.size(),
 606                         startFaceI
 607                     ).ptr()
 608                 );
 609                 startFaceI += pp.size();
 610             }
 611         }
 612
 613         forAll(patchSources, addedI)
 614         {
 615             const dictionary& dict = patchSources[addedI];
 616
 617             word patchName(dict.lookup("name"));
 618
 619             label destPatchI = patches.findPatchID(patchName);
 620
 621             if (destPatchI == -1)
 622             {
 623                 dictionary patchDict(dict.subDict("dictionary"));
 624
 625                 destPatchI = allPatches.size();
 626
 627                 Info<< "Adding new patch " << patchName
 628                     << " as patch " << destPatchI
 629                     << " from " << patchDict << endl;
 630
 631                 patchDict.set("nFaces", 0);
 632                 patchDict.set("startFace", startFaceI);
 633
 634                 // Add an empty patch.
 635                 allPatches.append
 636                 (
 637                     polyPatch::New
 638                     (
 639                         patchName,
 640                         patchDict,
 641                         destPatchI,
 642                         patches
 643                     ).ptr()
 644                 );
 645             }
 646         }
 647
 648         // Copy old patches.
 649         forAll(patches, patchI)
 650         {
 651             const polyPatch& pp = patches[patchI];
 652
 653             if (isA<processorPolyPatch>(pp))
 654             {
 655                 allPatches.append
 656                 (
 657                     pp.clone
 658                     (
 659                         patches,
 660                         patchI,
 661                         pp.size(),
 662                         startFaceI
 663                     ).ptr()
 664                 );
 665                 startFaceI += pp.size();
 666             }
 667         }
 668
 669         allPatches.shrink();
 670         mesh.removeBoundary();
 671         mesh.addPatches(allPatches);
 672
 673         Info<< endl;
 674     }
 675
 676
 677
 678     // 2. Repatch faces
 679     // ~~~~~~~~~~~~~~~~
 680
 681     directTopoChange meshMod(mesh);
 682
 683
 684     forAll(patchSources, addedI)
 685     {
 686         const dictionary& dict = patchSources[addedI];
 687
 688         word patchName(dict.lookup("name"));
 689
 690         label destPatchI = patches.findPatchID(patchName);
 691
 692         if (destPatchI == -1)
 693         {
 694             FatalErrorIn(args.executable()) << "patch " << patchName
 695                 << " not added. Problem." << abort(FatalError);
 696         }
 697
 698         word sourceType(dict.lookup("constructFrom"));
 699
 700         if (sourceType == "patches")
 701         {
 702             labelHashSet patchSources(patches.patchSet(dict.lookup("patches")));
 703
 704             // Repatch faces of the patches.
 705             forAllConstIter(labelHashSet, patchSources, iter)
 706             {
 707                 const polyPatch& pp = patches[iter.key()];
 708
 709                 Info<< "Moving faces from patch " << pp.name()
 710                     << " to patch " << destPatchI << endl;
 711
 712                 forAll(pp, i)
 713                 {
 714                     changePatchID
 715                     (
 716                         mesh,
 717                         pp.start() + i,
 718                         destPatchI,
 719                         meshMod
 720                     );
 721                 }
 722             }
 723         }
 724         else if (sourceType == "set")
 725         {
 726             word setName(dict.lookup("set"));
 727
 728             faceSet faces(mesh, setName);
 729
 730             Info<< "Read " << returnReduce(faces.size(), sumOp<label>())
 731                 << " faces from faceSet " << faces.name() << endl;
 732
 733             // Sort (since faceSet contains faces in arbitrary order)
 734             labelList faceLabels(faces.toc());
 735
 736             SortableList<label> patchFaces(faceLabels);
 737
 738             forAll(patchFaces, i)
 739             {
 740                 label faceI = patchFaces[i];
 741
 742                 if (mesh.isInternalFace(faceI))
 743                 {
 744                     FatalErrorIn(args.executable())
 745                         << "Face " << faceI << " specified in set "
 746                         << faces.name()
 747                         << " is not an external face of the mesh." << endl
 748                         << "This application can only repatch existing boundary"
 749                         << " faces." << exit(FatalError);
 750                 }
 751
 752                 changePatchID
 753                 (
 754                     mesh,
 755                     faceI,
 756                     destPatchI,
 757                     meshMod
 758                 );
 759             }
 760         }
 761         else
 762         {
 763             FatalErrorIn(args.executable())
 764                 << "Invalid source type " << sourceType << endl
 765                 << "Valid source types are 'patches' 'set'" << exit(FatalError);
 766         }
 767     }
 768     Info<< endl;
 769
 770
 771     // Change mesh, use inflation to reforce calculation of transformation
 772     // tensors.
 773     Info<< "Doing topology modification to order faces." << nl << endl;
 774     autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh, true);
 775     mesh.movePoints(map().preMotionPoints());
 776
 777     dumpCyclicMatch("coupled_", mesh);
 778
 779     // Synchronise points.
 780     if (!pointSync)
 781     {
 782         Info<< "Not synchronising points." << nl << endl;
 783     }
 784     else
 785     {
 786         Info<< "Synchronising points." << nl << endl;
 787
 788         // This is a bit tricky. Both normal and position might be out and
 789         // current separation also includes the normal
 790         // ( separation_ = (nf&(Cr - Cf))*nf ).
 791
 792         // For processor patches:
 793         // - disallow multiple separation/transformation. This basically
 794         //   excludes decomposed cyclics. Use the (probably 0) separation
 795         //   to align the points.
 796         // For cyclic patches:
 797         // - for separated ones use our own recalculated offset vector
 798         // - for rotational ones use current one.
 799
 800         forAll(mesh.boundaryMesh(), patchI)
 801         {
 802             const polyPatch& pp = mesh.boundaryMesh()[patchI];
 803
 804             if (pp.size() && isA<coupledPolyPatch>(pp))
 805             {
 806                 const coupledPolyPatch& cpp =
 807                     refCast<const coupledPolyPatch>(pp);
 808
 809                 if (cpp.separated())
 810                 {
 811                     Info<< "On coupled patch " << pp.name()
 812                         << " separation[0] was "
 813                         << cpp.separation()[0] << endl;
 814
 815                     if (isA<cyclicPolyPatch>(pp))
 816                     {
 817                         const cyclicPolyPatch& cycpp =
 818                             refCast<const cyclicPolyPatch>(pp);
 819
 820                         if (cycpp.transform() == cyclicPolyPatch::TRANSLATIONAL)
 821                         {
 822                             Info<< "On cyclic translation patch " << pp.name()
 823                                 << " forcing uniform separation of "
 824                                 << cycpp.separationVector() << endl;
 825                             const_cast<vectorField&>(cpp.separation()) =
 826                                 pointField(1, cycpp.separationVector());
 827                         }
 828                         else
 829                         {
 830                             const_cast<vectorField&>(cpp.separation()) =
 831                                 pointField
 832                                 (
 833                                     1,
 834                                     pp[pp.size()/2].centre(mesh.points())
 835                                   - pp[0].centre(mesh.points())
 836                                 );
 837                         }
 838                     }
 839                     else
 840                     {
 841                         const_cast<vectorField&>(cpp.separation())
 842                         .setSize(1);
 843                     }
 844                     Info<< "On coupled patch " << pp.name()
 845                         << " forcing uniform separation of "
 846                         << cpp.separation() << endl;
 847                 }
 848                 else if (!cpp.parallel())
 849                 {
 850                     Info<< "On coupled patch " << pp.name()
 851                         << " forcing uniform rotation of "
 852                         << cpp.forwardT()[0] << endl;
 853
 854                     const_cast<tensorField&>
 855                     (
 856                         cpp.forwardT()
 857                     ).setSize(1);
 858                     const_cast<tensorField&>
 859                     (
 860                         cpp.reverseT()
 861                     ).setSize(1);
 862
 863                     Info<< "On coupled patch " << pp.name()
 864                         << " forcing uniform rotation of "
 865                         << cpp.forwardT() << endl;
 866                 }
 867             }
 868         }
 869
 870         Info<< "Synchronising points." << endl;
 871
 872         pointField newPoints(mesh.points());
 873
 874         syncPoints
 875         (
 876             mesh,
 877             newPoints,
 878             nearestEqOp(),
 879             point(GREAT, GREAT, GREAT)
 880         );
 881
 882         scalarField diff(mag(newPoints-mesh.points()));
 883         Info<< "Points changed by average:" << gAverage(diff)
 884             << " max:" << gMax(diff) << nl << endl;
 885
 886         mesh.movePoints(newPoints);
 887     }
 888
 889     // 3. Remove zeros-sized patches
 890     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 891
 892     Info<< "Removing patches with no faces in them." << nl<< endl;
 893     filterPatches(mesh);
 894
 895
 896     dumpCyclicMatch("final_", mesh);
 897
 898
 899     // Set the precision of the points data to 10
 900     IOstream::defaultPrecision(10);
 901
 902     if (!overwrite)
 903     {
 904         runTime++;
 905     }
 906     else
 907     {
 908         mesh.setInstance(oldInstance);
 909     }
 910
 911     // Write resulting mesh
 912     Info<< "Writing repatched mesh to " << runTime.timeName() << nl << endl;
 913     mesh.write();
 914
 915     Info<< "End\n" << endl;
 916
 917     return 0;
 918 }
 919
 920
 921 // ************************************************************************* //