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

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