src/dynamicMesh/polyMeshAdder/faceCoupleInfo.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
   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 "faceCoupleInfo.H"
  27 #include "polyMesh.H"
  28 #include "matchPoints.H"
  29 #include "indirectPrimitivePatch.H"
  30 #include "meshTools.H"
  31 #include "treeDataFace.H"
  32 #include "indexedOctree.H"
  33 #include "OFstream.H"
  34
  35 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  36
  37 defineTypeNameAndDebug(Foam::faceCoupleInfo, 0);
  38
  39 const Foam::scalar Foam::faceCoupleInfo::angleTol_ = 1E-3;
  40
  41
  42 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  43
  44 //- Write edges
  45 void Foam::faceCoupleInfo::writeOBJ
  46 (
  47     const fileName& fName,
  48     const edgeList& edges,
  49     const pointField& points,
  50     const bool compact
  51 )
  52 {
  53     OFstream str(fName);
  54
  55     labelList pointMap(points.size(), -1);
  56
  57     if (compact)
  58     {
  59         label newPointI = 0;
  60
  61         forAll(edges, edgeI)
  62         {
  63             const edge& e = edges[edgeI];
  64
  65             forAll(e, eI)
  66             {
  67                 label pointI = e[eI];
  68
  69                 if (pointMap[pointI] == -1)
  70                 {
  71                     pointMap[pointI] = newPointI++;
  72
  73                     meshTools::writeOBJ(str, points[pointI]);
  74                 }
  75             }
  76         }
  77     }
  78     else
  79     {
  80         forAll(points, pointI)
  81         {
  82             meshTools::writeOBJ(str, points[pointI]);
  83         }
  84
  85         pointMap = identity(points.size());
  86     }
  87
  88     forAll(edges, edgeI)
  89     {
  90         const edge& e = edges[edgeI];
  91
  92         str<< "l " << pointMap[e[0]]+1 << ' ' << pointMap[e[1]]+1 << nl;
  93     }
  94 }
  95
  96
  97 //- Writes edges.
  98 void Foam::faceCoupleInfo::writeOBJ
  99 (
 100     const fileName& fName,
 101     const pointField& points0,
 102     const pointField& points1
 103 )
 104 {
 105     Pout<< "Writing connections as edges to " << fName << endl;
 106
 107     OFstream str(fName);
 108
 109     label vertI = 0;
 110
 111     forAll(points0, i)
 112     {
 113         meshTools::writeOBJ(str, points0[i]);
 114         vertI++;
 115         meshTools::writeOBJ(str, points1[i]);
 116         vertI++;
 117         str << "l " << vertI-1 << ' ' << vertI << nl;
 118     }
 119 }
 120
 121
 122 //- Writes face and point connectivity as .obj files.
 123 void Foam::faceCoupleInfo::writePointsFaces() const
 124 {
 125     const indirectPrimitivePatch& m = masterPatch();
 126     const indirectPrimitivePatch& s = slavePatch();
 127     const primitiveFacePatch& c = cutFaces();
 128
 129     // Patches
 130     {
 131         OFstream str("masterPatch.obj");
 132         Pout<< "Writing masterPatch to " << str.name() << endl;
 133         meshTools::writeOBJ(str, m.localFaces(), m.localPoints());
 134     }
 135     {
 136         OFstream str("slavePatch.obj");
 137         Pout<< "Writing slavePatch to " << str.name() << endl;
 138         meshTools::writeOBJ(str, s.localFaces(), s.localPoints());
 139     }
 140     {
 141         OFstream str("cutFaces.obj");
 142         Pout<< "Writing cutFaces to " << str.name() << endl;
 143         meshTools::writeOBJ(str, c.localFaces(), c.localPoints());
 144     }
 145
 146     // Point connectivity
 147     {
 148         Pout<< "Writing cutToMasterPoints to cutToMasterPoints.obj" << endl;
 149
 150         writeOBJ
 151         (
 152             "cutToMasterPoints.obj",
 153             m.localPoints(),
 154             pointField(c.localPoints(), masterToCutPoints_));
 155     }
 156     {
 157         Pout<< "Writing cutToSlavePoints to cutToSlavePoints.obj" << endl;
 158
 159         writeOBJ
 160         (
 161             "cutToSlavePoints.obj",
 162             s.localPoints(),
 163             pointField(c.localPoints(), slaveToCutPoints_)
 164         );
 165     }
 166
 167     // Face connectivity
 168     {
 169         Pout<< "Writing cutToMasterFaces to cutToMasterFaces.obj" << endl;
 170
 171         pointField equivMasterFaces(c.size());
 172
 173         forAll(cutToMasterFaces(), cutFaceI)
 174         {
 175             label masterFaceI = cutToMasterFaces()[cutFaceI];
 176
 177             if (masterFaceI != -1)
 178             {
 179                 equivMasterFaces[cutFaceI] = m[masterFaceI].centre(m.points());
 180             }
 181             else
 182             {
 183                 WarningIn("writePointsFaces()")
 184                     << "No master face for cut face " << cutFaceI
 185                     << " at position " << c[cutFaceI].centre(c.points())
 186                     << endl;
 187
 188                 equivMasterFaces[cutFaceI] = vector::zero;
 189             }
 190         }
 191
 192         writeOBJ
 193         (
 194             "cutToMasterFaces.obj",
 195             calcFaceCentres<List>(c, cutPoints(), 0, c.size()),
 196             equivMasterFaces
 197         );
 198     }
 199
 200     {
 201         Pout<< "Writing cutToSlaveFaces to cutToSlaveFaces.obj" << endl;
 202
 203         pointField equivSlaveFaces(c.size());
 204
 205         forAll(cutToSlaveFaces(), cutFaceI)
 206         {
 207             label slaveFaceI = cutToSlaveFaces()[cutFaceI];
 208
 209             equivSlaveFaces[cutFaceI] = s[slaveFaceI].centre(s.points());
 210         }
 211
 212         writeOBJ
 213         (
 214             "cutToSlaveFaces.obj",
 215             calcFaceCentres<List>(c, cutPoints(), 0, c.size()),
 216             equivSlaveFaces
 217         );
 218     }
 219
 220     Pout<< endl;
 221 }
 222
 223
 224 void Foam::faceCoupleInfo::writeEdges
 225 (
 226     const labelList& cutToMasterEdges,
 227     const labelList& cutToSlaveEdges
 228 ) const
 229 {
 230     const indirectPrimitivePatch& m = masterPatch();
 231     const indirectPrimitivePatch& s = slavePatch();
 232     const primitiveFacePatch& c = cutFaces();
 233
 234     // Edge connectivity
 235     {
 236         OFstream str("cutToMasterEdges.obj");
 237         Pout<< "Writing cutToMasterEdges to " << str.name() << endl;
 238
 239         label vertI = 0;
 240
 241         forAll(cutToMasterEdges, cutEdgeI)
 242         {
 243             if (cutToMasterEdges[cutEdgeI] != -1)
 244             {
 245                 const edge& masterEdge =
 246                     m.edges()[cutToMasterEdges[cutEdgeI]];
 247                 const edge& cutEdge = c.edges()[cutEdgeI];
 248
 249                 meshTools::writeOBJ(str, m.localPoints()[masterEdge[0]]);
 250                 vertI++;
 251                 meshTools::writeOBJ(str, m.localPoints()[masterEdge[1]]);
 252                 vertI++;
 253                 meshTools::writeOBJ(str, c.localPoints()[cutEdge[0]]);
 254                 vertI++;
 255                 meshTools::writeOBJ(str, c.localPoints()[cutEdge[1]]);
 256                 vertI++;
 257                 str << "l " << vertI-3 << ' ' << vertI-2 << nl;
 258                 str << "l " << vertI-3 << ' ' << vertI-1 << nl;
 259                 str << "l " << vertI-3 << ' ' << vertI << nl;
 260                 str << "l " << vertI-2 << ' ' << vertI-1 << nl;
 261                 str << "l " << vertI-2 << ' ' << vertI << nl;
 262                 str << "l " << vertI-1 << ' ' << vertI << nl;
 263             }
 264         }
 265     }
 266     {
 267         OFstream str("cutToSlaveEdges.obj");
 268         Pout<< "Writing cutToSlaveEdges to " << str.name() << endl;
 269
 270         label vertI = 0;
 271
 272         labelList slaveToCut(invert(s.nEdges(), cutToSlaveEdges));
 273
 274         forAll(slaveToCut, edgeI)
 275         {
 276             if (slaveToCut[edgeI] != -1)
 277             {
 278                 const edge& slaveEdge = s.edges()[edgeI];
 279                 const edge& cutEdge = c.edges()[slaveToCut[edgeI]];
 280
 281                 meshTools::writeOBJ(str, s.localPoints()[slaveEdge[0]]);
 282                 vertI++;
 283                 meshTools::writeOBJ(str, s.localPoints()[slaveEdge[1]]);
 284                 vertI++;
 285                 meshTools::writeOBJ(str, c.localPoints()[cutEdge[0]]);
 286                 vertI++;
 287                 meshTools::writeOBJ(str, c.localPoints()[cutEdge[1]]);
 288                 vertI++;
 289                 str << "l " << vertI-3 << ' ' << vertI-2 << nl;
 290                 str << "l " << vertI-3 << ' ' << vertI-1 << nl;
 291                 str << "l " << vertI-3 << ' ' << vertI << nl;
 292                 str << "l " << vertI-2 << ' ' << vertI-1 << nl;
 293                 str << "l " << vertI-2 << ' ' << vertI << nl;
 294                 str << "l " << vertI-1 << ' ' << vertI << nl;
 295             }
 296         }
 297     }
 298
 299     Pout<< endl;
 300 }
 301
 302
 303 // Given an edgelist and a map for the points on the edges it tries to find
 304 // the corresponding patch edges.
 305 Foam::labelList Foam::faceCoupleInfo::findMappedEdges
 306 (
 307     const edgeList& edges,
 308     const labelList& pointMap,
 309     const indirectPrimitivePatch& patch
 310 )
 311 {
 312     labelList toPatchEdges(edges.size());
 313
 314     forAll(toPatchEdges, edgeI)
 315     {
 316         const edge& e = edges[edgeI];
 317
 318         label v0 = pointMap[e[0]];
 319         label v1 = pointMap[e[1]];
 320
 321         toPatchEdges[edgeI] =
 322             meshTools::findEdge
 323             (
 324                 patch.edges(),
 325                 patch.pointEdges()[v0],
 326                 v0,
 327                 v1
 328             );
 329     }
 330     return toPatchEdges;
 331 }
 332
 333
 334 // Detect a cut edge which originates from two boundary faces having different
 335 // polyPatches.
 336 bool Foam::faceCoupleInfo::regionEdge
 337 (
 338     const polyMesh& slaveMesh,
 339     const label slaveEdgeI
 340 ) const
 341 {
 342     const labelList& eFaces = slavePatch().edgeFaces()[slaveEdgeI];
 343
 344     if (eFaces.size() == 1)
 345     {
 346         return true;
 347     }
 348     else
 349     {
 350         // Count how many different patches connected to this edge.
 351
 352         label patch0 = -1;
 353
 354         forAll(eFaces, i)
 355         {
 356             label faceI = eFaces[i];
 357
 358             label meshFaceI = slavePatch().addressing()[faceI];
 359
 360             label patchI = slaveMesh.boundaryMesh().whichPatch(meshFaceI);
 361
 362             if (patch0 == -1)
 363             {
 364                 patch0 = patchI;
 365             }
 366             else if (patchI != patch0)
 367             {
 368                 // Found two different patches connected to this edge.
 369                 return true;
 370             }
 371         }
 372         return false;
 373     }
 374 }
 375
 376
 377 // Find edge using pointI that is most aligned with vector between
 378 // master points. Patchdivision tells us whether or not to use
 379 // patch information to match edges.
 380 Foam::label Foam::faceCoupleInfo::mostAlignedCutEdge
 381 (
 382     const bool report,
 383     const polyMesh& slaveMesh,
 384     const bool patchDivision,
 385     const labelList& cutToMasterEdges,
 386     const labelList& cutToSlaveEdges,
 387     const label pointI,
 388     const label edgeStart,
 389     const label edgeEnd
 390 ) const
 391 {
 392     const pointField& localPoints = cutFaces().localPoints();
 393
 394     const labelList& pEdges = cutFaces().pointEdges()[pointI];
 395
 396     if (report)
 397     {
 398         Pout<< "mostAlignedEdge : finding nearest edge among "
 399             << UIndirectList<edge>(cutFaces().edges(), pEdges)()
 400             << " connected to point " << pointI
 401             << " coord:" << localPoints[pointI]
 402             << " running between " << edgeStart << " coord:"
 403             << localPoints[edgeStart]
 404             << " and " << edgeEnd << " coord:"
 405             << localPoints[edgeEnd]
 406             << endl;
 407     }
 408
 409     // Find the edge that gets us nearest end.
 410
 411     label maxEdgeI = -1;
 412     scalar maxCos = -GREAT;
 413
 414     forAll(pEdges, i)
 415     {
 416         label edgeI = pEdges[i];
 417
 418         if
 419         (
 420            !(
 421                 patchDivision
 422              && cutToMasterEdges[edgeI] == -1
 423             )
 424          || (
 425                 patchDivision
 426              && regionEdge(slaveMesh, cutToSlaveEdges[edgeI])
 427             )
 428         )
 429         {
 430             const edge& e = cutFaces().edges()[edgeI];
 431
 432             label otherPointI = e.otherVertex(pointI);
 433
 434             if (otherPointI == edgeEnd)
 435             {
 436                 // Shortcut: found edge end point.
 437                 if (report)
 438                 {
 439                     Pout<< "    mostAlignedEdge : found end point " << edgeEnd
 440                         << endl;
 441                 }
 442                 return edgeI;
 443             }
 444
 445             // Get angle between edge and edge to masterEnd
 446
 447             vector eVec(localPoints[otherPointI] - localPoints[pointI]);
 448
 449             scalar magEVec = mag(eVec);
 450
 451             if (magEVec < VSMALL)
 452             {
 453                 WarningIn("faceCoupleInfo::mostAlignedEdge")
 454                     << "Crossing zero sized edge " << edgeI
 455                     << " coords:" << localPoints[otherPointI]
 456                     << localPoints[pointI]
 457                     << " when walking from " << localPoints[edgeStart]
 458                     << " to " << localPoints[edgeEnd]
 459                     << endl;
 460                 return edgeI;
 461             }
 462
 463             eVec /= magEVec;
 464
 465             vector eToEndPoint(localPoints[edgeEnd] - localPoints[otherPointI]);
 466             eToEndPoint /= mag(eToEndPoint);
 467
 468             scalar cosAngle = eVec & eToEndPoint;
 469
 470             if (report)
 471             {
 472                 Pout<< "    edge:" << e << " points:" << localPoints[pointI]
 473                     << localPoints[otherPointI]
 474                     << "  vec:" << eVec
 475                     << "  vecToEnd:" << eToEndPoint
 476                     << " cosAngle:" << cosAngle
 477                     << endl;
 478             }
 479
 480             if (cosAngle > maxCos)
 481             {
 482                 maxCos = cosAngle;
 483                 maxEdgeI = edgeI;
 484             }
 485         }
 486     }
 487
 488     if (maxCos > 1 - angleTol_)
 489     {
 490         return maxEdgeI;
 491     }
 492     else
 493     {
 494         return -1;
 495     }
 496 }
 497
 498
 499 // Construct points to split points map (in cut addressing)
 500 void Foam::faceCoupleInfo::setCutEdgeToPoints(const labelList& cutToMasterEdges)
 501 {
 502     labelListList masterToCutEdges
 503     (
 504         invertOneToMany
 505         (
 506             masterPatch().nEdges(),
 507             cutToMasterEdges
 508         )
 509     );
 510
 511     const edgeList& cutEdges = cutFaces().edges();
 512
 513     // Size extra big so searching is faster
 514     cutEdgeToPoints_.resize
 515     (
 516         masterPatch().nEdges()
 517       + slavePatch().nEdges()
 518       + cutEdges.size()
 519     );
 520
 521     forAll(masterToCutEdges, masterEdgeI)
 522     {
 523         const edge& masterE = masterPatch().edges()[masterEdgeI];
 524
 525         //Pout<< "Master:" << masterPatch().localPoints()[masterE[0]] << ' '
 526         //    << masterPatch().localPoints()[masterE[1]] << endl;
 527
 528         const labelList& stringedEdges = masterToCutEdges[masterEdgeI];
 529
 530         if (stringedEdges.empty())
 531         {
 532             FatalErrorIn
 533             (
 534                 "faceCoupleInfo::setCutEdgeToPoints"
 535                 "(const labelList&)"
 536             )   << "Did not match all of master edges to cutFace edges"
 537                 << nl
 538                 << "First unmatched edge:" << masterEdgeI << " endPoints:"
 539                 << masterPatch().localPoints()[masterE[0]]
 540                 << masterPatch().localPoints()[masterE[1]]
 541                 << endl
 542                 << "This usually means that the slave patch is not a"
 543                 << " subdivision of the master patch"
 544                 << abort(FatalError);
 545         }
 546         else if (stringedEdges.size() > 1)
 547         {
 548             // String up the edges between e[0] and e[1]. Store the points
 549             // inbetween e[0] and e[1] (all in cutFaces() labels)
 550
 551             DynamicList<label> splitPoints(stringedEdges.size()-1);
 552
 553             // Unsplit edge endpoints
 554             const edge unsplitEdge
 555             (
 556                 masterToCutPoints_[masterE[0]],
 557                 masterToCutPoints_[masterE[1]]
 558             );
 559
 560             label startVertI = unsplitEdge[0];
 561             label startEdgeI = -1;
 562
 563             while (startVertI != unsplitEdge[1])
 564             {
 565                 // Loop over all string of edges. Update
 566                 // - startVertI : previous vertex
 567                 // - startEdgeI : previous edge
 568                 // and insert any points into splitPoints
 569
 570                 // For checking
 571                 label oldStart = startVertI;
 572
 573                 forAll(stringedEdges, i)
 574                 {
 575                     label edgeI = stringedEdges[i];
 576
 577                     if (edgeI != startEdgeI)
 578                     {
 579                         const edge& e = cutEdges[edgeI];
 580
 581                         //Pout<< "    cut:" << e << " at:"
 582                         //    << cutFaces().localPoints()[e[0]]
 583                         //    << ' ' << cutFaces().localPoints()[e[1]] << endl;
 584
 585                         if (e[0] == startVertI)
 586                         {
 587                             startEdgeI = edgeI;
 588                             startVertI = e[1];
 589                             if (e[1] != unsplitEdge[1])
 590                             {
 591                                 splitPoints.append(e[1]);
 592                             }
 593                             break;
 594                         }
 595                         else if (e[1] == startVertI)
 596                         {
 597                             startEdgeI = edgeI;
 598                             startVertI = e[0];
 599                             if (e[0] != unsplitEdge[1])
 600                             {
 601                                 splitPoints.append(e[0]);
 602                             }
 603                             break;
 604                         }
 605                     }
 606                 }
 607
 608                 // Check
 609                 if (oldStart == startVertI)
 610                 {
 611                     FatalErrorIn
 612                     (
 613                         "faceCoupleInfo::setCutEdgeToPoints"
 614                         "(const labelList&)"
 615                     )   << " unsplitEdge:" << unsplitEdge
 616                         << " does not correspond to split edges "
 617                         << UIndirectList<edge>(cutEdges, stringedEdges)()
 618                         << abort(FatalError);
 619                 }
 620             }
 621
 622             //Pout<< "For master edge:"
 623             //    << unsplitEdge
 624             //    << " Found stringed points "
 625             //    <<  UIndirectList<point>
 626             //        (
 627             //            cutFaces().localPoints(),
 628             //            splitPoints.shrink()
 629             //        )()
 630             //    << endl;
 631
 632             cutEdgeToPoints_.insert(unsplitEdge, splitPoints.shrink());
 633         }
 634     }
 635 }
 636
 637
 638 // Determines rotation for f1 to match up with f0, i.e. the index in f0 of
 639 // the first point of f1.
 640 Foam::label Foam::faceCoupleInfo::matchFaces
 641 (
 642     const scalar absTol,
 643     const pointField& points0,
 644     const face& f0,
 645     const pointField& points1,
 646     const face& f1,
 647     const bool sameOrientation
 648 )
 649 {
 650     if (f0.size() != f1.size())
 651     {
 652         FatalErrorIn
 653         (
 654             "faceCoupleInfo::matchFaces"
 655             "(const scalar, const face&, const pointField&"
 656             ", const face&, const pointField&)"
 657         )   << "Different sizes for supposedly matching faces." << nl
 658             << "f0:" << f0 << " coords:" << UIndirectList<point>(points0, f0)()
 659             << nl
 660             << "f1:" << f1 << " coords:" << UIndirectList<point>(points1, f1)()
 661             << abort(FatalError);
 662     }
 663
 664     const scalar absTolSqr = sqr(absTol);
 665
 666
 667     label matchFp = -1;
 668
 669     forAll(f0, startFp)
 670     {
 671         // See -if starting from startFp on f0- the two faces match.
 672         bool fullMatch = true;
 673
 674         label fp0 = startFp;
 675         label fp1 = 0;
 676
 677         forAll(f1, i)
 678         {
 679             scalar distSqr = Foam::magSqr(points0[f0[fp0]] - points1[f1[fp1]]);
 680
 681             if (distSqr > absTolSqr)
 682             {
 683                 fullMatch = false;
 684                 break;
 685             }
 686
 687             fp0 = f0.fcIndex(fp0);  // walk forward
 688
 689             if (sameOrientation)
 690             {
 691                 fp1 = f1.fcIndex(fp1);
 692             }
 693             else
 694             {
 695                 fp1 = f1.rcIndex(fp1);
 696             }
 697         }
 698
 699         if (fullMatch)
 700         {
 701             matchFp = startFp;
 702             break;
 703         }
 704     }
 705
 706     if (matchFp == -1)
 707     {
 708         FatalErrorIn
 709         (
 710             "faceCoupleInfo::matchFaces"
 711             "(const scalar, const face&, const pointField&"
 712             ", const face&, const pointField&)"
 713         )   << "No unique match between two faces" << nl
 714             << "Face " << f0 << " coords "
 715             << UIndirectList<point>(points0, f0)() << nl
 716             << "Face " << f1 << " coords "
 717             << UIndirectList<point>(points1, f1)()
 718             << "when using tolerance " << absTol
 719             << " and forwardMatching:" << sameOrientation
 720             << abort(FatalError);
 721     }
 722
 723     return matchFp;
 724 }
 725
 726
 727 // Find correspondence from patch points to cut points. This might
 728 // detect shared points so the output is a patch-to-cut point list
 729 // and a compaction list for the cut points (which will always be equal or more
 730 // connected than the patch).
 731 // Returns true if there are any duplicates.
 732 bool Foam::faceCoupleInfo::matchPointsThroughFaces
 733 (
 734     const scalar absTol,
 735     const pointField& cutPoints,
 736     const faceList& cutFaces,
 737     const pointField& patchPoints,
 738     const faceList& patchFaces,
 739     const bool sameOrientation,
 740
 741     labelList& patchToCutPoints,    // patch to (uncompacted) cut points
 742     labelList& cutToCompact,        // compaction list for cut points
 743     labelList& compactToCut         // inverse ,,
 744 )
 745 {
 746
 747     // From slave to cut point
 748     patchToCutPoints.setSize(patchPoints.size());
 749     patchToCutPoints = -1;
 750
 751     // Compaction list for cut points: either -1 or index into master which
 752     // gives the point to compact to.
 753     labelList cutPointRegion(cutPoints.size(), -1);
 754     DynamicList<label> cutPointRegionMaster;
 755
 756     forAll(patchFaces, patchFaceI)
 757     {
 758         const face& patchF = patchFaces[patchFaceI];
 759
 760         //const face& cutF = cutFaces[patchToCutFaces[patchFaceI]];
 761         const face& cutF = cutFaces[patchFaceI];
 762
 763         // Do geometric matching to get position of cutF[0] in patchF
 764         label patchFp = matchFaces
 765         (
 766             absTol,
 767             patchPoints,
 768             patchF,
 769             cutPoints,
 770             cutF,
 771             sameOrientation        // orientation
 772         );
 773
 774         forAll(cutF, cutFp)
 775         {
 776             label cutPointI = cutF[cutFp];
 777             label patchPointI = patchF[patchFp];
 778
 779             //const point& cutPt = cutPoints[cutPointI];
 780             //const point& patchPt = patchPoints[patchPointI];
 781             //if (mag(cutPt - patchPt) > SMALL)
 782             //{
 783             //    FatalErrorIn("matchPointsThroughFaces")
 784             //    << "cutP:" << cutPt
 785             //    << " patchP:" << patchPt
 786             //    << abort(FatalError);
 787             //}
 788
 789             if (patchToCutPoints[patchPointI] == -1)
 790             {
 791                 patchToCutPoints[patchPointI] = cutPointI;
 792             }
 793             else if (patchToCutPoints[patchPointI] != cutPointI)
 794             {
 795                 // Multiple cut points connecting to same patch.
 796                 // Check if already have region & region master for this set
 797                 label otherCutPointI = patchToCutPoints[patchPointI];
 798
 799                 //Pout<< "PatchPoint:" << patchPt
 800                 //    << " matches to:" << cutPointI
 801                 //    << " coord:" << cutPoints[cutPointI]
 802                 //    << " and to:" << otherCutPointI
 803                 //    << " coord:" << cutPoints[otherCutPointI]
 804                 //    << endl;
 805
 806                 if (cutPointRegion[otherCutPointI] != -1)
 807                 {
 808                     // Have region for this set. Copy.
 809                     label region = cutPointRegion[otherCutPointI];
 810                     cutPointRegion[cutPointI] = region;
 811
 812                     // Update region master with min point label
 813                     cutPointRegionMaster[region] = min
 814                     (
 815                         cutPointRegionMaster[region],
 816                         cutPointI
 817                     );
 818                 }
 819                 else
 820                 {
 821                     // Create new region.
 822                     label region = cutPointRegionMaster.size();
 823                     cutPointRegionMaster.append
 824                     (
 825                         min(cutPointI, otherCutPointI)
 826                     );
 827                     cutPointRegion[cutPointI] = region;
 828                     cutPointRegion[otherCutPointI] = region;
 829                 }
 830             }
 831
 832             if (sameOrientation)
 833             {
 834                 patchFp = patchF.fcIndex(patchFp);
 835             }
 836             else
 837             {
 838                 patchFp = patchF.rcIndex(patchFp);
 839             }
 840         }
 841     }
 842
 843     // Rework region&master into compaction array
 844     compactToCut.setSize(cutPointRegion.size());
 845     cutToCompact.setSize(cutPointRegion.size());
 846     cutToCompact = -1;
 847     label compactPointI = 0;
 848
 849     forAll(cutPointRegion, i)
 850     {
 851         if (cutPointRegion[i] == -1)
 852         {
 853             // Unduplicated point. Allocate new compacted point.
 854             cutToCompact[i] = compactPointI;
 855             compactToCut[compactPointI] = i;
 856             compactPointI++;
 857         }
 858         else
 859         {
 860             // Duplicate point. Get master.
 861
 862             label masterPointI = cutPointRegionMaster[cutPointRegion[i]];
 863
 864             if (cutToCompact[masterPointI] == -1)
 865             {
 866                 cutToCompact[masterPointI] = compactPointI;
 867                 compactToCut[compactPointI] = masterPointI;
 868                 compactPointI++;
 869             }
 870             cutToCompact[i] = cutToCompact[masterPointI];
 871         }
 872     }
 873     compactToCut.setSize(compactPointI);
 874
 875     return compactToCut.size() != cutToCompact.size();
 876 }
 877
 878
 879 // Return max distance from any point on cutF to masterF
 880 Foam::scalar Foam::faceCoupleInfo::maxDistance
 881 (
 882     const face& cutF,
 883     const pointField& cutPoints,
 884     const face& masterF,
 885     const pointField& masterPoints
 886 )
 887 {
 888     scalar maxDist = -GREAT;
 889
 890     forAll(cutF, fp)
 891     {
 892         const point& cutPt = cutPoints[cutF[fp]];
 893
 894         pointHit pHit = masterF.nearestPoint(cutPt, masterPoints);
 895
 896         maxDist = max(maxDist, pHit.distance());
 897     }
 898     return maxDist;
 899 }
 900
 901
 902 void Foam::faceCoupleInfo::findPerfectMatchingFaces
 903 (
 904     const primitiveMesh& mesh0,
 905     const primitiveMesh& mesh1,
 906     const scalar absTol,
 907
 908     labelList& mesh0Faces,
 909     labelList& mesh1Faces
 910 )
 911 {
 912     // Face centres of external faces (without invoking
 913     // mesh.faceCentres since mesh might have been clearedOut)
 914
 915     pointField fc0
 916     (
 917         calcFaceCentres<List>
 918         (
 919             mesh0.faces(),
 920             mesh0.points(),
 921             mesh0.nInternalFaces(),
 922             mesh0.nFaces() - mesh0.nInternalFaces()
 923         )
 924     );
 925
 926     pointField fc1
 927     (
 928         calcFaceCentres<List>
 929         (
 930             mesh1.faces(),
 931             mesh1.points(),
 932             mesh1.nInternalFaces(),
 933             mesh1.nFaces() - mesh1.nInternalFaces()
 934         )
 935     );
 936
 937
 938     if (debug)
 939     {
 940         Pout<< "Face matching tolerance : " << absTol << endl;
 941     }
 942
 943
 944     // Match geometrically
 945     labelList from1To0;
 946     bool matchedAllFaces = matchPoints
 947     (
 948         fc1,
 949         fc0,
 950         scalarField(fc1.size(), absTol),
 951         false,
 952         from1To0
 953     );
 954
 955     if (matchedAllFaces)
 956     {
 957         Warning
 958             << "faceCoupleInfo::faceCoupleInfo : "
 959             << "Matched ALL " << fc1.size()
 960             << " boundary faces of mesh0 to boundary faces of mesh1." << endl
 961             << "This is only valid if the mesh to add is fully"
 962             << " enclosed by the mesh it is added to." << endl;
 963     }
 964
 965
 966     // Collect matches.
 967     label nMatched = 0;
 968
 969     mesh0Faces.setSize(fc0.size());
 970     mesh1Faces.setSize(fc1.size());
 971
 972     forAll(from1To0, i)
 973     {
 974         if (from1To0[i] != -1)
 975         {
 976             mesh1Faces[nMatched] = i + mesh1.nInternalFaces();
 977             mesh0Faces[nMatched] = from1To0[i] + mesh0.nInternalFaces();
 978
 979             nMatched++;
 980         }
 981     }
 982
 983     mesh0Faces.setSize(nMatched);
 984     mesh1Faces.setSize(nMatched);
 985 }
 986
 987
 988 void Foam::faceCoupleInfo::findSlavesCoveringMaster
 989 (
 990     const primitiveMesh& mesh0,
 991     const primitiveMesh& mesh1,
 992     const scalar absTol,
 993
 994     labelList& mesh0Faces,
 995     labelList& mesh1Faces
 996 )
 997 {
 998     // Construct octree from all mesh0 boundary faces
 999     labelList bndFaces(mesh0.nFaces()-mesh0.nInternalFaces());
1000     forAll(bndFaces, i)
1001     {
1002         bndFaces[i] = mesh0.nInternalFaces() + i;
1003     }
1004
1005     treeBoundBox overallBb(mesh0.points());
1006
1007     Random rndGen(123456);
1008
1009     indexedOctree<treeDataFace> tree
1010     (
1011         treeDataFace    // all information needed to search faces
1012         (
1013             false,                      // do not cache bb
1014             mesh0,
1015             bndFaces                    // boundary faces only
1016         ),
1017         overallBb.extend(rndGen, 1E-4), // overall search domain
1018         8,                              // maxLevel
1019         10,                             // leafsize
1020         3.0                             // duplicity
1021     );
1022
1023     if (debug)
1024     {
1025         Pout<< "findSlavesCoveringMaster :"
1026             << " constructed octree for mesh0 boundary faces" << endl;
1027     }
1028
1029     // Search nearest face for every mesh1 boundary face.
1030
1031     labelHashSet mesh0Set(mesh0.nFaces() - mesh0.nInternalFaces());
1032     labelHashSet mesh1Set(mesh1.nFaces() - mesh1.nInternalFaces());
1033
1034     for
1035     (
1036         label mesh1FaceI = mesh1.nInternalFaces();
1037         mesh1FaceI < mesh1.nFaces();
1038         mesh1FaceI++
1039     )
1040     {
1041         const face& f1 = mesh1.faces()[mesh1FaceI];
1042
1043         // Generate face centre (prevent cellCentres() reconstruction)
1044         point fc(f1.centre(mesh1.points()));
1045
1046         pointIndexHit nearInfo = tree.findNearest(fc, Foam::sqr(absTol));
1047
1048         if (nearInfo.hit())
1049         {
1050             label mesh0FaceI = tree.shapes().faceLabels()[nearInfo.index()];
1051
1052             // Check if points of f1 actually lie on top of mesh0 face
1053             // This is the bit that might fail since if f0 is severely warped
1054             // and f1's points are not present on f0 (since subdivision)
1055             // then the distance of the points to f0 might be quite large
1056             // and the test will fail. This all should in fact be some kind
1057             // of walk from known corresponding points/edges/faces.
1058             scalar dist =
1059                 maxDistance
1060                 (
1061                     f1,
1062                     mesh1.points(),
1063                     mesh0.faces()[mesh0FaceI],
1064                     mesh0.points()
1065                 );
1066
1067             if (dist < absTol)
1068             {
1069                 mesh0Set.insert(mesh0FaceI);
1070                 mesh1Set.insert(mesh1FaceI);
1071             }
1072         }
1073     }
1074
1075     if (debug)
1076     {
1077         Pout<< "findSlavesCoveringMaster :"
1078             << " matched " << mesh1Set.size() << " mesh1 faces to "
1079             << mesh0Set.size() << " mesh0 faces" << endl;
1080     }
1081
1082     mesh0Faces = mesh0Set.toc();
1083     mesh1Faces = mesh1Set.toc();
1084 }
1085
1086
1087 // Grow cutToMasterFace across 'internal' edges.
1088 Foam::label Foam::faceCoupleInfo::growCutFaces
1089 (
1090     const labelList& cutToMasterEdges,
1091     Map<labelList>& candidates
1092 )
1093 {
1094     if (debug)
1095     {
1096         Pout<< "growCutFaces :"
1097             << " growing cut faces to masterPatch" << endl;
1098     }
1099
1100     label nTotChanged = 0;
1101
1102     while (true)
1103     {
1104         const labelListList& cutFaceEdges = cutFaces().faceEdges();
1105         const labelListList& cutEdgeFaces = cutFaces().edgeFaces();
1106
1107         label nChanged = 0;
1108
1109         forAll(cutToMasterFaces_, cutFaceI)
1110         {
1111             const label masterFaceI = cutToMasterFaces_[cutFaceI];
1112
1113             if (masterFaceI != -1)
1114             {
1115                 // We now have a cutFace for which we have already found a
1116                 // master face. Grow this masterface across any internal edge
1117                 // (internal: no corresponding master edge)
1118
1119                 const labelList& fEdges = cutFaceEdges[cutFaceI];
1120
1121                 forAll(fEdges, i)
1122                 {
1123                     const label cutEdgeI = fEdges[i];
1124
1125                     if (cutToMasterEdges[cutEdgeI] == -1)
1126                     {
1127                         // So this edge:
1128                         // - internal to the cutPatch (since all region edges
1129                         //   marked before)
1130                         // - on cutFaceI which corresponds to masterFace.
1131                         // Mark all connected faces with this masterFace.
1132
1133                         const labelList& eFaces = cutEdgeFaces[cutEdgeI];
1134
1135                         forAll(eFaces, j)
1136                         {
1137                             const label faceI = eFaces[j];
1138
1139                             if (cutToMasterFaces_[faceI] == -1)
1140                             {
1141                                 cutToMasterFaces_[faceI] = masterFaceI;
1142                                 candidates.erase(faceI);
1143                                 nChanged++;
1144                             }
1145                             else if (cutToMasterFaces_[faceI] != masterFaceI)
1146                             {
1147                                 const pointField& cutPoints =
1148                                     cutFaces().points();
1149                                 const pointField& masterPoints =
1150                                     masterPatch().points();
1151
1152                                 const edge& e = cutFaces().edges()[cutEdgeI];
1153
1154                                 label myMaster = cutToMasterFaces_[faceI];
1155                                 const face& myF = masterPatch()[myMaster];
1156
1157                                 const face& nbrF = masterPatch()[masterFaceI];
1158
1159                                 FatalErrorIn
1160                                 (
1161                                     "faceCoupleInfo::growCutFaces"
1162                                     "(const labelList&, Map<labelList>&)"
1163                                 )   << "Cut face "
1164                                     << cutFaces()[faceI].points(cutPoints)
1165                                     << " has master "
1166                                     << myMaster
1167                                     << " but also connects to nbr face "
1168                                     << cutFaces()[cutFaceI].points(cutPoints)
1169                                     << " with master " << masterFaceI
1170                                     << nl
1171                                     << "myMasterFace:"
1172                                     << myF.points(masterPoints)
1173                                     << "  nbrMasterFace:"
1174                                     << nbrF.points(masterPoints) << nl
1175                                     << "Across cut edge " << cutEdgeI
1176                                     << " coord:"
1177                                     << cutFaces().localPoints()[e[0]]
1178                                     << cutFaces().localPoints()[e[1]]
1179                                     << abort(FatalError);
1180                             }
1181                         }
1182                     }
1183                 }
1184             }
1185         }
1186
1187         if (debug)
1188         {
1189             Pout<< "growCutFaces : Grown an additional " << nChanged
1190                 << " cut to master face" << " correspondence" << endl;
1191         }
1192
1193         nTotChanged += nChanged;
1194
1195         if (nChanged == 0)
1196         {
1197             break;
1198         }
1199     }
1200
1201     return nTotChanged;
1202 }
1203
1204
1205 void Foam::faceCoupleInfo::checkMatch(const labelList& cutToMasterEdges) const
1206 {
1207     const pointField& cutLocalPoints = cutFaces().localPoints();
1208
1209     const pointField& masterLocalPoints = masterPatch().localPoints();
1210     const faceList& masterLocalFaces = masterPatch().localFaces();
1211
1212     forAll(cutToMasterEdges, cutEdgeI)
1213     {
1214         const edge& e = cutFaces().edges()[cutEdgeI];
1215
1216         if (cutToMasterEdges[cutEdgeI] == -1)
1217         {
1218             // Internal edge. Check that master face is same on both sides.
1219             const labelList& cutEFaces = cutFaces().edgeFaces()[cutEdgeI];
1220
1221             label masterFaceI = -1;
1222
1223             forAll(cutEFaces, i)
1224             {
1225                 label cutFaceI = cutEFaces[i];
1226
1227                 if (cutToMasterFaces_[cutFaceI] != -1)
1228                 {
1229                     if (masterFaceI == -1)
1230                     {
1231                         masterFaceI = cutToMasterFaces_[cutFaceI];
1232                     }
1233                     else if (masterFaceI != cutToMasterFaces_[cutFaceI])
1234                     {
1235                         label myMaster = cutToMasterFaces_[cutFaceI];
1236                         const face& myF = masterLocalFaces[myMaster];
1237
1238                         const face& nbrF = masterLocalFaces[masterFaceI];
1239
1240                         FatalErrorIn
1241                         (
1242                             "faceCoupleInfo::checkMatch(const labelList&) const"
1243                         )
1244                             << "Internal CutEdge " << e
1245                             << " coord:"
1246                             << cutLocalPoints[e[0]]
1247                             << cutLocalPoints[e[1]]
1248                             << " connects to master " << myMaster
1249                             << " and to master " << masterFaceI << nl
1250                             << "myMasterFace:"
1251                             << myF.points(masterLocalPoints)
1252                             << "  nbrMasterFace:"
1253                             << nbrF.points(masterLocalPoints)
1254                             << abort(FatalError);
1255                     }
1256                 }
1257             }
1258         }
1259     }
1260 }
1261
1262
1263 // Extends matching information by elimination across cutFaces using more
1264 // than one region edge. Updates cutToMasterFaces_ and sets candidates
1265 // which is for every cutface on a region edge the possible master faces.
1266 Foam::label Foam::faceCoupleInfo::matchEdgeFaces
1267 (
1268     const labelList& cutToMasterEdges,
1269     Map<labelList>& candidates
1270 )
1271 {
1272     // For every unassigned cutFaceI the possible list of master faces.
1273     candidates.clear();
1274     candidates.resize(cutFaces().size());
1275
1276     label nChanged = 0;
1277
1278     forAll(cutToMasterEdges, cutEdgeI)
1279     {
1280         label masterEdgeI = cutToMasterEdges[cutEdgeI];
1281
1282         if (masterEdgeI != -1)
1283         {
1284             // So cutEdgeI is matched to masterEdgeI. For all cut faces using
1285             // the cut edge store the master face as a candidate match.
1286
1287             const labelList& cutEFaces = cutFaces().edgeFaces()[cutEdgeI];
1288             const labelList& masterEFaces =
1289                 masterPatch().edgeFaces()[masterEdgeI];
1290
1291             forAll(cutEFaces, i)
1292             {
1293                 label cutFaceI = cutEFaces[i];
1294
1295                 if (cutToMasterFaces_[cutFaceI] == -1)
1296                 {
1297                     // So this face (cutFaceI) is on an edge (cutEdgeI) that
1298                     // is used by some master faces (masterEFaces). Check if
1299                     // the cutFaceI and some of these masterEFaces share more
1300                     // than one edge (which uniquely defines face)
1301
1302                     // Combine master faces with current set of candidate
1303                     // master faces.
1304                     Map<labelList>::iterator fnd = candidates.find(cutFaceI);
1305
1306                     if (fnd == candidates.end())
1307                     {
1308                         // No info yet for cutFaceI. Add all master faces as
1309                         // candidates
1310                         candidates.insert(cutFaceI, masterEFaces);
1311                     }
1312                     else
1313                     {
1314                         // From some other cutEdgeI there are already some
1315                         // candidate master faces. Check the overlap with
1316                         // the current set of master faces.
1317                         const labelList& masterFaces = fnd();
1318
1319                         DynamicList<label> newCandidates(masterFaces.size());
1320
1321                         forAll(masterEFaces, j)
1322                         {
1323                             if (findIndex(masterFaces, masterEFaces[j]) != -1)
1324                             {
1325                                 newCandidates.append(masterEFaces[j]);
1326                             }
1327                         }
1328
1329                         if (newCandidates.size() == 1)
1330                         {
1331                             // We found a perfect match. Delete entry from
1332                             // candidates map.
1333                             cutToMasterFaces_[cutFaceI] = newCandidates[0];
1334                             candidates.erase(cutFaceI);
1335                             nChanged++;
1336                         }
1337                         else
1338                         {
1339                             // Should not happen?
1340                             //Pout<< "On edge:" << cutEdgeI
1341                             //    << " have connected masterFaces:"
1342                             //    << masterEFaces
1343                             //    << " and from previous edge we have"
1344                             //    << " connected masterFaces" << masterFaces
1345                             //    << " . Overlap:" << newCandidates << endl;
1346
1347                             fnd() = newCandidates.shrink();
1348                         }
1349                     }
1350                 }
1351
1352             }
1353         }
1354     }
1355
1356     if (debug)
1357     {
1358         Pout<< "matchEdgeFaces : Found " << nChanged
1359             << " faces where there was"
1360             << " only one remaining choice for cut-master correspondence"
1361             << endl;
1362     }
1363
1364     return nChanged;
1365 }
1366
1367
1368 // Gets a list of cutFaces (that use a master edge) and the candidate
1369 // master faces.
1370 // Finds most aligned master face.
1371 Foam::label Foam::faceCoupleInfo::geometricMatchEdgeFaces
1372 (
1373     Map<labelList>& candidates
1374 )
1375 {
1376     const pointField& cutPoints = cutFaces().points();
1377
1378     label nChanged = 0;
1379
1380     // Mark all master faces that have been matched so far.
1381
1382     labelListList masterToCutFaces
1383     (
1384         invertOneToMany
1385         (
1386             masterPatch().size(),
1387             cutToMasterFaces_
1388         )
1389     );
1390
1391     forAllConstIter(Map<labelList>, candidates, iter)
1392     {
1393         label cutFaceI = iter.key();
1394
1395         const face& cutF = cutFaces()[cutFaceI];
1396
1397         if (cutToMasterFaces_[cutFaceI] == -1)
1398         {
1399             const labelList& masterFaces = iter();
1400
1401             // Find the best matching master face.
1402             scalar minDist = GREAT;
1403             label minMasterFaceI = -1;
1404
1405             forAll(masterFaces, i)
1406             {
1407                 label masterFaceI = masterFaces[i];
1408
1409                 if (masterToCutFaces[masterFaces[i]].empty())
1410                 {
1411                     scalar dist = maxDistance
1412                     (
1413                         cutF,
1414                         cutPoints,
1415                         masterPatch()[masterFaceI],
1416                         masterPatch().points()
1417                     );
1418
1419                     if (dist < minDist)
1420                     {
1421                         minDist = dist;
1422                         minMasterFaceI = masterFaceI;
1423                     }
1424                 }
1425             }
1426
1427             if (minMasterFaceI != -1)
1428             {
1429                 cutToMasterFaces_[cutFaceI] = minMasterFaceI;
1430                 masterToCutFaces[minMasterFaceI] = cutFaceI;
1431                 nChanged++;
1432             }
1433         }
1434     }
1435
1436     // (inefficiently) force candidates to be uptodate.
1437     forAll(cutToMasterFaces_, cutFaceI)
1438     {
1439         if (cutToMasterFaces_[cutFaceI] != -1)
1440         {
1441             candidates.erase(cutFaceI);
1442         }
1443     }
1444
1445
1446     if (debug)
1447     {
1448         Pout<< "geometricMatchEdgeFaces : Found " << nChanged
1449             << " faces where there was"
1450             << " only one remaining choice for cut-master correspondence"
1451             << endl;
1452     }
1453
1454     return nChanged;
1455 }
1456
1457
1458 // Calculate the set of cut faces inbetween master and slave patch
1459 // assuming perfect match (and optional face ordering on slave)
1460 void Foam::faceCoupleInfo::perfectPointMatch
1461 (
1462     const scalar absTol,
1463     const bool slaveFacesOrdered
1464 )
1465 {
1466     if (debug)
1467     {
1468         Pout<< "perfectPointMatch :"
1469             << " Matching master and slave to cut."
1470             << " Master and slave faces are identical" << nl;
1471
1472         if (slaveFacesOrdered)
1473         {
1474             Pout<< "and master and slave faces are ordered"
1475                 << " (on coupled patches)" << endl;
1476         }
1477         else
1478         {
1479             Pout<< "and master and slave faces are not ordered"
1480                 << " (on coupled patches)" << endl;
1481         }
1482     }
1483
1484     cutToMasterFaces_ = identity(masterPatch().size());
1485     cutPoints_ = masterPatch().localPoints();
1486     cutFacesPtr_.reset
1487     (
1488         new primitiveFacePatch
1489         (
1490             masterPatch().localFaces(),
1491             cutPoints_
1492         )
1493     );
1494     masterToCutPoints_ = identity(cutPoints_.size());
1495
1496
1497     // Cut faces to slave patch.
1498     bool matchedAllFaces = false;
1499
1500     if (slaveFacesOrdered)
1501     {
1502         cutToSlaveFaces_ = identity(cutFaces().size());
1503         matchedAllFaces = (cutFaces().size() == slavePatch().size());
1504     }
1505     else
1506     {
1507         // Faces do not have to be ordered (but all have
1508         // to match). Note: Faces will be already ordered if we enter here from
1509         // construct from meshes.
1510         matchedAllFaces = matchPoints
1511         (
1512             calcFaceCentres<List>
1513             (
1514                 cutFaces(),
1515                 cutPoints_,
1516                 0,
1517                 cutFaces().size()
1518             ),
1519             calcFaceCentres<IndirectList>
1520             (
1521                 slavePatch(),
1522                 slavePatch().points(),
1523                 0,
1524                 slavePatch().size()
1525             ),
1526             scalarField(slavePatch().size(), absTol),
1527             true,
1528             cutToSlaveFaces_
1529         );
1530     }
1531
1532
1533     if (!matchedAllFaces)
1534     {
1535         FatalErrorIn
1536         (
1537             "faceCoupleInfo::perfectPointMatch"
1538             "(const scalar, const bool)"
1539         )   << "Did not match all of the master faces to the slave faces"
1540             << endl
1541             << "This usually means that the slave patch and master patch"
1542             << " do not align to within " << absTol << " meter."
1543             << abort(FatalError);
1544     }
1545
1546
1547     // Find correspondence from slave points to cut points. This might
1548     // detect shared points so the output is a slave-to-cut point list
1549     // and a compaction list.
1550
1551     labelList cutToCompact, compactToCut;
1552     matchPointsThroughFaces
1553     (
1554         absTol,
1555         cutFaces().localPoints(),
1556         reorder(cutToSlaveFaces_, cutFaces().localFaces()),
1557         slavePatch().localPoints(),
1558         slavePatch().localFaces(),
1559         false,                      // slave and cut have opposite orientation
1560
1561         slaveToCutPoints_,          // slave to (uncompacted) cut points
1562         cutToCompact,               // compaction map: from cut to compacted
1563         compactToCut                // compaction map: from compacted to cut
1564     );
1565
1566
1567     // Use compaction lists to renumber cutPoints.
1568     cutPoints_ = UIndirectList<point>(cutPoints_, compactToCut)();
1569     {
1570         const faceList& cutLocalFaces = cutFaces().localFaces();
1571
1572         faceList compactFaces(cutLocalFaces.size());
1573         forAll(cutLocalFaces, i)
1574         {
1575             compactFaces[i] = renumber(cutToCompact, cutLocalFaces[i]);
1576         }
1577         cutFacesPtr_.reset
1578         (
1579             new primitiveFacePatch
1580             (
1581                 compactFaces,
1582                 cutPoints_
1583             )
1584         );
1585     }
1586     inplaceRenumber(cutToCompact, slaveToCutPoints_);
1587     inplaceRenumber(cutToCompact, masterToCutPoints_);
1588 }
1589
1590
1591 // Calculate the set of cut faces inbetween master and slave patch
1592 // assuming that slave patch is subdivision of masterPatch.
1593 void Foam::faceCoupleInfo::subDivisionMatch
1594 (
1595     const polyMesh& slaveMesh,
1596     const bool patchDivision,
1597     const scalar absTol
1598 )
1599 {
1600     if (debug)
1601     {
1602         Pout<< "subDivisionMatch :"
1603             << " Matching master and slave to cut."
1604             << " Slave can be subdivision of master but all master edges"
1605             << " have to be covered by slave edges." << endl;
1606     }
1607
1608     // Assume slave patch is subdivision of the master patch so cutFaces is a
1609     // copy of the slave (but reversed since orientation has to be like master).
1610     cutPoints_ = slavePatch().localPoints();
1611     {
1612         faceList cutFaces(slavePatch().size());
1613
1614         forAll(cutFaces, i)
1615         {
1616             cutFaces[i] = slavePatch().localFaces()[i].reverseFace();
1617         }
1618         cutFacesPtr_.reset(new primitiveFacePatch(cutFaces, cutPoints_));
1619     }
1620
1621     // Cut is copy of slave so addressing to slave is trivial.
1622     cutToSlaveFaces_ = identity(cutFaces().size());
1623     slaveToCutPoints_ = identity(slavePatch().nPoints());
1624
1625     // Cut edges to slave patch
1626     labelList cutToSlaveEdges
1627     (
1628         findMappedEdges
1629         (
1630             cutFaces().edges(),
1631             slaveToCutPoints_,  //note:should be cutToSlavePoints but since iden
1632             slavePatch()
1633         )
1634     );
1635
1636
1637     if (debug)
1638     {
1639         OFstream str("regionEdges.obj");
1640
1641         Pout<< "subDivisionMatch :"
1642             << " addressing for slave patch fully done."
1643             << " Dumping region edges to " << str.name() << endl;
1644
1645         label vertI = 0;
1646
1647         forAll(slavePatch().edges(), slaveEdgeI)
1648         {
1649             if (regionEdge(slaveMesh, slaveEdgeI))
1650             {
1651                 const edge& e = slavePatch().edges()[slaveEdgeI];
1652
1653                 meshTools::writeOBJ(str, slavePatch().localPoints()[e[0]]);
1654                 vertI++;
1655                 meshTools::writeOBJ(str, slavePatch().localPoints()[e[1]]);
1656                 vertI++;
1657                 str<< "l " << vertI-1 << ' ' << vertI << nl;
1658             }
1659         }
1660     }
1661
1662
1663     // Addressing from cut to master.
1664
1665     // Cut points to master patch
1666     // - determine master points to cut points. Has to be full!
1667     // - invert to get cut to master
1668     if (debug)
1669     {
1670         Pout<< "subDivisionMatch :"
1671             << " matching master points to cut points" << endl;
1672     }
1673
1674     bool matchedAllPoints = matchPoints
1675     (
1676         masterPatch().localPoints(),
1677         cutPoints_,
1678         scalarField(masterPatch().nPoints(), absTol),
1679         true,
1680         masterToCutPoints_
1681     );
1682
1683     if (!matchedAllPoints)
1684     {
1685         FatalErrorIn
1686         (
1687             "faceCoupleInfo::subDivisionMatch"
1688             "(const polyMesh&, const bool, const scalar)"
1689         )   << "Did not match all of the master points to the slave points"
1690             << endl
1691             << "This usually means that the slave patch is not a"
1692             << " subdivision of the master patch"
1693             << abort(FatalError);
1694     }
1695
1696
1697     // Do masterEdges to cutEdges :
1698     // - mark all edges between two masterEdge endpoints. (geometric test since
1699     //   this is the only distinction)
1700     // - this gives the borders inbetween which all cutfaces come from
1701     //   a single master face.
1702     if (debug)
1703     {
1704         Pout<< "subDivisionMatch :"
1705             << " matching cut edges to master edges" << endl;
1706     }
1707
1708     const edgeList& masterEdges = masterPatch().edges();
1709     const pointField& masterPoints = masterPatch().localPoints();
1710
1711     const edgeList& cutEdges = cutFaces().edges();
1712
1713     labelList cutToMasterEdges(cutFaces().nEdges(), -1);
1714
1715     forAll(masterEdges, masterEdgeI)
1716     {
1717         const edge& masterEdge = masterEdges[masterEdgeI];
1718
1719         label cutPoint0 = masterToCutPoints_[masterEdge[0]];
1720         label cutPoint1 = masterToCutPoints_[masterEdge[1]];
1721
1722         // Find edges between cutPoint0 and cutPoint1.
1723
1724         label cutPointI = cutPoint0;
1725         do
1726         {
1727             // Find edge (starting at pointI on cut), aligned with master
1728             // edge.
1729             label cutEdgeI =
1730                 mostAlignedCutEdge
1731                 (
1732                     false,
1733                     slaveMesh,
1734                     patchDivision,
1735                     cutToMasterEdges,
1736                     cutToSlaveEdges,
1737                     cutPointI,
1738                     cutPoint0,
1739                     cutPoint1
1740                 );
1741
1742             if (cutEdgeI == -1)
1743             {
1744                 // Problem. Report while matching to produce nice error message
1745                 mostAlignedCutEdge
1746                 (
1747                     true,
1748                     slaveMesh,
1749                     patchDivision,
1750                     cutToMasterEdges,
1751                     cutToSlaveEdges,
1752                     cutPointI,
1753                     cutPoint0,
1754                     cutPoint1
1755                 );
1756
1757                 Pout<< "Dumping unmatched pointEdges to errorEdges.obj"
1758                     << endl;
1759
1760                 writeOBJ
1761                 (
1762                     "errorEdges.obj",
1763                     edgeList
1764                     (
1765                         UIndirectList<edge>
1766                         (
1767                             cutFaces().edges(),
1768                             cutFaces().pointEdges()[cutPointI]
1769                         )
1770                     ),
1771                     cutFaces().localPoints(),
1772                     false
1773                 );
1774
1775                 FatalErrorIn
1776                 (
1777                     "faceCoupleInfo::subDivisionMatch"
1778                     "(const polyMesh&, const bool, const scalar)"
1779                 )   << "Problem in finding cut edges corresponding to"
1780                     << " master edge " << masterEdge
1781                     << " points:" << masterPoints[masterEdge[0]]
1782                     << ' ' << masterPoints[masterEdge[1]]
1783                     << " corresponding cut edge: (" << cutPoint0
1784                     << ' ' << cutPoint1
1785                     << abort(FatalError);
1786             }
1787
1788             cutToMasterEdges[cutEdgeI] = masterEdgeI;
1789
1790             cutPointI = cutEdges[cutEdgeI].otherVertex(cutPointI);
1791
1792         } while (cutPointI != cutPoint1);
1793     }
1794
1795     if (debug)
1796     {
1797         // Write all matched edges.
1798         writeEdges(cutToMasterEdges, cutToSlaveEdges);
1799     }
1800
1801     // Rework cutToMasterEdges into list of points inbetween two endpoints
1802     // (cutEdgeToPoints_)
1803     setCutEdgeToPoints(cutToMasterEdges);
1804
1805
1806     // Now that we have marked the cut edges that lie on top of master edges
1807     // we can find cut faces that have two (or more) of these edges.
1808     // Note that there can be multiple faces having two or more matched edges
1809     // since the cut faces can form a non-manifold surface(?)
1810     // So the matching is done as an elimination process where for every
1811     // cutFace on a matched edge we store the possible master faces and
1812     // eliminate more and more until we only have one possible master face
1813     // left.
1814
1815     if (debug)
1816     {
1817         Pout<< "subDivisionMatch :"
1818             << " matching (topological) cut faces to masterPatch" << endl;
1819     }
1820
1821     // For every unassigned cutFaceI the possible list of master faces.
1822     Map<labelList> candidates(cutFaces().size());
1823
1824     cutToMasterFaces_.setSize(cutFaces().size());
1825     cutToMasterFaces_ = -1;
1826
1827     while (true)
1828     {
1829         // See if there are any edges left where there is only one remaining
1830         // candidate.
1831         label nChanged = matchEdgeFaces(cutToMasterEdges, candidates);
1832
1833         checkMatch(cutToMasterEdges);
1834
1835
1836         // Now we should have found a face correspondence for every cutFace
1837         // that uses two or more cutEdges. Floodfill from these across
1838         // 'internal' cutedges (i.e. edges that do not have a master
1839         // equivalent) to determine some more cutToMasterFaces
1840         nChanged += growCutFaces(cutToMasterEdges, candidates);
1841
1842         checkMatch(cutToMasterEdges);
1843
1844         if (nChanged == 0)
1845         {
1846             break;
1847         }
1848     }
1849
1850     if (debug)
1851     {
1852         Pout<< "subDivisionMatch :"
1853             << " matching (geometric) cut faces to masterPatch" << endl;
1854     }
1855
1856     // Above should have matched all cases fully. Cannot prove this yet so
1857     // do any remaining unmatched edges by a geometric test.
1858     while (true)
1859     {
1860         // See if there are any edges left where there is only one remaining
1861         // candidate.
1862         label nChanged = geometricMatchEdgeFaces(candidates);
1863
1864         checkMatch(cutToMasterEdges);
1865
1866         nChanged += growCutFaces(cutToMasterEdges, candidates);
1867
1868         checkMatch(cutToMasterEdges);
1869
1870         if (nChanged == 0)
1871         {
1872             break;
1873         }
1874     }
1875
1876
1877     // All cut faces matched?
1878     forAll(cutToMasterFaces_, cutFaceI)
1879     {
1880         if (cutToMasterFaces_[cutFaceI] == -1)
1881         {
1882             const face& cutF = cutFaces()[cutFaceI];
1883
1884             FatalErrorIn
1885             (
1886                 "faceCoupleInfo::subDivisionMatch"
1887                 "(const polyMesh&, const bool, const scalar)"
1888             )   << "Did not match all of cutFaces to a master face" << nl
1889                 << "First unmatched cut face:" << cutFaceI << " with points:"
1890                 << UIndirectList<point>(cutFaces().points(), cutF)()
1891                 << nl
1892                 << "This usually means that the slave patch is not a"
1893                 << " subdivision of the master patch"
1894                 << abort(FatalError);
1895         }
1896     }
1897
1898     if (debug)
1899     {
1900         Pout<< "subDivisionMatch :"
1901             << " finished matching master and slave to cut" << endl;
1902     }
1903
1904     if (debug)
1905     {
1906         writePointsFaces();
1907     }
1908 }
1909
1910
1911 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
1912
1913 // Construct from mesh data
1914 Foam::faceCoupleInfo::faceCoupleInfo
1915 (
1916     const polyMesh& masterMesh,
1917     const polyMesh& slaveMesh,
1918     const scalar absTol,
1919     const bool perfectMatch
1920 )
1921 :
1922     masterPatchPtr_(NULL),
1923     slavePatchPtr_(NULL),
1924     cutPoints_(0),
1925     cutFacesPtr_(NULL),
1926     cutToMasterFaces_(0),
1927     masterToCutPoints_(0),
1928     cutToSlaveFaces_(0),
1929     slaveToCutPoints_(0),
1930     cutEdgeToPoints_(0)
1931 {
1932     // Get faces on both meshes that are aligned.
1933     // (not ordered i.e. masterToMesh[0] does
1934     // not couple to slaveToMesh[0]. This ordering is done later on)
1935     labelList masterToMesh;
1936     labelList slaveToMesh;
1937
1938     if (perfectMatch)
1939     {
1940         // Identical faces so use tight face-centre comparison.
1941         findPerfectMatchingFaces
1942         (
1943             masterMesh,
1944             slaveMesh,
1945             absTol,
1946             masterToMesh,
1947             slaveToMesh
1948         );
1949     }
1950     else
1951     {
1952         // Slave subdivision of master so use 'nearest'. Bit dodgy, especially
1953         // with using absTol (which is quite small)
1954         findSlavesCoveringMaster
1955         (
1956             masterMesh,
1957             slaveMesh,
1958             absTol,
1959             masterToMesh,
1960             slaveToMesh
1961         );
1962     }
1963
1964     // Construct addressing engines for both sides
1965     masterPatchPtr_.reset
1966     (
1967         new indirectPrimitivePatch
1968         (
1969             IndirectList<face>(masterMesh.faces(), masterToMesh),
1970             masterMesh.points()
1971         )
1972     );
1973
1974     slavePatchPtr_.reset
1975     (
1976         new indirectPrimitivePatch
1977         (
1978             IndirectList<face>(slaveMesh.faces(), slaveToMesh),
1979             slaveMesh.points()
1980         )
1981     );
1982
1983
1984     if (perfectMatch)
1985     {
1986         // Faces are perfectly aligned but probably not ordered.
1987         perfectPointMatch(absTol, false);
1988     }
1989     else
1990     {
1991         // Slave faces are subdivision of master face. Faces not ordered.
1992         subDivisionMatch(slaveMesh, false, absTol);
1993     }
1994
1995     if (debug)
1996     {
1997         writePointsFaces();
1998     }
1999 }
2000
2001
2002 // Slave is subdivision of master patch.
2003 // (so -both cover the same area -all of master points are present in slave)
2004 Foam::faceCoupleInfo::faceCoupleInfo
2005 (
2006     const polyMesh& masterMesh,
2007     const labelList& masterAddressing,
2008     const polyMesh& slaveMesh,
2009     const labelList& slaveAddressing,
2010     const scalar absTol,
2011     const bool perfectMatch,
2012     const bool orderedFaces,
2013     const bool patchDivision
2014 )
2015 :
2016     masterPatchPtr_
2017     (
2018         new indirectPrimitivePatch
2019         (
2020             IndirectList<face>(masterMesh.faces(), masterAddressing),
2021             masterMesh.points()
2022         )
2023     ),
2024     slavePatchPtr_
2025     (
2026         new indirectPrimitivePatch
2027         (
2028             IndirectList<face>(slaveMesh.faces(), slaveAddressing),
2029             slaveMesh.points()
2030         )
2031     ),
2032     cutPoints_(0),
2033     cutFacesPtr_(NULL),
2034     cutToMasterFaces_(0),
2035     masterToCutPoints_(0),
2036     cutToSlaveFaces_(0),
2037     slaveToCutPoints_(0),
2038     cutEdgeToPoints_(0)
2039 {
2040     if (perfectMatch && (masterAddressing.size() != slaveAddressing.size()))
2041     {
2042         FatalErrorIn
2043         (
2044             "faceCoupleInfo::faceCoupleInfo(const primitiveMesh&"
2045             ", const primitiveMesh&, const scalar, const bool"
2046         )   << "Perfect match specified but number of master and slave faces"
2047             << " differ." << endl
2048             << "master:" << masterAddressing.size()
2049             << "  slave:" << slaveAddressing.size()
2050             << abort(FatalError);
2051     }
2052
2053     if
2054     (
2055         masterAddressing.size()
2056      && min(masterAddressing) < masterMesh.nInternalFaces()
2057     )
2058     {
2059         FatalErrorIn
2060         (
2061             "faceCoupleInfo::faceCoupleInfo(const primitiveMesh&"
2062             ", const primitiveMesh&, const scalar, const bool"
2063         )   << "Supplied internal face on master mesh to couple." << nl
2064             << "Faces to be coupled have to be boundary faces."
2065             << abort(FatalError);
2066     }
2067     if
2068     (
2069         slaveAddressing.size()
2070      && min(slaveAddressing) < slaveMesh.nInternalFaces()
2071     )
2072     {
2073         FatalErrorIn
2074         (
2075             "faceCoupleInfo::faceCoupleInfo(const primitiveMesh&"
2076             ", const primitiveMesh&, const scalar, const bool"
2077         )   << "Supplied internal face on slave mesh to couple." << nl
2078             << "Faces to be coupled have to be boundary faces."
2079             << abort(FatalError);
2080     }
2081
2082
2083     if (perfectMatch)
2084     {
2085         perfectPointMatch(absTol, orderedFaces);
2086     }
2087     else
2088     {
2089         // Slave faces are subdivision of master face. Faces not ordered.
2090         subDivisionMatch(slaveMesh, patchDivision, absTol);
2091     }
2092
2093     if (debug)
2094     {
2095         writePointsFaces();
2096     }
2097 }
2098
2099
2100 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
2101
2102 Foam::faceCoupleInfo::~faceCoupleInfo()
2103 {}
2104
2105
2106 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
2107
2108 Foam::labelList Foam::faceCoupleInfo::faceLabels(const polyPatch& pp)
2109 {
2110     labelList faces(pp.size());
2111
2112     label faceI = pp.start();
2113
2114     forAll(pp, i)
2115     {
2116         faces[i] = faceI++;
2117     }
2118     return faces;
2119 }
2120
2121
2122 Foam::Map<Foam::label> Foam::faceCoupleInfo::makeMap(const labelList& lst)
2123 {
2124     Map<label> map(lst.size());
2125
2126     forAll(lst, i)
2127     {
2128         if (lst[i] != -1)
2129         {
2130             map.insert(i, lst[i]);
2131         }
2132     }
2133     return map;
2134 }
2135
2136
2137 Foam::Map<Foam::labelList> Foam::faceCoupleInfo::makeMap
2138 (
2139     const labelListList& lst
2140 )
2141 {
2142     Map<labelList> map(lst.size());
2143
2144     forAll(lst, i)
2145     {
2146         if (lst[i].size())
2147         {
2148             map.insert(i, lst[i]);
2149         }
2150     }
2151     return map;
2152 }
2153
2154
2155 // ************************************************************************* //