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 (C) 1991-2009 OpenCFD Ltd.
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of OpenFOAM.
  10
  11     OpenFOAM is free software; you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by the
  13     Free Software Foundation; either version 2 of the License, or (at your
  14     option) any later version.
  15
  16     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with OpenFOAM; if not, write to the Free Software Foundation,
  23     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
  24
  25 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 "polyTopoChange.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     polyTopoChange& 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 #   include "createNamedPolyMesh.H"
 526
 527     const bool overwrite = args.optionFound("overwrite");
 528
 529     Info<< "Reading createPatchDict." << nl << endl;
 530
 531     IOdictionary dict
 532     (
 533         IOobject
 534         (
 535             "createPatchDict",
 536             runTime.system(),
 537             (
 538                 regionName != polyMesh::defaultRegion
 539               ? regionName
 540               : word::null
 541             ),
 542             runTime,
 543             IOobject::MUST_READ,
 544             IOobject::NO_WRITE,
 545             false
 546         )
 547     );
 548
 549
 550     // Whether to synchronise points
 551     const Switch pointSync(dict.lookup("pointSync"));
 552
 553
 554     // Set the matching tolerance so we can read illegal meshes
 555     scalar tol = readScalar(dict.lookup("matchTolerance"));
 556     Info<< "Using relative tolerance " << tol
 557         << " to match up faces and points" << nl << endl;
 558     coupledPolyPatch::matchTol = tol;
 559
 560
 561     const word oldInstance = mesh.pointsInstance();
 562
 563     const polyBoundaryMesh& patches = mesh.boundaryMesh();
 564
 565     // If running parallel check same patches everywhere
 566     patches.checkParallelSync(true);
 567
 568
 569     dumpCyclicMatch("initial_", mesh);
 570
 571     // Read patch construct info from dictionary
 572     PtrList<dictionary> patchSources(dict.lookup("patchInfo"));
 573
 574
 575
 576     // 1. Add all new patches
 577     // ~~~~~~~~~~~~~~~~~~~~~~
 578
 579     if (patchSources.size())
 580     {
 581         // Old and new patches.
 582         DynamicList<polyPatch*> allPatches(patches.size()+patchSources.size());
 583
 584         label startFaceI = mesh.nInternalFaces();
 585
 586         // Copy old patches.
 587         forAll(patches, patchI)
 588         {
 589             const polyPatch& pp = patches[patchI];
 590
 591             if (!isA<processorPolyPatch>(pp))
 592             {
 593                 allPatches.append
 594                 (
 595                     pp.clone
 596                     (
 597                         patches,
 598                         patchI,
 599                         pp.size(),
 600                         startFaceI
 601                     ).ptr()
 602                 );
 603                 startFaceI += pp.size();
 604             }
 605         }
 606
 607         forAll(patchSources, addedI)
 608         {
 609             const dictionary& dict = patchSources[addedI];
 610
 611             word patchName(dict.lookup("name"));
 612
 613             label destPatchI = patches.findPatchID(patchName);
 614
 615             if (destPatchI == -1)
 616             {
 617                 dictionary patchDict(dict.subDict("dictionary"));
 618
 619                 destPatchI = allPatches.size();
 620
 621                 Info<< "Adding new patch " << patchName
 622                     << " as patch " << destPatchI
 623                     << " from " << patchDict << endl;
 624
 625                 patchDict.set("nFaces", 0);
 626                 patchDict.set("startFace", startFaceI);
 627
 628                 // Add an empty patch.
 629                 allPatches.append
 630                 (
 631                     polyPatch::New
 632                     (
 633                         patchName,
 634                         patchDict,
 635                         destPatchI,
 636                         patches
 637                     ).ptr()
 638                 );
 639             }
 640         }
 641
 642         // Copy old patches.
 643         forAll(patches, patchI)
 644         {
 645             const polyPatch& pp = patches[patchI];
 646
 647             if (isA<processorPolyPatch>(pp))
 648             {
 649                 allPatches.append
 650                 (
 651                     pp.clone
 652                     (
 653                         patches,
 654                         patchI,
 655                         pp.size(),
 656                         startFaceI
 657                     ).ptr()
 658                 );
 659                 startFaceI += pp.size();
 660             }
 661         }
 662
 663         allPatches.shrink();
 664         mesh.removeBoundary();
 665         mesh.addPatches(allPatches);
 666
 667         Info<< endl;
 668     }
 669
 670
 671
 672     // 2. Repatch faces
 673     // ~~~~~~~~~~~~~~~~
 674
 675     polyTopoChange meshMod(mesh);
 676
 677
 678     forAll(patchSources, addedI)
 679     {
 680         const dictionary& dict = patchSources[addedI];
 681
 682         word patchName(dict.lookup("name"));
 683
 684         label destPatchI = patches.findPatchID(patchName);
 685
 686         if (destPatchI == -1)
 687         {
 688             FatalErrorIn(args.executable()) << "patch " << patchName
 689                 << " not added. Problem." << abort(FatalError);
 690         }
 691
 692         word sourceType(dict.lookup("constructFrom"));
 693
 694         if (sourceType == "patches")
 695         {
 696             labelHashSet patchSources(patches.patchSet(dict.lookup("patches")));
 697
 698             // Repatch faces of the patches.
 699             forAllConstIter(labelHashSet, patchSources, iter)
 700             {
 701                 const polyPatch& pp = patches[iter.key()];
 702
 703                 Info<< "Moving faces from patch " << pp.name()
 704                     << " to patch " << destPatchI << endl;
 705
 706                 forAll(pp, i)
 707                 {
 708                     changePatchID
 709                     (
 710                         mesh,
 711                         pp.start() + i,
 712                         destPatchI,
 713                         meshMod
 714                     );
 715                 }
 716             }
 717         }
 718         else if (sourceType == "set")
 719         {
 720             word setName(dict.lookup("set"));
 721
 722             faceSet faces(mesh, setName);
 723
 724             Info<< "Read " << returnReduce(faces.size(), sumOp<label>())
 725                 << " faces from faceSet " << faces.name() << endl;
 726
 727             // Sort (since faceSet contains faces in arbitrary order)
 728             labelList faceLabels(faces.toc());
 729
 730             SortableList<label> patchFaces(faceLabels);
 731
 732             forAll(patchFaces, i)
 733             {
 734                 label faceI = patchFaces[i];
 735
 736                 if (mesh.isInternalFace(faceI))
 737                 {
 738                     FatalErrorIn(args.executable())
 739                         << "Face " << faceI << " specified in set "
 740                         << faces.name()
 741                         << " is not an external face of the mesh." << endl
 742                         << "This application can only repatch existing boundary"
 743                         << " faces." << exit(FatalError);
 744                 }
 745
 746                 changePatchID
 747                 (
 748                     mesh,
 749                     faceI,
 750                     destPatchI,
 751                     meshMod
 752                 );
 753             }
 754         }
 755         else
 756         {
 757             FatalErrorIn(args.executable())
 758                 << "Invalid source type " << sourceType << endl
 759                 << "Valid source types are 'patches' 'set'" << exit(FatalError);
 760         }
 761     }
 762     Info<< endl;
 763
 764
 765     // Change mesh, use inflation to reforce calculation of transformation
 766     // tensors.
 767     Info<< "Doing topology modification to order faces." << nl << endl;
 768     autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh, true);
 769     mesh.movePoints(map().preMotionPoints());
 770
 771     dumpCyclicMatch("coupled_", mesh);
 772
 773     // Synchronise points.
 774     if (!pointSync)
 775     {
 776         Info<< "Not synchronising points." << nl << endl;
 777     }
 778     else
 779     {
 780         Info<< "Synchronising points." << nl << endl;
 781
 782         // This is a bit tricky. Both normal and position might be out and
 783         // current separation also includes the normal
 784         // ( separation_ = (nf&(Cr - Cf))*nf ).
 785
 786         // For processor patches:
 787         // - disallow multiple separation/transformation. This basically
 788         //   excludes decomposed cyclics. Use the (probably 0) separation
 789         //   to align the points.
 790         // For cyclic patches:
 791         // - for separated ones use our own recalculated offset vector
 792         // - for rotational ones use current one.
 793
 794         forAll(mesh.boundaryMesh(), patchI)
 795         {
 796             const polyPatch& pp = mesh.boundaryMesh()[patchI];
 797
 798             if (pp.size() && isA<coupledPolyPatch>(pp))
 799             {
 800                 const coupledPolyPatch& cpp =
 801                     refCast<const coupledPolyPatch>(pp);
 802
 803                 if (cpp.separated())
 804                 {
 805                     Info<< "On coupled patch " << pp.name()
 806                         << " separation[0] was "
 807                         << cpp.separation()[0] << endl;
 808
 809                     if (isA<cyclicPolyPatch>(pp))
 810                     {
 811                         const cyclicPolyPatch& cycpp =
 812                             refCast<const cyclicPolyPatch>(pp);
 813
 814                         if (cycpp.transform() == cyclicPolyPatch::TRANSLATIONAL)
 815                         {
 816                             Info<< "On cyclic translation patch " << pp.name()
 817                                 << " forcing uniform separation of "
 818                                 << cycpp.separationVector() << endl;
 819                             const_cast<vectorField&>(cpp.separation()) =
 820                                 pointField(1, cycpp.separationVector());
 821                         }
 822                         else
 823                         {
 824                             const_cast<vectorField&>(cpp.separation()) =
 825                                 pointField
 826                                 (
 827                                     1,
 828                                     pp[pp.size()/2].centre(mesh.points())
 829                                   - pp[0].centre(mesh.points())
 830                                 );
 831                         }
 832                     }
 833                     else
 834                     {
 835                         const_cast<vectorField&>(cpp.separation())
 836                         .setSize(1);
 837                     }
 838                     Info<< "On coupled patch " << pp.name()
 839                         << " forcing uniform separation of "
 840                         << cpp.separation() << endl;
 841                 }
 842                 else if (!cpp.parallel())
 843                 {
 844                     Info<< "On coupled patch " << pp.name()
 845                         << " forcing uniform rotation of "
 846                         << cpp.forwardT()[0] << endl;
 847
 848                     const_cast<tensorField&>
 849                     (
 850                         cpp.forwardT()
 851                     ).setSize(1);
 852                     const_cast<tensorField&>
 853                     (
 854                         cpp.reverseT()
 855                     ).setSize(1);
 856
 857                     Info<< "On coupled patch " << pp.name()
 858                         << " forcing uniform rotation of "
 859                         << cpp.forwardT() << endl;
 860                 }
 861             }
 862         }
 863
 864         Info<< "Synchronising points." << endl;
 865
 866         pointField newPoints(mesh.points());
 867
 868         syncPoints
 869         (
 870             mesh,
 871             newPoints,
 872             nearestEqOp(),
 873             point(GREAT, GREAT, GREAT)
 874         );
 875
 876         scalarField diff(mag(newPoints-mesh.points()));
 877         Info<< "Points changed by average:" << gAverage(diff)
 878             << " max:" << gMax(diff) << nl << endl;
 879
 880         mesh.movePoints(newPoints);
 881     }
 882
 883     // 3. Remove zeros-sized patches
 884     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 885
 886     Info<< "Removing patches with no faces in them." << nl<< endl;
 887     filterPatches(mesh);
 888
 889
 890     dumpCyclicMatch("final_", mesh);
 891
 892
 893     // Set the precision of the points data to 10
 894     IOstream::defaultPrecision(10);
 895
 896     if (!overwrite)
 897     {
 898         runTime++;
 899     }
 900     else
 901     {
 902         mesh.setInstance(oldInstance);
 903     }
 904
 905     // Write resulting mesh
 906     Info<< "Writing repatched mesh to " << runTime.timeName() << nl << endl;
 907     mesh.write();
 908
 909     Info<< "End\n" << endl;
 910
 911     return 0;
 912 }
 913
 914
 915 // ************************************************************************* //