applications/utilities/mesh/advanced/splitCells/splitCells.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend 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 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 Description
  25     Utility to split cells with flat faces.
  26
  27     Uses a geometric cut with a plane dividing the edge angle into two so
  28     might produce funny cells. For hexes it will use by default a cut from
  29     edge onto opposite edge (i.e. purely topological).
  30
  31     Options:
  32     - split cells from cellSet only
  33     - use geometric cut for hexes as well
  34
  35     The angle is the angle between two faces sharing an edge as seen from
  36     inside each cell. So a cube will have all angles 90. If you want
  37     to split cells with cell centre outside use e.g. angle 200
  38
  39 \*---------------------------------------------------------------------------*/
  40
  41 #include "argList.H"
  42 #include "objectRegistry.H"
  43 #include "foamTime.H"
  44 #include "directTopoChange.H"
  45 #include "mapPolyMesh.H"
  46 #include "polyMesh.H"
  47 #include "cellCuts.H"
  48 #include "cellSet.H"
  49 #include "cellModeller.H"
  50 #include "meshCutter.H"
  51 #include "mathematicalConstants.H"
  52 #include "geomCellLooper.H"
  53 #include "plane.H"
  54 #include "edgeVertex.H"
  55 #include "meshTools.H"
  56 #include "ListOps.H"
  57
  58 using namespace Foam;
  59
  60 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  61
  62
  63 labelList pack(const boolList& lst)
  64 {
  65     labelList packedLst(lst.size());
  66     label packedI = 0;
  67
  68     forAll(lst, i)
  69     {
  70         if (lst[i])
  71         {
  72             packedLst[packedI++] = i;
  73         }
  74     }
  75     packedLst.setSize(packedI);
  76
  77     return packedLst;
  78 }
  79
  80
  81 scalarField pack(const boolList& lst, const scalarField& elems)
  82 {
  83     scalarField packedElems(lst.size());
  84     label packedI = 0;
  85
  86     forAll(lst, i)
  87     {
  88         if (lst[i])
  89         {
  90             packedElems[packedI++] = elems[i];
  91         }
  92     }
  93     packedElems.setSize(packedI);
  94
  95     return packedElems;
  96 }
  97
  98
  99 // Given sin and cos of max angle between normals calculate whether f0 and f1
 100 // on cellI make larger angle. Uses sinAngle only for quadrant detection.
 101 bool largerAngle
 102 (
 103     const primitiveMesh& mesh,
 104     const scalar cosAngle,
 105     const scalar sinAngle,
 106
 107     const label cellI,
 108     const label f0,             // face label
 109     const label f1,
 110
 111     const vector& n0,           // normal at f0
 112     const vector& n1
 113 )
 114 {
 115     const labelList& own = mesh.faceOwner();
 116
 117     bool sameFaceOrder = !((own[f0] == cellI) ^ (own[f1] == cellI));
 118
 119     // Get cos between faceArea vectors. Correct so flat angle (180 degrees)
 120     // gives -1.
 121     scalar normalCosAngle = n0 & n1;
 122
 123     if (sameFaceOrder)
 124     {
 125         normalCosAngle = -normalCosAngle;
 126     }
 127
 128
 129     // Get cos between faceCentre and normal vector to determine in
 130     // which quadrant angle is. (Is correct for unwarped faces only!)
 131     // Correct for non-outwards pointing normal.
 132     vector c1c0(mesh.faceCentres()[f1] - mesh.faceCentres()[f0]);
 133     c1c0 /= mag(c1c0) + VSMALL;
 134
 135     scalar fcCosAngle = n0 & c1c0;
 136
 137     if (own[f0] != cellI)
 138     {
 139         fcCosAngle = -fcCosAngle;
 140     }
 141
 142     if (sinAngle < 0.0)
 143     {
 144         // Looking for concave angles (quadrant 3 or 4)
 145         if (fcCosAngle <= 0)
 146         {
 147             // Angle is convex so smaller.
 148             return false;
 149         }
 150         else
 151         {
 152             if (normalCosAngle < cosAngle)
 153             {
 154                 return false;
 155             }
 156             else
 157             {
 158                 return true;
 159             }
 160         }
 161     }
 162     else
 163     {
 164         // Looking for convex angles (quadrant 1 or 2)
 165         if (fcCosAngle > 0)
 166         {
 167             // Concave angle
 168             return true;
 169         }
 170         else
 171         {
 172             // Convex. Check cos of normal vectors.
 173             if (normalCosAngle > cosAngle)
 174             {
 175                 return false;
 176             }
 177             else
 178             {
 179                 return true;
 180             }
 181         }
 182     }
 183 }
 184
 185
 186 // Split hex (and hex only) along edgeI creating two prisms
 187 bool splitHex
 188 (
 189     const polyMesh& mesh,
 190     const label cellI,
 191     const label edgeI,
 192
 193     DynamicList<label>& cutCells,
 194     DynamicList<labelList>& cellLoops,
 195     DynamicList<scalarField>& cellEdgeWeights
 196 )
 197 {
 198     // cut handling functions
 199     edgeVertex ev(mesh);
 200
 201     const edgeList& edges = mesh.edges();
 202     const faceList& faces = mesh.faces();
 203
 204     const edge& e = edges[edgeI];
 205
 206     // Get faces on the side, i.e. faces not using edge but still using one of
 207     // the edge endpoints.
 208
 209     label leftI = -1;
 210     label rightI = -1;
 211     label leftFp = -1;
 212     label rightFp = -1;
 213
 214     const cell& cFaces = mesh.cells()[cellI];
 215
 216     forAll(cFaces, i)
 217     {
 218         label faceI = cFaces[i];
 219
 220         const face& f = faces[faceI];
 221
 222         label fp0 = findIndex(f, e[0]);
 223         label fp1 = findIndex(f, e[1]);
 224
 225         if (fp0 == -1)
 226         {
 227             if (fp1 != -1)
 228             {
 229                 // Face uses e[1] but not e[0]
 230                 rightI = faceI;
 231                 rightFp = fp1;
 232
 233                 if (leftI != -1)
 234                 {
 235                     // Have both faces so exit
 236                     break;
 237                 }
 238             }
 239         }
 240         else
 241         {
 242             if (fp1 != -1)
 243             {
 244                 // Face uses both e[1] and e[0]
 245             }
 246             else
 247             {
 248                 leftI = faceI;
 249                 leftFp = fp0;
 250
 251                 if (rightI != -1)
 252                 {
 253                     break;
 254                 }
 255             }
 256         }
 257     }
 258
 259     if (leftI == -1 || rightI == -1)
 260     {
 261         FatalErrorIn("splitHex") << "Problem : leftI:" << leftI
 262             << " rightI:" << rightI << abort(FatalError);
 263     }
 264
 265     // Walk two vertices further on faces.
 266
 267     const face& leftF = faces[leftI];
 268
 269     label leftV = leftF[(leftFp + 2) % leftF.size()];
 270
 271     const face& rightF = faces[rightI];
 272
 273     label rightV = rightF[(rightFp + 2) % rightF.size()];
 274
 275     labelList loop(4);
 276     loop[0] = ev.vertToEVert(e[0]);
 277     loop[1] = ev.vertToEVert(leftV);
 278     loop[2] = ev.vertToEVert(rightV);
 279     loop[3] = ev.vertToEVert(e[1]);
 280
 281     scalarField loopWeights(4);
 282     loopWeights[0] = -GREAT;
 283     loopWeights[1] = -GREAT;
 284     loopWeights[2] = -GREAT;
 285     loopWeights[3] = -GREAT;
 286
 287     cutCells.append(cellI);
 288     cellLoops.append(loop);
 289     cellEdgeWeights.append(loopWeights);
 290
 291     return true;
 292 }
 293
 294
 295 // Split cellI along edgeI with a plane along halfNorm direction.
 296 bool splitCell
 297 (
 298     const polyMesh& mesh,
 299     const geomCellLooper& cellCutter,
 300
 301     const label cellI,
 302     const label edgeI,
 303     const vector& halfNorm,
 304
 305     const boolList& vertIsCut,
 306     const boolList& edgeIsCut,
 307     const scalarField& edgeWeight,
 308
 309     DynamicList<label>& cutCells,
 310     DynamicList<labelList>& cellLoops,
 311     DynamicList<scalarField>& cellEdgeWeights
 312 )
 313 {
 314     const edge& e = mesh.edges()[edgeI];
 315
 316     vector eVec = e.vec(mesh.points());
 317     eVec /= mag(eVec);
 318
 319     vector planeN = eVec ^ halfNorm;
 320
 321     // Slightly tilt plane to make it not cut edges exactly
 322     // halfway on fully regular meshes (since we want cuts
 323     // to be snapped to vertices)
 324     planeN += 0.01*halfNorm;
 325
 326     planeN /= mag(planeN);
 327
 328     // Define plane through edge
 329     plane cutPlane(mesh.points()[e.start()], planeN);
 330
 331     labelList loop;
 332     scalarField loopWeights;
 333
 334     if
 335     (
 336         cellCutter.cut
 337         (
 338             cutPlane,
 339             cellI,
 340             vertIsCut,
 341             edgeIsCut,
 342             edgeWeight,
 343             loop,
 344             loopWeights
 345         )
 346     )
 347     {
 348         // Did manage to cut cell. Copy into overall list.
 349         cutCells.append(cellI);
 350         cellLoops.append(loop);
 351         cellEdgeWeights.append(loopWeights);
 352
 353         return true;
 354     }
 355     else
 356     {
 357         return false;
 358     }
 359 }
 360
 361
 362 // Collects cuts for all cells in cellSet
 363 void collectCuts
 364 (
 365     const polyMesh& mesh,
 366     const geomCellLooper& cellCutter,
 367     const bool geometry,
 368     const scalar minCos,
 369     const scalar minSin,
 370     const cellSet& cellsToCut,
 371
 372     DynamicList<label>& cutCells,
 373     DynamicList<labelList>& cellLoops,
 374     DynamicList<scalarField>& cellEdgeWeights
 375 )
 376 {
 377     // Get data from mesh
 378     const cellShapeList& cellShapes = mesh.cellShapes();
 379     const labelList& own = mesh.faceOwner();
 380     const labelListList& cellEdges = mesh.cellEdges();
 381     const vectorField& faceAreas = mesh.faceAreas();
 382
 383     // Hex shape
 384     const cellModel& hex = *(cellModeller::lookup("hex"));
 385
 386     // cut handling functions
 387     edgeVertex ev(mesh);
 388
 389
 390     // Cut information per mesh entity
 391     boolList vertIsCut(mesh.nPoints(), false);
 392     boolList edgeIsCut(mesh.nEdges(), false);
 393     scalarField edgeWeight(mesh.nEdges(), -GREAT);
 394
 395     for
 396     (
 397         cellSet::const_iterator iter = cellsToCut.begin();
 398         iter != cellsToCut.end();
 399         ++iter
 400     )
 401     {
 402         label cellI = iter.key();
 403
 404         const labelList& cEdges = cellEdges[cellI];
 405
 406         forAll(cEdges, i)
 407         {
 408             label edgeI = cEdges[i];
 409
 410             label f0, f1;
 411             meshTools::getEdgeFaces(mesh, cellI, edgeI, f0, f1);
 412
 413             vector n0 = faceAreas[f0];
 414             n0 /= mag(n0);
 415
 416             vector n1 = faceAreas[f1];
 417             n1 /= mag(n1);
 418
 419             if
 420             (
 421                 largerAngle
 422                 (
 423                     mesh,
 424                     minCos,
 425                     minSin,
 426
 427                     cellI,
 428                     f0,
 429                     f1,
 430                     n0,
 431                     n1
 432                 )
 433             )
 434             {
 435                 bool splitOk = false;
 436
 437                 if (!geometry && cellShapes[cellI].model() == hex)
 438                 {
 439                     splitOk =
 440                         splitHex
 441                         (
 442                             mesh,
 443                             cellI,
 444                             edgeI,
 445
 446                             cutCells,
 447                             cellLoops,
 448                             cellEdgeWeights
 449                         );
 450                 }
 451                 else
 452                 {
 453                     vector halfNorm;
 454
 455                     if ((own[f0] == cellI) ^ (own[f1] == cellI))
 456                     {
 457                         // Opposite owner orientation
 458                         halfNorm = 0.5*(n0 - n1);
 459                     }
 460                     else
 461                     {
 462                         // Faces have same owner or same neighbour so
 463                         // normals point in same direction
 464                         halfNorm = 0.5*(n0 + n1);
 465                     }
 466
 467                     splitOk =
 468                         splitCell
 469                         (
 470                             mesh,
 471                             cellCutter,
 472                             cellI,
 473                             edgeI,
 474                             halfNorm,
 475
 476                             vertIsCut,
 477                             edgeIsCut,
 478                             edgeWeight,
 479
 480                             cutCells,
 481                             cellLoops,
 482                             cellEdgeWeights
 483                         );
 484                 }
 485
 486
 487                 if (splitOk)
 488                 {
 489                     // Update cell/edge/vertex wise info.
 490                     label index = cellLoops.size() - 1;
 491                     const labelList& loop = cellLoops[index];
 492                     const scalarField& loopWeights = cellEdgeWeights[index];
 493
 494                     forAll(loop, i)
 495                     {
 496                         label cut = loop[i];
 497
 498                         if (ev.isEdge(cut))
 499                         {
 500                             edgeIsCut[ev.getEdge(cut)] = true;
 501                             edgeWeight[ev.getEdge(cut)] = loopWeights[i];
 502                         }
 503                         else
 504                         {
 505                             vertIsCut[ev.getVertex(cut)] = true;
 506                         }
 507                     }
 508
 509                     // Stop checking edges for this cell.
 510                     break;
 511                 }
 512             }
 513         }
 514     }
 515
 516     cutCells.shrink();
 517     cellLoops.shrink();
 518     cellEdgeWeights.shrink();
 519 }
 520
 521
 522 // Main program:
 523
 524 int main(int argc, char *argv[])
 525 {
 526     argList::noParallel();
 527     argList::validOptions.insert("set", "cellSet name");
 528     argList::validOptions.insert("geometry", "");
 529     argList::validOptions.insert("tol", "edge snap tolerance");
 530     argList::validOptions.insert("overwrite", "");
 531     argList::validArgs.append("edge angle [0..360]");
 532
 533 #   include "setRootCase.H"
 534 #   include "createTime.H"
 535     runTime.functionObjects().off();
 536 #   include "createPolyMesh.H"
 537     const word oldInstance = mesh.pointsInstance();
 538
 539     scalar featureAngle(readScalar(IStringStream(args.additionalArgs()[0])()));
 540
 541     scalar radAngle = featureAngle * mathematicalConstant::pi/180.0;
 542     scalar minCos = Foam::cos(radAngle);
 543     scalar minSin = Foam::sin(radAngle);
 544
 545     bool readSet   = args.optionFound("set");
 546     bool geometry  = args.optionFound("geometry");
 547     bool overwrite = args.optionFound("overwrite");
 548
 549     scalar edgeTol = 0.2;
 550     args.optionReadIfPresent("tol", edgeTol);
 551
 552     Info<< "Trying to split cells with internal angles > feature angle\n" << nl
 553         << "featureAngle      : " << featureAngle << nl
 554         << "edge snapping tol : " << edgeTol << nl;
 555     if (readSet)
 556     {
 557         Info<< "candidate cells   : cellSet " << args.option("set") << nl;
 558     }
 559     else
 560     {
 561         Info<< "candidate cells   : all cells" << nl;
 562     }
 563     if (geometry)
 564     {
 565         Info<< "hex cuts          : geometric; using edge tolerance" << nl;
 566     }
 567     else
 568     {
 569         Info<< "hex cuts          : topological; cut to opposite edge" << nl;
 570     }
 571     Info<< endl;
 572
 573
 574     // Cell circumference cutter
 575     geomCellLooper cellCutter(mesh);
 576     // Snap all edge cuts close to endpoints to vertices.
 577     geomCellLooper::setSnapTol(edgeTol);
 578
 579     // Candidate cells to cut
 580     cellSet cellsToCut(mesh, "cellsToCut", mesh.nCells()/100);
 581
 582     if (readSet)
 583     {
 584         // Read cells to cut from cellSet
 585         cellSet cells(mesh, args.option("set"));
 586
 587         cellsToCut = cells;
 588     }
 589
 590     while (true)
 591     {
 592         if (!readSet)
 593         {
 594             // Try all cells for cutting
 595             for (label cellI = 0; cellI < mesh.nCells(); cellI++)
 596             {
 597                 cellsToCut.insert(cellI);
 598             }
 599         }
 600
 601
 602         // Cut information per cut cell
 603         DynamicList<label> cutCells(mesh.nCells()/10 + 10);
 604         DynamicList<labelList> cellLoops(mesh.nCells()/10 + 10);
 605         DynamicList<scalarField> cellEdgeWeights(mesh.nCells()/10 + 10);
 606
 607         collectCuts
 608         (
 609             mesh,
 610             cellCutter,
 611             geometry,
 612             minCos,
 613             minSin,
 614             cellsToCut,
 615
 616             cutCells,
 617             cellLoops,
 618             cellEdgeWeights
 619         );
 620
 621         cellSet cutSet(mesh, "cutSet", cutCells.size());
 622         forAll(cutCells, i)
 623         {
 624             cutSet.insert(cutCells[i]);
 625         }
 626
 627         // Gets cuts across cells from cuts through edges.
 628         Info<< "Writing " << cutSet.size() << " cells to cut to cellSet "
 629             << cutSet.instance()/cutSet.local()/cutSet.name()
 630             << endl << endl;
 631         cutSet.write();
 632
 633         // Analyze cuts for clashes.
 634         cellCuts cuts
 635         (
 636             mesh,
 637             cutCells,       // cells candidate for cutting
 638             cellLoops,
 639             cellEdgeWeights
 640         );
 641
 642         Info<< "Actually cut cells:" << cuts.nLoops() << nl << endl;
 643
 644         if (cuts.nLoops() == 0)
 645         {
 646             break;
 647         }
 648
 649         // Remove cut cells from cellsToCut  (Note:only relevant if -readSet)
 650         forAll(cuts.cellLoops(), cellI)
 651         {
 652             if (cuts.cellLoops()[cellI].size())
 653             {
 654                 //Info<< "Removing cut cell " << cellI << " from wishlist"
 655                 //    << endl;
 656                 cellsToCut.erase(cellI);
 657             }
 658         }
 659
 660         // At least some cells are cut.
 661         directTopoChange meshMod(mesh);
 662
 663         // Cutting engine
 664         meshCutter cutter(mesh);
 665
 666         // Insert mesh refinement into directTopoChange.
 667         cutter.setRefinement(cuts, meshMod);
 668
 669         // Do all changes
 670         Info<< "Morphing ..." << endl;
 671
 672         if (!overwrite)
 673         {
 674             runTime++;
 675         }
 676
 677         autoPtr<mapPolyMesh> morphMap = meshMod.changeMesh(mesh, false);
 678
 679         if (morphMap().hasMotionPoints())
 680         {
 681             mesh.movePoints(morphMap().preMotionPoints());
 682         }
 683
 684         // Update stored labels on meshCutter
 685         cutter.updateMesh(morphMap());
 686
 687         // Update cellSet
 688         cellsToCut.updateMesh(morphMap());
 689
 690         Info<< "Remaining:" << cellsToCut.size() << endl;
 691
 692         // Write resulting mesh
 693         if (overwrite)
 694         {
 695             mesh.setInstance(oldInstance);
 696         }
 697
 698         Info<< "Writing refined morphMesh to time " << runTime.timeName()
 699             << endl;
 700
 701         mesh.write();
 702     }
 703
 704     Info << "End\n" << endl;
 705
 706     return 0;
 707 }
 708
 709
 710 // ************************************************************************* //