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