src/meshTools/triSurface/booleanOps/intersectedSurface/intersectedSurface.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of OpenFOAM.
  10
  11     OpenFOAM is free software: you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by
  13     the Free Software Foundation, either version 3 of the License, or
  14     (at your option) any later version.
  15
  16     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "intersectedSurface.H"
  27 #include "surfaceIntersection.H"
  28 #include "faceList.H"
  29 #include "faceTriangulation.H"
  30 #include "treeBoundBox.H"
  31 #include "OFstream.H"
  32 #include "error.H"
  33 #include "meshTools.H"
  34 #include "edgeSurface.H"
  35 #include "DynamicList.H"
  36 #include "transform.H"
  37
  38 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  39
  40 defineTypeNameAndDebug(Foam::intersectedSurface, 0);
  41
  42 const Foam::label Foam::intersectedSurface::UNVISITED = 0;
  43 const Foam::label Foam::intersectedSurface::STARTTOEND = 1;
  44 const Foam::label Foam::intersectedSurface::ENDTOSTART = 2;
  45 const Foam::label Foam::intersectedSurface::BOTH = STARTTOEND | ENDTOSTART;
  46
  47 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  48
  49 // Write whole pointField and edges to stream
  50 void Foam::intersectedSurface::writeOBJ
  51 (
  52     const pointField& points,
  53     const edgeList& edges,
  54     Ostream& os
  55 )
  56 {
  57     forAll(points, pointI)
  58     {
  59         const point& pt = points[pointI];
  60
  61         os << "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << nl;
  62     }
  63     forAll(edges, edgeI)
  64     {
  65         const edge& e = edges[edgeI];
  66
  67         os << "l " << e.start()+1 << ' ' << e.end()+1 << nl;
  68     }
  69 }
  70
  71
  72 // Write whole pointField and selected edges to stream
  73 void Foam::intersectedSurface::writeOBJ
  74 (
  75     const pointField& points,
  76     const edgeList& edges,
  77     const labelList& faceEdges,
  78     Ostream& os
  79 )
  80 {
  81     forAll(points, pointI)
  82     {
  83         const point& pt = points[pointI];
  84
  85         os << "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << nl;
  86     }
  87     forAll(faceEdges, i)
  88     {
  89         const edge& e = edges[faceEdges[i]];
  90
  91         os << "l " << e.start()+1 << ' ' << e.end()+1 << nl;
  92     }
  93 }
  94
  95
  96 // write local points and edges to stream
  97 void Foam::intersectedSurface::writeLocalOBJ
  98 (
  99     const pointField& points,
 100     const edgeList& edges,
 101     const labelList& faceEdges,
 102     const fileName& fName
 103 )
 104 {
 105     OFstream os(fName);
 106
 107     labelList pointMap(points.size(), -1);
 108
 109     label maxVertI = 0;
 110
 111     forAll(faceEdges, i)
 112     {
 113         const edge& e = edges[faceEdges[i]];
 114
 115         forAll(e, i)
 116         {
 117             label pointI = e[i];
 118
 119             if (pointMap[pointI] == -1)
 120             {
 121                 const point& pt = points[pointI];
 122
 123                 os << "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << nl;
 124
 125                 pointMap[pointI] = maxVertI++;
 126             }
 127         }
 128     }
 129
 130     forAll(faceEdges, i)
 131     {
 132         const edge& e = edges[faceEdges[i]];
 133
 134         os << "l " << pointMap[e.start()]+1 << ' ' << pointMap[e.end()]+1
 135             << nl;
 136     }
 137 }
 138
 139
 140 // Write whole pointField and face to stream
 141 void Foam::intersectedSurface::writeOBJ
 142 (
 143     const pointField& points,
 144     const face& f,
 145     Ostream& os
 146 )
 147 {
 148     forAll(points, pointI)
 149     {
 150         const point& pt = points[pointI];
 151
 152         os << "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << nl;
 153     }
 154
 155     os << 'f';
 156
 157     forAll(f, fp)
 158     {
 159         os << ' ' << f[fp]+1;
 160     }
 161     os << nl;
 162 }
 163
 164
 165 // Print current visited state.
 166 void Foam::intersectedSurface::printVisit
 167 (
 168     const edgeList& edges,
 169     const labelList& edgeLabels,
 170     const Map<label>& visited
 171 )
 172 {
 173     Pout<< "Visited:" << nl;
 174     forAll(edgeLabels, i)
 175     {
 176         label edgeI = edgeLabels[i];
 177
 178         const edge& e = edges[edgeI];
 179
 180         label stat = visited[edgeI];
 181
 182         if (stat == UNVISITED)
 183         {
 184             Pout<< "    edge:" << edgeI << "  verts:" << e
 185                 << "  unvisited" << nl;
 186         }
 187         else if (stat == STARTTOEND)
 188         {
 189             Pout<< "    edge:" << edgeI << "  from " << e[0]
 190                 << " to " << e[1] << nl;
 191         }
 192         else if (stat == ENDTOSTART)
 193         {
 194             Pout<< "    edge:" << edgeI << "  from " << e[1]
 195                 << " to " << e[0] << nl;
 196         }
 197         else
 198         {
 199             Pout<< "    edge:" << edgeI << "  both " << e
 200                 << nl;
 201         }
 202     }
 203     Pout<< endl;
 204 }
 205
 206
 207 // Check if the two vertices that f0 and f1 share are in the same order on
 208 // both faces.
 209 bool Foam::intersectedSurface::sameEdgeOrder
 210 (
 211     const labelledTri& fA,
 212     const labelledTri& fB
 213 )
 214 {
 215     forAll(fA, fpA)
 216     {
 217         label fpB = findIndex(fB, fA[fpA]);
 218
 219         if (fpB != -1)
 220         {
 221             // Get prev/next vertex on fA
 222             label vA1 = fA[fA.fcIndex(fpA)];
 223             label vAMin1 = fA[fA.rcIndex(fpA)];
 224
 225             // Get prev/next vertex on fB
 226             label vB1 = fB[fB.fcIndex(fpB)];
 227             label vBMin1 = fB[fB.rcIndex(fpB)];
 228
 229             if (vA1 == vB1 || vAMin1 == vBMin1)
 230             {
 231                 return true;
 232             }
 233             else if (vA1 == vBMin1 || vAMin1 == vB1)
 234             {
 235                 // shared vertices in opposite order.
 236                 return false;
 237             }
 238             else
 239             {
 240                 FatalErrorIn("intersectedSurface::sameEdgeOrder")
 241                     << "Triangle:" << fA << " and triangle:" << fB
 242                     << " share a point but not an edge"
 243                     << abort(FatalError);
 244             }
 245         }
 246     }
 247
 248     FatalErrorIn("intersectedSurface::sameEdgeOrder")
 249         << "Triangle:" << fA << " and triangle:" << fB
 250         << " do not share an edge"
 251         << abort(FatalError);
 252
 253     return false;
 254 }
 255
 256
 257 void Foam::intersectedSurface::incCount
 258 (
 259     Map<label>& visited,
 260     const label key,
 261     const label offset
 262 )
 263 {
 264     Map<label>::iterator iter = visited.find(key);
 265
 266     if (iter == visited.end())
 267     {
 268         visited.insert(key, offset);
 269     }
 270     else
 271     {
 272         iter() += offset;
 273     }
 274 }
 275
 276
 277 // Calculate point to edge addressing for the face given by the edge
 278 // subset faceEdges. Constructs facePointEdges which for every point
 279 // gives a list of edge labels connected to it.
 280 Foam::Map<Foam::DynamicList<Foam::label> >
 281 Foam::intersectedSurface::calcPointEdgeAddressing
 282 (
 283     const edgeSurface& eSurf,
 284     const label faceI
 285 )
 286 {
 287     const pointField& points = eSurf.points();
 288     const edgeList& edges = eSurf.edges();
 289
 290     const labelList& fEdges = eSurf.faceEdges()[faceI];
 291
 292     Map<DynamicList<label> > facePointEdges(4*fEdges.size());
 293
 294     forAll(fEdges, i)
 295     {
 296         label edgeI = fEdges[i];
 297
 298         const edge& e = edges[edgeI];
 299
 300         // Add e.start to point-edges
 301         Map<DynamicList<label> >::iterator iter =
 302             facePointEdges.find(e.start());
 303
 304         if (iter == facePointEdges.end())
 305         {
 306             DynamicList<label> oneEdge;
 307             oneEdge.append(edgeI);
 308             facePointEdges.insert(e.start(), oneEdge);
 309         }
 310         else
 311         {
 312             iter().append(edgeI);
 313         }
 314
 315         // Add e.end to point-edges
 316         Map<DynamicList<label> >::iterator iter2 =
 317             facePointEdges.find(e.end());
 318
 319         if (iter2 == facePointEdges.end())
 320         {
 321             DynamicList<label> oneEdge;
 322             oneEdge.append(edgeI);
 323             facePointEdges.insert(e.end(), oneEdge);
 324         }
 325         else
 326         {
 327             iter2().append(edgeI);
 328         }
 329     }
 330
 331     // Shrink it
 332     forAllIter(Map< DynamicList<label> >, facePointEdges, iter)
 333     {
 334         iter().shrink();
 335
 336         // Check on dangling points.
 337         if (iter().empty())
 338         {
 339             FatalErrorIn
 340             (
 341                 "intersectedSurface::calcPointEdgeAddressing"
 342                 "(const edgeSurface&, const label)"
 343             )   << "Point:" << iter.key() << " used by too few edges:"
 344                 << iter() << abort(FatalError);
 345         }
 346     }
 347
 348     if (debug & 2)
 349     {
 350         // Print facePointEdges
 351         Pout<< "calcPointEdgeAddressing: face consisting of edges:" << endl;
 352         forAll(fEdges, i)
 353         {
 354             label edgeI = fEdges[i];
 355             const edge& e = edges[edgeI];
 356             Pout<< "    " << edgeI << ' ' << e
 357                 << points[e.start()]
 358                 << points[e.end()] << endl;
 359         }
 360
 361         Pout<< "    Constructed point-edge adressing:" << endl;
 362         forAllConstIter(Map< DynamicList<label> >, facePointEdges, iter)
 363         {
 364             Pout<< "    vertex " << iter.key() << " is connected to edges "
 365                 << iter() << endl;
 366         }
 367         Pout<<endl;
 368     }
 369
 370     return facePointEdges;
 371 }
 372
 373
 374 // Find next (triangle or cut) edge label coming from point prevVertI on
 375 // prevEdgeI doing a right handed walk (i.e. following right wall).
 376 // Note: normal is provided externally. Could be deducted from angle between
 377 // two triangle edges but these could be in line.
 378 Foam::label Foam::intersectedSurface::nextEdge
 379 (
 380     const edgeSurface& eSurf,
 381     const Map<label>& visited,
 382     const label faceI,
 383     const vector& n,
 384     const Map<DynamicList<label> >& facePointEdges,
 385     const label prevEdgeI,
 386     const label prevVertI
 387 )
 388 {
 389     const pointField& points = eSurf.points();
 390     const edgeList& edges = eSurf.edges();
 391     const labelList& fEdges = eSurf.faceEdges()[faceI];
 392
 393
 394     // Edges connected to prevVertI
 395     const DynamicList<label>& connectedEdges = facePointEdges[prevVertI];
 396
 397     if (connectedEdges.size() <= 1)
 398     {
 399         // Problem. Point not connected.
 400         {
 401             Pout<< "Writing face:" << faceI << " to face.obj" << endl;
 402             OFstream str("face.obj");
 403             writeOBJ(points, edges, fEdges, str);
 404
 405             Pout<< "Writing connectedEdges edges to faceEdges.obj" << endl;
 406             writeLocalOBJ(points, edges, connectedEdges, "faceEdges.obj");
 407         }
 408
 409         FatalErrorIn
 410         (
 411             "intersectedSurface::nextEdge(const pointField&, const edgeList&"
 412             ", const vector&, Map<DynamicList<label> >, const label"
 413             ", const label)"
 414         )   << "Problem: prevVertI:" << prevVertI << " on edge " << prevEdgeI
 415             << " has less than 2 connected edges."
 416             << " connectedEdges:" << connectedEdges << abort(FatalError);
 417
 418         return -1;
 419     }
 420
 421     if (connectedEdges.size() == 2)
 422     {
 423         // Simple case. Take other edge
 424         if (connectedEdges[0] == prevEdgeI)
 425         {
 426             if (debug & 2)
 427             {
 428                 Pout<< "Stepped from edge:" << edges[prevEdgeI]
 429                     << " to edge:" << edges[connectedEdges[1]]
 430                     << " chosen from candidates:" << connectedEdges << endl;
 431             }
 432             return connectedEdges[1];
 433         }
 434         else
 435         {
 436             if (debug & 2)
 437             {
 438                Pout<< "Stepped from edge:" << edges[prevEdgeI]
 439                    << " to edge:" << edges[connectedEdges[0]]
 440                    << " chosen from candidates:" << connectedEdges << endl;
 441             }
 442             return connectedEdges[0];
 443         }
 444     }
 445
 446
 447     // Multiple choices. Look at angle between edges.
 448
 449     const edge& prevE = edges[prevEdgeI];
 450
 451     // x-axis of coordinate system
 452     vector e0 = n ^ (points[prevE.otherVertex(prevVertI)] - points[prevVertI]);
 453     e0 /= mag(e0) + VSMALL;
 454
 455     // Get y-axis of coordinate system
 456     vector e1 = n ^ e0;
 457
 458     if (mag(mag(e1) - 1) > SMALL)
 459     {
 460         {
 461             Pout<< "Writing face:" << faceI << " to face.obj" << endl;
 462             OFstream str("face.obj");
 463             writeOBJ(points, edges, fEdges, str);
 464
 465             Pout<< "Writing connectedEdges edges to faceEdges.obj" << endl;
 466             writeLocalOBJ(points, edges, connectedEdges, "faceEdges.obj");
 467         }
 468
 469         FatalErrorIn("intersectedSurface::nextEdge")
 470             << "Unnormalized normal e1:" << e1
 471             << " formed from cross product of e0:" << e0 << " n:" << n
 472             << abort(FatalError);
 473     }
 474
 475
 476     //
 477     // Determine maximum angle over all connected (unvisited) edges.
 478     //
 479
 480     scalar maxAngle = -GREAT;
 481     label maxEdgeI = -1;
 482
 483     forAll(connectedEdges, connI)
 484     {
 485         label edgeI = connectedEdges[connI];
 486
 487         if (edgeI != prevEdgeI)
 488         {
 489             label stat = visited[edgeI];
 490
 491             const edge& e = edges[edgeI];
 492
 493             // Find out whether walk of edge from prevVert would be acceptible.
 494             if
 495             (
 496                 stat == UNVISITED
 497              || (
 498                     stat == STARTTOEND
 499                  && prevVertI == e[1]
 500                 )
 501              || (
 502                     stat == ENDTOSTART
 503                  && prevVertI == e[0]
 504                 )
 505             )
 506             {
 507                 // Calculate angle of edge with respect to base e0, e1
 508                 vector vec =
 509                     n ^ (points[e.otherVertex(prevVertI)] - points[prevVertI]);
 510
 511                 vec /= mag(vec) + VSMALL;
 512
 513                 scalar angle = pseudoAngle(e0, e1, vec);
 514
 515                 if (angle > maxAngle)
 516                 {
 517                     maxAngle = angle;
 518                     maxEdgeI = edgeI;
 519                 }
 520             }
 521         }
 522     }
 523
 524
 525     if (maxEdgeI == -1)
 526     {
 527         // No unvisited edge found
 528         {
 529             Pout<< "Writing face:" << faceI << " to face.obj" << endl;
 530             OFstream str("face.obj");
 531             writeOBJ(points, edges, fEdges, str);
 532
 533             Pout<< "Writing connectedEdges edges to faceEdges.obj" << endl;
 534             writeLocalOBJ(points, edges, connectedEdges, "faceEdges.obj");
 535         }
 536
 537         FatalErrorIn
 538         (
 539             "intersectedSurface::nextEdge(const pointField&, const edgeList&"
 540             ", const Map<label>&, const vector&"
 541             ", const Map<DynamicList<label> >&"
 542             ", const label, const label"
 543         )   << "Trying to step from edge " << edges[prevEdgeI]
 544             << ", vertex " << prevVertI
 545             << " but cannot find 'unvisited' edges among candidates:"
 546             << connectedEdges
 547             << abort(FatalError);
 548     }
 549
 550     if (debug & 2)
 551     {
 552         Pout<< "Stepped from edge:" << edges[prevEdgeI]
 553             << " to edge:" << maxEdgeI << " angle:" << edges[maxEdgeI]
 554             << " with angle:" << maxAngle
 555             << endl;
 556     }
 557
 558     return maxEdgeI;
 559 }
 560
 561
 562 // Create (polygonal) face by walking full circle starting from startVertI on
 563 // startEdgeI.
 564 // Uses nextEdge(..) to do the walking. Returns face. Updates visited with
 565 // the labels of visited edges.
 566 Foam::face Foam::intersectedSurface::walkFace
 567 (
 568     const edgeSurface& eSurf,
 569     const label faceI,
 570     const vector& n,
 571     const Map<DynamicList<label> >& facePointEdges,
 572
 573     const label startEdgeI,
 574     const label startVertI,
 575
 576     Map<label>& visited
 577 )
 578 {
 579     const pointField& points = eSurf.points();
 580     const edgeList& edges = eSurf.edges();
 581
 582     // Overestimate size of face
 583     face f(eSurf.faceEdges()[faceI].size());
 584
 585     label fp = 0;
 586
 587     label vertI = startVertI;
 588     label edgeI = startEdgeI;
 589
 590     while (true)
 591     {
 592         const edge& e = edges[edgeI];
 593
 594         if (debug & 2)
 595         {
 596             Pout<< "Now at:" << endl
 597                 << "    edge:" << edgeI << " vertices:" << e
 598                 << " positions:" << points[e.start()] << ' ' << points[e.end()]
 599                 << "    vertex:" << vertI << endl;
 600         }
 601
 602         // Mark edge as visited
 603         if (e[0] == vertI)
 604         {
 605             visited[edgeI] |= STARTTOEND;
 606         }
 607         else
 608         {
 609             visited[edgeI] |= ENDTOSTART;
 610         }
 611
 612
 613         // Store face vertex
 614         f[fp++] = vertI;
 615
 616         // step to other vertex
 617         vertI = e.otherVertex(vertI);
 618
 619         if (vertI == startVertI)
 620         {
 621             break;
 622         }
 623
 624         // step from vertex to next edge
 625         edgeI = nextEdge
 626         (
 627             eSurf,
 628             visited,
 629             faceI,
 630             n,
 631             facePointEdges,
 632             edgeI,
 633             vertI
 634         );
 635     }
 636
 637     f.setSize(fp);
 638
 639     return f;
 640 }
 641
 642
 643 void Foam::intersectedSurface::findNearestVisited
 644 (
 645     const edgeSurface& eSurf,
 646     const label faceI,
 647     const Map<DynamicList<label> >& facePointEdges,
 648     const Map<label>& pointVisited,
 649     const point& pt,
 650     const label excludePointI,
 651
 652     label& minVertI,
 653     scalar& minDist
 654 )
 655 {
 656     minVertI = -1;
 657     minDist = GREAT;
 658
 659     forAllConstIter(Map<label>, pointVisited, iter)
 660     {
 661         label pointI = iter.key();
 662
 663         if (pointI != excludePointI)
 664         {
 665             label nVisits = iter();
 666
 667             if (nVisits == 2*facePointEdges[pointI].size())
 668             {
 669                 // Fully visited (i.e. both sides of all edges)
 670                 scalar dist = mag(eSurf.points()[pointI] - pt);
 671
 672                 if (dist < minDist)
 673                 {
 674                     minDist = dist;
 675                     minVertI = pointI;
 676                 }
 677             }
 678         }
 679     }
 680
 681     if (minVertI == -1)
 682     {
 683         const labelList& fEdges = eSurf.faceEdges()[faceI];
 684
 685         SeriousErrorIn("intersectedSurface::findNearestVisited")
 686             << "Dumping face edges to faceEdges.obj" << endl;
 687
 688         writeLocalOBJ(eSurf.points(), eSurf.edges(), fEdges, "faceEdges.obj");
 689
 690         FatalErrorIn("intersectedSurface::findNearestVisited")
 691             << "No fully visited edge found for pt " << pt
 692             << abort(FatalError);
 693     }
 694 }
 695
 696
 697 // Sometimes there are bunches of edges that are not connected to the
 698 // other edges in the face. This routine tries to connect the loose
 699 // edges up to the 'proper' edges by adding two extra edges between a
 700 // properly connected edge and an unconnected one. Since at this level the
 701 // adressing is purely in form of points and a cloud of edges this can
 702 // be easily done.
 703 // (edges are to existing points. Don't want to introduce new vertices here
 704 // since then also neighbouring face would have to be split)
 705 Foam::faceList Foam::intersectedSurface::resplitFace
 706 (
 707     const triSurface& surf,
 708     const label faceI,
 709     const Map<DynamicList<label> >& facePointEdges,
 710     const Map<label>& visited,
 711     edgeSurface& eSurf
 712 )
 713 {
 714     {
 715         // Dump face for debugging.
 716         Pout<< "Writing face:" << faceI << " to face.obj" << endl;
 717         OFstream str("face.obj");
 718         writeOBJ(eSurf.points(), eSurf.edges(), eSurf.faceEdges()[faceI], str);
 719     }
 720
 721
 722     // Count the number of times point has been visited so we
 723     // can compare number to facePointEdges.
 724     Map<label> pointVisited(2*facePointEdges.size());
 725
 726
 727     {
 728         OFstream str0("visitedNone.obj");
 729         OFstream str1("visitedOnce.obj");
 730         OFstream str2("visitedTwice.obj");
 731         forAll(eSurf.points(), pointI)
 732         {
 733             const point& pt = eSurf.points()[pointI];
 734
 735             str0 << "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << nl;
 736             str1 << "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << nl;
 737             str2 << "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << nl;
 738         }
 739
 740
 741     forAllConstIter(Map<label>, visited, iter)
 742     {
 743         label edgeI = iter.key();
 744
 745         const edge& e = eSurf.edges()[edgeI];
 746
 747         label stat = iter();
 748
 749         if (stat == STARTTOEND || stat == ENDTOSTART)
 750         {
 751             incCount(pointVisited, e[0], 1);
 752             incCount(pointVisited, e[1], 1);
 753
 754             str1 << "l " << e[0]+1 << ' ' << e[1]+1 << nl;
 755         }
 756         else if (stat == BOTH)
 757         {
 758             incCount(pointVisited, e[0], 2);
 759             incCount(pointVisited, e[1], 2);
 760
 761             str2 << "l " << e[0]+1 << ' ' << e[1]+1 << nl;
 762         }
 763         else if (stat == UNVISITED)
 764         {
 765             incCount(pointVisited, e[0], 0);
 766             incCount(pointVisited, e[1], 0);
 767
 768             str0 << "l " << e[0]+1 << ' ' << e[1]+1 << nl;
 769         }
 770     }
 771     }
 772
 773
 774     {
 775         forAllConstIter(Map<label>, pointVisited, iter)
 776         {
 777             label pointI = iter.key();
 778
 779             label nVisits = iter();
 780
 781             Pout<< "point:" << pointI << "  nVisited:" << nVisits
 782                 << "  pointEdges:" << facePointEdges[pointI].size() << endl;
 783         }
 784     }
 785
 786
 787     // Find nearest point pair where one is not fully visited and
 788     // the other is.
 789
 790     label visitedVert0 = -1;
 791     label unvisitedVert0 = -1;
 792
 793     {
 794         scalar minDist = GREAT;
 795
 796         forAllConstIter(Map<DynamicList<label> >, facePointEdges, iter)
 797         {
 798             label pointI = iter.key();
 799
 800             label nVisits = pointVisited[pointI];
 801
 802             const DynamicList<label>& pEdges = iter();
 803
 804             if (nVisits < 2*pEdges.size())
 805             {
 806                 // Not fully visited. Find nearest fully visited.
 807
 808                 scalar nearDist;
 809                 label nearVertI;
 810
 811                 findNearestVisited
 812                 (
 813                     eSurf,
 814                     faceI,
 815                     facePointEdges,
 816                     pointVisited,
 817                     eSurf.points()[pointI],
 818                     -1,                         // Do not exclude vertex
 819                     nearVertI,
 820                     nearDist
 821                 );
 822
 823
 824                 if (nearDist < minDist)
 825                 {
 826                     minDist = nearDist;
 827                     visitedVert0 = nearVertI;
 828                     unvisitedVert0 = pointI;
 829                 }
 830             }
 831         }
 832     }
 833
 834
 835     // Find second intersection.
 836     label visitedVert1 = -1;
 837     label unvisitedVert1 = -1;
 838     {
 839         scalar minDist = GREAT;
 840
 841         forAllConstIter(Map<DynamicList<label> >, facePointEdges, iter)
 842         {
 843             label pointI = iter.key();
 844
 845             if (pointI != unvisitedVert0)
 846             {
 847                 label nVisits = pointVisited[pointI];
 848
 849                 const DynamicList<label>& pEdges = iter();
 850
 851                 if (nVisits < 2*pEdges.size())
 852                 {
 853                     // Not fully visited. Find nearest fully visited.
 854
 855                     scalar nearDist;
 856                     label nearVertI;
 857
 858                     findNearestVisited
 859                     (
 860                         eSurf,
 861                         faceI,
 862                         facePointEdges,
 863                         pointVisited,
 864                         eSurf.points()[pointI],
 865                         visitedVert0,           // vertex to exclude
 866                         nearVertI,
 867                         nearDist
 868                     );
 869
 870
 871                     if (nearDist < minDist)
 872                     {
 873                         minDist = nearDist;
 874                         visitedVert1 = nearVertI;
 875                         unvisitedVert1 = pointI;
 876                     }
 877                 }
 878             }
 879         }
 880     }
 881
 882
 883     Pout<< "resplitFace : adding intersection from " << visitedVert0
 884         << " to " << unvisitedVert0 << endl
 885         << " and from " << visitedVert1
 886         << " to " << unvisitedVert1 << endl;
 887
 888
 889 //    // Add the new intersection edges to the edgeSurface
 890 //    edgeList additionalEdges(2);
 891 //    additionalEdges[0] = edge(visitedVert0, unvisitedVert0);
 892 //    additionalEdges[1] = edge(visitedVert1, unvisitedVert1);
 893
 894     // Add the new intersection edges to the edgeSurface
 895     edgeList additionalEdges(1);
 896     additionalEdges[0] = edge(visitedVert0, unvisitedVert0);
 897
 898     eSurf.addIntersectionEdges(faceI, additionalEdges);
 899
 900     fileName newFName("face_" + Foam::name(faceI) + "_newEdges.obj");
 901     Pout<< "Dumping face:" << faceI << " to " << newFName << endl;
 902     writeLocalOBJ
 903     (
 904         eSurf.points(),
 905         eSurf.edges(),
 906         eSurf.faceEdges()[faceI],
 907         newFName
 908     );
 909
 910     // Retry splitFace. Use recursion since is rare situation.
 911     return splitFace(surf, faceI, eSurf);
 912 }
 913
 914
 915 Foam::faceList Foam::intersectedSurface::splitFace
 916 (
 917     const triSurface& surf,
 918     const label faceI,
 919     edgeSurface& eSurf
 920 )
 921 {
 922     // Alias
 923     const pointField& points = eSurf.points();
 924     const edgeList& edges = eSurf.edges();
 925     const labelList& fEdges = eSurf.faceEdges()[faceI];
 926
 927     // Create local (for the face only) point-edge connectivity.
 928     Map<DynamicList<label> > facePointEdges
 929     (
 930         calcPointEdgeAddressing
 931         (
 932             eSurf,
 933             faceI
 934         )
 935     );
 936
 937     // Order in which edges have been walked. Initialize outside edges.
 938     Map<label> visited(fEdges.size()*2);
 939
 940     forAll(fEdges, i)
 941     {
 942         label edgeI = fEdges[i];
 943
 944         if (eSurf.isSurfaceEdge(edgeI))
 945         {
 946             // Edge is edge from original surface so an outside edge for
 947             // the current face.
 948             label surfEdgeI = eSurf.parentEdge(edgeI);
 949
 950             label owner = surf.edgeOwner()[surfEdgeI];
 951
 952             if
 953             (
 954                 owner == faceI
 955              || sameEdgeOrder
 956                 (
 957                     surf.localFaces()[owner],
 958                     surf.localFaces()[faceI]
 959                 )
 960             )
 961             {
 962                 // Edge is in same order as current face.
 963                 // Mark off the opposite order.
 964                 visited.insert(edgeI, ENDTOSTART);
 965             }
 966             else
 967             {
 968                 // Edge is in same order as current face.
 969                 // Mark off the opposite order.
 970                 visited.insert(edgeI, STARTTOEND);
 971             }
 972         }
 973         else
 974         {
 975             visited.insert(edgeI, UNVISITED);
 976         }
 977     }
 978
 979
 980
 981     // Storage for faces
 982     DynamicList<face> faces(fEdges.size());
 983
 984     while (true)
 985     {
 986         // Find starting edge:
 987         // - unvisted triangle edge
 988         // - once visited intersection edge
 989         // Give priority to triangle edges.
 990         label startEdgeI = -1;
 991         label startVertI = -1;
 992
 993         forAll(fEdges, i)
 994         {
 995             label edgeI = fEdges[i];
 996
 997             const edge& e = edges[edgeI];
 998
 999             label stat = visited[edgeI];
1000
1001             if (stat == STARTTOEND)
1002             {
1003                 startEdgeI = edgeI;
1004                 startVertI = e[1];
1005
1006                 if (eSurf.isSurfaceEdge(edgeI))
1007                 {
1008                     break;
1009                 }
1010             }
1011             else if (stat == ENDTOSTART)
1012             {
1013                 startEdgeI = edgeI;
1014                 startVertI = e[0];
1015
1016                 if (eSurf.isSurfaceEdge(edgeI))
1017                 {
1018                     break;
1019                 }
1020             }
1021         }
1022
1023         if (startEdgeI == -1)
1024         {
1025             break;
1026         }
1027
1028         //Pout<< "splitFace: starting walk from edge:" << startEdgeI
1029         //    << ' ' << edges[startEdgeI] << " vertex:" << startVertI << endl;
1030
1031         //// Print current visited state.
1032         //printVisit(eSurf.edges(), fEdges, visited);
1033
1034         //{
1035         //    Pout<< "Writing face:" << faceI << " to face.obj" << endl;
1036         //    OFstream str("face.obj");
1037         //    writeOBJ(eSurf.points(), eSurf.edges(), fEdges, str);
1038         //}
1039
1040         faces.append
1041         (
1042             walkFace
1043             (
1044                 eSurf,
1045                 faceI,
1046                 surf.faceNormals()[faceI],
1047                 facePointEdges,
1048
1049                 startEdgeI,
1050                 startVertI,
1051
1052                 visited
1053             )
1054         );
1055     }
1056
1057
1058     // Check if any unvisited edges left.
1059     forAll(fEdges, i)
1060     {
1061         label edgeI = fEdges[i];
1062
1063         label stat = visited[edgeI];
1064
1065         if (eSurf.isSurfaceEdge(edgeI) && stat != BOTH)
1066         {
1067             SeriousErrorIn("Foam::intersectedSurface::splitFace")
1068                 << "Dumping face edges to faceEdges.obj" << endl;
1069
1070             writeLocalOBJ(points, edges, fEdges, "faceEdges.obj");
1071
1072             FatalErrorIn("intersectedSurface::splitFace")
1073                << "Problem: edge " << edgeI << " vertices "
1074                 << edges[edgeI] << " on face " << faceI
1075                 << " has visited status " << stat << " from a "
1076                  << "righthanded walk along all"
1077                 << " of the triangle edges. Are the original surfaces"
1078                 << " closed and non-intersecting?"
1079                 << abort(FatalError);
1080         }
1081         else if (stat != BOTH)
1082         {
1083             //{
1084             //    Pout<< "Dumping faces so far to faces.obj" << nl
1085             //        << faces << endl;
1086             //
1087             //    OFstream str("faces.obj");
1088             //
1089             //    forAll(faces, i)
1090             //    {
1091             //        writeOBJ(points, faces[i], str);
1092             //    }
1093             //}
1094
1095             Pout<< "** Resplitting **" << endl;
1096
1097             // Redo face after introducing extra edge. Edge introduced
1098             // should be one nearest to any fully visited edge.
1099             return resplitFace
1100             (
1101                 surf,
1102                 faceI,
1103                 facePointEdges,
1104                 visited,
1105                 eSurf
1106             );
1107         }
1108     }
1109
1110
1111     // See if normal needs flipping.
1112     faces.shrink();
1113
1114     vector n = faces[0].normal(eSurf.points());
1115
1116     if ((n & surf.faceNormals()[faceI]) < 0)
1117     {
1118         forAll(faces, i)
1119         {
1120             reverse(faces[i]);
1121         }
1122     }
1123
1124     return faces;
1125 }
1126
1127
1128 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
1129
1130 // Null constructor
1131 Foam::intersectedSurface::intersectedSurface()
1132 :
1133     triSurface(),
1134     intersectionEdges_(0),
1135     faceMap_(0),
1136     nSurfacePoints_(0)
1137 {}
1138
1139
1140 // Construct from components
1141 Foam::intersectedSurface::intersectedSurface(const triSurface& surf)
1142 :
1143     triSurface(surf),
1144     intersectionEdges_(0),
1145     faceMap_(0),
1146     nSurfacePoints_(0)
1147 {}
1148
1149
1150 // Construct from surface and intersection
1151 Foam::intersectedSurface::intersectedSurface
1152 (
1153     const triSurface& surf,
1154     const bool isFirstSurface,
1155     const surfaceIntersection& inter
1156 )
1157 :
1158     triSurface(),
1159     intersectionEdges_(0),
1160     faceMap_(0),
1161     nSurfacePoints_(surf.nPoints())
1162 {
1163     if (inter.cutPoints().empty() && inter.cutEdges().empty())
1164     {
1165         // No intersection. Make straight copy.
1166         triSurface::operator=(surf);
1167
1168         // Identity for face map
1169         faceMap_.setSize(size());
1170
1171         forAll(faceMap_, faceI)
1172         {
1173             faceMap_[faceI] = faceI;
1174         }
1175         return;
1176     }
1177
1178
1179     // Calculate face-edge addressing for surface + intersection.
1180     edgeSurface eSurf(surf, isFirstSurface, inter);
1181
1182     // Update point information for any removed points. (when are they removed?
1183     // - but better make sure)
1184     nSurfacePoints_ = eSurf.nSurfacePoints();
1185
1186     // Now we have full points, edges and edge addressing for surf. Start
1187     // extracting faces and triangulate them.
1188
1189     // Storage for new triangles of surface.
1190     DynamicList<labelledTri> newTris(eSurf.edges().size()/2);
1191
1192     // Start in newTris for decomposed face.
1193     labelList startTriI(surf.size(), 0);
1194
1195     forAll(surf, faceI)
1196     {
1197         startTriI[faceI] = newTris.size();
1198
1199         if (eSurf.faceEdges()[faceI].size() != surf.faceEdges()[faceI].size())
1200         {
1201             // Face has been cut by intersection.
1202             // Cut face into multiple subfaces. Use faceEdge information
1203             // from edgeSurface only. (original triSurface 'surf' is used for
1204             // faceNormals and regionnumber only)
1205             faceList newFaces
1206             (
1207                 splitFace
1208                 (
1209                     surf,
1210                     faceI,              // current triangle
1211                     eSurf               // face-edge description of surface
1212                                         // + intersection
1213                 )
1214             );
1215             forAll(newFaces, newFaceI)
1216             {
1217                 const face& newF = newFaces[newFaceI];
1218
1219 //                {
1220 //                    fileName fName
1221 //                    (
1222 //                        "face_"
1223 //                      + Foam::name(faceI)
1224 //                      + "_subFace_"
1225 //                      + Foam::name(newFaceI)
1226 //                      + ".obj"
1227 //                    );
1228 //                    Pout<< "Writing original face:" << faceI << " subFace:"
1229 //                        << newFaceI << " to " << fName << endl;
1230 //
1231 //                    OFstream str(fName);
1232 //
1233 //                    forAll(newF, fp)
1234 //                    {
1235 //                        meshTools::writeOBJ(str, eSurf.points()[newF[fp]]);
1236 //                    }
1237 //                    str << 'l';
1238 //                    forAll(newF, fp)
1239 //                    {
1240 //                        str << ' ' << fp+1;
1241 //                    }
1242 //                    str<< " 1" << nl;
1243 //                }
1244
1245
1246                 const vector& n = surf.faceNormals()[faceI];
1247                 const label region = surf[faceI].region();
1248
1249                 faceTriangulation tris(eSurf.points(), newF, n);
1250
1251                 forAll(tris, triI)
1252                 {
1253                     const triFace& t = tris[triI];
1254
1255                     forAll(t, i)
1256                     {
1257                         if (t[i] < 0 || t[i] >= eSurf.points().size())
1258                         {
1259                             FatalErrorIn
1260                             (
1261                                 "intersectedSurface::intersectedSurface"
1262                             )   << "Face triangulation of face " << faceI
1263                                 << " uses points outside range 0.."
1264                                 << eSurf.points().size()-1 << endl
1265                                 << "Triangulation:"
1266                                 << tris << abort(FatalError);
1267                         }
1268                     }
1269
1270                     newTris.append(labelledTri(t[0], t[1], t[2], region));
1271                 }
1272             }
1273         }
1274         else
1275         {
1276             // Face has not been cut at all. No need to renumber vertices since
1277             // eSurf keeps surface vertices first.
1278             newTris.append(surf.localFaces()[faceI]);
1279         }
1280     }
1281
1282     newTris.shrink();
1283
1284
1285     // Construct triSurface. Note that addressing will be compact since
1286     // edgeSurface is compact.
1287     triSurface::operator=
1288     (
1289         triSurface
1290         (
1291             newTris,
1292             surf.patches(),
1293             eSurf.points()
1294         )
1295     );
1296
1297     // Construct mapping back into original surface
1298     faceMap_.setSize(size());
1299
1300     for (label faceI = 0; faceI < surf.size()-1; faceI++)
1301     {
1302         for (label triI = startTriI[faceI]; triI < startTriI[faceI+1]; triI++)
1303         {
1304             faceMap_[triI] = faceI;
1305         }
1306     }
1307     for (label triI = startTriI[surf.size()-1]; triI < size(); triI++)
1308     {
1309         faceMap_[triI] = surf.size()-1;
1310     }
1311
1312
1313     // Find edges on *this which originate from 'cuts'. (i.e. newEdgeI >=
1314     // nSurfaceEdges) Renumber edges into local triSurface numbering.
1315
1316     intersectionEdges_.setSize(eSurf.edges().size() - eSurf.nSurfaceEdges());
1317
1318     label intersectionEdgeI = 0;
1319
1320     for
1321     (
1322         label edgeI = eSurf.nSurfaceEdges();
1323         edgeI < eSurf.edges().size();
1324         edgeI++
1325     )
1326     {
1327         // Get edge vertices in triSurface local numbering
1328         const edge& e = eSurf.edges()[edgeI];
1329         label surfStartI = meshPointMap()[e.start()];
1330         label surfEndI = meshPointMap()[e.end()];
1331
1332         // Find edge connected to surfStartI which also uses surfEndI.
1333         label surfEdgeI = -1;
1334
1335         const labelList& pEdges = pointEdges()[surfStartI];
1336
1337         forAll(pEdges, i)
1338         {
1339             const edge& surfE = edges()[pEdges[i]];
1340
1341             // Edge already connected to surfStart for sure. See if also
1342             // connects to surfEnd
1343             if (surfE.start() == surfEndI || surfE.end() == surfEndI)
1344             {
1345                 surfEdgeI = pEdges[i];
1346
1347                 break;
1348             }
1349         }
1350
1351         if (surfEdgeI != -1)
1352         {
1353             intersectionEdges_[intersectionEdgeI++] = surfEdgeI;
1354         }
1355         else
1356         {
1357             FatalErrorIn("intersectedSurface::intersectedSurface")
1358                 << "Cannot find edge among candidates " << pEdges
1359                 << " which uses points " << surfStartI
1360                 << " and " << surfEndI
1361                 << abort(FatalError);
1362         }
1363     }
1364 }
1365
1366
1367 // ************************************************************************* //