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