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