src/dynamicMesh/polyMeshAdder/faceCoupleInfo.C

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