applications/utilities/mesh/advanced/collapseEdges/collapseEdges.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     Collapse short edges and combines edges that are in line.
  27
  28     - collapse short edges. Length of edges to collapse provided as argument.
  29     - merge two edges if they are in line. Maximum angle provided as argument.
  30     - remove unused points.
  31
  32     Cannot remove cells. Can remove faces and points but does not check
  33     for nonsense resulting topology.
  34
  35     When collapsing an edge with one point on the boundary it will leave
  36     the boundary point intact. When both points inside it chooses random. When
  37     both points on boundary random again. Note: it should in fact use features
  38     where if one point is on a feature it collapses to that one. Alas we don't
  39     have features on a polyMesh.
  40
  41 \*---------------------------------------------------------------------------*/
  42
  43 #include "argList.H"
  44 #include "Time.H"
  45 #include "edgeCollapser.H"
  46 #include "polyTopoChange.H"
  47 #include "polyTopoChanger.H"
  48 #include "polyMesh.H"
  49 #include "mapPolyMesh.H"
  50 #include "mathematicalConstants.H"
  51 #include "PackedBoolList.H"
  52 #include "SortableList.H"
  53
  54 using namespace Foam;
  55
  56 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  57
  58 // Get faceEdges in order of face points, i.e. faceEdges[0] is between
  59 // f[0] and f[1]
  60 labelList getSortedEdges
  61 (
  62     const edgeList& edges,
  63     const labelList& f,
  64     const labelList& edgeLabels
  65 )
  66 {
  67     labelList faceEdges(edgeLabels.size(), -1);
  68
  69     // Find starting pos in f for every edgeLabels
  70     forAll(edgeLabels, i)
  71     {
  72         label edgeI = edgeLabels[i];
  73
  74         const edge& e = edges[edgeI];
  75
  76         label fp = findIndex(f, e[0]);
  77         label fp1 = f.fcIndex(fp);
  78
  79         if (f[fp1] == e[1])
  80         {
  81             // EdgeI between fp -> fp1
  82             faceEdges[fp] = edgeI;
  83         }
  84         else
  85         {
  86             // EdgeI between fp-1 -> fp
  87             faceEdges[f.rcIndex(fp)] = edgeI;
  88         }
  89     }
  90
  91     return faceEdges;
  92 }
  93
  94
  95 // Merges edges which are in straight line. I.e. edge split by point.
  96 label mergeEdges
  97 (
  98     const polyMesh& mesh,
  99     const scalar maxCos,
 100     edgeCollapser& collapser
 101 )
 102 {
 103     const pointField& points = mesh.points();
 104     const edgeList& edges = mesh.edges();
 105     const labelListList& pointEdges = mesh.pointEdges();
 106     const labelList& region = collapser.pointRegion();
 107     const labelList& master = collapser.pointRegionMaster();
 108
 109     label nCollapsed = 0;
 110
 111     forAll(pointEdges, pointI)
 112     {
 113         const labelList& pEdges = pointEdges[pointI];
 114
 115         if (pEdges.size() == 2)
 116         {
 117             const edge& leftE = edges[pEdges[0]];
 118             const edge& rightE = edges[pEdges[1]];
 119
 120             // Get the two vertices on both sides of the point
 121             label leftV = leftE.otherVertex(pointI);
 122             label rightV = rightE.otherVertex(pointI);
 123
 124             // Collapse only if none of the points part of merge network
 125             // or all of networks with different masters.
 126             label midMaster = -1;
 127             if (region[pointI] != -1)
 128             {
 129                 midMaster = master[region[pointI]];
 130             }
 131
 132             label leftMaster = -2;
 133             if (region[leftV] != -1)
 134             {
 135                 leftMaster = master[region[leftV]];
 136             }
 137
 138             label rightMaster = -3;
 139             if (region[rightV] != -1)
 140             {
 141                 rightMaster = master[region[rightV]];
 142             }
 143
 144             if
 145             (
 146                 midMaster != leftMaster
 147              && midMaster != rightMaster
 148              && leftMaster != rightMaster
 149             )
 150             {
 151                 // Check if the two edge are in line
 152                 vector leftVec = points[pointI] - points[leftV];
 153                 leftVec /= mag(leftVec) + VSMALL;
 154
 155                 vector rightVec = points[rightV] - points[pointI];
 156                 rightVec /= mag(rightVec) + VSMALL;
 157
 158                 if ((leftVec & rightVec) > maxCos)
 159                 {
 160                     // Collapse one (left) side of the edge. Make left vertex
 161                     // the master.
 162                     //if (collapser.unaffectedEdge(pEdges[0]))
 163                     {
 164                         collapser.collapseEdge(pEdges[0], leftV);
 165                         nCollapsed++;
 166                     }
 167                 }
 168             }
 169         }
 170     }
 171
 172     return nCollapsed;
 173 }
 174
 175
 176 // Return master point edge needs to be collapsed to (or -1)
 177 label edgeMaster(const PackedBoolList& boundaryPoint, const edge& e)
 178 {
 179     label masterPoint = -1;
 180
 181     // Collapse edge to boundary point.
 182     if (boundaryPoint.get(e[0]))
 183     {
 184         if (boundaryPoint.get(e[1]))
 185         {
 186             // Both points on boundary. Choose one to collapse to.
 187             // Note: should look at feature edges/points!
 188             masterPoint = e[0];
 189         }
 190         else
 191         {
 192             masterPoint = e[0];
 193         }
 194     }
 195     else
 196     {
 197         if (boundaryPoint.get(e[1]))
 198         {
 199             masterPoint = e[1];
 200         }
 201         else
 202         {
 203             // None on boundary. Choose arbitrary.
 204             // Note: should look at geometry?
 205             masterPoint = e[0];
 206         }
 207     }
 208     return masterPoint;
 209 }
 210
 211
 212 label collapseSmallEdges
 213 (
 214     const polyMesh& mesh,
 215     const PackedBoolList& boundaryPoint,
 216     const scalar minLen,
 217     edgeCollapser& collapser
 218 )
 219 {
 220     const pointField& points = mesh.points();
 221     const edgeList& edges = mesh.edges();
 222
 223     // Collapse all edges that are too small. Choose intelligently which
 224     // point to collapse edge to.
 225
 226     label nCollapsed = 0;
 227
 228     forAll(edges, edgeI)
 229     {
 230         const edge& e = edges[edgeI];
 231
 232         if (e.mag(points) < minLen)
 233         {
 234             label master = edgeMaster(boundaryPoint, e);
 235
 236             if (master != -1) // && collapser.unaffectedEdge(edgeI))
 237             {
 238                 collapser.collapseEdge(edgeI, master);
 239                 nCollapsed++;
 240             }
 241         }
 242     }
 243     return nCollapsed;
 244 }
 245
 246
 247 // Faces which have edges just larger than collapse length but faces which
 248 // are very small. This one tries to collapse them if it can be done with
 249 // edge collapse. For faces where a face gets replace by two edges use
 250 // collapseFaces
 251 label collapseHighAspectFaces
 252 (
 253     const polyMesh& mesh,
 254     const PackedBoolList& boundaryPoint,
 255     const scalar areaFac,
 256     const scalar edgeRatio,
 257     edgeCollapser& collapser
 258 )
 259 {
 260     const pointField& points = mesh.points();
 261     const edgeList& edges = mesh.edges();
 262     const faceList& faces = mesh.faces();
 263     const labelListList& faceEdges = mesh.faceEdges();
 264
 265     scalarField magArea(mag(mesh.faceAreas()));
 266
 267     label maxIndex = findMax(magArea);
 268
 269     scalar minArea = areaFac * magArea[maxIndex];
 270
 271     Info<< "Max face area:" << magArea[maxIndex] << endl
 272         << "Collapse area factor:" << areaFac << endl
 273         << "Collapse area:" << minArea << endl;
 274
 275     label nCollapsed = 0;
 276
 277     forAll(faces, faceI)
 278     {
 279         if (magArea[faceI] < minArea)
 280         {
 281             const face& f = faces[faceI];
 282
 283             // Get the edges in face point order
 284             labelList fEdges(getSortedEdges(edges, f, faceEdges[faceI]));
 285
 286             SortableList<scalar> lengths(fEdges.size());
 287             forAll(fEdges, i)
 288             {
 289                 lengths[i] = edges[fEdges[i]].mag(points);
 290             }
 291             lengths.sort();
 292
 293
 294             label edgeI = -1;
 295
 296             if (f.size() == 4)
 297             {
 298                 // Compare second largest to smallest
 299                 if (lengths[2] > edgeRatio*lengths[0])
 300                 {
 301                     // Collapse smallest only. Triangle should be cleared
 302                     // next time around.
 303                     edgeI = fEdges[lengths.indices()[0]];
 304                 }
 305             }
 306             else if (f.size() == 3)
 307             {
 308                 // Compare second largest to smallest
 309                 if (lengths[1] > edgeRatio*lengths[0])
 310                 {
 311                     edgeI = fEdges[lengths.indices()[0]];
 312                 }
 313             }
 314
 315
 316             if (edgeI != -1)
 317             {
 318                 label master = edgeMaster(boundaryPoint, edges[edgeI]);
 319
 320                 if (master != -1)// && collapser.unaffectedEdge(edgeI))
 321                 {
 322                     collapser.collapseEdge(edgeI, master);
 323                     nCollapsed++;
 324                 }
 325             }
 326         }
 327     }
 328
 329     return nCollapsed;
 330 }
 331
 332
 333 void set(const labelList& elems, const bool val, boolList& status)
 334 {
 335     forAll(elems, i)
 336     {
 337         status[elems[i]] = val;
 338     }
 339 }
 340
 341
 342 // Tries to simplify polygons to face of minSize (4=quad, 3=triangle)
 343 label simplifyFaces
 344 (
 345     const polyMesh& mesh,
 346     const PackedBoolList& boundaryPoint,
 347     const label minSize,
 348     const scalar lenGap,
 349     edgeCollapser& collapser
 350 )
 351 {
 352     const pointField& points = mesh.points();
 353     const edgeList& edges = mesh.edges();
 354     const faceList& faces = mesh.faces();
 355     const cellList& cells = mesh.cells();
 356     const labelListList& faceEdges = mesh.faceEdges();
 357     const labelList& faceOwner = mesh.faceOwner();
 358     const labelList& faceNeighbour = mesh.faceNeighbour();
 359     const labelListList& pointCells = mesh.pointCells();
 360     const labelListList& cellEdges = mesh.cellEdges();
 361
 362     label nCollapsed = 0;
 363
 364     boolList protectedEdge(mesh.nEdges(), false);
 365
 366     forAll(faces, faceI)
 367     {
 368         const face& f = faces[faceI];
 369
 370         if
 371         (
 372             f.size() > minSize
 373          && cells[faceOwner[faceI]].size() >= 6
 374          && (
 375                 mesh.isInternalFace(faceI)
 376              && cells[faceNeighbour[faceI]].size() >= 6
 377             )
 378         )
 379         {
 380             // Get the edges in face point order
 381             labelList fEdges(getSortedEdges(edges, f, faceEdges[faceI]));
 382
 383             SortableList<scalar> lengths(fEdges.size());
 384             forAll(fEdges, i)
 385             {
 386                 lengths[i] = edges[fEdges[i]].mag(points);
 387             }
 388             lengths.sort();
 389
 390
 391             // Now find a gap in length between consecutive elements greater
 392             // than lenGap.
 393
 394             label gapPos = -1;
 395
 396             for (label i = f.size()-1-minSize; i >= 0; --i)
 397             {
 398                 if (lengths[i+1] > lenGap*lengths[i])
 399                 {
 400                     gapPos = i;
 401
 402                     break;
 403                 }
 404             }
 405
 406             if (gapPos != -1)
 407             {
 408                 //for (label i = gapPos; i >= 0; --i)
 409                 label i = 0;  // Hack: collapse smallest edge only.
 410                 {
 411                     label edgeI = fEdges[lengths.indices()[i]];
 412
 413                     if (!protectedEdge[edgeI])
 414                     {
 415                         const edge& e = edges[edgeI];
 416
 417                         label master = edgeMaster(boundaryPoint, e);
 418
 419                         if (master != -1)
 420                         {
 421                             collapser.collapseEdge(edgeI, master);
 422
 423                             // Protect all other edges on all cells using edge
 424                             // points.
 425
 426                             const labelList& pCells0 = pointCells[e[0]];
 427
 428                             forAll(pCells0, i)
 429                             {
 430                                 set(cellEdges[pCells0[i]], true, protectedEdge);
 431                             }
 432                             const labelList& pCells1 = pointCells[e[1]];
 433
 434                             forAll(pCells1, i)
 435                             {
 436                                 set(cellEdges[pCells1[i]], true, protectedEdge);
 437                             }
 438
 439                             nCollapsed++;
 440                         }
 441                     }
 442                 }
 443             }
 444         }
 445     }
 446
 447     return nCollapsed;
 448 }
 449
 450
 451 // Main program:
 452
 453 int main(int argc, char *argv[])
 454 {
 455     argList::noParallel();
 456     argList::validOptions.insert("overwrite", "");
 457     argList::validArgs.append("edge length [m]");
 458     argList::validArgs.append("merge angle (degrees)");
 459
 460 #   include "setRootCase.H"
 461 #   include "createTime.H"
 462     runTime.functionObjects().off();
 463 #   include "createPolyMesh.H"
 464     const word oldInstance = mesh.pointsInstance();
 465
 466     scalar minLen(readScalar(IStringStream(args.additionalArgs()[0])()));
 467     scalar angle(readScalar(IStringStream(args.additionalArgs()[1])()));
 468     bool overwrite = args.optionFound("overwrite");
 469
 470     scalar maxCos = Foam::cos(angle*180/mathematicalConstant::pi);
 471
 472     Info<< "Merging:" << nl
 473         << "    edges with length less than " << minLen << " meters" << nl
 474         << "    edges split by a point with edges in line to within " << angle
 475         << " degrees" << nl
 476         << endl;
 477
 478
 479     bool meshChanged = false;
 480
 481     while (true)
 482     {
 483         const faceList& faces = mesh.faces();
 484
 485         // Get all points on the boundary
 486         PackedBoolList boundaryPoint(mesh.nPoints());
 487
 488         label nIntFaces = mesh.nInternalFaces();
 489         for (label faceI = nIntFaces; faceI < mesh.nFaces(); faceI++)
 490         {
 491             const face& f = faces[faceI];
 492
 493             forAll(f, fp)
 494             {
 495                 boundaryPoint.set(f[fp], 1);
 496             }
 497         }
 498
 499         // Edge collapsing engine
 500         edgeCollapser collapser(mesh);
 501
 502
 503         // Collapse all edges that are too small.
 504         label nCollapsed =
 505             collapseSmallEdges
 506             (
 507                 mesh,
 508                 boundaryPoint,
 509                 minLen,
 510                 collapser
 511             );
 512         Info<< "Collapsing " << nCollapsed << " small edges" << endl;
 513
 514
 515         // Remove midpoints on straight edges.
 516         if (nCollapsed == 0)
 517         {
 518             nCollapsed = mergeEdges(mesh, maxCos, collapser);
 519             Info<< "Collapsing " << nCollapsed << " in line edges" << endl;
 520         }
 521
 522
 523         // Remove small sliver faces that can be collapsed to single edge
 524         if (nCollapsed == 0)
 525         {
 526             nCollapsed =
 527                 collapseHighAspectFaces
 528                 (
 529                     mesh,
 530                     boundaryPoint,
 531                     1E-9,       // factor of largest face area
 532                     5,          // factor between smallest and largest edge on
 533                                 // face
 534                     collapser
 535                 );
 536             Info<< "Collapsing " << nCollapsed
 537                 << " small high aspect ratio faces" << endl;
 538         }
 539
 540         // Simplify faces to quads wherever possible
 541         //if (nCollapsed == 0)
 542         //{
 543         //    nCollapsed =
 544         //        simplifyFaces
 545         //        (
 546         //            mesh,
 547         //            boundaryPoint,
 548         //            4,              // minimum size of face
 549         //            0.2,            // gap in edge lengths on face
 550         //            collapser
 551         //        );
 552         //    Info<< "Collapsing " << nCollapsed << " polygonal faces" << endl;
 553         //}
 554
 555
 556         if (nCollapsed == 0)
 557         {
 558             break;
 559         }
 560
 561         polyTopoChange meshMod(mesh);
 562
 563         // Insert mesh refinement into polyTopoChange.
 564         collapser.setRefinement(meshMod);
 565
 566         // Do all changes
 567         Info<< "Morphing ..." << endl;
 568
 569         autoPtr<mapPolyMesh> morphMap = meshMod.changeMesh(mesh, false);
 570
 571         collapser.updateMesh(morphMap());
 572
 573         if (morphMap().hasMotionPoints())
 574         {
 575             mesh.movePoints(morphMap().preMotionPoints());
 576         }
 577
 578         meshChanged = true;
 579     }
 580
 581     if (meshChanged)
 582     {
 583         // Write resulting mesh
 584         if (!overwrite)
 585         {
 586             runTime++;
 587         }
 588         else
 589         {
 590             mesh.setInstance(oldInstance);
 591         }
 592
 593         Info<< "Writing collapsed mesh to time " << runTime.timeName() << endl;
 594
 595         mesh.write();
 596     }
 597
 598     Info << "End\n" << endl;
 599
 600     return 0;
 601 }
 602
 603
 604 // ************************************************************************* //