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