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

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of OpenFOAM.
  10
  11     OpenFOAM is free software; you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by the
  13     Free Software Foundation; either version 2 of the License, or (at your
  14     option) any later version.
  15
  16     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with OpenFOAM; if not, write to the Free Software Foundation,
  23     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
  24
  25 Description
  26     Utility to split cells with flat faces. Uses a geometric cut with a plane
  27     dividing the edge angle into two so might produce funny cells. For hexes
  28     it will use by default a cut from edge onto opposite edge (i.e. purely
  29     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 "Time.H"
  43 #include "polyTopoChange.H"
  44 #include "polyTopoChanger.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 #   include "createPolyMesh.H"
 536
 537     scalar featureAngle(readScalar(IStringStream(args.additionalArgs()[0])()));
 538
 539     scalar radAngle = featureAngle * mathematicalConstant::pi/180.0;
 540     scalar minCos = Foam::cos(radAngle);
 541     scalar minSin = Foam::sin(radAngle);
 542
 543     bool readSet = args.options().found("set");
 544     bool geometry = args.options().found("geometry");
 545     bool overwrite = args.options().found("overwrite");
 546
 547     scalar edgeTol = 0.2;
 548
 549     if (args.options().found("tol"))
 550     {
 551         edgeTol = readScalar(IStringStream(args.options()["tol"])());
 552     }
 553
 554     Info<< "Trying to split cells with internal angles > feature angle\n" << nl
 555         << "featureAngle      : " << featureAngle << nl
 556         << "edge snapping tol : " << edgeTol << nl;
 557     if (readSet)
 558     {
 559         Info<< "candidate cells   : cellSet " << args.options()["set"] << nl;
 560     }
 561     else
 562     {
 563         Info<< "candidate cells   : all cells" << nl;
 564     }
 565     if (geometry)
 566     {
 567         Info<< "hex cuts          : geometric; using edge tolerance" << nl;
 568     }
 569     else
 570     {
 571         Info<< "hex cuts          : topological; cut to opposite edge" << nl;
 572     }
 573     Info<< endl;
 574
 575
 576     // Cell circumference cutter
 577     geomCellLooper cellCutter(mesh);
 578     // Snap all edge cuts close to endpoints to vertices.
 579     geomCellLooper::setSnapTol(edgeTol);
 580
 581     // Candidate cells to cut
 582     cellSet cellsToCut(mesh, "cellsToCut", mesh.nCells()/100);
 583
 584     if (readSet)
 585     {
 586         // Read cells to cut from cellSet
 587         cellSet cells(mesh, args.options()["set"]);
 588
 589         cellsToCut = cells;
 590     }
 591
 592     while (true)
 593     {
 594         if (!readSet)
 595         {
 596             // Try all cells for cutting
 597             for (label cellI = 0; cellI < mesh.nCells(); cellI++)
 598             {
 599                 cellsToCut.insert(cellI);
 600             }
 601         }
 602
 603
 604         // Cut information per cut cell
 605         DynamicList<label> cutCells(mesh.nCells()/10 + 10);
 606         DynamicList<labelList> cellLoops(mesh.nCells()/10 + 10);
 607         DynamicList<scalarField> cellEdgeWeights(mesh.nCells()/10 + 10);
 608
 609         collectCuts
 610         (
 611             mesh,
 612             cellCutter,
 613             geometry,
 614             minCos,
 615             minSin,
 616             cellsToCut,
 617
 618             cutCells,
 619             cellLoops,
 620             cellEdgeWeights
 621         );
 622
 623         cellSet cutSet(mesh, "cutSet", cutCells.size());
 624         forAll(cutCells, i)
 625         {
 626             cutSet.insert(cutCells[i]);
 627         }
 628
 629         // Gets cuts across cells from cuts through edges.
 630         Info<< "Writing " << cutSet.size() << " cells to cut to cellSet "
 631             << cutSet.instance()/cutSet.local()/cutSet.name()
 632             << endl << endl;
 633         cutSet.write();
 634
 635         // Analyze cuts for clashes.
 636         cellCuts cuts
 637         (
 638             mesh,
 639             cutCells,       // cells candidate for cutting
 640             cellLoops,
 641             cellEdgeWeights
 642         );
 643
 644         Info<< "Actually cut cells:" << cuts.nLoops() << nl << endl;
 645
 646         if (cuts.nLoops() == 0)
 647         {
 648             break;
 649         }
 650
 651         // Remove cut cells from cellsToCut  (Note:only relevant if -readSet)
 652         forAll(cuts.cellLoops(), cellI)
 653         {
 654             if (cuts.cellLoops()[cellI].size() > 0)
 655             {
 656                 //Info<< "Removing cut cell " << cellI << " from wishlist"
 657                 //    << endl;
 658                 cellsToCut.erase(cellI);
 659             }
 660         }
 661
 662         // At least some cells are cut.
 663         polyTopoChange meshMod(mesh);
 664
 665         // Cutting engine
 666         meshCutter cutter(mesh);
 667
 668         // Insert mesh refinement into polyTopoChange.
 669         cutter.setRefinement(cuts, meshMod);
 670
 671         // Do all changes
 672         Info<< "Morphing ..." << endl;
 673
 674         if (!overwrite)
 675         {
 676             runTime++;
 677         }
 678
 679         autoPtr<mapPolyMesh> morphMap = meshMod.changeMesh(mesh, false);
 680
 681         if (morphMap().hasMotionPoints())
 682         {
 683             mesh.movePoints(morphMap().preMotionPoints());
 684         }
 685
 686         // Update stored labels on meshCutter
 687         cutter.updateMesh(morphMap());
 688
 689         // Update cellSet
 690         cellsToCut.updateMesh(morphMap());
 691
 692         Info<< "Remaining:" << cellsToCut.size() << endl;
 693
 694         // Write resulting mesh
 695         Info << "Writing refined morphMesh to time " << runTime.value() << endl;
 696
 697         mesh.write();
 698     }
 699
 700     Info << "End\n" << endl;
 701
 702     return 0;
 703 }
 704
 705
 706 // ************************************************************************* //