src/dynamicMesh/meshCut/cellCuts/cellCuts.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2004-2010 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
  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 "cellCuts.H"
  27 #include "polyMesh.H"
  28 #include "Time.H"
  29 #include "ListOps.H"
  30 #include "cellLooper.H"
  31 #include "refineCell.H"
  32 #include "meshTools.H"
  33 #include "geomCellLooper.H"
  34 #include "OFstream.H"
  35 #include "plane.H"
  36
  37 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  38
  39 defineTypeNameAndDebug(Foam::cellCuts, 0);
  40
  41
  42 // * * * * * * * * * * * * * Private Static Functions  * * * * * * * * * * * //
  43
  44 // Find val in first nElems elements of list.
  45 Foam::label Foam::cellCuts::findPartIndex
  46 (
  47     const labelList& elems,
  48     const label nElems,
  49     const label val
  50 )
  51 {
  52     for (label i = 0; i < nElems; i++)
  53     {
  54         if (elems[i] == val)
  55         {
  56             return i;
  57         }
  58     }
  59     return -1;
  60 }
  61
  62
  63 Foam::boolList Foam::cellCuts::expand
  64 (
  65     const label size,
  66     const labelList& labels
  67 )
  68 {
  69     boolList result(size, false);
  70
  71     forAll(labels, labelI)
  72     {
  73         result[labels[labelI]] = true;
  74     }
  75     return result;
  76 }
  77
  78
  79 Foam::scalarField Foam::cellCuts::expand
  80 (
  81     const label size,
  82     const labelList& labels,
  83     const scalarField& weights
  84 )
  85 {
  86     scalarField result(size, -GREAT);
  87
  88     forAll(labels, labelI)
  89     {
  90         result[labels[labelI]] = weights[labelI];
  91     }
  92     return result;
  93 }
  94
  95
  96 // Find first point in lst not in map.
  97 Foam::label Foam::cellCuts::firstUnique
  98 (
  99     const labelList& lst,
 100     const Map<label>& map
 101 )
 102 {
 103     forAll(lst, i)
 104     {
 105         if (!map.found(lst[i]))
 106         {
 107             return i;
 108         }
 109     }
 110     return -1;
 111 }
 112
 113
 114 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 115
 116 // Write cell and raw cuts on any of the elements
 117 void Foam::cellCuts::writeUncutOBJ
 118 (
 119     const fileName& dir,
 120     const label cellI
 121 ) const
 122 {
 123     //- Cell edges
 124     OFstream cutsStream(dir / "cell_" + name(cellI) + ".obj");
 125
 126     Pout<< "Writing cell for time " <<  mesh().time().timeName()
 127         << " to " << cutsStream.name() << nl;
 128
 129     meshTools::writeOBJ
 130     (
 131         cutsStream,
 132         mesh().cells(),
 133         mesh().faces(),
 134         mesh().points(),
 135         labelList(1, cellI)
 136     );
 137
 138     //- Loop cutting cell in two
 139     OFstream cutStream(dir / "cellCuts_" + name(cellI) + ".obj");
 140
 141     Pout<< "Writing raw cuts on cell for time " <<  mesh().time().timeName()
 142         << " to " << cutStream.name() << nl;
 143
 144     const labelList& cPoints = mesh().cellPoints()[cellI];
 145
 146     forAll(cPoints, i)
 147     {
 148         label pointI = cPoints[i];
 149         if (pointIsCut_[pointI])
 150         {
 151             meshTools::writeOBJ(cutStream, mesh().points()[pointI]);
 152         }
 153     }
 154
 155     const pointField& pts = mesh().points();
 156
 157     const labelList& cEdges = mesh().cellEdges()[cellI];
 158
 159     forAll(cEdges, i)
 160     {
 161         label edgeI = cEdges[i];
 162
 163         if (edgeIsCut_[edgeI])
 164         {
 165             const edge& e = mesh().edges()[edgeI];
 166
 167             const scalar w = edgeWeight_[edgeI];
 168
 169             meshTools::writeOBJ(cutStream, w*pts[e[1]] + (1-w)*pts[e[0]]);
 170         }
 171     }
 172 }
 173
 174
 175 void Foam::cellCuts::writeOBJ
 176 (
 177     const fileName& dir,
 178     const label cellI,
 179     const pointField& loopPoints,
 180     const labelList& anchors
 181 ) const
 182 {
 183     //- Cell edges
 184     OFstream cutsStream(dir / "cell_" + name(cellI) + ".obj");
 185
 186     Pout<< "Writing cell for time " <<  mesh().time().timeName()
 187         << " to " << cutsStream.name() << nl;
 188
 189     meshTools::writeOBJ
 190     (
 191         cutsStream,
 192         mesh().cells(),
 193         mesh().faces(),
 194         mesh().points(),
 195         labelList(1, cellI)
 196     );
 197
 198
 199     //- Loop cutting cell in two
 200     OFstream loopStream(dir / "cellLoop_" + name(cellI) + ".obj");
 201
 202     Pout<< "Writing loop for time " <<  mesh().time().timeName()
 203         << " to " << loopStream.name() << nl;
 204
 205     label vertI = 0;
 206
 207     writeOBJ(loopStream, loopPoints, vertI);
 208
 209
 210     //- Anchors for cell
 211     OFstream anchorStream(dir / "anchors_" + name(cellI) + ".obj");
 212
 213     Pout<< "Writing anchors for time " <<  mesh().time().timeName()
 214         << " to " << anchorStream.name() << endl;
 215
 216     forAll(anchors, i)
 217     {
 218         meshTools::writeOBJ(anchorStream, mesh().points()[anchors[i]]);
 219     }
 220 }
 221
 222
 223 // Find face on cell using the two edges.
 224 Foam::label Foam::cellCuts::edgeEdgeToFace
 225 (
 226     const label cellI,
 227     const label edgeA,
 228     const label edgeB
 229 ) const
 230 {
 231     const labelList& cFaces = mesh().cells()[cellI];
 232
 233     forAll(cFaces, cFaceI)
 234     {
 235         label faceI = cFaces[cFaceI];
 236
 237         const labelList& fEdges = mesh().faceEdges()[faceI];
 238
 239         if
 240         (
 241             findIndex(fEdges, edgeA) != -1
 242          && findIndex(fEdges, edgeB) != -1
 243         )
 244         {
 245            return faceI;
 246         }
 247     }
 248
 249     // Coming here means the loop is illegal since the two edges
 250     // are not shared by a face. We just mark loop as invalid and continue.
 251
 252     WarningIn
 253     (
 254         "Foam::cellCuts::edgeEdgeToFace"
 255         "(const label cellI, const label edgeA,"
 256         "const label edgeB) const"
 257     )   << "cellCuts : Cannot find face on cell "
 258         << cellI << " that has both edges " << edgeA << ' ' << edgeB << endl
 259         << "faces : " << cFaces << endl
 260         << "edgeA : " << mesh().edges()[edgeA] << endl
 261         << "edgeB : " << mesh().edges()[edgeB] << endl
 262         << "Marking the loop across this cell as invalid" << endl;
 263
 264     return -1;
 265 }
 266
 267
 268 // Find face on cell using an edge and a vertex.
 269 Foam::label Foam::cellCuts::edgeVertexToFace
 270 (
 271     const label cellI,
 272     const label edgeI,
 273     const label vertI
 274 ) const
 275 {
 276     const labelList& cFaces = mesh().cells()[cellI];
 277
 278     forAll(cFaces, cFaceI)
 279     {
 280         label faceI = cFaces[cFaceI];
 281
 282         const face& f = mesh().faces()[faceI];
 283
 284         const labelList& fEdges = mesh().faceEdges()[faceI];
 285
 286         if
 287         (
 288             findIndex(fEdges, edgeI) != -1
 289          && findIndex(f, vertI) != -1
 290         )
 291         {
 292            return faceI;
 293         }
 294     }
 295
 296     WarningIn
 297     (
 298         "Foam::cellCuts::edgeVertexToFace"
 299         "(const label cellI, const label edgeI, "
 300         "const label vertI) const"
 301     )   << "cellCuts : Cannot find face on cell "
 302         << cellI << " that has both edge " << edgeI << " and vertex "
 303         << vertI << endl
 304         << "faces : " << cFaces << endl
 305         << "edge : " << mesh().edges()[edgeI] << endl
 306         << "Marking the loop across this cell as invalid" << endl;
 307
 308     return -1;
 309 }
 310
 311
 312 // Find face using two vertices (guaranteed not to be along edge)
 313 Foam::label Foam::cellCuts::vertexVertexToFace
 314 (
 315     const label cellI,
 316     const label vertA,
 317     const label vertB
 318 ) const
 319 {
 320     const labelList& cFaces = mesh().cells()[cellI];
 321
 322     forAll(cFaces, cFaceI)
 323     {
 324         label faceI = cFaces[cFaceI];
 325
 326         const face& f = mesh().faces()[faceI];
 327
 328         if
 329         (
 330             findIndex(f, vertA) != -1
 331          && findIndex(f, vertB) != -1
 332         )
 333         {
 334            return faceI;
 335         }
 336     }
 337
 338     WarningIn("Foam::cellCuts::vertexVertexToFace")
 339         << "cellCuts : Cannot find face on cell "
 340         << cellI << " that has vertex " << vertA << " and vertex "
 341         << vertB << endl
 342         << "faces : " << cFaces << endl
 343         << "Marking the loop across this cell as invalid" << endl;
 344
 345     return -1;
 346 }
 347
 348
 349 void Foam::cellCuts::calcFaceCuts() const
 350 {
 351     if (faceCutsPtr_)
 352     {
 353         FatalErrorIn("cellCuts::calcFaceCuts()")
 354             << "faceCuts already calculated" << abort(FatalError);
 355     }
 356
 357     const faceList& faces = mesh().faces();
 358
 359
 360     faceCutsPtr_ = new labelListList(mesh().nFaces());
 361     labelListList& faceCuts = *faceCutsPtr_;
 362
 363     for (label faceI = 0; faceI < mesh().nFaces(); faceI++)
 364     {
 365         const face& f = faces[faceI];
 366
 367         // Big enough storage (possibly all points and all edges cut). Shrink
 368         // later on.
 369         labelList& cuts = faceCuts[faceI];
 370
 371         cuts.setSize(2*f.size());
 372
 373         label cutI = 0;
 374
 375         // Do point-edge-point walk over face and collect all cuts.
 376         // Problem is that we want to start from one of the endpoints of a
 377         // string of connected cuts; we don't want to start somewhere in the
 378         // middle.
 379
 380         // Pass1: find first point cut not preceeded by a cut.
 381         label startFp = -1;
 382
 383         forAll(f, fp)
 384         {
 385             label v0 = f[fp];
 386
 387             if (pointIsCut_[v0])
 388             {
 389                 label vMin1 = f[f.rcIndex(fp)];
 390
 391                 if
 392                 (
 393                     !pointIsCut_[vMin1]
 394                  && !edgeIsCut_[findEdge(faceI, v0, vMin1)]
 395                 )
 396                 {
 397                     cuts[cutI++] = vertToEVert(v0);
 398                     startFp = f.fcIndex(fp);
 399                     break;
 400                 }
 401             }
 402         }
 403
 404         // Pass2: first edge cut not preceeded by point cut
 405         if (startFp == -1)
 406         {
 407             forAll(f, fp)
 408             {
 409                 label fp1 = f.fcIndex(fp);
 410
 411                 label v0 = f[fp];
 412                 label v1 = f[fp1];
 413
 414                 label edgeI = findEdge(faceI, v0, v1);
 415
 416                 if (edgeIsCut_[edgeI] && !pointIsCut_[v0])
 417                 {
 418                     cuts[cutI++] = edgeToEVert(edgeI);
 419                     startFp = fp1;
 420                     break;
 421                 }
 422             }
 423         }
 424
 425         // Pass3: nothing found so far. Either face is not cut at all or all
 426         // vertices are cut. Start from 0.
 427         if (startFp == -1)
 428         {
 429             startFp = 0;
 430         }
 431
 432         // Store all cuts starting from startFp;
 433         label fp = startFp;
 434
 435         bool allVerticesCut = true;
 436
 437         forAll(f, i)
 438         {
 439             label fp1 = f.fcIndex(fp);
 440
 441             // Get the three items: current vertex, next vertex and edge
 442             // inbetween
 443             label v0 = f[fp];
 444             label v1 = f[fp1];
 445             label edgeI = findEdge(faceI, v0, v1);
 446
 447             if (pointIsCut_[v0])
 448             {
 449                 cuts[cutI++] = vertToEVert(v0);
 450             }
 451             else
 452             {
 453                 // For check if all vertices have been cut (= illegal)
 454                 allVerticesCut = false;
 455             }
 456
 457             if (edgeIsCut_[edgeI])
 458             {
 459                 cuts[cutI++] = edgeToEVert(edgeI);
 460             }
 461
 462             fp = fp1;
 463         }
 464
 465
 466         if (allVerticesCut)
 467         {
 468             WarningIn("Foam::cellCuts::calcFaceCuts() const")
 469                 << "Face " << faceI << " vertices " << f
 470                 << " has all its vertices cut. Not cutting face." << endl;
 471
 472             cutI = 0;
 473         }
 474
 475         // Remove duplicate starting point
 476         if (cutI > 1 && cuts[cutI-1] == cuts[0])
 477         {
 478             cutI--;
 479         }
 480         cuts.setSize(cutI);
 481     }
 482 }
 483
 484
 485 // Find edge on face using two vertices
 486 Foam::label Foam::cellCuts::findEdge
 487 (
 488     const label faceI,
 489     const label v0,
 490     const label v1
 491 ) const
 492 {
 493     const edgeList& edges = mesh().edges();
 494
 495     const labelList& fEdges = mesh().faceEdges()[faceI];
 496
 497     forAll(fEdges, i)
 498     {
 499         const edge& e = edges[fEdges[i]];
 500
 501         if
 502         (
 503             (e[0] == v0 && e[1] == v1)
 504          || (e[0] == v1 && e[1] == v0)
 505         )
 506         {
 507             return fEdges[i];
 508         }
 509     }
 510
 511     return -1;
 512 }
 513
 514
 515 // Check if there is a face on the cell on which all cuts are.
 516 Foam::label Foam::cellCuts::loopFace
 517 (
 518     const label cellI,
 519     const labelList& loop
 520 ) const
 521 {
 522     const cell& cFaces = mesh().cells()[cellI];
 523
 524     forAll(cFaces, cFaceI)
 525     {
 526         label faceI = cFaces[cFaceI];
 527
 528         const labelList& fEdges = mesh().faceEdges()[faceI];
 529         const face& f = mesh().faces()[faceI];
 530
 531         bool allOnFace = true;
 532
 533         forAll(loop, i)
 534         {
 535             label cut = loop[i];
 536
 537             if (isEdge(cut))
 538             {
 539                 if (findIndex(fEdges, getEdge(cut)) == -1)
 540                 {
 541                     // Edge not on face. Skip face.
 542                     allOnFace = false;
 543                     break;
 544                 }
 545             }
 546             else
 547             {
 548                 if (findIndex(f, getVertex(cut)) == -1)
 549                 {
 550                     // Vertex not on face. Skip face.
 551                     allOnFace = false;
 552                     break;
 553                 }
 554             }
 555         }
 556
 557         if (allOnFace)
 558         {
 559             // Found face where all elements of loop are on the face.
 560             return faceI;
 561         }
 562     }
 563     return -1;
 564 }
 565
 566
 567 // From point go into connected face
 568 bool Foam::cellCuts::walkPoint
 569 (
 570     const label cellI,
 571     const label startCut,
 572
 573     const label exclude0,
 574     const label exclude1,
 575
 576     const label otherCut,
 577
 578     label& nVisited,
 579     labelList& visited
 580 ) const
 581 {
 582     label vertI = getVertex(otherCut);
 583
 584     const labelList& pFaces = mesh().pointFaces()[vertI];
 585
 586     forAll(pFaces, pFaceI)
 587     {
 588         label otherFaceI = pFaces[pFaceI];
 589
 590         if
 591         (
 592             otherFaceI != exclude0
 593          && otherFaceI != exclude1
 594          && meshTools::faceOnCell(mesh(), cellI, otherFaceI)
 595         )
 596         {
 597             label oldNVisited = nVisited;
 598
 599             bool foundLoop =
 600                 walkCell
 601                 (
 602                     cellI,
 603                     startCut,
 604                     otherFaceI,
 605                     otherCut,
 606                     nVisited,
 607                     visited
 608                 );
 609
 610             if (foundLoop)
 611             {
 612                 return true;
 613             }
 614
 615             // No success. Restore state and continue
 616             nVisited = oldNVisited;
 617         }
 618     }
 619     return false;
 620 }
 621
 622
 623 // Cross cut (which is edge on faceI) onto next face
 624 bool Foam::cellCuts::crossEdge
 625 (
 626     const label cellI,
 627     const label startCut,
 628     const label faceI,
 629     const label otherCut,
 630
 631     label& nVisited,
 632     labelList& visited
 633 ) const
 634 {
 635     // Cross edge to other face
 636     label edgeI = getEdge(otherCut);
 637
 638     label otherFaceI = meshTools::otherFace(mesh(), cellI, faceI, edgeI);
 639
 640     // Store old state
 641     label oldNVisited = nVisited;
 642
 643     // Recurse to otherCut
 644     bool foundLoop =
 645         walkCell
 646         (
 647             cellI,
 648             startCut,
 649             otherFaceI,
 650             otherCut,
 651             nVisited,
 652             visited
 653         );
 654
 655     if (foundLoop)
 656     {
 657         return true;
 658     }
 659     else
 660     {
 661         // No success. Restore state (i.e. backtrack)
 662         nVisited = oldNVisited;
 663
 664         return false;
 665     }
 666 }
 667
 668
 669 bool Foam::cellCuts::addCut
 670 (
 671     const label cellI,
 672     const label cut,
 673     label& nVisited,
 674     labelList& visited
 675 ) const
 676 {
 677     // Check for duplicate cuts.
 678     if (findPartIndex(visited, nVisited, cut) != -1)
 679     {
 680         // Truncate (copy of) visited for ease of printing.
 681         labelList truncVisited(visited);
 682         truncVisited.setSize(nVisited);
 683
 684         Pout<< "For cell " << cellI << " : trying to add duplicate cut " << cut;
 685         labelList cuts(1, cut);
 686         writeCuts(Pout, cuts, loopWeights(cuts));
 687
 688         Pout<< " to path:";
 689         writeCuts(Pout, truncVisited, loopWeights(truncVisited));
 690         Pout<< endl;
 691
 692         return false;
 693     }
 694     else
 695     {
 696         visited[nVisited++] = cut;
 697
 698         return true;
 699     }
 700 }
 701
 702
 703 // Walk across faceI, storing cuts as you go. Returns last two cuts visisted.
 704 // Returns true if valid walk.
 705 bool Foam::cellCuts::walkFace
 706 (
 707     const label cellI,
 708     const label startCut,
 709     const label faceI,
 710     const label cut,
 711
 712     label& lastCut,
 713     label& beforeLastCut,
 714     label& nVisited,
 715     labelList& visited
 716 ) const
 717 {
 718     const labelList& fCuts = faceCuts()[faceI];
 719
 720     if (fCuts.size() < 2)
 721     {
 722         return false;
 723     }
 724
 725     // Easy case : two cuts.
 726     if (fCuts.size() == 2)
 727     {
 728         if (fCuts[0] == cut)
 729         {
 730             if (!addCut(cellI, cut, nVisited, visited))
 731             {
 732                 return false;
 733             }
 734
 735             beforeLastCut = cut;
 736             lastCut = fCuts[1];
 737
 738             return true;
 739         }
 740         else
 741         {
 742             if (!addCut(cellI, cut, nVisited, visited))
 743             {
 744                 return false;
 745             }
 746
 747             beforeLastCut = cut;
 748             lastCut = fCuts[0];
 749
 750             return true;
 751         }
 752     }
 753
 754     // Harder case: more than 2 cuts on face.
 755     // Should be start or end of string of cuts. Store all but last cut.
 756     if (fCuts[0] == cut)
 757     {
 758         // Walk forward
 759         for (label i = 0; i < fCuts.size()-1; i++)
 760         {
 761             if (!addCut(cellI, fCuts[i], nVisited, visited))
 762             {
 763                 return false;
 764             }
 765         }
 766         beforeLastCut = fCuts[fCuts.size()-2];
 767         lastCut = fCuts[fCuts.size()-1];
 768     }
 769     else if (fCuts[fCuts.size()-1] == cut)
 770     {
 771         for (label i = fCuts.size()-1; i >= 1; --i)
 772         {
 773             if (!addCut(cellI, fCuts[i], nVisited, visited))
 774             {
 775                 return false;
 776             }
 777         }
 778         beforeLastCut = fCuts[1];
 779         lastCut = fCuts[0];
 780     }
 781     else
 782     {
 783         WarningIn("Foam::cellCuts::walkFace")
 784             << "In middle of cut. cell:" << cellI << " face:" << faceI
 785             << " cuts:" << fCuts << " current cut:" << cut << endl;
 786
 787         return false;
 788     }
 789
 790     return true;
 791 }
 792
 793
 794
 795 // Walk across cuts (cut edges or cut vertices) of cell. Stops when hit cut
 796 // already visited. Returns true when loop of 3 or more vertices found.
 797 bool Foam::cellCuts::walkCell
 798 (
 799     const label cellI,
 800     const label startCut,
 801     const label faceI,
 802     const label cut,
 803
 804     label& nVisited,
 805     labelList& visited
 806 ) const
 807 {
 808     // Walk across face, storing cuts. Return last two cuts visited.
 809     label lastCut = -1;
 810     label beforeLastCut = -1;
 811
 812
 813     if (debug & 2)
 814     {
 815         Pout<< "For cell:" << cellI << " walked across face " << faceI
 816             << " from cut ";
 817         labelList cuts(1, cut);
 818         writeCuts(Pout, cuts, loopWeights(cuts));
 819         Pout<< endl;
 820     }
 821
 822     bool validWalk = walkFace
 823     (
 824         cellI,
 825         startCut,
 826         faceI,
 827         cut,
 828
 829         lastCut,
 830         beforeLastCut,
 831         nVisited,
 832         visited
 833     );
 834
 835     if (!validWalk)
 836     {
 837         return false;
 838     }
 839
 840     if (debug & 2)
 841     {
 842         Pout<< "    to last cut ";
 843         labelList cuts(1, lastCut);
 844         writeCuts(Pout, cuts, loopWeights(cuts));
 845         Pout<< endl;
 846     }
 847
 848     // Check if starting point reached.
 849     if (lastCut == startCut)
 850     {
 851         if (nVisited >= 3)
 852         {
 853             if (debug & 2)
 854             {
 855                 // Truncate (copy of) visited for ease of printing.
 856                 labelList truncVisited(visited);
 857                 truncVisited.setSize(nVisited);
 858
 859                 Pout<< "For cell " << cellI << " : found closed path:";
 860                 writeCuts(Pout, truncVisited, loopWeights(truncVisited));
 861                 Pout<< " closed by " << lastCut << endl;
 862             }
 863
 864             return true;
 865         }
 866         else
 867         {
 868             // Loop but too short. Probably folded back on itself.
 869             return false;
 870         }
 871     }
 872
 873
 874     // Check last cut and one before that to determine type.
 875
 876     // From beforeLastCut to lastCut:
 877     //  - from edge to edge
 878     //      (- always not along existing edge)
 879     //  - from edge to vertex
 880     //      - not along existing edge
 881     //      - along existing edge: not allowed?
 882     //  - from vertex to edge
 883     //      - not along existing edge
 884     //      - along existing edge. Not allowed. See above.
 885     //  - from vertex to vertex
 886     //      - not along existing edge
 887     //      - along existing edge
 888
 889     if (isEdge(beforeLastCut))
 890     {
 891         if (isEdge(lastCut))
 892         {
 893             // beforeLastCut=edge, lastCut=edge.
 894
 895             // Cross lastCut (=edge) into face not faceI
 896             return crossEdge
 897             (
 898                 cellI,
 899                 startCut,
 900                 faceI,
 901                 lastCut,
 902                 nVisited,
 903                 visited
 904             );
 905         }
 906         else
 907         {
 908             // beforeLastCut=edge, lastCut=vertex.
 909
 910             // Go from lastCut to all connected faces (except faceI)
 911             return walkPoint
 912             (
 913                 cellI,
 914                 startCut,
 915                 faceI,          // exclude0: don't cross back on current face
 916                 -1,             // exclude1
 917                 lastCut,
 918                 nVisited,
 919                 visited
 920             );
 921         }
 922     }
 923     else
 924     {
 925         if (isEdge(lastCut))
 926         {
 927             // beforeLastCut=vertex, lastCut=edge.
 928             return crossEdge
 929             (
 930                 cellI,
 931                 startCut,
 932                 faceI,
 933                 lastCut,
 934                 nVisited,
 935                 visited
 936             );
 937         }
 938         else
 939         {
 940             // beforeLastCut=vertex, lastCut=vertex. Check if along existing
 941             // edge.
 942             label edgeI =
 943                 findEdge
 944                 (
 945                     faceI,
 946                     getVertex(beforeLastCut),
 947                     getVertex(lastCut)
 948                 );
 949
 950             if (edgeI != -1)
 951             {
 952                 // Cut along existing edge. So is in fact on two faces.
 953                 // Get faces on both sides of the edge to make
 954                 // sure we dont fold back on to those.
 955
 956                 label f0, f1;
 957                 meshTools::getEdgeFaces(mesh(), cellI, edgeI, f0, f1);
 958
 959                 return walkPoint
 960                 (
 961                     cellI,
 962                     startCut,
 963                     f0,
 964                     f1,
 965                     lastCut,
 966                     nVisited,
 967                     visited
 968                 );
 969
 970             }
 971             else
 972             {
 973                 // Cross cut across face.
 974                 return walkPoint
 975                 (
 976                     cellI,
 977                     startCut,
 978                     faceI,      // face to exclude
 979                     -1,         // face to exclude
 980                     lastCut,
 981                     nVisited,
 982                     visited
 983                 );
 984             }
 985         }
 986     }
 987 }
 988
 989
 990 // Determine for every cut cell the loop (= face) it is cut by. Done by starting
 991 // from a cut edge or cut vertex and walking across faces, from cut to cut,
 992 // until starting cut hit.
 993 // If multiple loops are possible across a cell circumference takes the first
 994 // one found.
 995 void Foam::cellCuts::calcCellLoops(const labelList& cutCells)
 996 {
 997     // Calculate cuts per face.
 998     const labelListList& allFaceCuts = faceCuts();
 999
1000     // Per cell the number of faces with valid cuts. Is used as quick
1001     // rejection to see if cell can be cut.
1002     labelList nCutFaces(mesh().nCells(), 0);
1003
1004     forAll(allFaceCuts, faceI)
1005     {
1006         const labelList& fCuts = allFaceCuts[faceI];
1007
1008         if (fCuts.size() == mesh().faces()[faceI].size())
1009         {
1010             // Too many cuts on face. WalkCell would get very upset so disable.
1011             nCutFaces[mesh().faceOwner()[faceI]] = labelMin;
1012
1013             if (mesh().isInternalFace(faceI))
1014             {
1015                 nCutFaces[mesh().faceNeighbour()[faceI]] = labelMin;
1016             }
1017         }
1018         else if (fCuts.size() >= 2)
1019         {
1020             // Could be valid cut. Update count for owner and neighbour.
1021             nCutFaces[mesh().faceOwner()[faceI]]++;
1022
1023             if (mesh().isInternalFace(faceI))
1024             {
1025                 nCutFaces[mesh().faceNeighbour()[faceI]]++;
1026             }
1027         }
1028     }
1029
1030
1031     // Stack of visited cuts (nVisited used as stack pointer)
1032     // Size big enough.
1033     labelList visited(mesh().nPoints());
1034
1035     forAll(cutCells, i)
1036     {
1037         label cellI = cutCells[i];
1038
1039         bool validLoop = false;
1040
1041         // Quick rejection: has enough faces that are cut?
1042         if (nCutFaces[cellI] >= 3)
1043         {
1044             const labelList& cFaces = mesh().cells()[cellI];
1045
1046             if (debug & 2)
1047             {
1048                 Pout<< "cell:" << cellI << " cut faces:" << endl;
1049                 forAll(cFaces, i)
1050                 {
1051                     label faceI = cFaces[i];
1052                     const labelList& fCuts = allFaceCuts[faceI];
1053
1054                     Pout<< "    face:" << faceI << " cuts:";
1055                     writeCuts(Pout, fCuts, loopWeights(fCuts));
1056                     Pout<< endl;
1057                 }
1058             }
1059
1060             label nVisited = 0;
1061
1062             // Determine the first cut face to start walking from.
1063             forAll(cFaces, cFaceI)
1064             {
1065                 label faceI = cFaces[cFaceI];
1066
1067                 const labelList& fCuts = allFaceCuts[faceI];
1068
1069                 // Take first or last cut of multiple on face.
1070                 // Note that in calcFaceCuts
1071                 // we have already made sure this is the start or end of a
1072                 // string of cuts.
1073                 if (fCuts.size() >= 2)
1074                 {
1075                     // Try walking from start of fCuts.
1076                     nVisited = 0;
1077
1078                     if (debug & 2)
1079                     {
1080                         Pout<< "cell:" << cellI
1081                             << " start walk at face:" << faceI
1082                             << " cut:";
1083                         labelList cuts(1, fCuts[0]);
1084                         writeCuts(Pout, cuts, loopWeights(cuts));
1085                         Pout<< endl;
1086                     }
1087
1088                     validLoop =
1089                         walkCell
1090                         (
1091                             cellI,
1092                             fCuts[0],
1093                             faceI,
1094                             fCuts[0],
1095
1096                             nVisited,
1097                             visited
1098                         );
1099
1100                     if (validLoop)
1101                     {
1102                         break;
1103                     }
1104
1105                     // No need to try and walk from end of cuts since
1106                     // all paths already tried by walkCell.
1107                 }
1108             }
1109
1110             if (validLoop)
1111             {
1112                 // Copy nVisited elements out of visited (since visited is
1113                 // never truncated for efficiency reasons)
1114
1115                 labelList& loop = cellLoops_[cellI];
1116
1117                 loop.setSize(nVisited);
1118
1119                 for (label i = 0; i < nVisited; i++)
1120                 {
1121                     loop[i] = visited[i];
1122                 }
1123             }
1124             else
1125             {
1126                 // Invalid loop. Leave cellLoops_[cellI] zero size which
1127                 // flags this.
1128                 Pout<< "calcCellLoops(const labelList&) : did not find valid"
1129                     << " loop for cell " << cellI << endl;
1130                 // Dump cell and cuts on cell.
1131                 writeUncutOBJ(".", cellI);
1132
1133                 cellLoops_[cellI].setSize(0);
1134             }
1135         }
1136         else
1137         {
1138             //Pout<< "calcCellLoops(const labelList&) : did not find valid"
1139             //    << " loop for cell " << cellI << " since not enough cut faces"
1140             //    << endl;
1141             cellLoops_[cellI].setSize(0);
1142         }
1143     }
1144 }
1145
1146
1147 // Walk unset edges of single cell from starting point and marks visited
1148 // edges and vertices with status.
1149 void Foam::cellCuts::walkEdges
1150 (
1151     const label cellI,
1152     const label pointI,
1153     const label status,
1154
1155     Map<label>& edgeStatus,
1156     Map<label>& pointStatus
1157 ) const
1158 {
1159     if (pointStatus.insert(pointI, status))
1160     {
1161         // First visit to pointI
1162
1163         const labelList& pEdges = mesh().pointEdges()[pointI];
1164
1165         forAll(pEdges, pEdgeI)
1166         {
1167             label edgeI = pEdges[pEdgeI];
1168
1169             if
1170             (
1171                 meshTools::edgeOnCell(mesh(), cellI, edgeI)
1172              && edgeStatus.insert(edgeI, status)
1173             )
1174             {
1175                 // First visit to edgeI so recurse.
1176
1177                 label v2 = mesh().edges()[edgeI].otherVertex(pointI);
1178
1179                 walkEdges(cellI, v2, status, edgeStatus, pointStatus);
1180             }
1181         }
1182     }
1183 }
1184
1185
1186 // Invert anchor point selection.
1187 Foam::labelList Foam::cellCuts::nonAnchorPoints
1188 (
1189     const labelList& cellPoints,
1190     const labelList& anchorPoints,
1191     const labelList& loop
1192 ) const
1193 {
1194     labelList newElems(cellPoints.size());
1195     label newElemI = 0;
1196
1197     forAll(cellPoints, i)
1198     {
1199         label pointI = cellPoints[i];
1200
1201         if
1202         (
1203             findIndex(anchorPoints, pointI) == -1
1204          && findIndex(loop, vertToEVert(pointI)) == -1
1205         )
1206         {
1207             newElems[newElemI++] = pointI;
1208         }
1209     }
1210
1211     newElems.setSize(newElemI);
1212
1213     return newElems;
1214 }
1215
1216
1217 //- Check anchor points on 'outside' of loop
1218 bool Foam::cellCuts::loopAnchorConsistent
1219 (
1220     const label cellI,
1221     const pointField& loopPts,
1222     const labelList& anchorPoints
1223 ) const
1224 {
1225     // Create identity face for ease of calculation of normal etc.
1226     face f(identity(loopPts.size()));
1227
1228     vector normal = f.normal(loopPts);
1229     point ctr = f.centre(loopPts);
1230
1231
1232     // Get average position of anchor points.
1233     vector avg(vector::zero);
1234
1235     forAll(anchorPoints, ptI)
1236     {
1237         avg += mesh().points()[anchorPoints[ptI]];
1238     }
1239     avg /= anchorPoints.size();
1240
1241
1242     if (((avg - ctr) & normal) > 0)
1243     {
1244         return true;
1245     }
1246     else
1247     {
1248         return false;
1249     }
1250 }
1251
1252
1253 // Determines set of anchor points given a loop. The loop should split the
1254 // cell into (one or) two sets of vertices. The set of vertices that is
1255 // on the 'normal' side of the loop is the anchor set.
1256 // Returns true if valid set, false otherwise.
1257 bool Foam::cellCuts::calcAnchors
1258 (
1259     const label cellI,
1260     const labelList& loop,
1261     const pointField& loopPts,
1262
1263     labelList& anchorPoints
1264 ) const
1265 {
1266     const edgeList& edges = mesh().edges();
1267
1268     const labelList& cPoints = mesh().cellPoints()[cellI];
1269     const labelList& cEdges = mesh().cellEdges()[cellI];
1270     const cell& cFaces = mesh().cells()[cellI];
1271
1272     // Points on loop
1273
1274     // Status of point:
1275     // 0 - on loop
1276     // 1 - point set 1
1277     // 2 - point set 2
1278     Map<label> pointStatus(2*cPoints.size());
1279     Map<label> edgeStatus(2*cEdges.size());
1280
1281     // Mark loop vertices
1282     forAll(loop, i)
1283     {
1284         label cut = loop[i];
1285
1286         if (isEdge(cut))
1287         {
1288             edgeStatus.insert(getEdge(cut), 0);
1289         }
1290         else
1291         {
1292             pointStatus.insert(getVertex(cut), 0);
1293         }
1294     }
1295     // Since edges between two cut vertices have not been marked, mark them
1296     // explicitly
1297     forAll(cEdges, i)
1298     {
1299         label edgeI = cEdges[i];
1300         const edge& e = edges[edgeI];
1301
1302         if (pointStatus.found(e[0]) && pointStatus.found(e[1]))
1303         {
1304             edgeStatus.insert(edgeI, 0);
1305         }
1306     }
1307
1308
1309     // Find uncut starting vertex
1310     label uncutIndex = firstUnique(cPoints, pointStatus);
1311
1312     if (uncutIndex == -1)
1313     {
1314         WarningIn("Foam::cellCuts::calcAnchors")
1315             << "Invalid loop " << loop << " for cell " << cellI << endl
1316             << "Can not find point on cell which is not cut by loop."
1317             << endl;
1318
1319         writeOBJ(".", cellI, loopPts, labelList(0));
1320
1321         return false;
1322     }
1323
1324     // Walk unset vertices and edges and mark with 1 in pointStatus, edgeStatus
1325     walkEdges(cellI, cPoints[uncutIndex], 1, edgeStatus, pointStatus);
1326
1327     // Find new uncut starting vertex
1328     uncutIndex = firstUnique(cPoints, pointStatus);
1329
1330     if (uncutIndex == -1)
1331     {
1332         // All vertices either in loop or in anchor. So split is along single
1333         // face.
1334         WarningIn("Foam::cellCuts::calcAnchors")
1335             << "Invalid loop " << loop << " for cell " << cellI << endl
1336             << "All vertices of cell are either in loop or in anchor set"
1337             << endl;
1338
1339         writeOBJ(".", cellI, loopPts, labelList(0));
1340
1341         return false;
1342     }
1343
1344     // Walk unset vertices and edges and mark with 2. These are the
1345     // pointset 2.
1346     walkEdges(cellI, cPoints[uncutIndex], 2, edgeStatus, pointStatus);
1347
1348     // Collect both sets in lists.
1349     DynamicList<label> connectedPoints(cPoints.size());
1350     DynamicList<label> otherPoints(cPoints.size());
1351
1352     forAllConstIter(Map<label>, pointStatus, iter)
1353     {
1354         if (iter() == 1)
1355         {
1356             connectedPoints.append(iter.key());
1357         }
1358         else if (iter() == 2)
1359         {
1360             otherPoints.append(iter.key());
1361         }
1362     }
1363     connectedPoints.shrink();
1364     otherPoints.shrink();
1365
1366     // Check that all points have been used.
1367     uncutIndex = firstUnique(cPoints, pointStatus);
1368
1369     if (uncutIndex != -1)
1370     {
1371         WarningIn("Foam::cellCuts::calcAnchors")
1372             << "Invalid loop " << loop << " for cell " << cellI
1373             << " since it splits the cell into more than two cells" << endl;
1374
1375         writeOBJ(".", cellI, loopPts, connectedPoints);
1376
1377         return false;
1378     }
1379
1380
1381     // Check that both parts (connectedPoints, otherPoints) have enough faces.
1382     labelHashSet connectedFaces(2*cFaces.size());
1383     labelHashSet otherFaces(2*cFaces.size());
1384
1385     forAllConstIter(Map<label>, pointStatus, iter)
1386     {
1387         label pointI = iter.key();
1388
1389         const labelList& pFaces = mesh().pointFaces()[pointI];
1390
1391         if (iter() == 1)
1392         {
1393             forAll(pFaces, pFaceI)
1394             {
1395                 if (meshTools::faceOnCell(mesh(), cellI, pFaces[pFaceI]))
1396                 {
1397                     connectedFaces.insert(pFaces[pFaceI]);
1398                 }
1399             }
1400         }
1401         else if (iter() == 2)
1402         {
1403             forAll(pFaces, pFaceI)
1404             {
1405                 if (meshTools::faceOnCell(mesh(), cellI, pFaces[pFaceI]))
1406                 {
1407                     otherFaces.insert(pFaces[pFaceI]);
1408                 }
1409             }
1410         }
1411     }
1412
1413     if (connectedFaces.size() < 3)
1414     {
1415         WarningIn("Foam::cellCuts::calcAnchors")
1416             << "Invalid loop " << loop << " for cell " << cellI
1417             << " since would have too few faces on one side." << nl
1418             << "All faces:" << cFaces << endl;
1419
1420         writeOBJ(".", cellI, loopPts, connectedPoints);
1421
1422         return false;
1423     }
1424
1425     if (otherFaces.size() < 3)
1426     {
1427         WarningIn("Foam::cellCuts::calcAnchors")
1428             << "Invalid loop " << loop << " for cell " << cellI
1429             << " since would have too few faces on one side." << nl
1430             << "All faces:" << cFaces << endl;
1431
1432         writeOBJ(".", cellI, loopPts, otherPoints);
1433
1434         return false;
1435     }
1436
1437
1438     // Check that faces are split into two regions and not more.
1439     // When walking across the face the only transition of pointStatus is
1440     // from set1 to loop to set2 or back. Not allowed is from set1 to loop to
1441     // set1.
1442     {
1443         forAll(cFaces, i)
1444         {
1445             label faceI = cFaces[i];
1446
1447             const face& f = mesh().faces()[faceI];
1448
1449             bool hasSet1 = false;
1450             bool hasSet2 = false;
1451
1452             label prevStat = pointStatus[f[0]];
1453
1454             forAll(f, fp)
1455             {
1456                 label v0 = f[fp];
1457                 label pStat = pointStatus[v0];
1458
1459                 if (pStat == prevStat)
1460                 {
1461                 }
1462                 else if (pStat == 0)
1463                 {
1464                     // Loop.
1465                 }
1466                 else if (pStat == 1)
1467                 {
1468                     if (hasSet1)
1469                     {
1470                         // Second occurence of set1.
1471                         WarningIn("Foam::cellCuts::calcAnchors")
1472                             << "Invalid loop " << loop << " for cell " << cellI
1473                             << " since face " << f << " would be split into"
1474                             << " more than two faces" << endl;
1475
1476                         writeOBJ(".", cellI, loopPts, otherPoints);
1477
1478                         return false;
1479                     }
1480
1481                     hasSet1 = true;
1482                 }
1483                 else if (pStat == 2)
1484                 {
1485                     if (hasSet2)
1486                     {
1487                         // Second occurence of set1.
1488                         WarningIn("Foam::cellCuts::calcAnchors")
1489                             << "Invalid loop " << loop << " for cell " << cellI
1490                             << " since face " << f << " would be split into"
1491                             << " more than two faces" << endl;
1492
1493                         writeOBJ(".", cellI, loopPts, otherPoints);
1494
1495                         return false;
1496                     }
1497
1498                     hasSet2 = true;
1499                 }
1500                 else
1501                 {
1502                     FatalErrorIn("Foam::cellCuts::calcAnchors")
1503                         << abort(FatalError);
1504                 }
1505
1506                 prevStat = pStat;
1507
1508
1509                 label v1 = f.nextLabel(fp);
1510                 label edgeI = findEdge(faceI, v0, v1);
1511
1512                 label eStat = edgeStatus[edgeI];
1513
1514                 if (eStat == prevStat)
1515                 {
1516                 }
1517                 else if (eStat == 0)
1518                 {
1519                     // Loop.
1520                 }
1521                 else if (eStat == 1)
1522                 {
1523                     if (hasSet1)
1524                     {
1525                         // Second occurence of set1.
1526                         WarningIn("Foam::cellCuts::calcAnchors")
1527                             << "Invalid loop " << loop << " for cell " << cellI
1528                             << " since face " << f << " would be split into"
1529                             << " more than two faces" << endl;
1530
1531                         writeOBJ(".", cellI, loopPts, otherPoints);
1532
1533                         return false;
1534                     }
1535
1536                     hasSet1 = true;
1537                 }
1538                 else if (eStat == 2)
1539                 {
1540                     if (hasSet2)
1541                     {
1542                         // Second occurence of set1.
1543                         WarningIn("Foam::cellCuts::calcAnchors")
1544                             << "Invalid loop " << loop << " for cell " << cellI
1545                             << " since face " << f << " would be split into"
1546                             << " more than two faces" << endl;
1547
1548                         writeOBJ(".", cellI, loopPts, otherPoints);
1549
1550                         return false;
1551                     }
1552
1553                     hasSet2 = true;
1554                 }
1555                 prevStat = eStat;
1556             }
1557         }
1558     }
1559
1560
1561
1562
1563     // Check which one of point sets to use.
1564     bool loopOk = loopAnchorConsistent(cellI, loopPts, connectedPoints);
1565
1566     //if (debug)
1567     {
1568         // Additional check: are non-anchor points on other side?
1569         bool otherLoopOk = loopAnchorConsistent(cellI, loopPts, otherPoints);
1570
1571         if (loopOk == otherLoopOk)
1572         {
1573             // Both sets of points are supposedly on the same side as the
1574             // loop normal. Oops.
1575
1576             WarningIn("Foam::cellCuts::calcAnchors")
1577                 << " For cell:" << cellI
1578                 << " achorpoints and nonanchorpoints are geometrically"
1579                 << " on same side!" << endl
1580                 << "cellPoints:" << cPoints << endl
1581                 << "loop:" << loop << endl
1582                 << "anchors:" << connectedPoints << endl
1583                 << "otherPoints:" << otherPoints << endl;
1584
1585             writeOBJ(".", cellI, loopPts, connectedPoints);
1586         }
1587     }
1588
1589     if (loopOk)
1590     {
1591         // connectedPoints on 'outside' of loop so these are anchor points
1592         anchorPoints.transfer(connectedPoints);
1593         connectedPoints.clear();
1594     }
1595     else
1596     {
1597         anchorPoints.transfer(otherPoints);
1598         otherPoints.clear();
1599     }
1600     return true;
1601 }
1602
1603
1604 Foam::pointField Foam::cellCuts::loopPoints
1605 (
1606     const labelList& loop,
1607     const scalarField& loopWeights
1608 ) const
1609 {
1610     pointField loopPts(loop.size());
1611
1612     forAll(loop, fp)
1613     {
1614         loopPts[fp] = coord(loop[fp], loopWeights[fp]);
1615     }
1616     return loopPts;
1617 }
1618
1619
1620 // Returns weights of loop. Inverse of loopPoints.
1621 Foam::scalarField Foam::cellCuts::loopWeights(const labelList& loop) const
1622 {
1623     scalarField weights(loop.size());
1624
1625     forAll(loop, fp)
1626     {
1627         label cut = loop[fp];
1628
1629         if (isEdge(cut))
1630         {
1631             label edgeI = getEdge(cut);
1632
1633             weights[fp] = edgeWeight_[edgeI];
1634         }
1635         else
1636         {
1637             weights[fp] = -GREAT;
1638         }
1639     }
1640     return weights;
1641 }
1642
1643
1644 // Check if cut edges in loop are compatible with ones in edgeIsCut_
1645 bool Foam::cellCuts::validEdgeLoop
1646 (
1647     const labelList& loop,
1648     const scalarField& loopWeights
1649 ) const
1650 {
1651     forAll(loop, fp)
1652     {
1653         label cut = loop[fp];
1654
1655         if (isEdge(cut))
1656         {
1657             label edgeI = getEdge(cut);
1658
1659             // Check: cut compatible only if can be snapped to existing one.
1660             if (edgeIsCut_[edgeI])
1661             {
1662                 scalar edgeLen = mesh().edges()[edgeI].mag(mesh().points());
1663
1664                 if
1665                 (
1666                     mag(loopWeights[fp] - edgeWeight_[edgeI])
1667                   > geomCellLooper::snapTol()*edgeLen
1668                 )
1669                 {
1670                     // Positions differ too much->would create two cuts.
1671                     return false;
1672                 }
1673             }
1674         }
1675     }
1676     return true;
1677 }
1678
1679
1680 // Counts cuts on face. Includes cuts through vertices and through edges.
1681 // Assumes that if edge is cut both in edgeIsCut and in loop that the position
1682 // of the cut is the same.
1683 Foam::label Foam::cellCuts::countFaceCuts
1684 (
1685     const label faceI,
1686     const labelList& loop
1687 ) const
1688 {
1689     label nCuts = 0;
1690
1691     // Count cut vertices
1692     const face& f = mesh().faces()[faceI];
1693
1694     forAll(f, fp)
1695     {
1696         label vertI = f[fp];
1697
1698         // Vertex already cut or mentioned in current loop.
1699         if
1700         (
1701             pointIsCut_[vertI]
1702          || (findIndex(loop, vertToEVert(vertI)) != -1)
1703         )
1704         {
1705             nCuts++;
1706         }
1707     }
1708
1709     // Count cut edges.
1710     const labelList& fEdges = mesh().faceEdges()[faceI];
1711
1712     forAll(fEdges, fEdgeI)
1713     {
1714         label edgeI = fEdges[fEdgeI];
1715
1716         // Edge already cut or mentioned in current loop.
1717         if
1718         (
1719             edgeIsCut_[edgeI]
1720          || (findIndex(loop, edgeToEVert(edgeI)) != -1)
1721         )
1722         {
1723             nCuts++;
1724         }
1725     }
1726
1727     return nCuts;
1728 }
1729
1730
1731 // Determine compatibility of loop with existing cut pattern. Does not use
1732 // cut-addressing (faceCuts_, cutCuts_)
1733 bool Foam::cellCuts::conservativeValidLoop
1734 (
1735     const label cellI,
1736     const labelList& loop
1737 ) const
1738 {
1739
1740     if (loop.size() < 2)
1741     {
1742         return false;
1743     }
1744
1745     forAll(loop, cutI)
1746     {
1747         if (isEdge(loop[cutI]))
1748         {
1749             label edgeI = getEdge(loop[cutI]);
1750
1751             if (edgeIsCut_[edgeI])
1752             {
1753                 // edge compatibility already checked.
1754             }
1755             else
1756             {
1757                 // Quick rejection: vertices of edge itself cannot be cut.
1758                 const edge& e = mesh().edges()[edgeI];
1759
1760                 if (pointIsCut_[e.start()] || pointIsCut_[e.end()])
1761                 {
1762                     return false;
1763                 }
1764
1765
1766                 // Check faces using this edge
1767                 const labelList& eFaces = mesh().edgeFaces()[edgeI];
1768
1769                 forAll(eFaces, eFaceI)
1770                 {
1771                     label nCuts = countFaceCuts(eFaces[eFaceI], loop);
1772
1773                     if (nCuts > 2)
1774                     {
1775                         return false;
1776                     }
1777                 }
1778             }
1779         }
1780         else
1781         {
1782             // Vertex cut
1783
1784             label vertI = getVertex(loop[cutI]);
1785
1786             if (!pointIsCut_[vertI])
1787             {
1788                 // New cut through vertex.
1789
1790                 // Check edges using vertex.
1791                 const labelList& pEdges = mesh().pointEdges()[vertI];
1792
1793                 forAll(pEdges, pEdgeI)
1794                 {
1795                     label edgeI = pEdges[pEdgeI];
1796
1797                     if (edgeIsCut_[edgeI])
1798                     {
1799                         return false;
1800                     }
1801                 }
1802
1803                 // Check faces using vertex.
1804                 const labelList& pFaces = mesh().pointFaces()[vertI];
1805
1806                 forAll(pFaces, pFaceI)
1807                 {
1808                     label nCuts = countFaceCuts(pFaces[pFaceI], loop);
1809
1810                     if (nCuts > 2)
1811                     {
1812                         return false;
1813                     }
1814                 }
1815             }
1816         }
1817     }
1818
1819     // All ok.
1820     return true;
1821 }
1822
1823
1824 // Determine compatibility of loop with existing cut pattern. Does not use
1825 // derived cut-addressing (faceCuts), only pointIsCut, edgeIsCut.
1826 // Adds any cross-cuts found to newFaceSplitCut and sets cell points on
1827 // one side of the loop in anchorPoints.
1828 bool Foam::cellCuts::validLoop
1829 (
1830     const label cellI,
1831     const labelList& loop,
1832     const scalarField& loopWeights,
1833
1834     Map<edge>& newFaceSplitCut,
1835     labelList& anchorPoints
1836 ) const
1837 {
1838     if (loop.size() < 2)
1839     {
1840         return false;
1841     }
1842
1843     if (debug & 4)
1844     {
1845         // Allow as fallback the 'old' loop checking where only a single
1846         // cut per face is allowed.
1847         if (!conservativeValidLoop(cellI, loop))
1848         {
1849             return  false;
1850         }
1851     }
1852
1853     forAll(loop, fp)
1854     {
1855         label cut = loop[fp];
1856         label nextCut = loop[(fp+1) % loop.size()];
1857
1858         // Label (if any) of face cut (so cut not along existing edge)
1859         label meshFaceI = -1;
1860
1861         if (isEdge(cut))
1862         {
1863             label edgeI = getEdge(cut);
1864
1865             // Look one cut ahead to find if it is along existing edge.
1866
1867             if (isEdge(nextCut))
1868             {
1869                 // From edge to edge -> cross cut
1870                 label nextEdgeI = getEdge(nextCut);
1871
1872                 // Find face and mark as to be split.
1873                 meshFaceI = edgeEdgeToFace(cellI, edgeI, nextEdgeI);
1874
1875                 if (meshFaceI == -1)
1876                 {
1877                     // Can't find face using both cut edges.
1878                     return false;
1879                 }
1880             }
1881             else
1882             {
1883                 // From edge to vertex -> cross cut only if no existing edge.
1884
1885                 label nextVertI = getVertex(nextCut);
1886                 const edge& e = mesh().edges()[edgeI];
1887
1888                 if (e.start() != nextVertI && e.end() != nextVertI)
1889                 {
1890                     // New edge. Find face and mark as to be split.
1891                     meshFaceI = edgeVertexToFace(cellI, edgeI, nextVertI);
1892
1893                     if (meshFaceI == -1)
1894                     {
1895                         // Can't find face. Ilegal.
1896                         return false;
1897                     }
1898                 }
1899             }
1900         }
1901         else
1902         {
1903             // Cut is vertex.
1904             label vertI = getVertex(cut);
1905
1906             if (isEdge(nextCut))
1907             {
1908                 // From vertex to edge -> cross cut only if no existing edge.
1909                 label nextEdgeI = getEdge(nextCut);
1910
1911                 const edge& nextE = mesh().edges()[nextEdgeI];
1912
1913                 if (nextE.start() != vertI && nextE.end() != vertI)
1914                 {
1915                     // Should be cross cut. Find face.
1916                     meshFaceI = edgeVertexToFace(cellI, nextEdgeI, vertI);
1917
1918                     if (meshFaceI == -1)
1919                     {
1920                         return false;
1921                     }
1922                 }
1923             }
1924             else
1925             {
1926                 // From vertex to vertex -> cross cut only if no existing edge.
1927                 label nextVertI = getVertex(nextCut);
1928
1929                 if (meshTools::findEdge(mesh(), vertI, nextVertI) == -1)
1930                 {
1931                     // New edge. Find face.
1932                     meshFaceI = vertexVertexToFace(cellI, vertI, nextVertI);
1933
1934                     if (meshFaceI == -1)
1935                     {
1936                         return false;
1937                     }
1938                 }
1939             }
1940         }
1941
1942         if (meshFaceI != -1)
1943         {
1944             // meshFaceI is cut across along cut-nextCut (so not along existing
1945             // edge). Check if this is compatible with existing pattern.
1946             edge cutEdge(cut, nextCut);
1947
1948             Map<edge>::const_iterator iter = faceSplitCut_.find(meshFaceI);
1949
1950             if (iter == faceSplitCut_.end())
1951             {
1952                 // Face not yet cut so insert.
1953                 newFaceSplitCut.insert(meshFaceI, cutEdge);
1954             }
1955             else
1956             {
1957                 // Face already cut. Ok if same edge.
1958                 if (iter() != cutEdge)
1959                 {
1960                     return false;
1961                 }
1962             }
1963         }
1964     }
1965
1966     // Is there a face on which all cuts are?
1967     label faceContainingLoop = loopFace(cellI, loop);
1968
1969     if (faceContainingLoop != -1)
1970     {
1971         WarningIn("Foam::cellCuts::validLoop")
1972             << "Found loop on cell " << cellI << " with all points"
1973             << " on face " << faceContainingLoop << endl;
1974
1975         //writeOBJ(".", cellI, loopPoints(loop, loopWeights), labelList(0));
1976
1977         return false;
1978     }
1979
1980     // Calculate anchor points
1981     // Final success is determined by whether anchor points can be determined.
1982     return calcAnchors
1983     (
1984         cellI,
1985         loop,
1986         loopPoints(loop, loopWeights),
1987         anchorPoints
1988     );
1989 }
1990
1991
1992 // Update basic cut information (pointIsCut, edgeIsCut) from cellLoops.
1993 // Assumes cellLoops_ and edgeWeight_ already set and consistent.
1994 // Does not use any other information.
1995 void Foam::cellCuts::setFromCellLoops()
1996 {
1997     // 'Uncut' edges/vertices that are not used in loops.
1998     pointIsCut_ = false;
1999
2000     edgeIsCut_ = false;
2001
2002     faceSplitCut_.clear();
2003
2004     forAll(cellLoops_, cellI)
2005     {
2006         const labelList& loop = cellLoops_[cellI];
2007
2008         if (loop.size())
2009         {
2010             // Storage for cross-face cuts
2011             Map<edge> faceSplitCuts(loop.size());
2012             // Storage for points on one side of cell.
2013             labelList anchorPoints;
2014
2015             if
2016             (
2017                !validLoop
2018                 (
2019                     cellI,
2020                     loop,
2021                     loopWeights(loop),
2022                     faceSplitCuts,
2023                     anchorPoints
2024                 )
2025             )
2026             {
2027                 //writeOBJ(".", cellI, loopPoints(cellI), anchorPoints);
2028
2029                 //FatalErrorIn("cellCuts::setFromCellLoops()")
2030                 WarningIn("cellCuts::setFromCellLoops")
2031                     << "Illegal loop " << loop
2032                     << " when recreating cut-addressing"
2033                     << " from existing cellLoops for cell " << cellI
2034                     << endl;
2035                     //<< abort(FatalError);
2036
2037                 cellLoops_[cellI].setSize(0);
2038                 cellAnchorPoints_[cellI].setSize(0);
2039             }
2040             else
2041             {
2042                 // Copy anchor points.
2043                 cellAnchorPoints_[cellI].transfer(anchorPoints);
2044
2045                 // Copy faceSplitCuts into overall faceSplit info.
2046                 forAllConstIter(Map<edge>, faceSplitCuts, iter)
2047                 {
2048                     faceSplitCut_.insert(iter.key(), iter());
2049                 }
2050
2051                 // Update edgeIsCut, pointIsCut information
2052                 forAll(loop, cutI)
2053                 {
2054                     label cut = loop[cutI];
2055
2056                     if (isEdge(cut))
2057                     {
2058                         edgeIsCut_[getEdge(cut)] = true;
2059                     }
2060                     else
2061                     {
2062                         pointIsCut_[getVertex(cut)] = true;
2063                     }
2064                 }
2065             }
2066         }
2067     }
2068
2069     // Reset edge weights
2070     forAll(edgeIsCut_, edgeI)
2071     {
2072         if (!edgeIsCut_[edgeI])
2073         {
2074             // Weight not used. Set to illegal value to make any use fall over.
2075             edgeWeight_[edgeI] = -GREAT;
2076         }
2077     }
2078 }
2079
2080
2081 // Upate basic cut information from single cellLoop. Returns true if loop
2082 // was valid.
2083 bool Foam::cellCuts::setFromCellLoop
2084 (
2085     const label cellI,
2086     const labelList& loop,
2087     const scalarField& loopWeights
2088 )
2089 {
2090     // Dump loop for debugging.
2091     if (debug)
2092     {
2093         OFstream str("last_cell.obj");
2094
2095         str<< "# edges of cell " << cellI << nl;
2096
2097         meshTools::writeOBJ
2098         (
2099             str,
2100             mesh().cells(),
2101             mesh().faces(),
2102             mesh().points(),
2103             labelList(1, cellI)
2104         );
2105
2106
2107         OFstream loopStr("last_loop.obj");
2108
2109         loopStr<< "# looppoints for cell " << cellI << nl;
2110
2111         pointField pointsOfLoop = loopPoints(loop, loopWeights);
2112
2113         forAll(pointsOfLoop, i)
2114         {
2115             meshTools::writeOBJ(loopStr, pointsOfLoop[i]);
2116         }
2117
2118         str << 'l';
2119
2120         forAll(pointsOfLoop, i)
2121         {
2122             str << ' ' << i + 1;
2123         }
2124         str << ' ' << 1 << nl;
2125     }
2126
2127     bool okLoop = false;
2128
2129     if (validEdgeLoop(loop, loopWeights))
2130     {
2131         // Storage for cuts across face
2132         Map<edge> faceSplitCuts(loop.size());
2133         // Storage for points on one side of cell
2134         labelList anchorPoints;
2135
2136         okLoop =
2137             validLoop(cellI, loop, loopWeights, faceSplitCuts, anchorPoints);
2138
2139         if (okLoop)
2140         {
2141             // Valid loop. Copy cellLoops and anchorPoints
2142             cellLoops_[cellI] = loop;
2143             cellAnchorPoints_[cellI].transfer(anchorPoints);
2144
2145             // Copy split cuts
2146             forAllConstIter(Map<edge>, faceSplitCuts, iter)
2147             {
2148                 faceSplitCut_.insert(iter.key(), iter());
2149             }
2150
2151
2152             // Update edgeIsCut, pointIsCut information
2153             forAll(loop, cutI)
2154             {
2155                 label cut = loop[cutI];
2156
2157                 if (isEdge(cut))
2158                 {
2159                     label edgeI = getEdge(cut);
2160
2161                     edgeIsCut_[edgeI] = true;
2162
2163                     edgeWeight_[edgeI] = loopWeights[cutI];
2164                 }
2165                 else
2166                 {
2167                     label vertI = getVertex(cut);
2168
2169                     pointIsCut_[vertI] = true;
2170                 }
2171             }
2172         }
2173     }
2174
2175     return okLoop;
2176 }
2177
2178
2179 // Update basic cut information from cellLoops. Checks for consistency with
2180 // existing cut pattern.
2181 void Foam::cellCuts::setFromCellLoops
2182 (
2183     const labelList& cellLabels,
2184     const labelListList& cellLoops,
2185     const List<scalarField>& cellLoopWeights
2186 )
2187 {
2188     // 'Uncut' edges/vertices that are not used in loops.
2189     pointIsCut_ = false;
2190
2191     edgeIsCut_ = false;
2192
2193     forAll(cellLabels, cellLabelI)
2194     {
2195         label cellI = cellLabels[cellLabelI];
2196
2197         const labelList& loop = cellLoops[cellLabelI];
2198
2199         if (loop.size())
2200         {
2201             const scalarField& loopWeights = cellLoopWeights[cellLabelI];
2202
2203             if (setFromCellLoop(cellI, loop, loopWeights))
2204             {
2205                 // Valid loop. Call above will have upated all already.
2206             }
2207             else
2208             {
2209                 // Clear cellLoops
2210                 cellLoops_[cellI].setSize(0);
2211             }
2212         }
2213     }
2214 }
2215
2216
2217 // Cut cells and update basic cut information from cellLoops. Checks each loop
2218 // for consistency with existing cut pattern.
2219 void Foam::cellCuts::setFromCellCutter
2220 (
2221     const cellLooper& cellCutter,
2222     const List<refineCell>& refCells
2223 )
2224 {
2225     // 'Uncut' edges/vertices that are not used in loops.
2226     pointIsCut_ = false;
2227
2228     edgeIsCut_ = false;
2229
2230     // storage for loop of cuts (cut vertices and/or cut edges)
2231     labelList cellLoop;
2232     scalarField cellLoopWeights;
2233
2234     // For debugging purposes
2235     DynamicList<label> invalidCutCells(2);
2236     DynamicList<labelList> invalidCutLoops(2);
2237     DynamicList<scalarField> invalidCutLoopWeights(2);
2238
2239
2240     forAll(refCells, refCellI)
2241     {
2242         const refineCell& refCell = refCells[refCellI];
2243
2244         label cellI = refCell.cellNo();
2245
2246         const vector& refDir = refCell.direction();
2247
2248
2249         // Cut cell. Determines cellLoop and cellLoopWeights
2250         bool goodCut =
2251             cellCutter.cut
2252             (
2253                 refDir,
2254                 cellI,
2255
2256                 pointIsCut_,
2257                 edgeIsCut_,
2258                 edgeWeight_,
2259
2260                 cellLoop,
2261                 cellLoopWeights
2262             );
2263
2264         // Check whether edge refinement is on a per face basis compatible with
2265         // current pattern.
2266         if (goodCut)
2267         {
2268             if (setFromCellLoop(cellI, cellLoop, cellLoopWeights))
2269             {
2270                 // Valid loop. Will have updated all info already.
2271             }
2272             else
2273             {
2274                 cellLoops_[cellI].setSize(0);
2275
2276                 // Discarded by validLoop
2277                 if (debug)
2278                 {
2279                     invalidCutCells.append(cellI);
2280                     invalidCutLoops.append(cellLoop);
2281                     invalidCutLoopWeights.append(cellLoopWeights);
2282                 }
2283             }
2284         }
2285         else
2286         {
2287             // Clear cellLoops
2288             cellLoops_[cellI].setSize(0);
2289         }
2290     }
2291
2292     if (debug && invalidCutCells.size())
2293     {
2294         invalidCutCells.shrink();
2295         invalidCutLoops.shrink();
2296         invalidCutLoopWeights.shrink();
2297
2298         fileName cutsFile("invalidLoopCells.obj");
2299
2300         Pout<< "cellCuts : writing inValidLoops cells to " << cutsFile << endl;
2301
2302         OFstream cutsStream(cutsFile);
2303
2304         meshTools::writeOBJ
2305         (
2306             cutsStream,
2307             mesh().cells(),
2308             mesh().faces(),
2309             mesh().points(),
2310             invalidCutCells
2311         );
2312
2313         fileName loopsFile("invalidLoops.obj");
2314
2315         Pout<< "cellCuts : writing inValidLoops loops to " << loopsFile << endl;
2316
2317         OFstream loopsStream(loopsFile);
2318
2319         label vertI = 0;
2320
2321         forAll(invalidCutLoops, i)
2322         {
2323             writeOBJ
2324             (
2325                 loopsStream,
2326                 loopPoints(invalidCutLoops[i], invalidCutLoopWeights[i]),
2327                 vertI
2328             );
2329         }
2330     }
2331 }
2332
2333
2334 // Same as one before but cut plane prescribed (instead of just normal)
2335 void Foam::cellCuts::setFromCellCutter
2336 (
2337     const cellLooper& cellCutter,
2338     const labelList& cellLabels,
2339     const List<plane>& cellCutPlanes
2340 )
2341 {
2342     // 'Uncut' edges/vertices that are not used in loops.
2343     pointIsCut_ = false;
2344
2345     edgeIsCut_ = false;
2346
2347     // storage for loop of cuts (cut vertices and/or cut edges)
2348     labelList cellLoop;
2349     scalarField cellLoopWeights;
2350
2351     // For debugging purposes
2352     DynamicList<label> invalidCutCells(2);
2353     DynamicList<labelList> invalidCutLoops(2);
2354     DynamicList<scalarField> invalidCutLoopWeights(2);
2355
2356
2357     forAll(cellLabels, i)
2358     {
2359         label cellI = cellLabels[i];
2360
2361         // Cut cell. Determines cellLoop and cellLoopWeights
2362         bool goodCut =
2363             cellCutter.cut
2364             (
2365                 cellCutPlanes[i],
2366                 cellI,
2367
2368                 pointIsCut_,
2369                 edgeIsCut_,
2370                 edgeWeight_,
2371
2372                 cellLoop,
2373                 cellLoopWeights
2374             );
2375
2376         // Check whether edge refinement is on a per face basis compatible with
2377         // current pattern.
2378         if (goodCut)
2379         {
2380             if (setFromCellLoop(cellI, cellLoop, cellLoopWeights))
2381             {
2382                 // Valid loop. Will have updated all info already.
2383             }
2384             else
2385             {
2386                 cellLoops_[cellI].setSize(0);
2387
2388                 // Discarded by validLoop
2389                 if (debug)
2390                 {
2391                     invalidCutCells.append(cellI);
2392                     invalidCutLoops.append(cellLoop);
2393                     invalidCutLoopWeights.append(cellLoopWeights);
2394                 }
2395             }
2396         }
2397         else
2398         {
2399             // Clear cellLoops
2400             cellLoops_[cellI].setSize(0);
2401         }
2402     }
2403
2404     if (debug && invalidCutCells.size())
2405     {
2406         invalidCutCells.shrink();
2407         invalidCutLoops.shrink();
2408         invalidCutLoopWeights.shrink();
2409
2410         fileName cutsFile("invalidLoopCells.obj");
2411
2412         Pout<< "cellCuts : writing inValidLoops cells to " << cutsFile << endl;
2413
2414         OFstream cutsStream(cutsFile);
2415
2416         meshTools::writeOBJ
2417         (
2418             cutsStream,
2419             mesh().cells(),
2420             mesh().faces(),
2421             mesh().points(),
2422             invalidCutCells
2423         );
2424
2425         fileName loopsFile("invalidLoops.obj");
2426
2427         Pout<< "cellCuts : writing inValidLoops loops to " << loopsFile << endl;
2428
2429         OFstream loopsStream(loopsFile);
2430
2431         label vertI = 0;
2432
2433         forAll(invalidCutLoops, i)
2434         {
2435             writeOBJ
2436             (
2437                 loopsStream,
2438                 loopPoints(invalidCutLoops[i], invalidCutLoopWeights[i]),
2439                 vertI
2440             );
2441         }
2442     }
2443 }
2444
2445
2446 // Set orientation of loops
2447 void Foam::cellCuts::orientPlanesAndLoops()
2448 {
2449     // Determine anchorPoints if not yet done by validLoop.
2450     forAll(cellLoops_, cellI)
2451     {
2452         const labelList& loop = cellLoops_[cellI];
2453
2454         if (loop.size() && cellAnchorPoints_[cellI].empty())
2455         {
2456             // Leave anchor points empty if illegal loop.
2457             calcAnchors
2458             (
2459                 cellI,
2460                 loop,
2461                 loopPoints(cellI),
2462                 cellAnchorPoints_[cellI]
2463             );
2464         }
2465     }
2466
2467     if (debug & 2)
2468     {
2469         Pout<< "cellAnchorPoints:" << endl;
2470     }
2471     forAll(cellAnchorPoints_, cellI)
2472     {
2473         if (cellLoops_[cellI].size())
2474         {
2475             if (cellAnchorPoints_[cellI].empty())
2476             {
2477                 FatalErrorIn("orientPlanesAndLoops()")
2478                     << "No anchor points for cut cell " << cellI << endl
2479                     << "cellLoop:" << cellLoops_[cellI] << abort(FatalError);
2480             }
2481
2482             if (debug & 2)
2483             {
2484                 Pout<< "    cell:" << cellI << " anchored at "
2485                     << cellAnchorPoints_[cellI] << endl;
2486             }
2487         }
2488     }
2489
2490     // Calculate number of valid cellLoops
2491     nLoops_ = 0;
2492
2493     forAll(cellLoops_, cellI)
2494     {
2495         if (cellLoops_[cellI].size())
2496         {
2497             nLoops_++;
2498         }
2499     }
2500 }
2501
2502
2503 // Do all: calculate addressing, calculate loops splitting cells
2504 void Foam::cellCuts::calcLoopsAndAddressing(const labelList& cutCells)
2505 {
2506     // Sanity check on weights
2507     forAll(edgeIsCut_, edgeI)
2508     {
2509         if (edgeIsCut_[edgeI])
2510         {
2511             scalar weight = edgeWeight_[edgeI];
2512
2513             if (weight < 0 || weight > 1)
2514             {
2515                 FatalErrorIn
2516                 (
2517                     "cellCuts::calcLoopsAndAddressing(const labelList&)"
2518                 )   << "Weight out of range [0,1]. Edge " << edgeI
2519                     << " verts:" << mesh().edges()[edgeI]
2520                     << " weight:" << weight << abort(FatalError);
2521             }
2522         }
2523         else
2524         {
2525             // Weight not used. Set to illegal value to make any use fall over.
2526             edgeWeight_[edgeI] = -GREAT;
2527         }
2528     }
2529
2530
2531     // Calculate faces that split cells in two
2532     calcCellLoops(cutCells);
2533
2534     if (debug & 2)
2535     {
2536         Pout<< "-- cellLoops --" << endl;
2537         forAll(cellLoops_, cellI)
2538         {
2539             const labelList& loop = cellLoops_[cellI];
2540
2541             if (loop.size())
2542             {
2543                 Pout<< "cell:" << cellI << "  ";
2544                 writeCuts(Pout, loop, loopWeights(loop));
2545                 Pout<< endl;
2546             }
2547         }
2548     }
2549
2550     // Redo basic cut information (pointIsCut, edgeIsCut, faceSplitCut)
2551     // using cellLoop only.
2552     setFromCellLoops();
2553 }
2554
2555
2556 void Foam::cellCuts::check() const
2557 {
2558     // Check weights for unsnapped values
2559     forAll(edgeIsCut_, edgeI)
2560     {
2561         if (edgeIsCut_[edgeI])
2562         {
2563             if
2564             (
2565                 edgeWeight_[edgeI] < geomCellLooper::snapTol()
2566              || edgeWeight_[edgeI] > (1 - geomCellLooper::snapTol())
2567             )
2568             {
2569                 // Should have been snapped.
2570                 //FatalErrorIn("cellCuts::check()")
2571                 WarningIn("cellCuts::check()")
2572                     << "edge:" << edgeI << " vertices:"
2573                     << mesh().edges()[edgeI]
2574                     << " weight:" << edgeWeight_[edgeI] << " should have been"
2575                     << " snapped to one of its endpoints"
2576                     //<< abort(FatalError);
2577                     << endl;
2578             }
2579         }
2580         else
2581         {
2582             if (edgeWeight_[edgeI] > - 1)
2583             {
2584                 FatalErrorIn("cellCuts::check()")
2585                     << "edge:" << edgeI << " vertices:"
2586                     << mesh().edges()[edgeI]
2587                     << " weight:" << edgeWeight_[edgeI] << " is not cut but"
2588                     << " its weight is not marked invalid"
2589                     << abort(FatalError);
2590             }
2591         }
2592     }
2593
2594     // Check that all elements of cellloop are registered
2595     forAll(cellLoops_, cellI)
2596     {
2597         const labelList& loop = cellLoops_[cellI];
2598
2599         forAll(loop, i)
2600         {
2601             label cut = loop[i];
2602
2603             if
2604             (
2605                 (isEdge(cut) && !edgeIsCut_[getEdge(cut)])
2606              || (!isEdge(cut) && !pointIsCut_[getVertex(cut)])
2607             )
2608             {
2609                 labelList cuts(1, cut);
2610                 writeCuts(Pout, cuts, loopWeights(cuts));
2611
2612                 FatalErrorIn("cellCuts::check()")
2613                     << "cell:" << cellI << " loop:"
2614                     << loop
2615                     << " cut:" << cut << " is not marked as cut"
2616                     << abort(FatalError);
2617             }
2618         }
2619     }
2620
2621     // Check that no elements of cell loop are anchor point.
2622     forAll(cellLoops_, cellI)
2623     {
2624         const labelList& anchors = cellAnchorPoints_[cellI];
2625
2626         const labelList& loop = cellLoops_[cellI];
2627
2628         if (loop.size() && anchors.empty())
2629         {
2630             FatalErrorIn("cellCuts::check()")
2631                 << "cell:" << cellI << " loop:" << loop
2632                 << " has no anchor points"
2633                 << abort(FatalError);
2634         }
2635
2636
2637         forAll(loop, i)
2638         {
2639             label cut = loop[i];
2640
2641             if
2642             (
2643                 !isEdge(cut)
2644              && findIndex(anchors, getVertex(cut)) != -1
2645             )
2646             {
2647                 FatalErrorIn("cellCuts::check()")
2648                     << "cell:" << cellI << " loop:" << loop
2649                     << " anchor points:" << anchors
2650                     << " anchor:" << getVertex(cut) << " is part of loop"
2651                     << abort(FatalError);
2652             }
2653         }
2654     }
2655
2656
2657     // Check that cut faces have a neighbour that is cut.
2658     forAllConstIter(Map<edge>, faceSplitCut_, iter)
2659     {
2660         label faceI = iter.key();
2661
2662         if (mesh().isInternalFace(faceI))
2663         {
2664             label own = mesh().faceOwner()[faceI];
2665             label nei = mesh().faceNeighbour()[faceI];
2666
2667             if (cellLoops_[own].empty() && cellLoops_[nei].empty())
2668             {
2669                 FatalErrorIn("cellCuts::check()")
2670                     << "Internal face:" << faceI << " cut by " << iter()
2671                     << " has owner:" << own
2672                     << " and neighbour:" << nei
2673                     << " that are both uncut"
2674                     << abort(FatalError);
2675             }
2676         }
2677         else
2678         {
2679             label own = mesh().faceOwner()[faceI];
2680
2681             if (cellLoops_[own].empty())
2682             {
2683                 FatalErrorIn("cellCuts::check()")
2684                     << "Boundary face:" << faceI << " cut by " << iter()
2685                     << " has owner:" << own
2686                     << " that is uncut"
2687                     << abort(FatalError);
2688             }
2689         }
2690     }
2691 }
2692
2693
2694 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
2695
2696 // Construct from cells to cut and pattern of cuts
2697 Foam::cellCuts::cellCuts
2698 (
2699     const polyMesh& mesh,
2700     const labelList& cutCells,
2701     const labelList& meshVerts,
2702     const labelList& meshEdges,
2703     const scalarField& meshEdgeWeights
2704 )
2705 :
2706     edgeVertex(mesh),
2707     pointIsCut_(expand(mesh.nPoints(), meshVerts)),
2708     edgeIsCut_(expand(mesh.nEdges(), meshEdges)),
2709     edgeWeight_(expand(mesh.nEdges(), meshEdges, meshEdgeWeights)),
2710     faceCutsPtr_(NULL),
2711     faceSplitCut_(cutCells.size()),
2712     cellLoops_(mesh.nCells()),
2713     nLoops_(-1),
2714     cellAnchorPoints_(mesh.nCells())
2715 {
2716     if (debug)
2717     {
2718         Pout<< "cellCuts : constructor from cut verts and edges" << endl;
2719     }
2720
2721     calcLoopsAndAddressing(cutCells);
2722
2723     // Calculate planes and flip cellLoops if nessecary
2724     orientPlanesAndLoops();
2725
2726     if (debug)
2727     {
2728         check();
2729     }
2730
2731     clearOut();
2732
2733     if (debug)
2734     {
2735         Pout<< "cellCuts : leaving constructor from cut verts and edges"
2736             << endl;
2737     }
2738 }
2739
2740
2741 // Construct from pattern of cuts. Finds out itself which cells are cut.
2742 // (can go wrong if e.g. all neighbours of cell are refined)
2743 Foam::cellCuts::cellCuts
2744 (
2745     const polyMesh& mesh,
2746     const labelList& meshVerts,
2747     const labelList& meshEdges,
2748     const scalarField& meshEdgeWeights
2749 )
2750 :
2751     edgeVertex(mesh),
2752     pointIsCut_(expand(mesh.nPoints(), meshVerts)),
2753     edgeIsCut_(expand(mesh.nEdges(), meshEdges)),
2754     edgeWeight_(expand(mesh.nEdges(), meshEdges, meshEdgeWeights)),
2755     faceCutsPtr_(NULL),
2756     faceSplitCut_(mesh.nFaces()/10 + 1),
2757     cellLoops_(mesh.nCells()),
2758     nLoops_(-1),
2759     cellAnchorPoints_(mesh.nCells())
2760 {
2761     if (debug)
2762     {
2763         Pout<< "cellCuts : constructor from cellLoops" << endl;
2764     }
2765
2766     calcLoopsAndAddressing(identity(mesh.nCells()));
2767
2768     // Calculate planes and flip cellLoops if nessecary
2769     orientPlanesAndLoops();
2770
2771     if (debug)
2772     {
2773         check();
2774     }
2775
2776     clearOut();
2777
2778     if (debug)
2779     {
2780         Pout<< "cellCuts : leaving constructor from cellLoops" << endl;
2781     }
2782 }
2783
2784
2785 // Construct from complete cellLoops. Assumes correct cut pattern.
2786 // Only constructs cut-cut addressing and cellAnchorPoints
2787 Foam::cellCuts::cellCuts
2788 (
2789     const polyMesh& mesh,
2790     const labelList& cellLabels,
2791     const labelListList& cellLoops,
2792     const List<scalarField>& cellEdgeWeights
2793 )
2794 :
2795     edgeVertex(mesh),
2796     pointIsCut_(mesh.nPoints(), false),
2797     edgeIsCut_(mesh.nEdges(), false),
2798     edgeWeight_(mesh.nEdges(), -GREAT),
2799     faceCutsPtr_(NULL),
2800     faceSplitCut_(cellLabels.size()),
2801     cellLoops_(mesh.nCells()),
2802     nLoops_(-1),
2803     cellAnchorPoints_(mesh.nCells())
2804 {
2805     if (debug)
2806     {
2807         Pout<< "cellCuts : constructor from cellLoops" << endl;
2808     }
2809
2810     // Update pointIsCut, edgeIsCut, faceSplitCut from cell loops.
2811     // Makes sure cuts are consistent
2812     setFromCellLoops(cellLabels, cellLoops, cellEdgeWeights);
2813
2814     // Calculate planes and flip cellLoops if nessecary
2815     orientPlanesAndLoops();
2816
2817     if (debug)
2818     {
2819         check();
2820     }
2821
2822     clearOut();
2823
2824     if (debug)
2825     {
2826         Pout<< "cellCuts : leaving constructor from cellLoops" << endl;
2827     }
2828 }
2829
2830
2831 // Construct from list of cells to cut and cell cutter.
2832 Foam::cellCuts::cellCuts
2833 (
2834     const polyMesh& mesh,
2835     const cellLooper& cellCutter,
2836     const List<refineCell>& refCells
2837 )
2838 :
2839     edgeVertex(mesh),
2840     pointIsCut_(mesh.nPoints(), false),
2841     edgeIsCut_(mesh.nEdges(), false),
2842     edgeWeight_(mesh.nEdges(), -GREAT),
2843     faceCutsPtr_(NULL),
2844     faceSplitCut_(refCells.size()),
2845     cellLoops_(mesh.nCells()),
2846     nLoops_(-1),
2847     cellAnchorPoints_(mesh.nCells())
2848 {
2849     if (debug)
2850     {
2851         Pout<< "cellCuts : constructor from cellCutter" << endl;
2852     }
2853
2854     // Update pointIsCut, edgeIsCut, faceSplitCut from cell loops.
2855     // Makes sure cuts are consistent
2856     setFromCellCutter(cellCutter, refCells);
2857
2858     // Calculate planes and flip cellLoops if nessecary
2859     orientPlanesAndLoops();
2860
2861     if (debug)
2862     {
2863         check();
2864     }
2865
2866     clearOut();
2867
2868     if (debug)
2869     {
2870         Pout<< "cellCuts : leaving constructor from cellCutter" << endl;
2871     }
2872 }
2873
2874
2875 // Construct from list of cells to cut, plane to cut with and cell cutter.
2876 Foam::cellCuts::cellCuts
2877 (
2878     const polyMesh& mesh,
2879     const cellLooper& cellCutter,
2880     const labelList& cellLabels,
2881     const List<plane>& cutPlanes
2882 )
2883 :
2884     edgeVertex(mesh),
2885     pointIsCut_(mesh.nPoints(), false),
2886     edgeIsCut_(mesh.nEdges(), false),
2887     edgeWeight_(mesh.nEdges(), -GREAT),
2888     faceCutsPtr_(NULL),
2889     faceSplitCut_(cellLabels.size()),
2890     cellLoops_(mesh.nCells()),
2891     nLoops_(-1),
2892     cellAnchorPoints_(mesh.nCells())
2893 {
2894     if (debug)
2895     {
2896         Pout<< "cellCuts : constructor from cellCutter with prescribed plane"
2897             << endl;
2898     }
2899
2900     // Update pointIsCut, edgeIsCut, faceSplitCut from cell loops.
2901     // Makes sure cuts are consistent
2902     setFromCellCutter(cellCutter, cellLabels, cutPlanes);
2903
2904     // Calculate planes and flip cellLoops if nessecary
2905     orientPlanesAndLoops();
2906
2907     if (debug)
2908     {
2909         check();
2910     }
2911
2912     clearOut();
2913
2914     if (debug)
2915     {
2916         Pout<< "cellCuts : leaving constructor from cellCutter with prescribed"
2917             << " plane" << endl;
2918     }
2919 }
2920
2921
2922 // Construct from components
2923 Foam::cellCuts::cellCuts
2924 (
2925     const polyMesh& mesh,
2926     const boolList& pointIsCut,
2927     const boolList& edgeIsCut,
2928     const scalarField& edgeWeight,
2929     const Map<edge>& faceSplitCut,
2930     const labelListList& cellLoops,
2931     const label nLoops,
2932     const labelListList& cellAnchorPoints
2933 )
2934 :
2935     edgeVertex(mesh),
2936     pointIsCut_(pointIsCut),
2937     edgeIsCut_(edgeIsCut),
2938     edgeWeight_(edgeWeight),
2939     faceCutsPtr_(NULL),
2940     faceSplitCut_(faceSplitCut),
2941     cellLoops_(cellLoops),
2942     nLoops_(nLoops),
2943     cellAnchorPoints_(cellAnchorPoints)
2944 {
2945     if (debug)
2946     {
2947         Pout<< "cellCuts : constructor from components" << endl;
2948         Pout<< "cellCuts : leaving constructor from components" << endl;
2949     }
2950 }
2951
2952
2953 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
2954
2955 Foam::cellCuts::~cellCuts()
2956 {
2957     clearOut();
2958 }
2959
2960
2961 void Foam::cellCuts::clearOut()
2962 {
2963     deleteDemandDrivenData(faceCutsPtr_);
2964 }
2965
2966
2967 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
2968
2969 Foam::pointField Foam::cellCuts::loopPoints(const label cellI) const
2970 {
2971     const labelList& loop = cellLoops_[cellI];
2972
2973     pointField loopPts(loop.size());
2974
2975     forAll(loop, fp)
2976     {
2977         label cut = loop[fp];
2978
2979         if (isEdge(cut))
2980         {
2981             label edgeI = getEdge(cut);
2982
2983             loopPts[fp] = coord(cut, edgeWeight_[edgeI]);
2984         }
2985         else
2986         {
2987             loopPts[fp] = coord(cut, -GREAT);
2988         }
2989     }
2990     return loopPts;
2991 }
2992
2993
2994 // Flip loop for cell
2995 void Foam::cellCuts::flip(const label cellI)
2996 {
2997     labelList& loop = cellLoops_[cellI];
2998
2999     reverse(loop);
3000
3001     // Reverse anchor point set.
3002     cellAnchorPoints_[cellI] =
3003         nonAnchorPoints
3004         (
3005             mesh().cellPoints()[cellI],
3006             cellAnchorPoints_[cellI],
3007             loop
3008         );
3009 }
3010
3011
3012 // Flip loop only
3013 void Foam::cellCuts::flipLoopOnly(const label cellI)
3014 {
3015     labelList& loop = cellLoops_[cellI];
3016
3017     reverse(loop);
3018 }
3019
3020
3021 void Foam::cellCuts::writeOBJ
3022 (
3023     Ostream& os,
3024     const pointField& loopPts,
3025     label& vertI
3026 ) const
3027 {
3028     label startVertI = vertI;
3029
3030     forAll(loopPts, fp)
3031     {
3032         const point& pt = loopPts[fp];
3033
3034         os << "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << endl;
3035
3036         vertI++;
3037     }
3038
3039     os  << 'f';
3040     forAll(loopPts, fp)
3041     {
3042         os  << ' ' << startVertI + fp + 1;
3043     }
3044     os  << endl;
3045 }
3046
3047
3048 void Foam::cellCuts::writeOBJ(Ostream& os) const
3049 {
3050     label vertI = 0;
3051
3052     forAll(cellLoops_, cellI)
3053     {
3054         writeOBJ(os, loopPoints(cellI), vertI);
3055     }
3056 }
3057
3058
3059 void Foam::cellCuts::writeCellOBJ(const fileName& dir, const label cellI) const
3060 {
3061     const labelList& anchors = cellAnchorPoints_[cellI];
3062
3063     writeOBJ(dir, cellI, loopPoints(cellI), anchors);
3064 }
3065
3066
3067 // ************************************************************************* //