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) 2011 OpenFOAM Foundation
   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 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 "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 "unitConversion.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     forAllConstIter(cellSet, cellsToCut, iter)
 396     {
 397         const label cellI = iter.key();
 398         const labelList& cEdges = cellEdges[cellI];
 399
 400         forAll(cEdges, i)
 401         {
 402             label edgeI = cEdges[i];
 403
 404             label f0, f1;
 405             meshTools::getEdgeFaces(mesh, cellI, edgeI, f0, f1);
 406
 407             vector n0 = faceAreas[f0];
 408             n0 /= mag(n0);
 409
 410             vector n1 = faceAreas[f1];
 411             n1 /= mag(n1);
 412
 413             if
 414             (
 415                 largerAngle
 416                 (
 417                     mesh,
 418                     minCos,
 419                     minSin,
 420
 421                     cellI,
 422                     f0,
 423                     f1,
 424                     n0,
 425                     n1
 426                 )
 427             )
 428             {
 429                 bool splitOk = false;
 430
 431                 if (!geometry && cellShapes[cellI].model() == hex)
 432                 {
 433                     splitOk =
 434                         splitHex
 435                         (
 436                             mesh,
 437                             cellI,
 438                             edgeI,
 439
 440                             cutCells,
 441                             cellLoops,
 442                             cellEdgeWeights
 443                         );
 444                 }
 445                 else
 446                 {
 447                     vector halfNorm;
 448
 449                     if ((own[f0] == cellI) ^ (own[f1] == cellI))
 450                     {
 451                         // Opposite owner orientation
 452                         halfNorm = 0.5*(n0 - n1);
 453                     }
 454                     else
 455                     {
 456                         // Faces have same owner or same neighbour so
 457                         // normals point in same direction
 458                         halfNorm = 0.5*(n0 + n1);
 459                     }
 460
 461                     splitOk =
 462                         splitCell
 463                         (
 464                             mesh,
 465                             cellCutter,
 466                             cellI,
 467                             edgeI,
 468                             halfNorm,
 469
 470                             vertIsCut,
 471                             edgeIsCut,
 472                             edgeWeight,
 473
 474                             cutCells,
 475                             cellLoops,
 476                             cellEdgeWeights
 477                         );
 478                 }
 479
 480
 481                 if (splitOk)
 482                 {
 483                     // Update cell/edge/vertex wise info.
 484                     label index = cellLoops.size() - 1;
 485                     const labelList& loop = cellLoops[index];
 486                     const scalarField& loopWeights = cellEdgeWeights[index];
 487
 488                     forAll(loop, i)
 489                     {
 490                         label cut = loop[i];
 491
 492                         if (ev.isEdge(cut))
 493                         {
 494                             edgeIsCut[ev.getEdge(cut)] = true;
 495                             edgeWeight[ev.getEdge(cut)] = loopWeights[i];
 496                         }
 497                         else
 498                         {
 499                             vertIsCut[ev.getVertex(cut)] = true;
 500                         }
 501                     }
 502
 503                     // Stop checking edges for this cell.
 504                     break;
 505                 }
 506             }
 507         }
 508     }
 509
 510     cutCells.shrink();
 511     cellLoops.shrink();
 512     cellEdgeWeights.shrink();
 513 }
 514
 515
 516 // Main program:
 517
 518 int main(int argc, char *argv[])
 519 {
 520     argList::addNote
 521     (
 522         "split cells with flat faces"
 523     );
 524     #include "addOverwriteOption.H"
 525     argList::noParallel();
 526     argList::validArgs.append("edgeAngle [0..360]");
 527
 528     argList::addOption
 529     (
 530         "set",
 531         "name",
 532         "split cells from specified cellSet only"
 533     );
 534     argList::addBoolOption
 535     (
 536         "geometry",
 537         "use geometric cut for hexes as well"
 538     );
 539     argList::addOption
 540     (
 541         "tol",
 542         "scalar", "edge snap tolerance (default 0.2)"
 543     );
 544
 545     #include "setRootCase.H"
 546     #include "createTime.H"
 547     runTime.functionObjects().off();
 548     #include "createPolyMesh.H"
 549     const word oldInstance = mesh.pointsInstance();
 550
 551     const scalar featureAngle = args.argRead<scalar>(1);
 552     const scalar minCos = Foam::cos(degToRad(featureAngle));
 553     const scalar minSin = Foam::sin(degToRad(featureAngle));
 554
 555     const bool readSet   = args.optionFound("set");
 556     const bool geometry  = args.optionFound("geometry");
 557     const bool overwrite = args.optionFound("overwrite");
 558
 559     const scalar edgeTol = args.optionLookupOrDefault("tol", 0.2);
 560
 561     Info<< "Trying to split cells with internal angles > feature angle\n" << nl
 562         << "featureAngle      : " << featureAngle << nl
 563         << "edge snapping tol : " << edgeTol << nl;
 564     if (readSet)
 565     {
 566         Info<< "candidate cells   : cellSet " << args["set"] << nl;
 567     }
 568     else
 569     {
 570         Info<< "candidate cells   : all cells" << nl;
 571     }
 572     if (geometry)
 573     {
 574         Info<< "hex cuts          : geometric; using edge tolerance" << nl;
 575     }
 576     else
 577     {
 578         Info<< "hex cuts          : topological; cut to opposite edge" << nl;
 579     }
 580     Info<< endl;
 581
 582
 583     // Cell circumference cutter
 584     geomCellLooper cellCutter(mesh);
 585     // Snap all edge cuts close to endpoints to vertices.
 586     geomCellLooper::setSnapTol(edgeTol);
 587
 588     // Candidate cells to cut
 589     cellSet cellsToCut(mesh, "cellsToCut", mesh.nCells()/100);
 590
 591     if (readSet)
 592     {
 593         // Read cells to cut from cellSet
 594         cellSet cells(mesh, args["set"]);
 595
 596         cellsToCut = cells;
 597     }
 598
 599     while (true)
 600     {
 601         if (!readSet)
 602         {
 603             // Try all cells for cutting
 604             for (label cellI = 0; cellI < mesh.nCells(); cellI++)
 605             {
 606                 cellsToCut.insert(cellI);
 607             }
 608         }
 609
 610
 611         // Cut information per cut cell
 612         DynamicList<label> cutCells(mesh.nCells()/10 + 10);
 613         DynamicList<labelList> cellLoops(mesh.nCells()/10 + 10);
 614         DynamicList<scalarField> cellEdgeWeights(mesh.nCells()/10 + 10);
 615
 616         collectCuts
 617         (
 618             mesh,
 619             cellCutter,
 620             geometry,
 621             minCos,
 622             minSin,
 623             cellsToCut,
 624
 625             cutCells,
 626             cellLoops,
 627             cellEdgeWeights
 628         );
 629
 630         cellSet cutSet(mesh, "cutSet", cutCells.size());
 631         forAll(cutCells, i)
 632         {
 633             cutSet.insert(cutCells[i]);
 634         }
 635
 636         // Gets cuts across cells from cuts through edges.
 637         Info<< "Writing " << cutSet.size() << " cells to cut to cellSet "
 638             << cutSet.instance()/cutSet.local()/cutSet.name()
 639             << endl << endl;
 640         cutSet.write();
 641
 642         // Analyze cuts for clashes.
 643         cellCuts cuts
 644         (
 645             mesh,
 646             cutCells,       // cells candidate for cutting
 647             cellLoops,
 648             cellEdgeWeights
 649         );
 650
 651         Info<< "Actually cut cells:" << cuts.nLoops() << nl << endl;
 652
 653         if (cuts.nLoops() == 0)
 654         {
 655             break;
 656         }
 657
 658         // Remove cut cells from cellsToCut  (Note:only relevant if -readSet)
 659         forAll(cuts.cellLoops(), cellI)
 660         {
 661             if (cuts.cellLoops()[cellI].size())
 662             {
 663                 //Info<< "Removing cut cell " << cellI << " from wishlist"
 664                 //    << endl;
 665                 cellsToCut.erase(cellI);
 666             }
 667         }
 668
 669         // At least some cells are cut.
 670         polyTopoChange meshMod(mesh);
 671
 672         // Cutting engine
 673         meshCutter cutter(mesh);
 674
 675         // Insert mesh refinement into polyTopoChange.
 676         cutter.setRefinement(cuts, meshMod);
 677
 678         // Do all changes
 679         Info<< "Morphing ..." << endl;
 680
 681         if (!overwrite)
 682         {
 683             runTime++;
 684         }
 685
 686         autoPtr<mapPolyMesh> morphMap = meshMod.changeMesh(mesh, false);
 687
 688         if (morphMap().hasMotionPoints())
 689         {
 690             mesh.movePoints(morphMap().preMotionPoints());
 691         }
 692
 693         // Update stored labels on meshCutter
 694         cutter.updateMesh(morphMap());
 695
 696         // Update cellSet
 697         cellsToCut.updateMesh(morphMap());
 698
 699         Info<< "Remaining:" << cellsToCut.size() << endl;
 700
 701         // Write resulting mesh
 702         if (overwrite)
 703         {
 704             mesh.setInstance(oldInstance);
 705         }
 706
 707         Info<< "Writing refined morphMesh to time " << runTime.timeName()
 708             << endl;
 709
 710         mesh.write();
 711     }
 712
 713     Info<< "End\n" << endl;
 714
 715     return 0;
 716 }
 717
 718
 719 // ************************************************************************* //