src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.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 \*---------------------------------------------------------------------------*/
  26
  27 #include "cyclicPolyPatch.H"
  28 #include "addToRunTimeSelectionTable.H"
  29 #include "polyBoundaryMesh.H"
  30 #include "polyMesh.H"
  31 #include "demandDrivenData.H"
  32 #include "OFstream.H"
  33 #include "patchZones.H"
  34 #include "matchPoints.H"
  35 #include "EdgeMap.H"
  36 #include "Time.H"
  37 #include "transformList.H"
  38
  39 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  40
  41 namespace Foam
  42 {
  43     defineTypeNameAndDebug(cyclicPolyPatch, 0);
  44
  45     addToRunTimeSelectionTable(polyPatch, cyclicPolyPatch, word);
  46     addToRunTimeSelectionTable(polyPatch, cyclicPolyPatch, dictionary);
  47
  48
  49 template<>
  50 const char* NamedEnum<cyclicPolyPatch::transformType, 3>::names[] =
  51 {
  52     "unknown",
  53     "rotational",
  54     "translational"
  55 };
  56
  57 const NamedEnum<cyclicPolyPatch::transformType, 3>
  58     cyclicPolyPatch::transformTypeNames;
  59 }
  60
  61
  62 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  63
  64 Foam::label Foam::cyclicPolyPatch::findMaxArea
  65 (
  66     const pointField& points,
  67     const faceList& faces
  68 )
  69 {
  70     label maxI = -1;
  71     scalar maxAreaSqr = -GREAT;
  72
  73     forAll(faces, faceI)
  74     {
  75         scalar areaSqr = magSqr(faces[faceI].normal(points));
  76
  77         if (areaSqr > maxAreaSqr)
  78         {
  79             maxAreaSqr = areaSqr;
  80             maxI = faceI;
  81         }
  82     }
  83     return maxI;
  84 }
  85
  86
  87 void Foam::cyclicPolyPatch::calcTransforms()
  88 {
  89     if (size())
  90     {
  91         const pointField& points = this->points();
  92
  93         // Determine geometric quantities on the two halves
  94         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  95
  96         primitivePatch half0
  97         (
  98             SubList<face>
  99             (
 100                 *this,
 101                 size()/2
 102             ),
 103             points
 104         );
 105
 106         pointField half0Ctrs(calcFaceCentres(half0, half0.points()));
 107
 108         scalarField half0Tols(calcFaceTol(half0, half0.points(), half0Ctrs));
 109
 110         primitivePatch half1
 111         (
 112             SubList<face>
 113             (
 114                 *this,
 115                 size()/2,
 116                 size()/2
 117             ),
 118             points
 119         );
 120         pointField half1Ctrs(calcFaceCentres(half1, half1.points()));
 121
 122         // Dump halves
 123         if (debug)
 124         {
 125             fileName casePath(boundaryMesh().mesh().time().path());
 126
 127             fileName nm0(casePath/name()+"_half0_faces.obj");
 128             Pout<< "cyclicPolyPatch::calcTransforms : Writing half0"
 129                 << " faces to OBJ file " << nm0 << endl;
 130             writeOBJ(nm0, half0, half0.points());
 131
 132             fileName nm1(casePath/name()+"_half1_faces.obj");
 133             Pout<< "cyclicPolyPatch::calcTransforms : Writing half1"
 134                 << " faces to OBJ file " << nm1 << endl;
 135             writeOBJ(nm1, half1, half1.points());
 136
 137             OFstream str(casePath/name()+"_half0_to_half1.obj");
 138             label vertI = 0;
 139             Pout<< "cyclicPolyPatch::calcTransforms :"
 140                 << " Writing coupled face centres as lines to " << str.name()
 141                 << endl;
 142             forAll(half0Ctrs, i)
 143             {
 144                 const point& p0 = half0Ctrs[i];
 145                 str << "v " << p0.x() << ' ' << p0.y() << ' ' << p0.z() << nl;
 146                 vertI++;
 147                 const point& p1 = half1Ctrs[i];
 148                 str << "v " << p1.x() << ' ' << p1.y() << ' ' << p1.z() << nl;
 149                 vertI++;
 150                 str << "l " << vertI-1 << ' ' << vertI << nl;
 151             }
 152         }
 153
 154         vectorField half0Normals(half0.size());
 155         vectorField half1Normals(half1.size());
 156
 157         for (label facei = 0; facei < size()/2; facei++)
 158         {
 159             half0Normals[facei] = operator[](facei).normal(points);
 160             label nbrFacei = facei+size()/2;
 161             half1Normals[facei] = operator[](nbrFacei).normal(points);
 162
 163             scalar magSf = mag(half0Normals[facei]);
 164             scalar nbrMagSf = mag(half1Normals[facei]);
 165             scalar avSf = (magSf + nbrMagSf)/2.0;
 166
 167             if (magSf < ROOTVSMALL && nbrMagSf < ROOTVSMALL)
 168             {
 169                 // Undetermined normal. Use dummy normal to force separation
 170                 // check. (note use of sqrt(VSMALL) since that is how mag
 171                 // scales)
 172                 half0Normals[facei] = point(1, 0, 0);
 173                 half1Normals[facei] = half0Normals[facei];
 174             }
 175             else if (mag(magSf - nbrMagSf)/avSf > coupledPolyPatch::matchTol)
 176             {
 177                 FatalErrorIn
 178                 (
 179                     "cyclicPolyPatch::calcTransforms()"
 180                 )   << "face " << facei << " area does not match neighbour "
 181                     << nbrFacei << " by "
 182                     << 100*mag(magSf - nbrMagSf)/avSf
 183                     << "% -- possible face ordering problem." << endl
 184                     << "patch:" << name()
 185                     << " my area:" << magSf
 186                     << " neighbour area:" << nbrMagSf
 187                     << " matching tolerance:" << coupledPolyPatch::matchTol
 188                      << endl
 189                     << "Mesh face:" << start()+facei
 190                     << " vertices:"
 191                     << UIndirectList<point>(points, operator[](facei))()
 192                     << endl
 193                     << "Neighbour face:" << start()+nbrFacei
 194                     << " vertices:"
 195                     << UIndirectList<point>(points, operator[](nbrFacei))()
 196                     << endl
 197                     << "Rerun with cyclic debug flag set"
 198                     << " for more information." << exit(FatalError);
 199             }
 200             else
 201             {
 202                 half0Normals[facei] /= magSf;
 203                 half1Normals[facei] /= nbrMagSf;
 204             }
 205         }
 206
 207
 208         // See if transformation is prescribed
 209         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 210
 211         switch (transform_)
 212         {
 213             case ROTATIONAL:
 214             {
 215                 // Specified single rotation tensor.
 216
 217                 // Get best fitting face and its opposite number
 218                 label face0 = getConsistentRotationFace(half0Ctrs);
 219                 label face1 = face0;
 220
 221                 vector n0 =
 222                     (
 223                         (half0Ctrs[face0] - rotationCentre_)
 224                       ^ rotationAxis_
 225                     );
 226                 vector n1 =
 227                     (
 228                         (half1Ctrs[face1] - rotationCentre_)
 229                       ^ -rotationAxis_
 230                     );
 231                 n0 /= mag(n0) + VSMALL;
 232                 n1 /= mag(n1) + VSMALL;
 233
 234                 if (debug)
 235                 {
 236                     Pout<< "cyclicPolyPatch::calcTransforms :"
 237                         << " Specified rotation :"
 238                         << " n0:" << n0 << " n1:" << n1 << endl;
 239                 }
 240
 241                 // Calculate transformation tensors from face0,1 only.
 242                 // Note: can use tight tolerance now.
 243                 calcTransformTensors
 244                 (
 245                     pointField(1, half0Ctrs[face0]),
 246                     pointField(1, half1Ctrs[face1]),
 247                     vectorField(1, n0),
 248                     vectorField(1, n1),
 249                     scalarField(1, half0Tols[face0]),
 250                     1E-4
 251                 );
 252
 253                 break;
 254             }
 255
 256             default:
 257             {
 258                 // Calculate transformation tensors from all faces.
 259                 calcTransformTensors
 260                 (
 261                     half0Ctrs,
 262                     half1Ctrs,
 263                     half0Normals,
 264                     half1Normals,
 265                     half0Tols
 266                 );
 267
 268                 break;
 269             }
 270         }
 271     }
 272 }
 273
 274
 275 // Get geometric zones of patch by looking at normals.
 276 // Method 1: any edge with sharpish angle is edge between two halves.
 277 //           (this will handle e.g. wedge geometries).
 278 //           Also two fully disconnected regions will be handled this way.
 279 // Method 2: sort faces into two halves based on face normal.
 280 bool Foam::cyclicPolyPatch::getGeometricHalves
 281 (
 282     const primitivePatch& pp,
 283     labelList& half0ToPatch,
 284     labelList& half1ToPatch
 285 ) const
 286 {
 287     // Calculate normals
 288     const vectorField& faceNormals = pp.faceNormals();
 289
 290     // Find edges with sharp angles.
 291     boolList regionEdge(pp.nEdges(), false);
 292
 293     const labelListList& edgeFaces = pp.edgeFaces();
 294
 295     label nRegionEdges = 0;
 296
 297     forAll(edgeFaces, edgeI)
 298     {
 299         const labelList& eFaces = edgeFaces[edgeI];
 300
 301         // Check manifold edges for sharp angle.
 302         // (Non-manifold already handled by patchZones)
 303         if (eFaces.size() == 2)
 304         {
 305             if ((faceNormals[eFaces[0]] & faceNormals[eFaces[1]])< featureCos_)
 306             {
 307                 regionEdge[edgeI] = true;
 308
 309                 nRegionEdges++;
 310             }
 311         }
 312     }
 313
 314
 315     // For every face determine zone it is connected to (without crossing
 316     // any regionEdge)
 317     patchZones ppZones(pp, regionEdge);
 318
 319     if (debug)
 320     {
 321         Pout<< "cyclicPolyPatch::getGeometricHalves : "
 322             << "Found " << nRegionEdges << " edges on patch " << name()
 323             << " where the cos of the angle between two connected faces"
 324             << " was less than " << featureCos_ << nl
 325             << "Patch divided by these and by single sides edges into "
 326             << ppZones.nZones() << " parts." << endl;
 327
 328
 329         // Dumping zones to obj files.
 330
 331         labelList nZoneFaces(ppZones.nZones());
 332
 333         for (label zoneI = 0; zoneI < ppZones.nZones(); zoneI++)
 334         {
 335             OFstream stream
 336             (
 337                 boundaryMesh().mesh().time().path()
 338                /name()+"_zone_"+Foam::name(zoneI)+".obj"
 339             );
 340             Pout<< "cyclicPolyPatch::getGeometricHalves : Writing zone "
 341                 << zoneI << " face centres to OBJ file " << stream.name()
 342                 << endl;
 343
 344             labelList zoneFaces(findIndices(ppZones, zoneI));
 345
 346             forAll(zoneFaces, i)
 347             {
 348                 writeOBJ(stream, pp[zoneFaces[i]].centre(pp.points()));
 349             }
 350
 351             nZoneFaces[zoneI] = zoneFaces.size();
 352         }
 353     }
 354
 355
 356     if (ppZones.nZones() == 2)
 357     {
 358         half0ToPatch = findIndices(ppZones, 0);
 359         half1ToPatch = findIndices(ppZones, 1);
 360     }
 361     else
 362     {
 363         if (debug)
 364         {
 365             Pout<< "cyclicPolyPatch::getGeometricHalves :"
 366                 << " falling back to face-normal comparison" << endl;
 367         }
 368         label n0Faces = 0;
 369         half0ToPatch.setSize(pp.size());
 370
 371         label n1Faces = 0;
 372         half1ToPatch.setSize(pp.size());
 373
 374         // Compare to face 0 normal.
 375         forAll(faceNormals, faceI)
 376         {
 377             if ((faceNormals[faceI] & faceNormals[0]) > 0)
 378             {
 379                 half0ToPatch[n0Faces++] = faceI;
 380             }
 381             else
 382             {
 383                 half1ToPatch[n1Faces++] = faceI;
 384             }
 385         }
 386         half0ToPatch.setSize(n0Faces);
 387         half1ToPatch.setSize(n1Faces);
 388
 389         if (debug)
 390         {
 391             Pout<< "cyclicPolyPatch::getGeometricHalves :"
 392                 << " Number of faces per zone:("
 393                 << n0Faces << ' ' << n1Faces << ')' << endl;
 394         }
 395     }
 396
 397     if (half0ToPatch.size() != half1ToPatch.size())
 398     {
 399         fileName casePath(boundaryMesh().mesh().time().path());
 400
 401         // Dump halves
 402         {
 403             fileName nm0(casePath/name()+"_half0_faces.obj");
 404             Pout<< "cyclicPolyPatch::getGeometricHalves : Writing half0"
 405                 << " faces to OBJ file " << nm0 << endl;
 406             writeOBJ(nm0, UIndirectList<face>(pp, half0ToPatch)(), pp.points());
 407
 408             fileName nm1(casePath/name()+"_half1_faces.obj");
 409             Pout<< "cyclicPolyPatch::getGeometricHalves : Writing half1"
 410                 << " faces to OBJ file " << nm1 << endl;
 411             writeOBJ(nm1, UIndirectList<face>(pp, half1ToPatch)(), pp.points());
 412         }
 413
 414         // Dump face centres
 415         {
 416             OFstream str0(casePath/name()+"_half0.obj");
 417             Pout<< "cyclicPolyPatch::getGeometricHalves : Writing half0"
 418                 << " face centres to OBJ file " << str0.name() << endl;
 419
 420             forAll(half0ToPatch, i)
 421             {
 422                 writeOBJ(str0, pp[half0ToPatch[i]].centre(pp.points()));
 423             }
 424
 425             OFstream str1(casePath/name()+"_half1.obj");
 426             Pout<< "cyclicPolyPatch::getGeometricHalves : Writing half1"
 427                 << " face centres to OBJ file " << str1.name() << endl;
 428             forAll(half1ToPatch, i)
 429             {
 430                 writeOBJ(str1, pp[half1ToPatch[i]].centre(pp.points()));
 431             }
 432         }
 433
 434         SeriousErrorIn
 435         (
 436             "cyclicPolyPatch::getGeometricHalves"
 437             "(const primitivePatch&, labelList&, labelList&) const"
 438         )   << "Patch " << name() << " gets decomposed in two zones of"
 439             << "inequal size: " << half0ToPatch.size()
 440             << " and " << half1ToPatch.size() << endl
 441             << "This means that the patch is either not two separate regions"
 442             << " or one region where the angle between the different regions"
 443             << " is not sufficiently sharp." << endl
 444             << "Please adapt the featureCos setting." << endl
 445             << "Continuing with incorrect face ordering from now on!" << endl;
 446
 447         return false;
 448     }
 449     else
 450     {
 451         return true;
 452     }
 453 }
 454
 455
 456 // Given a split of faces into left and right half calculate the centres
 457 // and anchor points. Transform the left points so they align with the
 458 // right ones.
 459 void Foam::cyclicPolyPatch::getCentresAndAnchors
 460 (
 461     const primitivePatch& pp,
 462     const faceList& half0Faces,
 463     const faceList& half1Faces,
 464
 465     pointField& ppPoints,
 466     pointField& half0Ctrs,
 467     pointField& half1Ctrs,
 468     pointField& anchors0,
 469     scalarField& tols
 470 ) const
 471 {
 472     // Get geometric data on both halves.
 473     half0Ctrs = calcFaceCentres(half0Faces, pp.points());
 474     anchors0 = getAnchorPoints(half0Faces, pp.points());
 475     half1Ctrs = calcFaceCentres(half1Faces, pp.points());
 476
 477     switch (transform_)
 478     {
 479         case ROTATIONAL:
 480         {
 481             label face0 = getConsistentRotationFace(half0Ctrs);
 482             label face1 = getConsistentRotationFace(half1Ctrs);
 483
 484             vector n0 = ((half0Ctrs[face0] - rotationCentre_) ^ rotationAxis_);
 485             vector n1 = ((half1Ctrs[face1] - rotationCentre_) ^ -rotationAxis_);
 486             n0 /= mag(n0) + VSMALL;
 487             n1 /= mag(n1) + VSMALL;
 488
 489             if (debug)
 490             {
 491                 Pout<< "cyclicPolyPatch::getCentresAndAnchors :"
 492                     << " Specified rotation :"
 493                     << " n0:" << n0 << " n1:" << n1 << endl;
 494             }
 495
 496             // Rotation (around origin)
 497             const tensor reverseT(rotationTensor(n0, -n1));
 498
 499             // Rotation
 500             forAll(half0Ctrs, faceI)
 501             {
 502                 half0Ctrs[faceI] = Foam::transform(reverseT, half0Ctrs[faceI]);
 503                 anchors0[faceI] = Foam::transform(reverseT, anchors0[faceI]);
 504             }
 505
 506             ppPoints = Foam::transform(reverseT, pp.points());
 507
 508             break;
 509         }
 510         //- Problem: usually specified translation is not accurate enough
 511         //- to get proper match so keep automatic determination over here.
 512         //case TRANSLATIONAL:
 513         //{
 514         //    // Transform 0 points.
 515         //
 516         //    if (debug)
 517         //    {
 518         //        Pout<< "cyclicPolyPatch::getCentresAndAnchors :"
 519         //            << "Specified translation : " << separationVector_ << endl;
 520         //    }
 521         //
 522         //    half0Ctrs += separationVector_;
 523         //    anchors0 += separationVector_;
 524         //    break;
 525         //}
 526         default:
 527         {
 528             // Assumes that cyclic is planar. This is also the initial
 529             // condition for patches without faces.
 530
 531             // Determine the face with max area on both halves. These
 532             // two faces are used to determine the transformation tensors
 533             label max0I = findMaxArea(pp.points(), half0Faces);
 534             vector n0 = half0Faces[max0I].normal(pp.points());
 535             n0 /= mag(n0) + VSMALL;
 536
 537             label max1I = findMaxArea(pp.points(), half1Faces);
 538             vector n1 = half1Faces[max1I].normal(pp.points());
 539             n1 /= mag(n1) + VSMALL;
 540
 541             if (mag(n0 & n1) < 1-coupledPolyPatch::matchTol)
 542             {
 543                 if (debug)
 544                 {
 545                     Pout<< "cyclicPolyPatch::getCentresAndAnchors :"
 546                         << " Detected rotation :"
 547                         << " n0:" << n0 << " n1:" << n1 << endl;
 548                 }
 549
 550                 // Rotation (around origin)
 551                 const tensor reverseT(rotationTensor(n0, -n1));
 552
 553                 // Rotation
 554                 forAll(half0Ctrs, faceI)
 555                 {
 556                     half0Ctrs[faceI] = Foam::transform
 557                     (
 558                         reverseT,
 559                         half0Ctrs[faceI]
 560                     );
 561                     anchors0[faceI] = Foam::transform
 562                     (
 563                         reverseT,
 564                         anchors0[faceI]
 565                     );
 566                 }
 567                 ppPoints = Foam::transform(reverseT, pp.points());
 568             }
 569             else
 570             {
 571                 // Parallel translation. Get average of all used points.
 572
 573                 primitiveFacePatch half0(half0Faces, pp.points());
 574                 const pointField& half0Pts = half0.localPoints();
 575                 const point ctr0(sum(half0Pts)/half0Pts.size());
 576
 577                 primitiveFacePatch half1(half1Faces, pp.points());
 578                 const pointField& half1Pts = half1.localPoints();
 579                 const point ctr1(sum(half1Pts)/half1Pts.size());
 580
 581                 if (debug)
 582                 {
 583                     Pout<< "cyclicPolyPatch::getCentresAndAnchors :"
 584                         << " Detected translation :"
 585                         << " n0:" << n0 << " n1:" << n1
 586                         << " ctr0:" << ctr0 << " ctr1:" << ctr1 << endl;
 587                 }
 588
 589                 half0Ctrs += ctr1 - ctr0;
 590                 anchors0 += ctr1 - ctr0;
 591                 ppPoints = pp.points() + ctr1 - ctr0;
 592             }
 593             break;
 594         }
 595     }
 596
 597
 598     // Calculate typical distance per face
 599     tols = calcFaceTol(half1Faces, pp.points(), half1Ctrs);
 600 }
 601
 602
 603 // Calculates faceMap and rotation. Returns true if all ok.
 604 bool Foam::cyclicPolyPatch::matchAnchors
 605 (
 606     const bool report,
 607     const primitivePatch& pp,
 608     const labelList& half0ToPatch,
 609     const pointField& anchors0,
 610
 611     const labelList& half1ToPatch,
 612     const faceList& half1Faces,
 613     const labelList& from1To0,
 614
 615     const scalarField& tols,
 616
 617     labelList& faceMap,
 618     labelList& rotation
 619 ) const
 620 {
 621     // Set faceMap such that half0 faces get first and corresponding half1
 622     // faces last.
 623
 624     forAll(half0ToPatch, half0FaceI)
 625     {
 626         // Label in original patch
 627         label patchFaceI = half0ToPatch[half0FaceI];
 628
 629         faceMap[patchFaceI] = half0FaceI;
 630
 631         // No rotation
 632         rotation[patchFaceI] = 0;
 633     }
 634
 635     bool fullMatch = true;
 636
 637     forAll(from1To0, half1FaceI)
 638     {
 639         label patchFaceI = half1ToPatch[half1FaceI];
 640
 641         // This face has to match the corresponding one on half0.
 642         label half0FaceI = from1To0[half1FaceI];
 643
 644         label newFaceI = half0FaceI + pp.size()/2;
 645
 646         faceMap[patchFaceI] = newFaceI;
 647
 648         // Rotate patchFaceI such that its f[0] aligns with that of
 649         // the corresponding face
 650         // (which after shuffling will be at position half0FaceI)
 651
 652         const point& wantedAnchor = anchors0[half0FaceI];
 653
 654         rotation[newFaceI] = getRotation
 655         (
 656             pp.points(),
 657             half1Faces[half1FaceI],
 658             wantedAnchor,
 659             tols[half1FaceI]
 660         );
 661
 662         if (rotation[newFaceI] == -1)
 663         {
 664             fullMatch = false;
 665
 666             if (report)
 667             {
 668                 const face& f = half1Faces[half1FaceI];
 669                 SeriousErrorIn
 670                 (
 671                     "cyclicPolyPatch::matchAnchors(..)"
 672                 )   << "Patch:" << name() << " : "
 673                     << "Cannot find point on face " << f
 674                     << " with vertices:"
 675                     << UIndirectList<point>(pp.points(), f)()
 676                     << " that matches point " << wantedAnchor
 677                     << " when matching the halves of cyclic patch " << name()
 678                     << endl
 679                     << "Continuing with incorrect face ordering from now on!"
 680                     << endl;
 681             }
 682         }
 683     }
 684     return fullMatch;
 685 }
 686
 687
 688 Foam::label Foam::cyclicPolyPatch::getConsistentRotationFace
 689 (
 690     const pointField& faceCentres
 691 ) const
 692 {
 693     const scalarField magRadSqr =
 694         magSqr((faceCentres - rotationCentre_) ^ rotationAxis_);
 695     scalarField axisLen = (faceCentres - rotationCentre_) & rotationAxis_;
 696     axisLen = axisLen - min(axisLen);
 697     const scalarField magLenSqr = magRadSqr + axisLen*axisLen;
 698
 699     label rotFace = -1;
 700     scalar maxMagLenSqr = -GREAT;
 701     scalar maxMagRadSqr = -GREAT;
 702     forAll(faceCentres, i)
 703     {
 704         if (magLenSqr[i] >= maxMagLenSqr)
 705         {
 706             if (magRadSqr[i] > maxMagRadSqr)
 707             {
 708                 rotFace = i;
 709                 maxMagLenSqr = magLenSqr[i];
 710                 maxMagRadSqr = magRadSqr[i];
 711             }
 712         }
 713     }
 714
 715     if (debug)
 716     {
 717         Info<< "getConsistentRotationFace(const pointField&)" << nl
 718             << "    rotFace = " << rotFace << nl
 719             << "    point =  " << faceCentres[rotFace] << endl;
 720     }
 721
 722     return rotFace;
 723 }
 724
 725
 726 // * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * * * * * //
 727
 728 Foam::cyclicPolyPatch::cyclicPolyPatch
 729 (
 730     const word& name,
 731     const label size,
 732     const label start,
 733     const label index,
 734     const polyBoundaryMesh& bm
 735 )
 736 :
 737     coupledPolyPatch(name, size, start, index, bm),
 738     coupledPointsPtr_(NULL),
 739     coupledEdgesPtr_(NULL),
 740     featureCos_(0.9),
 741     transform_(UNKNOWN),
 742     rotationAxis_(vector::zero),
 743     rotationCentre_(point::zero),
 744     separationVector_(vector::zero)
 745 {}
 746
 747
 748 Foam::cyclicPolyPatch::cyclicPolyPatch
 749 (
 750     const word& name,
 751     const dictionary& dict,
 752     const label index,
 753     const polyBoundaryMesh& bm
 754 )
 755 :
 756     coupledPolyPatch(name, dict, index, bm),
 757     coupledPointsPtr_(NULL),
 758     coupledEdgesPtr_(NULL),
 759     featureCos_(0.9),
 760     transform_(UNKNOWN),
 761     rotationAxis_(vector::zero),
 762     rotationCentre_(point::zero),
 763     separationVector_(vector::zero)
 764 {
 765     dict.readIfPresent("featureCos", featureCos_);
 766
 767     if (dict.found("transform"))
 768     {
 769         transform_ = transformTypeNames.read(dict.lookup("transform"));
 770         switch (transform_)
 771         {
 772             case ROTATIONAL:
 773             {
 774                 dict.lookup("rotationAxis") >> rotationAxis_;
 775                 dict.lookup("rotationCentre") >> rotationCentre_;
 776                 break;
 777             }
 778             case TRANSLATIONAL:
 779             {
 780                 dict.lookup("separationVector") >> separationVector_;
 781                 break;
 782             }
 783             default:
 784             {
 785                 // no additional info required
 786             }
 787         }
 788     }
 789 }
 790
 791
 792 Foam::cyclicPolyPatch::cyclicPolyPatch
 793 (
 794     const cyclicPolyPatch& pp,
 795     const polyBoundaryMesh& bm
 796 )
 797 :
 798     coupledPolyPatch(pp, bm),
 799     coupledPointsPtr_(NULL),
 800     coupledEdgesPtr_(NULL),
 801     featureCos_(pp.featureCos_),
 802     transform_(pp.transform_),
 803     rotationAxis_(pp.rotationAxis_),
 804     rotationCentre_(pp.rotationCentre_),
 805     separationVector_(pp.separationVector_)
 806 {}
 807
 808
 809 Foam::cyclicPolyPatch::cyclicPolyPatch
 810 (
 811     const cyclicPolyPatch& pp,
 812     const polyBoundaryMesh& bm,
 813     const label index,
 814     const label newSize,
 815     const label newStart
 816 )
 817 :
 818     coupledPolyPatch(pp, bm, index, newSize, newStart),
 819     coupledPointsPtr_(NULL),
 820     coupledEdgesPtr_(NULL),
 821     featureCos_(pp.featureCos_),
 822     transform_(pp.transform_),
 823     rotationAxis_(pp.rotationAxis_),
 824     rotationCentre_(pp.rotationCentre_),
 825     separationVector_(pp.separationVector_)
 826 {}
 827
 828
 829 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 830
 831 Foam::cyclicPolyPatch::~cyclicPolyPatch()
 832 {
 833     deleteDemandDrivenData(coupledPointsPtr_);
 834     deleteDemandDrivenData(coupledEdgesPtr_);
 835 }
 836
 837
 838
 839 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 840
 841 void Foam::cyclicPolyPatch::initGeometry()
 842 {
 843     polyPatch::initGeometry();
 844 }
 845
 846 void Foam::cyclicPolyPatch::calcGeometry()
 847 {
 848     polyPatch::calcGeometry();
 849     calcTransforms();
 850 }
 851
 852 void Foam::cyclicPolyPatch::initMovePoints(const pointField& p)
 853 {
 854     polyPatch::initMovePoints(p);
 855 }
 856
 857 void Foam::cyclicPolyPatch::movePoints(const pointField& p)
 858 {
 859     polyPatch::movePoints(p);
 860     calcTransforms();
 861 }
 862
 863 void Foam::cyclicPolyPatch::initUpdateMesh()
 864 {
 865     polyPatch::initUpdateMesh();
 866 }
 867
 868 void Foam::cyclicPolyPatch::updateMesh()
 869 {
 870     polyPatch::updateMesh();
 871     deleteDemandDrivenData(coupledPointsPtr_);
 872     deleteDemandDrivenData(coupledEdgesPtr_);
 873 }
 874
 875
 876 const Foam::edgeList& Foam::cyclicPolyPatch::coupledPoints() const
 877 {
 878     if (!coupledPointsPtr_)
 879     {
 880         // Look at cyclic patch as two halves, A and B.
 881         // Now all we know is that relative face index in halfA is same
 882         // as coupled face in halfB and also that the 0th vertex
 883         // corresponds.
 884
 885         // From halfA point to halfB or -1.
 886         labelList coupledPoint(nPoints(), -1);
 887
 888         for (label patchFaceA = 0; patchFaceA < size()/2; patchFaceA++)
 889         {
 890             const face& fA = localFaces()[patchFaceA];
 891
 892             forAll(fA, indexA)
 893             {
 894                 label patchPointA = fA[indexA];
 895
 896                 if (coupledPoint[patchPointA] == -1)
 897                 {
 898                     const face& fB = localFaces()[patchFaceA + size()/2];
 899
 900                     label indexB = (fB.size() - indexA) % fB.size();
 901
 902                     // Filter out points on wedge axis
 903                     if (patchPointA != fB[indexB])
 904                     {
 905                         coupledPoint[patchPointA] = fB[indexB];
 906                     }
 907                 }
 908             }
 909         }
 910
 911         coupledPointsPtr_ = new edgeList(nPoints());
 912         edgeList& connected = *coupledPointsPtr_;
 913
 914         // Extract coupled points.
 915         label connectedI = 0;
 916
 917         forAll(coupledPoint, i)
 918         {
 919             if (coupledPoint[i] != -1)
 920             {
 921                 connected[connectedI++] = edge(i, coupledPoint[i]);
 922             }
 923         }
 924
 925         connected.setSize(connectedI);
 926
 927         if (debug)
 928         {
 929             OFstream str
 930             (
 931                 boundaryMesh().mesh().time().path()
 932                /"coupledPoints.obj"
 933             );
 934             label vertI = 0;
 935
 936             Pout<< "Writing file " << str.name() << " with coordinates of "
 937                 << "coupled points" << endl;
 938
 939             forAll(connected, i)
 940             {
 941                 const point& a = points()[meshPoints()[connected[i][0]]];
 942                 const point& b = points()[meshPoints()[connected[i][1]]];
 943
 944                 str<< "v " << a.x() << ' ' << a.y() << ' ' << a.z() << nl;
 945                 str<< "v " << b.x() << ' ' << b.y() << ' ' << b.z() << nl;
 946                 vertI += 2;
 947
 948                 str<< "l " << vertI-1 << ' ' << vertI << nl;
 949             }
 950         }
 951
 952         // Remove any addressing calculated for the coupled edges calculation
 953         const_cast<primitivePatch&>(static_cast<const primitivePatch&>(*this))
 954             .clearOut();
 955     }
 956     return *coupledPointsPtr_;
 957 }
 958
 959
 960 const Foam::edgeList& Foam::cyclicPolyPatch::coupledEdges() const
 961 {
 962     if (!coupledEdgesPtr_)
 963     {
 964         // Build map from points on halfA to points on halfB.
 965         const edgeList& pointCouples = coupledPoints();
 966
 967         Map<label> aToB(2*pointCouples.size());
 968
 969         forAll(pointCouples, i)
 970         {
 971             const edge& e = pointCouples[i];
 972
 973             aToB.insert(e[0], e[1]);
 974         }
 975
 976         // Map from edge on half A to points (in halfB indices)
 977         EdgeMap<label> edgeMap(nEdges());
 978
 979         for (label patchFaceA = 0; patchFaceA < size()/2; patchFaceA++)
 980         {
 981             const labelList& fEdges = faceEdges()[patchFaceA];
 982
 983             forAll(fEdges, i)
 984             {
 985                 label edgeI = fEdges[i];
 986
 987                 const edge& e = edges()[edgeI];
 988
 989                 // Convert edge end points to corresponding points on halfB
 990                 // side.
 991                 Map<label>::const_iterator fnd0 = aToB.find(e[0]);
 992                 if (fnd0 != aToB.end())
 993                 {
 994                     Map<label>::const_iterator fnd1 = aToB.find(e[1]);
 995                     if (fnd1 != aToB.end())
 996                     {
 997                         edgeMap.insert(edge(fnd0(), fnd1()), edgeI);
 998                     }
 999                 }
1000             }
1001         }
1002
1003         coupledEdgesPtr_ = new edgeList(nEdges()/2);
1004         edgeList& coupledEdges = *coupledEdgesPtr_;
1005         label coupleI = 0;
1006
1007         for (label patchFaceB = size()/2; patchFaceB < size(); patchFaceB++)
1008         {
1009             const labelList& fEdges = faceEdges()[patchFaceB];
1010
1011             forAll(fEdges, i)
1012             {
1013                 label edgeI = fEdges[i];
1014
1015                 const edge& e = edges()[edgeI];
1016
1017                 // Look up halfA edge from HashTable.
1018                 EdgeMap<label>::iterator iter = edgeMap.find(e);
1019
1020                 if (iter != edgeMap.end())
1021                 {
1022                     label halfAEdgeI = iter();
1023
1024                     // Store correspondence. Filter out edges on wedge axis.
1025                     if (halfAEdgeI != edgeI)
1026                     {
1027                         coupledEdges[coupleI++] = edge(halfAEdgeI, edgeI);
1028                     }
1029
1030                     // Remove so we build unique list
1031                     edgeMap.erase(iter);
1032                 }
1033             }
1034         }
1035         coupledEdges.setSize(coupleI);
1036
1037
1038         // Some checks
1039
1040         forAll(coupledEdges, i)
1041         {
1042             const edge& e = coupledEdges[i];
1043
1044             if (e[0] == e[1] || e[0] < 0 || e[1] < 0)
1045             {
1046                 FatalErrorIn("cyclicPolyPatch::coupledEdges() const")
1047                     << "Problem : at position " << i
1048                     << " illegal couple:" << e
1049                     << abort(FatalError);
1050             }
1051         }
1052
1053         if (debug)
1054         {
1055             OFstream str
1056             (
1057                 boundaryMesh().mesh().time().path()
1058                /"coupledEdges.obj"
1059             );
1060             label vertI = 0;
1061
1062             Pout<< "Writing file " << str.name() << " with centres of "
1063                 << "coupled edges" << endl;
1064
1065             forAll(coupledEdges, i)
1066             {
1067                 const edge& e = coupledEdges[i];
1068
1069                 const point& a = edges()[e[0]].centre(localPoints());
1070                 const point& b = edges()[e[1]].centre(localPoints());
1071
1072                 str<< "v " << a.x() << ' ' << a.y() << ' ' << a.z() << nl;
1073                 str<< "v " << b.x() << ' ' << b.y() << ' ' << b.z() << nl;
1074                 vertI += 2;
1075
1076                 str<< "l " << vertI-1 << ' ' << vertI << nl;
1077             }
1078         }
1079
1080         // Remove any addressing calculated for the coupled edges calculation
1081         const_cast<primitivePatch&>(static_cast<const primitivePatch&>(*this))
1082             .clearOut();
1083     }
1084     return *coupledEdgesPtr_;
1085 }
1086
1087
1088 void Foam::cyclicPolyPatch::initOrder(const primitivePatch& pp) const
1089 {}
1090
1091
1092 //  Return new ordering. Ordering is -faceMap: for every face index
1093 //  the new face -rotation:for every new face the clockwise shift
1094 //  of the original face. Return false if nothing changes (faceMap
1095 //  is identity, rotation is 0)
1096 bool Foam::cyclicPolyPatch::order
1097 (
1098     const primitivePatch& pp,
1099     labelList& faceMap,
1100     labelList& rotation
1101 ) const
1102 {
1103     faceMap.setSize(pp.size());
1104     faceMap = -1;
1105
1106     rotation.setSize(pp.size());
1107     rotation = 0;
1108
1109     if (pp.empty())
1110     {
1111         // No faces, nothing to change.
1112         return false;
1113     }
1114
1115     if (pp.size()&1)
1116     {
1117         FatalErrorIn("cyclicPolyPatch::order(..)")
1118             << "Size of cyclic " << name() << " should be a multiple of 2"
1119             << ". It is " << pp.size() << abort(FatalError);
1120     }
1121
1122     label halfSize = pp.size()/2;
1123
1124     // Supplied primitivePatch already with new points.
1125     // Cyclics are limited to one transformation tensor
1126     // currently anyway (i.e. straight plane) so should not be too big a
1127     // problem.
1128
1129
1130     // Indices of faces on half0
1131     labelList half0ToPatch;
1132     // Indices of faces on half1
1133     labelList half1ToPatch;
1134
1135
1136     // 1. Test if already correctly oriented by starting from trivial ordering.
1137     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1138
1139     half0ToPatch = identity(halfSize);
1140     half1ToPatch = half0ToPatch + halfSize;
1141
1142     // Get faces
1143     faceList half0Faces(UIndirectList<face>(pp, half0ToPatch));
1144     faceList half1Faces(UIndirectList<face>(pp, half1ToPatch));
1145
1146     // Get geometric quantities
1147     pointField half0Ctrs, half1Ctrs, anchors0, ppPoints;
1148     scalarField tols;
1149     getCentresAndAnchors
1150     (
1151         pp,
1152         half0Faces,
1153         half1Faces,
1154
1155         ppPoints,
1156         half0Ctrs,
1157         half1Ctrs,
1158         anchors0,
1159         tols
1160     );
1161
1162     // Geometric match of face centre vectors
1163     labelList from1To0;
1164     bool matchedAll = matchPoints
1165     (
1166         half1Ctrs,
1167         half0Ctrs,
1168         tols,
1169         false,
1170         from1To0
1171     );
1172
1173     if (debug)
1174     {
1175         Pout<< "cyclicPolyPatch::order : test if already ordered:"
1176             << matchedAll << endl;
1177
1178         // Dump halves
1179         fileName nm0("match1_"+name()+"_half0_faces.obj");
1180         Pout<< "cyclicPolyPatch::order : Writing half0"
1181             << " faces to OBJ file " << nm0 << endl;
1182         writeOBJ(nm0, half0Faces, ppPoints);
1183
1184         fileName nm1("match1_"+name()+"_half1_faces.obj");
1185         Pout<< "cyclicPolyPatch::order : Writing half1"
1186             << " faces to OBJ file " << nm1 << endl;
1187         writeOBJ(nm1, half1Faces, pp.points());
1188
1189         OFstream ccStr
1190         (
1191             boundaryMesh().mesh().time().path()
1192            /"match1_"+ name() + "_faceCentres.obj"
1193         );
1194         Pout<< "cyclicPolyPatch::order : "
1195             << "Dumping currently found cyclic match as lines between"
1196             << " corresponding face centres to file " << ccStr.name()
1197             << endl;
1198
1199         // Recalculate untransformed face centres
1200         //pointField rawHalf0Ctrs = calcFaceCentres(half0Faces, pp.points());
1201         label vertI = 0;
1202
1203         forAll(half1Ctrs, i)
1204         {
1205             //if (from1To0[i] != -1)
1206             {
1207                 // Write edge between c1 and c0
1208                 //const point& c0 = rawHalf0Ctrs[from1To0[i]];
1209                 //const point& c0 = half0Ctrs[from1To0[i]];
1210                 const point& c0 = half0Ctrs[i];
1211                 const point& c1 = half1Ctrs[i];
1212                 writeOBJ(ccStr, c0, c1, vertI);
1213             }
1214         }
1215     }
1216
1217
1218     // 2. Ordered in pairs (so 0,1 coupled and 2,3 etc.)
1219     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1220
1221     if (!matchedAll)
1222     {
1223         label faceI = 0;
1224         for (label i = 0; i < halfSize; i++)
1225         {
1226             half0ToPatch[i] = faceI++;
1227             half1ToPatch[i] = faceI++;
1228         }
1229
1230         // And redo all matching
1231         half0Faces = UIndirectList<face>(pp, half0ToPatch);
1232         half1Faces = UIndirectList<face>(pp, half1ToPatch);
1233
1234         getCentresAndAnchors
1235         (
1236             pp,
1237             half0Faces,
1238             half1Faces,
1239
1240             ppPoints,
1241             half0Ctrs,
1242             half1Ctrs,
1243             anchors0,
1244             tols
1245         );
1246
1247         // Geometric match of face centre vectors
1248         matchedAll = matchPoints
1249         (
1250             half1Ctrs,
1251             half0Ctrs,
1252             tols,
1253             false,
1254             from1To0
1255         );
1256
1257         if (debug)
1258         {
1259             Pout<< "cyclicPolyPatch::order : test if pairwise ordered:"
1260                 << matchedAll << endl;
1261             // Dump halves
1262             fileName nm0("match2_"+name()+"_half0_faces.obj");
1263             Pout<< "cyclicPolyPatch::order : Writing half0"
1264                 << " faces to OBJ file " << nm0 << endl;
1265             writeOBJ(nm0, half0Faces, ppPoints);
1266
1267             fileName nm1("match2_"+name()+"_half1_faces.obj");
1268             Pout<< "cyclicPolyPatch::order : Writing half1"
1269                 << " faces to OBJ file " << nm1 << endl;
1270             writeOBJ(nm1, half1Faces, pp.points());
1271
1272             OFstream ccStr
1273             (
1274                 boundaryMesh().mesh().time().path()
1275                /"match2_"+name()+"_faceCentres.obj"
1276             );
1277             Pout<< "cyclicPolyPatch::order : "
1278                 << "Dumping currently found cyclic match as lines between"
1279                 << " corresponding face centres to file " << ccStr.name()
1280                 << endl;
1281
1282             // Recalculate untransformed face centres
1283             //pointField rawHalf0Ctrs = calcFaceCentres(half0Faces, pp.points());
1284             label vertI = 0;
1285
1286             forAll(half1Ctrs, i)
1287             {
1288                 if (from1To0[i] != -1)
1289                 {
1290                     // Write edge between c1 and c0
1291                     //const point& c0 = rawHalf0Ctrs[from1To0[i]];
1292                     const point& c0 = half0Ctrs[from1To0[i]];
1293                     const point& c1 = half1Ctrs[i];
1294                     writeOBJ(ccStr, c0, c1, vertI);
1295                 }
1296             }
1297         }
1298     }
1299
1300
1301     // 3. Baffles(coincident faces) converted into cyclics (e.g. jump)
1302     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1303
1304     if (!matchedAll)
1305     {
1306         label baffleI = 0;
1307
1308         forAll(*this, faceI)
1309         {
1310             const face& f = pp.localFaces()[faceI];
1311             const labelList& pFaces = pp.pointFaces()[f[0]];
1312
1313             label matchedFaceI = -1;
1314
1315             forAll(pFaces, i)
1316             {
1317                 label otherFaceI = pFaces[i];
1318
1319                 if (otherFaceI > faceI)
1320                 {
1321                     const face& otherF = pp.localFaces()[otherFaceI];
1322
1323                     // Note: might pick up two similar oriented faces
1324                     //       (but that is illegal anyway)
1325                     if (f == otherF)
1326                     {
1327                         matchedFaceI = otherFaceI;
1328                         break;
1329                     }
1330                 }
1331             }
1332
1333             if (matchedFaceI != -1)
1334             {
1335                 half0ToPatch[baffleI] = faceI;
1336                 half1ToPatch[baffleI] = matchedFaceI;
1337                 baffleI++;
1338             }
1339         }
1340
1341         if (baffleI == halfSize)
1342         {
1343             // And redo all matching
1344             half0Faces = UIndirectList<face>(pp, half0ToPatch);
1345             half1Faces = UIndirectList<face>(pp, half1ToPatch);
1346
1347             getCentresAndAnchors
1348             (
1349                 pp,
1350                 half0Faces,
1351                 half1Faces,
1352
1353                 ppPoints,
1354                 half0Ctrs,
1355                 half1Ctrs,
1356                 anchors0,
1357                 tols
1358             );
1359
1360             // Geometric match of face centre vectors
1361             matchedAll = matchPoints
1362             (
1363                 half1Ctrs,
1364                 half0Ctrs,
1365                 tols,
1366                 false,
1367                 from1To0
1368             );
1369
1370             if (debug)
1371             {
1372                 Pout<< "cyclicPolyPatch::order : test if baffles:"
1373                     << matchedAll << endl;
1374                 // Dump halves
1375                 fileName nm0("match3_"+name()+"_half0_faces.obj");
1376                 Pout<< "cyclicPolyPatch::order : Writing half0"
1377                     << " faces to OBJ file " << nm0 << endl;
1378                 writeOBJ(nm0, half0Faces, ppPoints);
1379
1380                 fileName nm1("match3_"+name()+"_half1_faces.obj");
1381                 Pout<< "cyclicPolyPatch::order : Writing half1"
1382                     << " faces to OBJ file " << nm1 << endl;
1383                 writeOBJ(nm1, half1Faces, pp.points());
1384
1385                 OFstream ccStr
1386                 (
1387                     boundaryMesh().mesh().time().path()
1388                    /"match3_"+ name() + "_faceCentres.obj"
1389                 );
1390                 Pout<< "cyclicPolyPatch::order : "
1391                     << "Dumping currently found cyclic match as lines between"
1392                     << " corresponding face centres to file " << ccStr.name()
1393                     << endl;
1394
1395                 // Recalculate untransformed face centres
1396                 //pointField rawHalf0Ctrs = calcFaceCentres(half0Faces, pp.points());
1397                 label vertI = 0;
1398
1399                 forAll(half1Ctrs, i)
1400                 {
1401                     if (from1To0[i] != -1)
1402                     {
1403                         // Write edge between c1 and c0
1404                         //const point& c0 = rawHalf0Ctrs[from1To0[i]];
1405                         const point& c0 = half0Ctrs[from1To0[i]];
1406                         const point& c1 = half1Ctrs[i];
1407                         writeOBJ(ccStr, c0, c1, vertI);
1408                     }
1409                 }
1410             }
1411         }
1412     }
1413
1414
1415     // 4. Automatic geometric ordering
1416     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1417
1418     if (!matchedAll)
1419     {
1420         // Split faces according to feature angle or topology
1421         bool okSplit = getGeometricHalves(pp, half0ToPatch, half1ToPatch);
1422
1423         if (!okSplit)
1424         {
1425             // Did not split into two equal parts.
1426             return false;
1427         }
1428
1429         // And redo all matching
1430         half0Faces = UIndirectList<face>(pp, half0ToPatch);
1431         half1Faces = UIndirectList<face>(pp, half1ToPatch);
1432
1433         getCentresAndAnchors
1434         (
1435             pp,
1436             half0Faces,
1437             half1Faces,
1438
1439             ppPoints,
1440             half0Ctrs,
1441             half1Ctrs,
1442             anchors0,
1443             tols
1444         );
1445
1446         // Geometric match of face centre vectors
1447         matchedAll = matchPoints
1448         (
1449             half1Ctrs,
1450             half0Ctrs,
1451             tols,
1452             false,
1453             from1To0
1454         );
1455
1456         if (debug)
1457         {
1458             Pout<< "cyclicPolyPatch::order : automatic ordering result:"
1459                 << matchedAll << endl;
1460             // Dump halves
1461             fileName nm0("match4_"+name()+"_half0_faces.obj");
1462             Pout<< "cyclicPolyPatch::order : Writing half0"
1463                 << " faces to OBJ file " << nm0 << endl;
1464             writeOBJ(nm0, half0Faces, ppPoints);
1465
1466             fileName nm1("match4_"+name()+"_half1_faces.obj");
1467             Pout<< "cyclicPolyPatch::order : Writing half1"
1468                 << " faces to OBJ file " << nm1 << endl;
1469             writeOBJ(nm1, half1Faces, pp.points());
1470
1471             OFstream ccStr
1472             (
1473                 boundaryMesh().mesh().time().path()
1474                /"match4_"+ name() + "_faceCentres.obj"
1475             );
1476             Pout<< "cyclicPolyPatch::order : "
1477                 << "Dumping currently found cyclic match as lines between"
1478                 << " corresponding face centres to file " << ccStr.name()
1479                 << endl;
1480
1481             // Recalculate untransformed face centres
1482             //pointField rawHalf0Ctrs = calcFaceCentres(half0Faces, pp.points());
1483             label vertI = 0;
1484
1485             forAll(half1Ctrs, i)
1486             {
1487                 if (from1To0[i] != -1)
1488                 {
1489                     // Write edge between c1 and c0
1490                     //const point& c0 = rawHalf0Ctrs[from1To0[i]];
1491                     const point& c0 = half0Ctrs[from1To0[i]];
1492                     const point& c1 = half1Ctrs[i];
1493                     writeOBJ(ccStr, c0, c1, vertI);
1494                 }
1495             }
1496         }
1497     }
1498
1499
1500     if (!matchedAll || debug)
1501     {
1502         // Dump halves
1503         fileName nm0(name()+"_half0_faces.obj");
1504         Pout<< "cyclicPolyPatch::order : Writing half0"
1505             << " faces to OBJ file " << nm0 << endl;
1506         writeOBJ(nm0, half0Faces, pp.points());
1507
1508         fileName nm1(name()+"_half1_faces.obj");
1509         Pout<< "cyclicPolyPatch::order : Writing half1"
1510             << " faces to OBJ file " << nm1 << endl;
1511         writeOBJ(nm1, half1Faces, pp.points());
1512
1513         OFstream ccStr
1514         (
1515             boundaryMesh().mesh().time().path()
1516            /name() + "_faceCentres.obj"
1517         );
1518         Pout<< "cyclicPolyPatch::order : "
1519             << "Dumping currently found cyclic match as lines between"
1520             << " corresponding face centres to file " << ccStr.name()
1521             << endl;
1522
1523         // Recalculate untransformed face centres
1524         //pointField rawHalf0Ctrs = calcFaceCentres(half0Faces, pp.points());
1525         label vertI = 0;
1526
1527         forAll(half1Ctrs, i)
1528         {
1529             if (from1To0[i] != -1)
1530             {
1531                 // Write edge between c1 and c0
1532                 //const point& c0 = rawHalf0Ctrs[from1To0[i]];
1533                 const point& c0 = half0Ctrs[from1To0[i]];
1534                 const point& c1 = half1Ctrs[i];
1535                 writeOBJ(ccStr, c0, c1, vertI);
1536             }
1537         }
1538     }
1539
1540
1541     if (!matchedAll)
1542     {
1543         SeriousErrorIn
1544         (
1545             "cyclicPolyPatch::order"
1546             "(const primitivePatch&, labelList&, labelList&) const"
1547         )   << "Patch:" << name() << " : "
1548             << "Cannot match vectors to faces on both sides of patch" << endl
1549             << "    Perhaps your faces do not match?"
1550             << " The obj files written contain the current match." << endl
1551             << "    Continuing with incorrect face ordering from now on!"
1552             << endl;
1553
1554             return false;
1555     }
1556
1557
1558     // Set faceMap such that half0 faces get first and corresponding half1
1559     // faces last.
1560     matchAnchors
1561     (
1562         true,                   // report if anchor matching error
1563         pp,
1564         half0ToPatch,
1565         anchors0,
1566         half1ToPatch,
1567         half1Faces,
1568         from1To0,
1569         tols,
1570         faceMap,
1571         rotation
1572     );
1573
1574     // Return false if no change neccesary, true otherwise.
1575
1576     forAll(faceMap, faceI)
1577     {
1578         if (faceMap[faceI] != faceI || rotation[faceI] != 0)
1579         {
1580             return true;
1581         }
1582     }
1583
1584     return false;
1585 }
1586
1587
1588 void Foam::cyclicPolyPatch::write(Ostream& os) const
1589 {
1590     polyPatch::write(os);
1591     os.writeKeyword("featureCos") << featureCos_ << token::END_STATEMENT << nl;
1592     switch (transform_)
1593     {
1594         case ROTATIONAL:
1595         {
1596             os.writeKeyword("transform") << transformTypeNames[ROTATIONAL]
1597                 << token::END_STATEMENT << nl;
1598             os.writeKeyword("rotationAxis") << rotationAxis_
1599                 << token::END_STATEMENT << nl;
1600             os.writeKeyword("rotationCentre") << rotationCentre_
1601                 << token::END_STATEMENT << nl;
1602             break;
1603         }
1604         case TRANSLATIONAL:
1605         {
1606             os.writeKeyword("transform") << transformTypeNames[TRANSLATIONAL]
1607                 << token::END_STATEMENT << nl;
1608             os.writeKeyword("separationVector") << separationVector_
1609                 << token::END_STATEMENT << nl;
1610             break;
1611         }
1612         default:
1613         {
1614             // no additional info to write
1615         }
1616     }
1617 }
1618
1619
1620 // ************************************************************************* //