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