src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/oldCyclic/oldCyclicPolyPatch.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2004-2011 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
  13     the Free Software Foundation, either version 3 of the License, or
  14     (at your 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, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "oldCyclicPolyPatch.H"
  27 #include "addToRunTimeSelectionTable.H"
  28 #include "polyBoundaryMesh.H"
  29 #include "polyMesh.H"
  30 #include "demandDrivenData.H"
  31 #include "OFstream.H"
  32 #include "patchZones.H"
  33 #include "matchPoints.H"
  34 #include "Time.H"
  35 #include "transformList.H"
  36 #include "cyclicPolyPatch.H"
  37
  38 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  39
  40 namespace Foam
  41 {
  42     defineTypeNameAndDebug(oldCyclicPolyPatch, 0);
  43
  44     addToRunTimeSelectionTable(polyPatch, oldCyclicPolyPatch, word);
  45     addToRunTimeSelectionTable(polyPatch, oldCyclicPolyPatch, dictionary);
  46
  47     template<>
  48     const char* Foam::NamedEnum
  49     <
  50         Foam::oldCyclicPolyPatch::transformType,
  51         3
  52     >::names[] =
  53     {
  54         "unknown",
  55         "rotational",
  56         "translational"
  57     };
  58 }
  59
  60 const Foam::NamedEnum<Foam::oldCyclicPolyPatch::transformType, 3>
  61     Foam::oldCyclicPolyPatch::transformTypeNames;
  62
  63 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  64
  65 Foam::pointField Foam::oldCyclicPolyPatch::calcFaceCentres
  66 (
  67     const UList<face>& faces,
  68     const pointField& points
  69 )
  70 {
  71     pointField ctrs(faces.size());
  72
  73     forAll(faces, faceI)
  74     {
  75         ctrs[faceI] = faces[faceI].centre(points);
  76     }
  77
  78     return ctrs;
  79 }
  80
  81
  82 Foam::pointField Foam::oldCyclicPolyPatch::getAnchorPoints
  83 (
  84     const UList<face>& faces,
  85     const pointField& points
  86 )
  87 {
  88     pointField anchors(faces.size());
  89
  90     forAll(faces, faceI)
  91     {
  92         anchors[faceI] = points[faces[faceI][0]];
  93     }
  94
  95     return anchors;
  96 }
  97
  98
  99 Foam::label Foam::oldCyclicPolyPatch::findMaxArea
 100 (
 101     const pointField& points,
 102     const faceList& faces
 103 )
 104 {
 105     label maxI = -1;
 106     scalar maxAreaSqr = -GREAT;
 107
 108     forAll(faces, faceI)
 109     {
 110         scalar areaSqr = magSqr(faces[faceI].normal(points));
 111
 112         if (areaSqr > maxAreaSqr)
 113         {
 114             maxAreaSqr = areaSqr;
 115             maxI = faceI;
 116         }
 117     }
 118     return maxI;
 119 }
 120
 121
 122 // Get geometric zones of patch by looking at normals.
 123 // Method 1: any edge with sharpish angle is edge between two halves.
 124 //           (this will handle e.g. wedge geometries).
 125 //           Also two fully disconnected regions will be handled this way.
 126 // Method 2: sort faces into two halves based on face normal.
 127 bool Foam::oldCyclicPolyPatch::getGeometricHalves
 128 (
 129     const primitivePatch& pp,
 130     labelList& half0ToPatch,
 131     labelList& half1ToPatch
 132 ) const
 133 {
 134     // Calculate normals
 135     const vectorField& faceNormals = pp.faceNormals();
 136
 137     // Find edges with sharp angles.
 138     boolList regionEdge(pp.nEdges(), false);
 139
 140     const labelListList& edgeFaces = pp.edgeFaces();
 141
 142     label nRegionEdges = 0;
 143
 144     forAll(edgeFaces, edgeI)
 145     {
 146         const labelList& eFaces = edgeFaces[edgeI];
 147
 148         // Check manifold edges for sharp angle.
 149         // (Non-manifold already handled by patchZones)
 150         if (eFaces.size() == 2)
 151         {
 152             if ((faceNormals[eFaces[0]] & faceNormals[eFaces[1]])< featureCos_)
 153             {
 154                 regionEdge[edgeI] = true;
 155
 156                 nRegionEdges++;
 157             }
 158         }
 159     }
 160
 161
 162     // For every face determine zone it is connected to (without crossing
 163     // any regionEdge)
 164     patchZones ppZones(pp, regionEdge);
 165
 166     if (debug)
 167     {
 168         Pout<< "oldCyclicPolyPatch::getGeometricHalves : "
 169             << "Found " << nRegionEdges << " edges on patch " << name()
 170             << " where the cos of the angle between two connected faces"
 171             << " was less than " << featureCos_ << nl
 172             << "Patch divided by these and by single sides edges into "
 173             << ppZones.nZones() << " parts." << endl;
 174
 175
 176         // Dumping zones to obj files.
 177
 178         labelList nZoneFaces(ppZones.nZones());
 179
 180         for (label zoneI = 0; zoneI < ppZones.nZones(); zoneI++)
 181         {
 182             OFstream stream
 183             (
 184                 boundaryMesh().mesh().time().path()
 185                /name()+"_zone_"+Foam::name(zoneI)+".obj"
 186             );
 187             Pout<< "oldCyclicPolyPatch::getGeometricHalves : Writing zone "
 188                 << zoneI << " face centres to OBJ file " << stream.name()
 189                 << endl;
 190
 191             labelList zoneFaces(findIndices(ppZones, zoneI));
 192
 193             forAll(zoneFaces, i)
 194             {
 195                 writeOBJ(stream, pp[zoneFaces[i]].centre(pp.points()));
 196             }
 197
 198             nZoneFaces[zoneI] = zoneFaces.size();
 199         }
 200     }
 201
 202
 203     if (ppZones.nZones() == 2)
 204     {
 205         half0ToPatch = findIndices(ppZones, 0);
 206         half1ToPatch = findIndices(ppZones, 1);
 207     }
 208     else
 209     {
 210         if (debug)
 211         {
 212             Pout<< "oldCyclicPolyPatch::getGeometricHalves :"
 213                 << " falling back to face-normal comparison" << endl;
 214         }
 215         label n0Faces = 0;
 216         half0ToPatch.setSize(pp.size());
 217
 218         label n1Faces = 0;
 219         half1ToPatch.setSize(pp.size());
 220
 221         // Compare to face 0 normal.
 222         forAll(faceNormals, faceI)
 223         {
 224             if ((faceNormals[faceI] & faceNormals[0]) > 0)
 225             {
 226                 half0ToPatch[n0Faces++] = faceI;
 227             }
 228             else
 229             {
 230                 half1ToPatch[n1Faces++] = faceI;
 231             }
 232         }
 233         half0ToPatch.setSize(n0Faces);
 234         half1ToPatch.setSize(n1Faces);
 235
 236         if (debug)
 237         {
 238             Pout<< "oldCyclicPolyPatch::getGeometricHalves :"
 239                 << " Number of faces per zone:("
 240                 << n0Faces << ' ' << n1Faces << ')' << endl;
 241         }
 242     }
 243
 244     if (half0ToPatch.size() != half1ToPatch.size())
 245     {
 246         fileName casePath(boundaryMesh().mesh().time().path());
 247
 248         // Dump halves
 249         {
 250             fileName nm0(casePath/name()+"_half0_faces.obj");
 251             Pout<< "oldCyclicPolyPatch::getGeometricHalves : Writing half0"
 252                 << " faces to OBJ file " << nm0 << endl;
 253             writeOBJ(nm0, UIndirectList<face>(pp, half0ToPatch)(), pp.points());
 254
 255             fileName nm1(casePath/name()+"_half1_faces.obj");
 256             Pout<< "oldCyclicPolyPatch::getGeometricHalves : Writing half1"
 257                 << " faces to OBJ file " << nm1 << endl;
 258             writeOBJ(nm1, UIndirectList<face>(pp, half1ToPatch)(), pp.points());
 259         }
 260
 261         // Dump face centres
 262         {
 263             OFstream str0(casePath/name()+"_half0.obj");
 264             Pout<< "oldCyclicPolyPatch::getGeometricHalves : Writing half0"
 265                 << " face centres to OBJ file " << str0.name() << endl;
 266
 267             forAll(half0ToPatch, i)
 268             {
 269                 writeOBJ(str0, pp[half0ToPatch[i]].centre(pp.points()));
 270             }
 271
 272             OFstream str1(casePath/name()+"_half1.obj");
 273             Pout<< "oldCyclicPolyPatch::getGeometricHalves : Writing half1"
 274                 << " face centres to OBJ file " << str1.name() << endl;
 275             forAll(half1ToPatch, i)
 276             {
 277                 writeOBJ(str1, pp[half1ToPatch[i]].centre(pp.points()));
 278             }
 279         }
 280
 281         SeriousErrorIn
 282         (
 283             "oldCyclicPolyPatch::getGeometricHalves"
 284             "(const primitivePatch&, labelList&, labelList&) const"
 285         )   << "Patch " << name() << " gets decomposed in two zones of"
 286             << "inequal size: " << half0ToPatch.size()
 287             << " and " << half1ToPatch.size() << endl
 288             << "This means that the patch is either not two separate regions"
 289             << " or one region where the angle between the different regions"
 290             << " is not sufficiently sharp." << endl
 291             << "Please adapt the featureCos setting." << endl
 292             << "Continuing with incorrect face ordering from now on!" << endl;
 293
 294         return false;
 295     }
 296     else
 297     {
 298         return true;
 299     }
 300 }
 301
 302
 303 // Given a split of faces into left and right half calculate the centres
 304 // and anchor points. Transform the left points so they align with the
 305 // right ones.
 306 void Foam::oldCyclicPolyPatch::getCentresAndAnchors
 307 (
 308     const primitivePatch& pp,
 309     const faceList& half0Faces,
 310     const faceList& half1Faces,
 311
 312     pointField& ppPoints,
 313     pointField& half0Ctrs,
 314     pointField& half1Ctrs,
 315     pointField& anchors0,
 316     scalarField& tols
 317 ) const
 318 {
 319     // Get geometric data on both halves.
 320     half0Ctrs = calcFaceCentres(half0Faces, pp.points());
 321     anchors0 = getAnchorPoints(half0Faces, pp.points());
 322     half1Ctrs = calcFaceCentres(half1Faces, pp.points());
 323
 324     switch (transform_)
 325     {
 326         case ROTATIONAL:
 327         {
 328             label face0 = getConsistentRotationFace(half0Ctrs);
 329             label face1 = getConsistentRotationFace(half1Ctrs);
 330
 331             vector n0 = ((half0Ctrs[face0] - rotationCentre_) ^ rotationAxis_);
 332             vector n1 = ((half1Ctrs[face1] - rotationCentre_) ^ -rotationAxis_);
 333             n0 /= mag(n0) + VSMALL;
 334             n1 /= mag(n1) + VSMALL;
 335
 336             if (debug)
 337             {
 338                 Pout<< "oldCyclicPolyPatch::getCentresAndAnchors :"
 339                     << " Specified rotation :"
 340                     << " n0:" << n0 << " n1:" << n1 << endl;
 341             }
 342
 343             // Rotation (around origin)
 344             const tensor reverseT(rotationTensor(n0, -n1));
 345
 346             // Rotation
 347             forAll(half0Ctrs, faceI)
 348             {
 349                 half0Ctrs[faceI] = Foam::transform(reverseT, half0Ctrs[faceI]);
 350                 anchors0[faceI] = Foam::transform(reverseT, anchors0[faceI]);
 351             }
 352
 353             ppPoints = Foam::transform(reverseT, pp.points());
 354
 355             break;
 356         }
 357         //- Problem: usually specified translation is not accurate enough
 358         //- to get proper match so keep automatic determination over here.
 359         //case TRANSLATIONAL:
 360         //{
 361         //    // Transform 0 points.
 362         //
 363         //    if (debug)
 364         //    {
 365         //        Pout<< "oldCyclicPolyPatch::getCentresAndAnchors :"
 366         //            << "Specified translation : " << separationVector_
 367         //            << endl;
 368         //    }
 369         //
 370         //    half0Ctrs += separationVector_;
 371         //    anchors0 += separationVector_;
 372         //    break;
 373         //}
 374         default:
 375         {
 376             // Assumes that cyclic is planar. This is also the initial
 377             // condition for patches without faces.
 378
 379             // Determine the face with max area on both halves. These
 380             // two faces are used to determine the transformation tensors
 381             label max0I = findMaxArea(pp.points(), half0Faces);
 382             vector n0 = half0Faces[max0I].normal(pp.points());
 383             n0 /= mag(n0) + VSMALL;
 384
 385             label max1I = findMaxArea(pp.points(), half1Faces);
 386             vector n1 = half1Faces[max1I].normal(pp.points());
 387             n1 /= mag(n1) + VSMALL;
 388
 389             if (mag(n0 & n1) < 1-matchTolerance())
 390             {
 391                 if (debug)
 392                 {
 393                     Pout<< "oldCyclicPolyPatch::getCentresAndAnchors :"
 394                         << " Detected rotation :"
 395                         << " n0:" << n0 << " n1:" << n1 << endl;
 396                 }
 397
 398                 // Rotation (around origin)
 399                 const tensor reverseT(rotationTensor(n0, -n1));
 400
 401                 // Rotation
 402                 forAll(half0Ctrs, faceI)
 403                 {
 404                     half0Ctrs[faceI] = Foam::transform
 405                     (
 406                         reverseT,
 407                         half0Ctrs[faceI]
 408                     );
 409                     anchors0[faceI] = Foam::transform
 410                     (
 411                         reverseT,
 412                         anchors0[faceI]
 413                     );
 414                 }
 415                 ppPoints = Foam::transform(reverseT, pp.points());
 416             }
 417             else
 418             {
 419                 // Parallel translation. Get average of all used points.
 420
 421                 primitiveFacePatch half0(half0Faces, pp.points());
 422                 const pointField& half0Pts = half0.localPoints();
 423                 const point ctr0(sum(half0Pts)/half0Pts.size());
 424
 425                 primitiveFacePatch half1(half1Faces, pp.points());
 426                 const pointField& half1Pts = half1.localPoints();
 427                 const point ctr1(sum(half1Pts)/half1Pts.size());
 428
 429                 if (debug)
 430                 {
 431                     Pout<< "oldCyclicPolyPatch::getCentresAndAnchors :"
 432                         << " Detected translation :"
 433                         << " n0:" << n0 << " n1:" << n1
 434                         << " ctr0:" << ctr0 << " ctr1:" << ctr1 << endl;
 435                 }
 436
 437                 half0Ctrs += ctr1 - ctr0;
 438                 anchors0 += ctr1 - ctr0;
 439                 ppPoints = pp.points() + ctr1 - ctr0;
 440             }
 441             break;
 442         }
 443     }
 444
 445
 446     // Calculate typical distance per face
 447     tols = calcFaceTol(matchTolerance(), half1Faces, pp.points(), half1Ctrs);
 448 }
 449
 450
 451 // Calculates faceMap and rotation. Returns true if all ok.
 452 bool Foam::oldCyclicPolyPatch::matchAnchors
 453 (
 454     const bool report,
 455     const primitivePatch& pp,
 456     const labelList& half0ToPatch,
 457     const pointField& anchors0,
 458
 459     const labelList& half1ToPatch,
 460     const faceList& half1Faces,
 461     const labelList& from1To0,
 462
 463     const scalarField& tols,
 464
 465     labelList& faceMap,
 466     labelList& rotation
 467 ) const
 468 {
 469     // Set faceMap such that half0 faces get first and corresponding half1
 470     // faces last.
 471
 472     forAll(half0ToPatch, half0FaceI)
 473     {
 474         // Label in original patch
 475         label patchFaceI = half0ToPatch[half0FaceI];
 476
 477         faceMap[patchFaceI] = half0FaceI;
 478
 479         // No rotation
 480         rotation[patchFaceI] = 0;
 481     }
 482
 483     bool fullMatch = true;
 484
 485     forAll(from1To0, half1FaceI)
 486     {
 487         label patchFaceI = half1ToPatch[half1FaceI];
 488
 489         // This face has to match the corresponding one on half0.
 490         label half0FaceI = from1To0[half1FaceI];
 491
 492         label newFaceI = half0FaceI + pp.size()/2;
 493
 494         faceMap[patchFaceI] = newFaceI;
 495
 496         // Rotate patchFaceI such that its f[0] aligns with that of
 497         // the corresponding face
 498         // (which after shuffling will be at position half0FaceI)
 499
 500         const point& wantedAnchor = anchors0[half0FaceI];
 501
 502         rotation[newFaceI] = getRotation
 503         (
 504             pp.points(),
 505             half1Faces[half1FaceI],
 506             wantedAnchor,
 507             tols[half1FaceI]
 508         );
 509
 510         if (rotation[newFaceI] == -1)
 511         {
 512             fullMatch = false;
 513
 514             if (report)
 515             {
 516                 const face& f = half1Faces[half1FaceI];
 517                 SeriousErrorIn
 518                 (
 519                     "oldCyclicPolyPatch::matchAnchors(..)"
 520                 )   << "Patch:" << name() << " : "
 521                     << "Cannot find point on face " << f
 522                     << " with vertices:"
 523                     << UIndirectList<point>(pp.points(), f)()
 524                     << " that matches point " << wantedAnchor
 525                     << " when matching the halves of cyclic patch " << name()
 526                     << endl
 527                     << "Continuing with incorrect face ordering from now on!"
 528                     << endl;
 529             }
 530         }
 531     }
 532     return fullMatch;
 533 }
 534
 535
 536 Foam::label Foam::oldCyclicPolyPatch::getConsistentRotationFace
 537 (
 538     const pointField& faceCentres
 539 ) const
 540 {
 541     const scalarField magRadSqr
 542     (
 543         magSqr((faceCentres - rotationCentre_) ^ rotationAxis_)
 544     );
 545     scalarField axisLen
 546     (
 547         (faceCentres - rotationCentre_) & rotationAxis_
 548     );
 549     axisLen = axisLen - min(axisLen);
 550     const scalarField magLenSqr
 551     (
 552         magRadSqr + axisLen*axisLen
 553     );
 554
 555     label rotFace = -1;
 556     scalar maxMagLenSqr = -GREAT;
 557     scalar maxMagRadSqr = -GREAT;
 558     forAll(faceCentres, i)
 559     {
 560         if (magLenSqr[i] >= maxMagLenSqr)
 561         {
 562             if (magRadSqr[i] > maxMagRadSqr)
 563             {
 564                 rotFace = i;
 565                 maxMagLenSqr = magLenSqr[i];
 566                 maxMagRadSqr = magRadSqr[i];
 567             }
 568         }
 569     }
 570
 571     if (debug)
 572     {
 573         Info<< "getConsistentRotationFace(const pointField&)" << nl
 574             << "    rotFace = " << rotFace << nl
 575             << "    point =  " << faceCentres[rotFace] << endl;
 576     }
 577
 578     return rotFace;
 579 }
 580
 581
 582 // * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * * * * * //
 583
 584 Foam::oldCyclicPolyPatch::oldCyclicPolyPatch
 585 (
 586     const word& name,
 587     const label size,
 588     const label start,
 589     const label index,
 590     const polyBoundaryMesh& bm
 591 )
 592 :
 593     coupledPolyPatch(name, size, start, index, bm),
 594     featureCos_(0.9),
 595     transform_(UNKNOWN),
 596     rotationAxis_(vector::zero),
 597     rotationCentre_(point::zero),
 598     separationVector_(vector::zero)
 599 {}
 600
 601
 602 Foam::oldCyclicPolyPatch::oldCyclicPolyPatch
 603 (
 604     const word& name,
 605     const dictionary& dict,
 606     const label index,
 607     const polyBoundaryMesh& bm
 608 )
 609 :
 610     coupledPolyPatch(name, dict, index, bm),
 611     featureCos_(0.9),
 612     transform_(UNKNOWN),
 613     rotationAxis_(vector::zero),
 614     rotationCentre_(point::zero),
 615     separationVector_(vector::zero)
 616 {
 617     if (dict.found("neighbourPatch"))
 618     {
 619         FatalIOErrorIn
 620         (
 621             "oldCyclicPolyPatch::oldCyclicPolyPatch\n"
 622             "(\n"
 623             "    const word& name,\n"
 624             "    const dictionary& dict,\n"
 625             "    const label index,\n"
 626             "    const polyBoundaryMesh& bm\n"
 627             ")",
 628             dict
 629         )   << "Found \"neighbourPatch\" entry when reading cyclic patch "
 630             << name << endl
 631             << "Is this mesh already with split cyclics?" << endl
 632             << "If so run a newer version that supports it"
 633             << ", if not comment out the \"neighbourPatch\" entry and re-run"
 634             << exit(FatalIOError);
 635     }
 636
 637     dict.readIfPresent("featureCos", featureCos_);
 638
 639     if (dict.found("transform"))
 640     {
 641         transform_ = transformTypeNames.read(dict.lookup("transform"));
 642         switch (transform_)
 643         {
 644             case ROTATIONAL:
 645             {
 646                 dict.lookup("rotationAxis") >> rotationAxis_;
 647                 dict.lookup("rotationCentre") >> rotationCentre_;
 648                 break;
 649             }
 650             case TRANSLATIONAL:
 651             {
 652                 dict.lookup("separationVector") >> separationVector_;
 653                 break;
 654             }
 655             default:
 656             {
 657                 // no additional info required
 658             }
 659         }
 660     }
 661 }
 662
 663
 664 Foam::oldCyclicPolyPatch::oldCyclicPolyPatch
 665 (
 666     const oldCyclicPolyPatch& pp,
 667     const polyBoundaryMesh& bm
 668 )
 669 :
 670     coupledPolyPatch(pp, bm),
 671     featureCos_(pp.featureCos_),
 672     transform_(pp.transform_),
 673     rotationAxis_(pp.rotationAxis_),
 674     rotationCentre_(pp.rotationCentre_),
 675     separationVector_(pp.separationVector_)
 676 {}
 677
 678
 679 Foam::oldCyclicPolyPatch::oldCyclicPolyPatch
 680 (
 681     const oldCyclicPolyPatch& pp,
 682     const polyBoundaryMesh& bm,
 683     const label index,
 684     const label newSize,
 685     const label newStart
 686 )
 687 :
 688     coupledPolyPatch(pp, bm, index, newSize, newStart),
 689     featureCos_(pp.featureCos_),
 690     transform_(pp.transform_),
 691     rotationAxis_(pp.rotationAxis_),
 692     rotationCentre_(pp.rotationCentre_),
 693     separationVector_(pp.separationVector_)
 694 {}
 695
 696
 697 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 698
 699 Foam::oldCyclicPolyPatch::~oldCyclicPolyPatch()
 700 {}
 701
 702
 703 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 704
 705 void Foam::oldCyclicPolyPatch::initGeometry(PstreamBuffers& pBufs)
 706 {
 707     polyPatch::initGeometry(pBufs);
 708 }
 709
 710
 711 void Foam::oldCyclicPolyPatch::calcGeometry
 712 (
 713     const primitivePatch& referPatch,
 714     const pointField& thisCtrs,
 715     const vectorField& thisAreas,
 716     const pointField& thisCc,
 717     const pointField& nbrCtrs,
 718     const vectorField& nbrAreas,
 719     const pointField& nbrCc
 720 )
 721 {}
 722
 723
 724 void Foam::oldCyclicPolyPatch::calcGeometry(PstreamBuffers& pBufs)
 725 {}
 726
 727
 728 void Foam::oldCyclicPolyPatch::initMovePoints
 729 (
 730     PstreamBuffers& pBufs,
 731     const pointField& p
 732 )
 733 {
 734     polyPatch::initMovePoints(pBufs, p);
 735 }
 736
 737
 738 void Foam::oldCyclicPolyPatch::movePoints
 739 (
 740     PstreamBuffers& pBufs,
 741     const pointField& p
 742 )
 743 {
 744     polyPatch::movePoints(pBufs, p);
 745 }
 746
 747
 748 void Foam::oldCyclicPolyPatch::initUpdateMesh(PstreamBuffers& pBufs)
 749 {
 750     polyPatch::initUpdateMesh(pBufs);
 751 }
 752
 753
 754 void Foam::oldCyclicPolyPatch::updateMesh(PstreamBuffers& pBufs)
 755 {
 756     polyPatch::updateMesh(pBufs);
 757 }
 758
 759
 760 void Foam::oldCyclicPolyPatch::initOrder
 761 (
 762     PstreamBuffers&,
 763     const primitivePatch& pp
 764 ) const
 765 {}
 766
 767
 768 //  Return new ordering. Ordering is -faceMap: for every face index
 769 //  the new face -rotation:for every new face the clockwise shift
 770 //  of the original face. Return false if nothing changes (faceMap
 771 //  is identity, rotation is 0)
 772 bool Foam::oldCyclicPolyPatch::order
 773 (
 774     PstreamBuffers&,
 775     const primitivePatch& pp,
 776     labelList& faceMap,
 777     labelList& rotation
 778 ) const
 779 {
 780     faceMap.setSize(pp.size());
 781     faceMap = -1;
 782
 783     rotation.setSize(pp.size());
 784     rotation = 0;
 785
 786     if (pp.empty())
 787     {
 788         // No faces, nothing to change.
 789         return false;
 790     }
 791
 792     if (pp.size()&1)
 793     {
 794         FatalErrorIn("oldCyclicPolyPatch::order(..)")
 795             << "Size of cyclic " << name() << " should be a multiple of 2"
 796             << ". It is " << pp.size() << abort(FatalError);
 797     }
 798
 799     label halfSize = pp.size()/2;
 800
 801     // Supplied primitivePatch already with new points.
 802     // Cyclics are limited to one transformation tensor
 803     // currently anyway (i.e. straight plane) so should not be too big a
 804     // problem.
 805
 806
 807     // Indices of faces on half0
 808     labelList half0ToPatch;
 809     // Indices of faces on half1
 810     labelList half1ToPatch;
 811
 812
 813     // 1. Test if already correctly oriented by starting from trivial ordering.
 814     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 815
 816     half0ToPatch = identity(halfSize);
 817     half1ToPatch = half0ToPatch + halfSize;
 818
 819     // Get faces
 820     faceList half0Faces(UIndirectList<face>(pp, half0ToPatch));
 821     faceList half1Faces(UIndirectList<face>(pp, half1ToPatch));
 822
 823     // Get geometric quantities
 824     pointField half0Ctrs, half1Ctrs, anchors0, ppPoints;
 825     scalarField tols;
 826     getCentresAndAnchors
 827     (
 828         pp,
 829         half0Faces,
 830         half1Faces,
 831
 832         ppPoints,
 833         half0Ctrs,
 834         half1Ctrs,
 835         anchors0,
 836         tols
 837     );
 838
 839     // Geometric match of face centre vectors
 840     labelList from1To0;
 841     bool matchedAll = matchPoints
 842     (
 843         half1Ctrs,
 844         half0Ctrs,
 845         tols,
 846         false,
 847         from1To0
 848     );
 849
 850     if (debug)
 851     {
 852         Pout<< "oldCyclicPolyPatch::order : test if already ordered:"
 853             << matchedAll << endl;
 854
 855         // Dump halves
 856         fileName nm0("match1_"+name()+"_half0_faces.obj");
 857         Pout<< "oldCyclicPolyPatch::order : Writing half0"
 858             << " faces to OBJ file " << nm0 << endl;
 859         writeOBJ(nm0, half0Faces, ppPoints);
 860
 861         fileName nm1("match1_"+name()+"_half1_faces.obj");
 862         Pout<< "oldCyclicPolyPatch::order : Writing half1"
 863             << " faces to OBJ file " << nm1 << endl;
 864         writeOBJ(nm1, half1Faces, pp.points());
 865
 866         OFstream ccStr
 867         (
 868             boundaryMesh().mesh().time().path()
 869            /"match1_"+ name() + "_faceCentres.obj"
 870         );
 871         Pout<< "oldCyclicPolyPatch::order : "
 872             << "Dumping currently found cyclic match as lines between"
 873             << " corresponding face centres to file " << ccStr.name()
 874             << endl;
 875
 876         // Recalculate untransformed face centres
 877         //pointField rawHalf0Ctrs = calcFaceCentres(half0Faces, pp.points());
 878         label vertI = 0;
 879
 880         forAll(half1Ctrs, i)
 881         {
 882             //if (from1To0[i] != -1)
 883             {
 884                 // Write edge between c1 and c0
 885                 //const point& c0 = rawHalf0Ctrs[from1To0[i]];
 886                 //const point& c0 = half0Ctrs[from1To0[i]];
 887                 const point& c0 = half0Ctrs[i];
 888                 const point& c1 = half1Ctrs[i];
 889                 writeOBJ(ccStr, c0, c1, vertI);
 890             }
 891         }
 892     }
 893
 894
 895     // 2. Ordered in pairs (so 0,1 coupled and 2,3 etc.)
 896     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 897
 898     if (!matchedAll)
 899     {
 900         label faceI = 0;
 901         for (label i = 0; i < halfSize; i++)
 902         {
 903             half0ToPatch[i] = faceI++;
 904             half1ToPatch[i] = faceI++;
 905         }
 906
 907         // And redo all matching
 908         half0Faces = UIndirectList<face>(pp, half0ToPatch);
 909         half1Faces = UIndirectList<face>(pp, half1ToPatch);
 910
 911         getCentresAndAnchors
 912         (
 913             pp,
 914             half0Faces,
 915             half1Faces,
 916
 917             ppPoints,
 918             half0Ctrs,
 919             half1Ctrs,
 920             anchors0,
 921             tols
 922         );
 923
 924         // Geometric match of face centre vectors
 925         matchedAll = matchPoints
 926         (
 927             half1Ctrs,
 928             half0Ctrs,
 929             tols,
 930             false,
 931             from1To0
 932         );
 933
 934         if (debug)
 935         {
 936             Pout<< "oldCyclicPolyPatch::order : test if pairwise ordered:"
 937                 << matchedAll << endl;
 938             // Dump halves
 939             fileName nm0("match2_"+name()+"_half0_faces.obj");
 940             Pout<< "oldCyclicPolyPatch::order : Writing half0"
 941                 << " faces to OBJ file " << nm0 << endl;
 942             writeOBJ(nm0, half0Faces, ppPoints);
 943
 944             fileName nm1("match2_"+name()+"_half1_faces.obj");
 945             Pout<< "oldCyclicPolyPatch::order : Writing half1"
 946                 << " faces to OBJ file " << nm1 << endl;
 947             writeOBJ(nm1, half1Faces, pp.points());
 948
 949             OFstream ccStr
 950             (
 951                 boundaryMesh().mesh().time().path()
 952                /"match2_"+name()+"_faceCentres.obj"
 953             );
 954             Pout<< "oldCyclicPolyPatch::order : "
 955                 << "Dumping currently found cyclic match as lines between"
 956                 << " corresponding face centres to file " << ccStr.name()
 957                 << endl;
 958
 959             // Recalculate untransformed face centres
 960             label vertI = 0;
 961
 962             forAll(half1Ctrs, i)
 963             {
 964                 if (from1To0[i] != -1)
 965                 {
 966                     // Write edge between c1 and c0
 967                     const point& c0 = half0Ctrs[from1To0[i]];
 968                     const point& c1 = half1Ctrs[i];
 969                     writeOBJ(ccStr, c0, c1, vertI);
 970                 }
 971             }
 972         }
 973     }
 974
 975
 976     // 3. Baffles(coincident faces) converted into cyclics (e.g. jump)
 977     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 978
 979     if (!matchedAll)
 980     {
 981         label baffleI = 0;
 982
 983         forAll(pp, faceI)
 984         {
 985             const face& f = pp.localFaces()[faceI];
 986             const labelList& pFaces = pp.pointFaces()[f[0]];
 987
 988             label matchedFaceI = -1;
 989
 990             forAll(pFaces, i)
 991             {
 992                 label otherFaceI = pFaces[i];
 993
 994                 if (otherFaceI > faceI)
 995                 {
 996                     const face& otherF = pp.localFaces()[otherFaceI];
 997
 998                     // Note: might pick up two similar oriented faces
 999                     //       (but that is illegal anyway)
1000                     if (f == otherF)
1001                     {
1002                         matchedFaceI = otherFaceI;
1003                         break;
1004                     }
1005                 }
1006             }
1007
1008             if (matchedFaceI != -1)
1009             {
1010                 half0ToPatch[baffleI] = faceI;
1011                 half1ToPatch[baffleI] = matchedFaceI;
1012                 baffleI++;
1013             }
1014         }
1015
1016         if (baffleI == halfSize)
1017         {
1018             // And redo all matching
1019             half0Faces = UIndirectList<face>(pp, half0ToPatch);
1020             half1Faces = UIndirectList<face>(pp, half1ToPatch);
1021
1022             getCentresAndAnchors
1023             (
1024                 pp,
1025                 half0Faces,
1026                 half1Faces,
1027
1028                 ppPoints,
1029                 half0Ctrs,
1030                 half1Ctrs,
1031                 anchors0,
1032                 tols
1033             );
1034
1035             // Geometric match of face centre vectors
1036             matchedAll = matchPoints
1037             (
1038                 half1Ctrs,
1039                 half0Ctrs,
1040                 tols,
1041                 false,
1042                 from1To0
1043             );
1044
1045             if (debug)
1046             {
1047                 Pout<< "oldCyclicPolyPatch::order : test if baffles:"
1048                     << matchedAll << endl;
1049                 // Dump halves
1050                 fileName nm0("match3_"+name()+"_half0_faces.obj");
1051                 Pout<< "oldCyclicPolyPatch::order : Writing half0"
1052                     << " faces to OBJ file " << nm0 << endl;
1053                 writeOBJ(nm0, half0Faces, ppPoints);
1054
1055                 fileName nm1("match3_"+name()+"_half1_faces.obj");
1056                 Pout<< "oldCyclicPolyPatch::order : Writing half1"
1057                     << " faces to OBJ file " << nm1 << endl;
1058                 writeOBJ(nm1, half1Faces, pp.points());
1059
1060                 OFstream ccStr
1061                 (
1062                     boundaryMesh().mesh().time().path()
1063                    /"match3_"+ name() + "_faceCentres.obj"
1064                 );
1065                 Pout<< "oldCyclicPolyPatch::order : "
1066                     << "Dumping currently found cyclic match as lines between"
1067                     << " corresponding face centres to file " << ccStr.name()
1068                     << endl;
1069
1070                 // Recalculate untransformed face centres
1071                 label vertI = 0;
1072
1073                 forAll(half1Ctrs, i)
1074                 {
1075                     if (from1To0[i] != -1)
1076                     {
1077                         // Write edge between c1 and c0
1078                         const point& c0 = half0Ctrs[from1To0[i]];
1079                         const point& c1 = half1Ctrs[i];
1080                         writeOBJ(ccStr, c0, c1, vertI);
1081                     }
1082                 }
1083             }
1084         }
1085     }
1086
1087
1088     // 4. Automatic geometric ordering
1089     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1090
1091     if (!matchedAll)
1092     {
1093         // Split faces according to feature angle or topology
1094         bool okSplit = getGeometricHalves(pp, half0ToPatch, half1ToPatch);
1095
1096         if (!okSplit)
1097         {
1098             // Did not split into two equal parts.
1099             return false;
1100         }
1101
1102         // And redo all matching
1103         half0Faces = UIndirectList<face>(pp, half0ToPatch);
1104         half1Faces = UIndirectList<face>(pp, half1ToPatch);
1105
1106         getCentresAndAnchors
1107         (
1108             pp,
1109             half0Faces,
1110             half1Faces,
1111
1112             ppPoints,
1113             half0Ctrs,
1114             half1Ctrs,
1115             anchors0,
1116             tols
1117         );
1118
1119         // Geometric match of face centre vectors
1120         matchedAll = matchPoints
1121         (
1122             half1Ctrs,
1123             half0Ctrs,
1124             tols,
1125             false,
1126             from1To0
1127         );
1128
1129         if (debug)
1130         {
1131             Pout<< "oldCyclicPolyPatch::order : automatic ordering result:"
1132                 << matchedAll << endl;
1133             // Dump halves
1134             fileName nm0("match4_"+name()+"_half0_faces.obj");
1135             Pout<< "oldCyclicPolyPatch::order : Writing half0"
1136                 << " faces to OBJ file " << nm0 << endl;
1137             writeOBJ(nm0, half0Faces, ppPoints);
1138
1139             fileName nm1("match4_"+name()+"_half1_faces.obj");
1140             Pout<< "oldCyclicPolyPatch::order : Writing half1"
1141                 << " faces to OBJ file " << nm1 << endl;
1142             writeOBJ(nm1, half1Faces, pp.points());
1143
1144             OFstream ccStr
1145             (
1146                 boundaryMesh().mesh().time().path()
1147                /"match4_"+ name() + "_faceCentres.obj"
1148             );
1149             Pout<< "oldCyclicPolyPatch::order : "
1150                 << "Dumping currently found cyclic match as lines between"
1151                 << " corresponding face centres to file " << ccStr.name()
1152                 << endl;
1153
1154             // Recalculate untransformed face centres
1155             label vertI = 0;
1156
1157             forAll(half1Ctrs, i)
1158             {
1159                 if (from1To0[i] != -1)
1160                 {
1161                     // Write edge between c1 and c0
1162                     const point& c0 = half0Ctrs[from1To0[i]];
1163                     const point& c1 = half1Ctrs[i];
1164                     writeOBJ(ccStr, c0, c1, vertI);
1165                 }
1166             }
1167         }
1168     }
1169
1170
1171     if (!matchedAll || debug)
1172     {
1173         // Dump halves
1174         fileName nm0(name()+"_half0_faces.obj");
1175         Pout<< "oldCyclicPolyPatch::order : Writing half0"
1176             << " faces to OBJ file " << nm0 << endl;
1177         writeOBJ(nm0, half0Faces, pp.points());
1178
1179         fileName nm1(name()+"_half1_faces.obj");
1180         Pout<< "oldCyclicPolyPatch::order : Writing half1"
1181             << " faces to OBJ file " << nm1 << endl;
1182         writeOBJ(nm1, half1Faces, pp.points());
1183
1184         OFstream ccStr
1185         (
1186             boundaryMesh().mesh().time().path()
1187            /name() + "_faceCentres.obj"
1188         );
1189         Pout<< "oldCyclicPolyPatch::order : "
1190             << "Dumping currently found cyclic match as lines between"
1191             << " corresponding face centres to file " << ccStr.name()
1192             << endl;
1193
1194         // Recalculate untransformed face centres
1195         //pointField rawHalf0Ctrs = calcFaceCentres(half0Faces, pp.points());
1196         label vertI = 0;
1197
1198         forAll(half1Ctrs, i)
1199         {
1200             if (from1To0[i] != -1)
1201             {
1202                 // Write edge between c1 and c0
1203                 //const point& c0 = rawHalf0Ctrs[from1To0[i]];
1204                 const point& c0 = half0Ctrs[from1To0[i]];
1205                 const point& c1 = half1Ctrs[i];
1206                 writeOBJ(ccStr, c0, c1, vertI);
1207             }
1208         }
1209     }
1210
1211
1212     if (!matchedAll)
1213     {
1214         SeriousErrorIn
1215         (
1216             "oldCyclicPolyPatch::order"
1217             "(const primitivePatch&, labelList&, labelList&) const"
1218         )   << "Patch:" << name() << " : "
1219             << "Cannot match vectors to faces on both sides of patch" << endl
1220             << "    Perhaps your faces do not match?"
1221             << " The obj files written contain the current match." << endl
1222             << "    Continuing with incorrect face ordering from now on!"
1223             << endl;
1224
1225             return false;
1226     }
1227
1228
1229     // Set faceMap such that half0 faces get first and corresponding half1
1230     // faces last.
1231     matchAnchors
1232     (
1233         true,                   // report if anchor matching error
1234         pp,
1235         half0ToPatch,
1236         anchors0,
1237         half1ToPatch,
1238         half1Faces,
1239         from1To0,
1240         tols,
1241         faceMap,
1242         rotation
1243     );
1244
1245     // Return false if no change neccesary, true otherwise.
1246
1247     forAll(faceMap, faceI)
1248     {
1249         if (faceMap[faceI] != faceI || rotation[faceI] != 0)
1250         {
1251             return true;
1252         }
1253     }
1254
1255     return false;
1256 }
1257
1258
1259 void Foam::oldCyclicPolyPatch::write(Ostream& os) const
1260 {
1261     // Replacement of polyPatch::write to write 'cyclic' instead of type():
1262     os.writeKeyword("type") << cyclicPolyPatch::typeName
1263         << token::END_STATEMENT << nl;
1264     patchIdentifier::write(os);
1265     os.writeKeyword("nFaces") << size() << token::END_STATEMENT << nl;
1266     os.writeKeyword("startFace") << start() << token::END_STATEMENT << nl;
1267
1268
1269     os.writeKeyword("featureCos") << featureCos_ << token::END_STATEMENT << nl;
1270     switch (transform_)
1271     {
1272         case ROTATIONAL:
1273         {
1274             os.writeKeyword("transform") << transformTypeNames[transform_]
1275                 << token::END_STATEMENT << nl;
1276             os.writeKeyword("rotationAxis") << rotationAxis_
1277                 << token::END_STATEMENT << nl;
1278             os.writeKeyword("rotationCentre") << rotationCentre_
1279                 << token::END_STATEMENT << nl;
1280             break;
1281         }
1282         case TRANSLATIONAL:
1283         {
1284             os.writeKeyword("transform") << transformTypeNames[transform_]
1285                 << token::END_STATEMENT << nl;
1286             os.writeKeyword("separationVector") << separationVector_
1287                 << token::END_STATEMENT << nl;
1288             break;
1289         }
1290         default:
1291         {
1292             // no additional info to write
1293         }
1294     }
1295
1296     WarningIn("oldCyclicPolyPatch::write(Ostream& os) const")
1297         << "Please run foamUpgradeCyclics to convert these old-style"
1298         << " cyclics into two separate cyclics patches."
1299         << endl;
1300 }
1301
1302
1303 // ************************************************************************* //