src/meshTools/triSurface/booleanOps/surfaceIntersection/surfaceIntersection.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of OpenFOAM.
  10
  11     OpenFOAM is free software; you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by the
  13     Free Software Foundation; either version 2 of the License, or (at your
  14     option) any later version.
  15
  16     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with OpenFOAM; if not, write to the Free Software Foundation,
  23     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
  24
  25 Description
  26
  27 \*---------------------------------------------------------------------------*/
  28
  29 #include "surfaceIntersection.H"
  30 #include "triSurfaceSearch.H"
  31 #include "labelPairLookup.H"
  32 #include "OFstream.H"
  33 #include "HashSet.H"
  34 #include "triSurface.H"
  35 #include "pointIndexHit.H"
  36 #include "octreeDataTriSurface.H"
  37 #include "octree.H"
  38 #include "mergePoints.H"
  39 #include "plane.H"
  40 #include "edgeIntersections.H"
  41
  42 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  43
  44 defineTypeNameAndDebug(Foam::surfaceIntersection, 0);
  45
  46
  47 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  48
  49 // Checks if there exists a special topological situation that causes
  50 // edge and the face it hit not to be recognized.
  51 //
  52 // For now if the face shares a point with the edge
  53 bool Foam::surfaceIntersection::excludeEdgeHit
  54 (
  55     const triSurface& surf,
  56     const label edgeI,
  57     const label faceI,
  58     const scalar
  59 )
  60 {
  61     const labelledTri& f = surf.localFaces()[faceI];
  62
  63     const edge& e = surf.edges()[edgeI];
  64
  65     if
  66     (
  67         (f[0] == e.start())
  68      || (f[0] == e.end())
  69      || (f[1] == e.start())
  70      || (f[1] == e.end())
  71      || (f[2] == e.start())
  72      || (f[2] == e.end())
  73     )
  74     {
  75         return true;
  76
  77 //        // Get edge vector
  78 //        vector eVec = e.vec(surf.localPoints());
  79 //        eVec /= mag(eVec) + VSMALL;
  80 //
  81 //        const labelList& eLabels = surf.faceEdges()[faceI];
  82 //
  83 //        // Get edge vector of 0th edge of face
  84 //        vector e0Vec = surf.edges()[eLabels[0]].vec(surf.localPoints());
  85 //        e0Vec /= mag(e0Vec) + VSMALL;
  86 //
  87 //        vector n = e0Vec ^ eVec;
  88 //
  89 //        if (mag(n) < SMALL)
  90 //        {
  91 //            // e0 is aligned with e. Choose next edge of face.
  92 //            vector e1Vec = surf.edges()[eLabels[1]].vec(surf.localPoints());
  93 //            e1Vec /= mag(e1Vec) + VSMALL;
  94 //
  95 //            n = e1Vec ^ eVec;
  96 //
  97 //            if (mag(n) < SMALL)
  98 //            {
  99 //                // Problematic triangle. Two edges aligned with edgeI. Give
 100 //                // up.
 101 //                return true;
 102 //            }
 103 //        }
 104 //
 105 //        // Check if same as faceNormal
 106 //        if (mag(n & surf.faceNormals()[faceI]) > 1-tol)
 107 //        {
 108 //
 109 //            Pout<< "edge:" << e << "  face:" << faceI
 110 //                << "  e0Vec:" << e0Vec << "  n:" << n
 111 //                << "  normalComponent:" << (n & surf.faceNormals()[faceI])
 112 //                << "  tol:" << tol << endl;
 113 //
 114 //            return true;
 115 //        }
 116 //        else
 117 //        {
 118 //            return false;
 119 //        }
 120     }
 121     else
 122     {
 123         return false;
 124     }
 125 }
 126
 127
 128 //// Find intersection of plane with edges of hitFaceI. Returns
 129 //// - edgeI
 130 //// - intersection point
 131 //Foam::pointIndexHit Foam::surfaceIntersection::faceEdgeIntersection
 132 //(
 133 //    const triSurface& surf,
 134 //    const label hitFaceI,
 135 //
 136 //    const vector& n,
 137 //    const point& eStart,
 138 //    const point& eEnd
 139 //)
 140 //{
 141 //    pointIndexHit pInter;
 142 //
 143 //    const pointField& points = surf.points();
 144 //
 145 //    const labelledTri& f = surf.localFaces()[hitFaceI];
 146 //
 147 //    // Plane for intersect test.
 148 //    plane pl(eStart, n);
 149 //
 150 //    forAll(f, fp)
 151 //    {
 152 //        label fp1 = (fp + 1) % 3;
 153 //
 154 //        const point& start = points[f[fp]];
 155 //        const point& end = points[f[fp1]];
 156 //
 157 //        vector eVec(end - start);
 158 //
 159 //        scalar s = pl.normalIntersect(start, eVec);
 160 //
 161 //        if (s < 0 || s > 1)
 162 //        {
 163 //            pInter.setPoint(start + s*eVec);
 164 //
 165 //            // Check if is correct one: orientation walking
 166 //            //  eStart - eEnd - hitPoint should be opposite n
 167 //            vector n2(triPointRef(start, end, pInter.hitPoint()).normal());
 168 //
 169 //            Pout<< "plane normal:" << n
 170 //                << "  start:" << start << "  end:" << end
 171 //                << "  hit at:" << pInter.hitPoint()
 172 //                << "  resulting normal:" << n2 << endl;
 173 //
 174 //            if ((n2 & n) < 0)
 175 //            {
 176 //                pInter.setHit();
 177 //
 178 //                // Find corresponding edge between f[fp] f[fp1]
 179 //                label edgeI =
 180 //                    meshTools::findEdge
 181 //                    (
 182 //                        surf.edges(),
 183 //                        surf.faceEdges()[hitFaceI],
 184 //                        f[fp],
 185 //                        f[fp1]
 186 //                    );
 187 //
 188 //                pInter.setIndex(edgeI);
 189 //
 190 //                return pInter;
 191 //            }
 192 //        }
 193 //    }
 194 //
 195 //    FatalErrorIn("surfaceIntersection::borderEdgeIntersection")
 196 //        << "Did not find intersection of plane " << pl
 197 //        << " with edges of face " << hitFaceI << " verts:" << f
 198 //        << abort(FatalError);
 199 //
 200 //    return pInter;
 201 //}
 202
 203
 204
 205 void Foam::surfaceIntersection::storeIntersection
 206 (
 207     const bool isFirstSurf,
 208     const labelList& facesA,
 209     const label faceB,
 210     DynamicList<edge>& allCutEdges,
 211     DynamicList<point>& allCutPoints
 212 )
 213 {
 214
 215     forAll(facesA, facesAI)
 216     {
 217         label faceA = facesA[facesAI];
 218
 219         // Combine two faces. Always make sure the face from the first surface
 220         // is element 0.
 221         FixedList<label, 2> twoFaces;
 222         if (isFirstSurf)
 223         {
 224             twoFaces[0] = faceA;
 225             twoFaces[1] = faceB;
 226         }
 227         else
 228         {
 229             twoFaces[0] = faceB;
 230             twoFaces[1] = faceA;
 231         }
 232
 233         labelPairLookup::const_iterator iter = facePairToVertex_.find(twoFaces);
 234
 235         if (iter == facePairToVertex_.end())
 236         {
 237             // New intersection. Store face-face intersection.
 238             facePairToVertex_.insert(twoFaces, allCutPoints.size()-1);
 239         }
 240         else
 241         {
 242             // Second occurrence of surf1-surf2 intersection.
 243             // Or rather the face on surf1 intersects a face on
 244             // surface2 twice -> we found edge.
 245
 246             // Check whether perhaps degenerate
 247             const point& prevHit = allCutPoints[*iter];
 248
 249             const point& thisHit = allCutPoints[allCutPoints.size()-1];
 250
 251             if (mag(prevHit - thisHit) < SMALL)
 252             {
 253                 WarningIn
 254                 (
 255                     "Foam::surfaceIntersection::storeIntersection"
 256                     "(const bool isFirstSurf, const labelList& facesA,"
 257                     "const label faceB, DynamicList<edge>& allCutEdges,"
 258                     "DynamicList<point>& allCutPoints)"
 259                 )   << "Encountered degenerate edge between face "
 260                     << twoFaces[0] << " on first surface"
 261                     << " and face " << twoFaces[1] << " on second surface"
 262                     << endl
 263                     << "Point on first surface:" << prevHit << endl
 264                     << "Point on second surface:" << thisHit << endl
 265                     << endl;
 266             }
 267             else
 268             {
 269                 allCutEdges.append(edge(*iter, allCutPoints.size()-1));
 270
 271                 // Remember face on surf
 272                 facePairToEdge_.insert(twoFaces, allCutEdges.size()-1);
 273             }
 274         }
 275     }
 276 }
 277
 278
 279 // Classify cut of edge of surface1 with surface2:
 280 // 1- point of edge hits point on surface2
 281 // 2- edge pierces point on surface2
 282 // 3- point of edge hits edge on surface2
 283 // 4- edge pierces edge on surface2
 284 // 5- point of edge hits face on surface2
 285 // 6- edge pierces face on surface2
 286 //
 287 // Note that handling of 2 and 4 should be the same but with surface1 and
 288 // surface2 reversed.
 289 void Foam::surfaceIntersection::classifyHit
 290 (
 291     const triSurface& surf1,
 292     const scalarField& surf1PointTol,
 293     const triSurface& surf2,
 294     const bool isFirstSurf,
 295     const label edgeI,
 296     const scalar tolDim,
 297     const pointIndexHit& pHit,
 298
 299     DynamicList<edge>& allCutEdges,
 300     DynamicList<point>& allCutPoints,
 301     List<DynamicList<label> >& surfEdgeCuts
 302 )
 303 {
 304     const edge& e = surf1.edges()[edgeI];
 305
 306     const labelList& facesA = surf1.edgeFaces()[edgeI];
 307
 308     // Label of face on surface2 edgeI intersected
 309     label surf2FaceI = pHit.index();
 310
 311     // Classify point on surface2
 312
 313     const labelledTri& f2 = surf2.localFaces()[surf2FaceI];
 314
 315     const pointField& surf2Pts = surf2.localPoints();
 316
 317     label nearType;
 318     label nearLabel;
 319
 320     (void)triPointRef
 321     (
 322         surf2Pts[f2[0]],
 323         surf2Pts[f2[1]],
 324         surf2Pts[f2[2]]
 325     ).classify(pHit.hitPoint(), tolDim, nearType, nearLabel);
 326
 327     // Classify points on edge of surface1
 328     label edgeEnd =
 329         classify
 330         (
 331             surf1PointTol[e.start()],
 332             surf1PointTol[e.end()],
 333             pHit.hitPoint(),
 334             e,
 335             surf1.localPoints()
 336         );
 337
 338     if (nearType == triPointRef::POINT)
 339     {
 340         if (edgeEnd >= 0)
 341         {
 342             // 1. Point hits point. Do nothing.
 343             if (debug&2)
 344             {
 345                 Pout<< pHit.hitPoint() << " is surf1:"
 346                     << " end point of edge " << e
 347                     << " surf2: vertex " << f2[nearLabel]
 348                     << " coord:" << surf2Pts[f2[nearLabel]] << endl;
 349             }
 350         }
 351         else
 352         {
 353             // 2. Edge hits point. Cut edge with new point.
 354             if (debug&2)
 355             {
 356                 Pout<< pHit.hitPoint() << " is surf1:"
 357                     << " somewhere on edge " << e
 358                     << " surf2: vertex " << f2[nearLabel]
 359                     << " coord:" << surf2Pts[f2[nearLabel]] << endl;
 360             }
 361
 362             allCutPoints.append(pHit.hitPoint());
 363             surfEdgeCuts[edgeI].append(allCutPoints.size()-1);
 364
 365             const labelList& facesB = surf2.pointFaces()[f2[nearLabel]];
 366
 367             forAll(facesB, faceBI)
 368             {
 369                 storeIntersection
 370                 (
 371                     isFirstSurf,
 372                     facesA,
 373                     facesB[faceBI],
 374                     allCutEdges,
 375                     allCutPoints
 376                 );
 377             }
 378         }
 379     }
 380     else if (nearType == triPointRef::EDGE)
 381     {
 382         if (edgeEnd >= 0)
 383         {
 384             // 3. Point hits edge. Do nothing on this side. Reverse
 385             // is handled by 2 (edge hits point)
 386             label edge2I = getEdge(surf2, surf2FaceI, nearLabel);
 387             const edge& e2 = surf2.edges()[edge2I];
 388
 389             if (debug&2)
 390             {
 391                 Pout<< pHit.hitPoint() << " is surf1:"
 392                     << " end point of edge " << e
 393                     << " surf2: edge " << e2
 394                     << " coords:" << surf2Pts[e2.start()]
 395                     << surf2Pts[e2.end()] << endl;
 396             }
 397         }
 398         else
 399         {
 400             // 4. Edge hits edge.
 401
 402             // Cut edge with new point (creates duplicates when
 403             // doing the surf2 with surf1 intersection but these
 404             // are merged later on)
 405
 406             label edge2I = getEdge(surf2, surf2FaceI, nearLabel);
 407             const edge& e2 = surf2.edges()[edge2I];
 408
 409             if (debug&2)
 410             {
 411                 Pout<< pHit.hitPoint() << " is surf1:"
 412                     << " somewhere on edge " << e
 413                     << " surf2: edge " << e2
 414                     << " coords:" << surf2Pts[e2.start()]
 415                     << surf2Pts[e2.end()] << endl;
 416             }
 417
 418             allCutPoints.append(pHit.hitPoint());
 419             surfEdgeCuts[edgeI].append(allCutPoints.size()-1);
 420
 421             // edge hits all faces on surf2 connected to the edge
 422
 423             if (isFirstSurf)
 424             {
 425                 // edge-edge intersection is symmetric, store only
 426                 // once.
 427                 // edge hits all faces on surf2 connected to the
 428                 // edge
 429
 430                 const labelList& facesB = surf2.edgeFaces()[edge2I];
 431
 432                 forAll(facesB, faceBI)
 433                 {
 434                     storeIntersection
 435                     (
 436                         isFirstSurf,
 437                         facesA,
 438                         facesB[faceBI],
 439                         allCutEdges,
 440                         allCutPoints
 441                     );
 442                 }
 443             }
 444         }
 445     }
 446     else
 447     {
 448         if (edgeEnd >= 0)
 449         {
 450             // 5. Point hits face. Do what? Introduce
 451             // point & triangulation in face?
 452             if (debug&2)
 453             {
 454                 Pout<< pHit.hitPoint() << " is surf1:"
 455                     << " end point of edge " << e
 456                     << " surf2: face " << surf2FaceI
 457                     << endl;
 458             }
 459
 460             //
 461             // Look exactly at what side (of surf2) edge is. Leave out ones on
 462             // inside of surf2 (i.e. on opposite side of normal)
 463             //
 464
 465             // Vertex on/near surf2
 466             label nearVert = -1;
 467
 468             if (edgeEnd == 0)
 469             {
 470                 nearVert = e.start();
 471             }
 472             else
 473             {
 474                 nearVert = e.end();
 475             }
 476
 477             const point& nearPt = surf1.localPoints()[nearVert];
 478
 479             // Vertex away from surf2
 480             label otherVert = e.otherVertex(nearVert);
 481
 482             const point& otherPt = surf1.localPoints()[otherVert];
 483
 484
 485             if (debug)
 486             {
 487                 Pout
 488                     << pHit.hitPoint() << " is surf1:"
 489                     << " end point of edge " << e << " coord:"
 490                     << surf1.localPoints()[nearVert]
 491                     << " surf2: face " << surf2FaceI << endl;
 492             }
 493
 494             vector eVec = otherPt - nearPt;
 495
 496             if ((surf2.faceNormals()[surf2FaceI] & eVec) > 0)
 497             {
 498                 // otherVert on outside of surf2
 499
 500                 // Shift hitPoint a bit along edge.
 501                 //point hitPt = nearPt + 0.1*eVec;
 502                 point hitPt = nearPt;
 503
 504                 if (debug&2)
 505                 {
 506                     Pout<< "Shifted " << pHit.hitPoint()
 507                         << " to " << hitPt
 508                         << " along edge:" << e
 509                         << " coords:" << surf1.localPoints()[e.start()]
 510                         << surf1.localPoints()[e.end()] << endl;
 511                 }
 512
 513                 // Reclassify as normal edge-face pierce (see below)
 514
 515                 allCutPoints.append(hitPt);
 516                 surfEdgeCuts[edgeI].append(allCutPoints.size()-1);
 517
 518                 // edge hits single face only
 519                 storeIntersection
 520                 (
 521                     isFirstSurf,
 522                     facesA,
 523                     surf2FaceI,
 524                     allCutEdges,
 525                     allCutPoints
 526                 );
 527             }
 528             else
 529             {
 530                 if (debug&2)
 531                 {
 532                     Pout<< "Discarding " << pHit.hitPoint()
 533                         << " since edge " << e << " on inside of surf2."
 534                         << " surf2 normal:" << surf2.faceNormals()[surf2FaceI]
 535                         << endl;
 536                 }
 537             }
 538         }
 539         else
 540         {
 541             // 6. Edge pierces face. 'Normal' situation.
 542             if (debug&2)
 543             {
 544                 Pout<< pHit.hitPoint() << " is surf1:"
 545                     << " somewhere on edge " << e
 546                     << " surf2: face " << surf2FaceI
 547                     << endl;
 548             }
 549
 550             // edgeI intersects surf2. Store point.
 551             allCutPoints.append(pHit.hitPoint());
 552             surfEdgeCuts[edgeI].append(allCutPoints.size()-1);
 553
 554             // edge hits single face only
 555             storeIntersection
 556             (
 557                 isFirstSurf,
 558                 facesA,
 559                 surf2FaceI,
 560                 allCutEdges,
 561                 allCutPoints
 562             );
 563         }
 564     }
 565     if (debug&2)
 566     {
 567         Pout<< endl;
 568     }
 569 }
 570
 571
 572 // Cut all edges of surf1 with surf2. Sets
 573 // - cutPoints          : coordinates of cutPoints
 574 // - cutEdges           : newly created edges between cutPoints
 575 // - facePairToVertex   : hash from face1I and face2I to cutPoint
 576 // - facePairToEdge     : hash from face1I and face2I to cutEdge
 577 // - surfEdgeCuts       : gives for each edge the cutPoints
 578 //                        (in order from start to end)
 579 //
 580 void Foam::surfaceIntersection::doCutEdges
 581 (
 582     const triSurface& surf1,
 583     const triSurfaceSearch& querySurf2,
 584     const bool isFirstSurf,
 585     const bool isSelfIntersection,
 586
 587     DynamicList<edge>& allCutEdges,
 588     DynamicList<point>& allCutPoints,
 589     List<DynamicList<label> >& surfEdgeCuts
 590 )
 591 {
 592     scalar oldTol = intersection::setPlanarTol(1E-3);
 593
 594     const pointField& surf1Pts = surf1.localPoints();
 595
 596     // Calculate local (to point) tolerance based on min edge length.
 597     scalarField surf1PointTol(surf1Pts.size());
 598
 599     forAll(surf1PointTol, pointI)
 600     {
 601         surf1PointTol[pointI] =
 602             intersection::planarTol()
 603           * minEdgeLen(surf1, pointI);
 604     }
 605
 606     const triSurface& surf2 = querySurf2.surface();
 607
 608     forAll(surf1.edges(), edgeI)
 609     {
 610         const edge& e = surf1.edges()[edgeI];
 611
 612         point pStart = surf1Pts[e.start()];
 613         const point& pEnd = surf1Pts[e.end()];
 614
 615         const point tolVec = intersection::planarTol()*(pEnd-pStart);
 616         const scalar tolDim = mag(tolVec);
 617
 618         bool doTrack = false;
 619         do
 620         {
 621             pointIndexHit pHit = querySurf2.tree().findLine(pStart, pEnd);
 622
 623             if (pHit.hit())
 624             {
 625                 if (isSelfIntersection)
 626                 {
 627                     // Skip all intersections which are hit at endpoints of
 628                     // edge.
 629                     // Problem is that if faces are almost coincident the
 630                     // intersection point will be calculated quite incorrectly
 631                     // The error might easily be larger than 1% of the edge
 632                     // length.
 633                     // So what we do here is to exclude hit faces if our edge
 634                     // is in their plane and they share a point with the edge.
 635
 636                     // Label of face on surface2 edgeI intersected
 637                     label hitFaceI = pHit.index();
 638
 639                     if
 640                     (
 641                         !excludeEdgeHit
 642                         (
 643                             surf1,
 644                             edgeI,
 645                             hitFaceI,
 646                             0.1         // 1-cos of angle between normals
 647                         )
 648                     )
 649                     {
 650                         // Classify point on surface1
 651                         label edgeEnd = classify
 652                         (
 653                             surf1PointTol[e.start()],
 654                             surf1PointTol[e.end()],
 655                             pHit.hitPoint(),
 656                             e,
 657                             surf1Pts
 658                         );
 659
 660                         if (edgeEnd < 0)
 661                         {
 662                             if (debug)
 663                             {
 664                                 Pout<< "edge:" << edgeI << " vertices:" << e
 665                                     << "  start:" << surf1Pts[e.start()]
 666                                     << "  end:" << surf1Pts[e.end()]
 667                                     << "  hit:" << pHit.hitPoint()
 668                                     << "  tolDim:" << tolDim
 669                                     << "  planarTol:"
 670                                     << intersection::planarTol()
 671                                     << endl;
 672                             }
 673                             allCutPoints.append(pHit.hitPoint());
 674                             surfEdgeCuts[edgeI].append(allCutPoints.size()-1);
 675                         }
 676                     }
 677                 }
 678                 else
 679                 {
 680                     classifyHit
 681                     (
 682                         surf1,
 683                         surf1PointTol,
 684                         surf2,
 685                         isFirstSurf,
 686                         edgeI,
 687                         tolDim,
 688                         pHit,
 689
 690                         allCutEdges,
 691                         allCutPoints,
 692                         surfEdgeCuts
 693                     );
 694                 }
 695
 696                 if (mag(pHit.hitPoint() - pEnd) < tolDim)
 697                 {
 698                     doTrack = false;
 699                 }
 700                 else
 701                 {
 702                     pStart = pHit.hitPoint() + tolVec;
 703
 704                     doTrack = true;
 705                 }
 706             }
 707             else
 708             {
 709                 doTrack = false;
 710             }
 711         }
 712         while(doTrack);
 713     }
 714     intersection::setPlanarTol(oldTol);
 715 }
 716
 717
 718 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 719
 720 // Null constructor
 721 Foam::surfaceIntersection::surfaceIntersection()
 722 :
 723     cutPoints_(0),
 724     cutEdges_(0),
 725     facePairToVertex_(0),
 726     facePairToEdge_(0),
 727     surf1EdgeCuts_(0),
 728     surf2EdgeCuts_(0)
 729 {}
 730
 731
 732 // Construct from two surfaces
 733 Foam::surfaceIntersection::surfaceIntersection
 734 (
 735     const triSurfaceSearch& query1,
 736     const triSurfaceSearch& query2
 737 )
 738 :
 739     cutPoints_(0),
 740     cutEdges_(0),
 741     facePairToVertex_(2*max(query1.surface().size(), query2.surface().size())),
 742     facePairToEdge_(2*max(query1.surface().size(), query2.surface().size())),
 743     surf1EdgeCuts_(0),
 744     surf2EdgeCuts_(0)
 745 {
 746     const triSurface& surf1 = query1.surface();
 747     const triSurface& surf2 = query2.surface();
 748
 749     //
 750     // Cut all edges of surf1 with surf2.
 751     //
 752     if (debug)
 753     {
 754         Pout<< "Cutting surf1 edges" << endl;
 755     }
 756
 757
 758     DynamicList<edge> allCutEdges(surf1.nEdges()/20);
 759     DynamicList<point> allCutPoints(surf1.nPoints()/20);
 760
 761
 762     // From edge to cut index on surface1
 763     List<DynamicList<label> > edgeCuts1(query1.surface().nEdges());
 764
 765     doCutEdges
 766     (
 767         surf1,
 768         query2,
 769         true,               // is first surface; construct labelPair in correct
 770                             // order
 771         false,              // not self intersection
 772
 773         allCutEdges,
 774         allCutPoints,
 775         edgeCuts1
 776     );
 777     // Transfer to straight labelListList
 778     transfer(edgeCuts1, surf1EdgeCuts_);
 779
 780
 781     //
 782     // Cut all edges of surf2 with surf1.
 783     //
 784     if (debug)
 785     {
 786         Pout<< "Cutting surf2 edges" << endl;
 787     }
 788
 789     // From edge to cut index
 790     List<DynamicList<label> > edgeCuts2(query2.surface().nEdges());
 791
 792     doCutEdges
 793     (
 794         surf2,
 795         query1,
 796         false,              // is second surface
 797         false,              // not self intersection
 798
 799         allCutEdges,
 800         allCutPoints,
 801         edgeCuts2
 802     );
 803
 804     // Transfer to straight label(List)List
 805     transfer(edgeCuts2, surf2EdgeCuts_);
 806     cutEdges_.transfer(allCutEdges);
 807     cutPoints_.transfer(allCutPoints);
 808
 809
 810     if (debug)
 811     {
 812         Pout<< "surfaceIntersection : Intersection generated:"
 813             << endl
 814             << "    points:" << cutPoints_.size() << endl
 815             << "    edges :" << cutEdges_.size() << endl;
 816
 817         Pout<< "surfaceIntersection : Writing intersection to intEdges.obj"
 818             << endl;
 819
 820         OFstream intStream("intEdges.obj");
 821         writeOBJ(cutPoints_, cutEdges_, intStream);
 822
 823         // Dump all cut edges to files
 824         Pout<< "Dumping cut edges of surface1 to surf1EdgeCuts.obj" << endl;
 825         OFstream edge1Stream("surf1EdgeCuts.obj");
 826         writeIntersectedEdges(surf1, surf1EdgeCuts_, edge1Stream);
 827
 828         Pout<< "Dumping cut edges of surface2 to surf2EdgeCuts.obj" << endl;
 829         OFstream edge2Stream("surf2EdgeCuts.obj");
 830         writeIntersectedEdges(surf2, surf2EdgeCuts_, edge2Stream);
 831     }
 832 }
 833
 834
 835 // Construct from full intersection Poutrmation
 836 Foam::surfaceIntersection::surfaceIntersection
 837 (
 838     const triSurface& surf1,
 839     const edgeIntersections& intersections1,
 840     const triSurface& surf2,
 841     const edgeIntersections& intersections2
 842 )
 843 :
 844     cutPoints_(0),
 845     cutEdges_(0),
 846     facePairToVertex_(2*max(surf1.size(), surf2.size())),
 847     facePairToEdge_(2*max(surf1.size(), surf2.size())),
 848     surf1EdgeCuts_(0),
 849     surf2EdgeCuts_(0)
 850 {
 851
 852     // All intersection Pout (so for both surfaces)
 853     DynamicList<edge> allCutEdges((surf1.nEdges() + surf2.nEdges())/20);
 854     DynamicList<point> allCutPoints((surf1.nPoints() + surf2.nPoints())/20);
 855
 856
 857     // Cut all edges of surf1 with surf2
 858     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 859
 860     if (debug)
 861     {
 862         Pout<< "Storing surf1 intersections" << endl;
 863     }
 864
 865     {
 866         // From edge to cut index on surface1
 867         List<DynamicList<label> > edgeCuts1(surf1.nEdges());
 868
 869         forAll(intersections1, edgeI)
 870         {
 871             const List<pointIndexHit>& intersections = intersections1[edgeI];
 872
 873             forAll(intersections, i)
 874             {
 875                 const pointIndexHit& pHit = intersections[i];
 876
 877                 // edgeI intersects surf2. Store point.
 878                 allCutPoints.append(pHit.hitPoint());
 879                 edgeCuts1[edgeI].append(allCutPoints.size()-1);
 880
 881                 storeIntersection
 882                 (
 883                     true,                       // is first surface
 884                     surf1.edgeFaces()[edgeI],
 885                     pHit.index(),               // surf2FaceI
 886                     allCutEdges,
 887                     allCutPoints
 888                 );
 889             }
 890         }
 891
 892         // Transfer to straight labelListList
 893         transfer(edgeCuts1, surf1EdgeCuts_);
 894     }
 895
 896
 897
 898     // Cut all edges of surf2 with surf1
 899     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 900
 901     if (debug)
 902     {
 903         Pout<< "Storing surf2 intersections" << endl;
 904     }
 905
 906     {
 907         // From edge to cut index on surface2
 908         List<DynamicList<label> > edgeCuts2(surf2.nEdges());
 909
 910         forAll(intersections2, edgeI)
 911         {
 912             const List<pointIndexHit>& intersections = intersections2[edgeI];
 913
 914             forAll(intersections, i)
 915             {
 916                 const pointIndexHit& pHit = intersections[i];
 917
 918                 // edgeI intersects surf1. Store point.
 919                 allCutPoints.append(pHit.hitPoint());
 920                 edgeCuts2[edgeI].append(allCutPoints.size()-1);
 921
 922                 storeIntersection
 923                 (
 924                     false,                      // is second surface
 925                     surf2.edgeFaces()[edgeI],
 926                     pHit.index(),               // surf2FaceI
 927                     allCutEdges,
 928                     allCutPoints
 929                 );
 930             }
 931         }
 932
 933         // Transfer to surf2EdgeCuts_ (straight labelListList)
 934         transfer(edgeCuts2, surf2EdgeCuts_);
 935     }
 936
 937
 938     // Transfer to straight label(List)List
 939     cutEdges_.transfer(allCutEdges);
 940     cutPoints_.transfer(allCutPoints);
 941
 942
 943     if (debug)
 944     {
 945         Pout<< "surfaceIntersection : Intersection generated:"
 946             << endl
 947             << "    points:" << cutPoints_.size() << endl
 948             << "    edges :" << cutEdges_.size() << endl;
 949
 950         Pout<< "surfaceIntersection : Writing intersection to intEdges.obj"
 951             << endl;
 952
 953         OFstream intStream("intEdges.obj");
 954         writeOBJ(cutPoints_, cutEdges_, intStream);
 955
 956         // Dump all cut edges to files
 957         Pout<< "Dumping cut edges of surface1 to surf1EdgeCuts.obj" << endl;
 958         OFstream edge1Stream("surf1EdgeCuts.obj");
 959         writeIntersectedEdges(surf1, surf1EdgeCuts_, edge1Stream);
 960
 961         Pout<< "Dumping cut edges of surface2 to surf2EdgeCuts.obj" << endl;
 962         OFstream edge2Stream("surf2EdgeCuts.obj");
 963         writeIntersectedEdges(surf2, surf2EdgeCuts_, edge2Stream);
 964     }
 965
 966
 967     // Debugging stuff
 968     {
 969         // Check all facePairToVertex is used.
 970         labelHashSet usedPoints;
 971
 972         forAllConstIter(labelPairLookup, facePairToEdge_, iter)
 973         {
 974             label edgeI = iter();
 975
 976             const edge& e = cutEdges_[edgeI];
 977
 978             usedPoints.insert(e[0]);
 979             usedPoints.insert(e[1]);
 980         }
 981
 982         forAllConstIter(labelPairLookup, facePairToVertex_, iter)
 983         {
 984             label pointI = iter();
 985
 986             if (!usedPoints.found(pointI))
 987             {
 988                 FatalErrorIn("surfaceIntersection::surfaceIntersection")
 989                     << "Problem: cut point:" << pointI
 990                     << " coord:" << cutPoints_[pointI]
 991                     << " not used by any edge" << abort(FatalError);
 992             }
 993         }
 994     }
 995 }
 996
 997
 998 // Construct from single surface. Used to test for self-intersection.
 999 Foam::surfaceIntersection::surfaceIntersection
1000 (
1001     const triSurfaceSearch& query1
1002 )
1003 :
1004     cutPoints_(0),
1005     cutEdges_(0),
1006     facePairToVertex_(2*query1.surface().size()),
1007     facePairToEdge_(2*query1.surface().size()),
1008     surf1EdgeCuts_(0),
1009     surf2EdgeCuts_(0)
1010 {
1011     const triSurface& surf1 = query1.surface();
1012
1013     //
1014     // Cut all edges of surf1 with surf1 itself.
1015     //
1016     if (debug)
1017     {
1018         Pout<< "Cutting surf1 edges" << endl;
1019     }
1020
1021     DynamicList<edge> allCutEdges;
1022     DynamicList<point> allCutPoints;
1023
1024     // From edge to cut index on surface1
1025     List<DynamicList<label> > edgeCuts1(query1.surface().nEdges());
1026
1027     doCutEdges
1028     (
1029         surf1,
1030         query1,
1031         true,               // is first surface; construct labelPair in correct
1032                             // order
1033         true,               // self intersection
1034
1035
1036         allCutEdges,
1037         allCutPoints,
1038         edgeCuts1
1039     );
1040
1041     // Transfer to straight label(List)List
1042     transfer(edgeCuts1, surf1EdgeCuts_);
1043     cutEdges_.transfer(allCutEdges);
1044     cutPoints_.transfer(allCutPoints);
1045
1046     // Shortcut.
1047     if (cutPoints_.empty() && cutEdges_.empty())
1048     {
1049         if (debug)
1050         {
1051             Pout<< "Empty intersection" << endl;
1052         }
1053         return;
1054     }
1055
1056     //
1057     // Remove duplicate points (from edge-point or edge-edge cutting)
1058     //
1059
1060     // Get typical dimension.
1061     scalar minEdgeLen = GREAT;
1062     forAll(surf1.edges(), edgeI)
1063     {
1064         minEdgeLen = min
1065         (
1066             minEdgeLen,
1067             surf1.edges()[edgeI].mag(surf1.localPoints())
1068         );
1069     }
1070
1071     // Merge points
1072     labelList pointMap;
1073     pointField newPoints;
1074
1075     bool hasMerged = mergePoints
1076     (
1077         cutPoints_,
1078         minEdgeLen*intersection::planarTol(),
1079         false,
1080         pointMap,
1081         newPoints
1082     );
1083
1084     if (hasMerged)
1085     {
1086         if (debug)
1087         {
1088             Pout<< "Merged:" << hasMerged
1089                 << "  mergeDist:" << minEdgeLen*intersection::planarTol()
1090                 << "  cutPoints:" << cutPoints_.size()
1091                 << "  newPoints:" << newPoints.size()
1092                 << endl;
1093         }
1094
1095         // Copy points
1096         cutPoints_.transfer(newPoints);
1097
1098         // Renumber vertices referenced by edges
1099         forAll(cutEdges_, edgeI)
1100         {
1101             edge& e = cutEdges_[edgeI];
1102
1103             e.start() = pointMap[e.start()];
1104             e.end() = pointMap[e.end()];
1105
1106             if (e.mag(cutPoints_) < minEdgeLen*intersection::planarTol())
1107             {
1108                 if (debug)
1109                 {
1110                     Pout<< "Degenerate cut:" << edgeI << " vertices:" << e
1111                         << " coords:" << cutPoints_[e.start()] << ' '
1112                         << cutPoints_[e.end()] << endl;
1113                 }
1114             }
1115         }
1116
1117         // Renumber vertices referenced by edgeCut lists. Remove duplicates.
1118         forAll(surf1EdgeCuts_, edgeI)
1119         {
1120             // Get indices of cutPoints this edge is cut by
1121             labelList& cutVerts = surf1EdgeCuts_[edgeI];
1122
1123             removeDuplicates(pointMap, cutVerts);
1124         }
1125     }
1126
1127     if (debug)
1128     {
1129         Pout<< "surfaceIntersection : Intersection generated and compressed:"
1130             << endl
1131             << "    points:" << cutPoints_.size() << endl
1132             << "    edges :" << cutEdges_.size() << endl;
1133
1134
1135         Pout<< "surfaceIntersection : Writing intersection to intEdges.obj"
1136             << endl;
1137
1138         OFstream intStream("intEdges.obj");
1139         writeOBJ(cutPoints_, cutEdges_, intStream);
1140     }
1141
1142     if (debug)
1143     {
1144         // Dump all cut edges to files
1145         Pout<< "Dumping cut edges of surface1 to surf1EdgeCuts.obj" << endl;
1146         OFstream edge1Stream("surf1EdgeCuts.obj");
1147         writeIntersectedEdges(surf1, surf1EdgeCuts_, edge1Stream);
1148     }
1149 }
1150
1151
1152 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
1153
1154 const Foam::pointField& Foam::surfaceIntersection::cutPoints() const
1155 {
1156     return cutPoints_;
1157 }
1158
1159
1160 const Foam::edgeList& Foam::surfaceIntersection::cutEdges() const
1161 {
1162     return cutEdges_;
1163 }
1164
1165
1166 const Foam::labelPairLookup& Foam::surfaceIntersection::facePairToEdge() const
1167 {
1168     return facePairToEdge_;
1169 }
1170
1171
1172 const Foam::labelListList& Foam::surfaceIntersection::edgeCuts
1173 (
1174     const bool isFirstSurf
1175 ) const
1176 {
1177     if (isFirstSurf)
1178     {
1179         return surf1EdgeCuts_;
1180     }
1181     else
1182     {
1183         return surf2EdgeCuts_;
1184     }
1185 }
1186
1187
1188 const Foam::labelListList& Foam::surfaceIntersection::surf1EdgeCuts() const
1189 {
1190     return surf1EdgeCuts_;
1191 }
1192
1193
1194 const Foam::labelListList& Foam::surfaceIntersection::surf2EdgeCuts() const
1195 {
1196     return surf2EdgeCuts_;
1197 }
1198
1199
1200 // ************************************************************************* //