src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.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 \*---------------------------------------------------------------------------*/
  26
  27 #include "edgeCollapser.H"
  28 #include "polyMesh.H"
  29 #include "polyTopoChange.H"
  30 #include "ListOps.H"
  31
  32 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  33
  34 Foam::label Foam::edgeCollapser::findIndex
  35 (
  36     const labelList& elems,
  37     const label start,
  38     const label size,
  39     const label val
  40 )
  41 {
  42     for (label i = start; i < size; i++)
  43     {
  44         if (elems[i] == val)
  45         {
  46             return i;
  47         }
  48     }
  49     return -1;
  50 }
  51
  52
  53 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  54
  55 // Changes region of connected set of points
  56 Foam::label Foam::edgeCollapser::changePointRegion
  57 (
  58     const label pointI,
  59     const label oldRegion,
  60     const label newRegion
  61 )
  62 {
  63     label nChanged = 0;
  64
  65     if (pointRegion_[pointI] == oldRegion)
  66     {
  67         pointRegion_[pointI] = newRegion;
  68         nChanged++;
  69
  70         // Step to neighbouring point across edges with same region number
  71
  72         const labelList& pEdges = mesh_.pointEdges()[pointI];
  73
  74         forAll(pEdges, i)
  75         {
  76             label otherPointI = mesh_.edges()[pEdges[i]].otherVertex(pointI);
  77
  78             nChanged += changePointRegion(otherPointI, oldRegion, newRegion);
  79         }
  80     }
  81     return nChanged;
  82 }
  83
  84
  85 bool Foam::edgeCollapser::pointRemoved(const label pointI) const
  86 {
  87     label region = pointRegion_[pointI];
  88
  89     if (region == -1 || pointRegionMaster_[region] == pointI)
  90     {
  91         return false;
  92     }
  93     else
  94     {
  95         return true;
  96     }
  97 }
  98
  99
 100 void Foam::edgeCollapser::filterFace(const label faceI, face& f) const
 101 {
 102     label newFp = 0;
 103
 104     forAll(f, fp)
 105     {
 106         label pointI = f[fp];
 107
 108         label region = pointRegion_[pointI];
 109
 110         if (region == -1)
 111         {
 112             f[newFp++] = pointI;
 113         }
 114         else
 115         {
 116             label master = pointRegionMaster_[region];
 117
 118             if (findIndex(f, 0, newFp, master) == -1)
 119             {
 120                 f[newFp++] = master;
 121             }
 122         }
 123     }
 124
 125     // Check for pinched face. Tries to correct
 126     // - consecutive duplicate vertex. Removes duplicate vertex.
 127     // - duplicate vertex with one other vertex in between (spike).
 128     // Both of these should not really occur! and should be checked before
 129     // collapsing edges.
 130
 131     const label size = newFp;
 132
 133     newFp = 2;
 134
 135     for (label fp = 2; fp < size; fp++)
 136     {
 137         label fp1 = fp-1;
 138         label fp2 = fp-2;
 139
 140         label pointI = f[fp];
 141
 142         // Search for previous occurrence.
 143         label index = findIndex(f, 0, fp, pointI);
 144
 145         if (index == fp1)
 146         {
 147             WarningIn
 148             (
 149                 "Foam::edgeCollapser::filterFace(const label faceI, "
 150                 "face& f) const"
 151             )   << "Removing consecutive duplicate vertex in face "
 152                 << f << endl;
 153             // Don't store current pointI
 154         }
 155         else if (index == fp2)
 156         {
 157             WarningIn
 158             (
 159                 "Foam::edgeCollapser::filterFace(const label faceI, "
 160                 "face& f) const"
 161             )   << "Removing non-consecutive duplicate vertex in face "
 162                 << f << endl;
 163             // Don't store current pointI and remove previous
 164             newFp--;
 165         }
 166         else if (index != -1)
 167         {
 168             WarningIn
 169             (
 170                 "Foam::edgeCollapser::filterFace(const label faceI, "
 171                 "face& f) const"
 172             )   << "Pinched face " << f << endl;
 173             f[newFp++] = pointI;
 174         }
 175         else
 176         {
 177             f[newFp++] = pointI;
 178         }
 179     }
 180
 181     f.setSize(newFp);
 182 }
 183
 184
 185 // Debugging.
 186 void Foam::edgeCollapser::printRegions() const
 187 {
 188     forAll(pointRegionMaster_, regionI)
 189     {
 190         label master = pointRegionMaster_[regionI];
 191
 192         if (master != -1)
 193         {
 194             Info<< "Region:" << regionI << nl
 195                 << "    master:" << master
 196                 << ' ' << mesh_.points()[master] << nl;
 197
 198             forAll(pointRegion_, pointI)
 199             {
 200                 if (pointRegion_[pointI] == regionI && pointI != master)
 201                 {
 202                     Info<< "    slave:" << pointI
 203                         << ' ' <<  mesh_.points()[pointI] << nl;
 204                 }
 205             }
 206         }
 207     }
 208 }
 209
 210
 211
 212 // Collapse list of edges
 213 void Foam::edgeCollapser::collapseEdges(const labelList& edgeLabels)
 214 {
 215     const edgeList& edges = mesh_.edges();
 216
 217     forAll(edgeLabels, i)
 218     {
 219         label edgeI = edgeLabels[i];
 220         const edge& e = edges[edgeI];
 221
 222         label region0 = pointRegion_[e[0]];
 223         label region1 = pointRegion_[e[1]];
 224
 225         if (region0 == -1)
 226         {
 227             if (region1 == -1)
 228             {
 229                 // Both unaffected. Choose ad lib.
 230                 collapseEdge(edgeI, e[0]);
 231             }
 232             else
 233             {
 234                 // Collapse to whatever e[1] collapses
 235                 collapseEdge(edgeI, e[1]);
 236             }
 237         }
 238         else
 239         {
 240             if (region1 == -1)
 241             {
 242                 // Collapse to whatever e[0] collapses
 243                 collapseEdge(edgeI, e[0]);
 244             }
 245             else
 246             {
 247                 // Both collapsed.
 248                 if (pointRegionMaster_[region0] == e[0])
 249                 {
 250                     // e[0] is a master
 251                     collapseEdge(edgeI, e[0]);
 252                 }
 253                 else if (pointRegionMaster_[region1] == e[1])
 254                 {
 255                     // e[1] is a master
 256                     collapseEdge(edgeI, e[1]);
 257                 }
 258                 else
 259                 {
 260                     // Dont know
 261                     collapseEdge(edgeI, e[0]);
 262                 }
 263             }
 264         }
 265     }
 266 }
 267
 268
 269 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 270
 271 // Construct from mesh
 272 Foam::edgeCollapser::edgeCollapser(const polyMesh& mesh)
 273 :
 274     mesh_(mesh),
 275     pointRegion_(mesh.nPoints(), -1),
 276     pointRegionMaster_(mesh.nPoints() / 100),
 277     freeRegions_()
 278 {}
 279
 280
 281 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 282
 283 bool Foam::edgeCollapser::unaffectedEdge(const label edgeI) const
 284 {
 285     const edge& e = mesh_.edges()[edgeI];
 286
 287     return (pointRegion_[e[0]] == -1) && (pointRegion_[e[1]] == -1);
 288 }
 289
 290
 291 bool Foam::edgeCollapser::collapseEdge(const label edgeI, const label master)
 292 {
 293     const edge& e = mesh_.edges()[edgeI];
 294
 295     label pointRegion0 = pointRegion_[e[0]];
 296     label pointRegion1 = pointRegion_[e[1]];
 297
 298     if (pointRegion0 == -1)
 299     {
 300         if (pointRegion1 == -1)
 301         {
 302             // Both endpoints not collapsed. Create new region.
 303
 304             label freeRegion = -1;
 305
 306             if (freeRegions_.size())
 307             {
 308                 freeRegion = freeRegions_.removeHead();
 309
 310                 if (pointRegionMaster_[freeRegion] != -1)
 311                 {
 312                     FatalErrorIn
 313                     ("edgeCollapser::collapseEdge(const label, const label)")
 314                         << "Problem : freeed region :" << freeRegion
 315                         << " has already master "
 316                         << pointRegionMaster_[freeRegion]
 317                         << abort(FatalError);
 318                 }
 319             }
 320             else
 321             {
 322                 // If no region found create one. This is the only place where
 323                 // new regions are created.
 324                 freeRegion = pointRegionMaster_.size();
 325             }
 326
 327             pointRegion_[e[0]] = freeRegion;
 328             pointRegion_[e[1]] = freeRegion;
 329
 330             pointRegionMaster_(freeRegion) = master;
 331         }
 332         else
 333         {
 334             // e[1] is part of collapse network, e[0] not. Add e0 to e1 region.
 335             pointRegion_[e[0]] = pointRegion1;
 336
 337             pointRegionMaster_[pointRegion1] = master;
 338         }
 339     }
 340     else
 341     {
 342         if (pointRegion1 == -1)
 343         {
 344             // e[0] is part of collapse network. Add e1 to e0 region
 345             pointRegion_[e[1]] = pointRegion0;
 346
 347             pointRegionMaster_[pointRegion0] = master;
 348         }
 349         else if (pointRegion0 != pointRegion1)
 350         {
 351             // Both part of collapse network. Merge the two regions.
 352
 353             // Use the smaller region number for the whole network.
 354             label minRegion = min(pointRegion0, pointRegion1);
 355             label maxRegion = max(pointRegion0, pointRegion1);
 356
 357             // Use minRegion as region for combined net, free maxRegion.
 358             pointRegionMaster_[minRegion] = master;
 359             pointRegionMaster_[maxRegion] = -1;
 360             freeRegions_.insert(maxRegion);
 361
 362             if (minRegion != pointRegion0)
 363             {
 364                 changePointRegion(e[0], pointRegion0, minRegion);
 365             }
 366             if (minRegion != pointRegion1)
 367             {
 368                 changePointRegion(e[1], pointRegion1, minRegion);
 369             }
 370         }
 371     }
 372
 373     return true;
 374 }
 375
 376
 377 bool Foam::edgeCollapser::setRefinement(polyTopoChange& meshMod)
 378 {
 379     const cellList& cells = mesh_.cells();
 380     const labelList& faceOwner = mesh_.faceOwner();
 381     const labelList& faceNeighbour = mesh_.faceNeighbour();
 382     const labelListList& pointFaces = mesh_.pointFaces();
 383     const labelListList& cellEdges = mesh_.cellEdges();
 384
 385     // Print regions:
 386     //printRegions()
 387
 388     bool meshChanged = false;
 389
 390
 391     // Current faces (is also collapseStatus: f.size() < 3)
 392     faceList newFaces(mesh_.faces());
 393
 394     // Current cellCollapse status
 395     boolList cellRemoved(mesh_.nCells(), false);
 396
 397
 398     do
 399     {
 400         // Update face collapse from edge collapses
 401         forAll(newFaces, faceI)
 402         {
 403             filterFace(faceI, newFaces[faceI]);
 404         }
 405
 406         // Check if faces to be collapsed cause cells to become collapsed.
 407         label nCellCollapsed = 0;
 408
 409         forAll(cells, cellI)
 410         {
 411             if (!cellRemoved[cellI])
 412             {
 413                 const cell& cFaces = cells[cellI];
 414
 415                 label nFaces = cFaces.size();
 416
 417                 forAll(cFaces, i)
 418                 {
 419                     label faceI = cFaces[i];
 420
 421                     if (newFaces[faceI].size() < 3)
 422                     {
 423                         --nFaces;
 424
 425                         if (nFaces < 4)
 426                         {
 427                             Info<< "Cell:" << cellI
 428                                 << " uses faces:" << cFaces
 429                                 << " of which too many are marked for removal:"
 430                                 << endl
 431                                 << "   ";
 432                             forAll(cFaces, j)
 433                             {
 434                                 if (newFaces[cFaces[j]].size() < 3)
 435                                 {
 436                                     Info<< ' '<< cFaces[j];
 437                                 }
 438                             }
 439                             Info<< endl;
 440
 441                             cellRemoved[cellI] = true;
 442
 443                             // Collapse all edges of cell to nothing
 444                             collapseEdges(cellEdges[cellI]);
 445
 446                             nCellCollapsed++;
 447
 448                             break;
 449                         }
 450                     }
 451                 }
 452             }
 453         }
 454
 455         if (nCellCollapsed == 0)
 456         {
 457             break;
 458         }
 459     } while(true);
 460
 461
 462     // Keep track of faces that have been done already.
 463     boolList doneFace(mesh_.nFaces(), false);
 464
 465     {
 466         // Mark points used.
 467         boolList usedPoint(mesh_.nPoints(), false);
 468
 469
 470         forAll(cellRemoved, cellI)
 471         {
 472             if (cellRemoved[cellI])
 473             {
 474                 meshMod.removeCell(cellI, -1);
 475             }
 476         }
 477
 478
 479         // Remove faces
 480         forAll(newFaces, faceI)
 481         {
 482             const face& f = newFaces[faceI];
 483
 484             if (f.size() < 3)
 485             {
 486                 meshMod.removeFace(faceI, -1);
 487                 meshChanged = true;
 488
 489                 // Mark face as been done.
 490                 doneFace[faceI] = true;
 491             }
 492             else
 493             {
 494                 // Kept face. Mark vertices
 495                 forAll(f, fp)
 496                 {
 497                     usedPoint[f[fp]] = true;
 498                 }
 499             }
 500         }
 501
 502         // Remove unused vertices that have not been marked for removal already
 503         forAll(usedPoint, pointI)
 504         {
 505             if (!usedPoint[pointI] && !pointRemoved(pointI))
 506             {
 507                 meshMod.removePoint(pointI, -1);
 508                 meshChanged = true;
 509             }
 510         }
 511     }
 512
 513
 514
 515     // Remove points.
 516     forAll(pointRegion_, pointI)
 517     {
 518         if (pointRemoved(pointI))
 519         {
 520             meshMod.removePoint(pointI, -1);
 521             meshChanged = true;
 522         }
 523     }
 524
 525
 526
 527     const polyBoundaryMesh& boundaryMesh = mesh_.boundaryMesh();
 528     const faceZoneMesh& faceZones = mesh_.faceZones();
 529
 530
 531     // Renumber faces that use points
 532     forAll(pointRegion_, pointI)
 533     {
 534         if (pointRemoved(pointI))
 535         {
 536             const labelList& changedFaces = pointFaces[pointI];
 537
 538             forAll(changedFaces, changedFaceI)
 539             {
 540                 label faceI = changedFaces[changedFaceI];
 541
 542                 if (!doneFace[faceI])
 543                 {
 544                     doneFace[faceI] = true;
 545
 546                     // Get current zone info
 547                     label zoneID = faceZones.whichZone(faceI);
 548
 549                     bool zoneFlip = false;
 550
 551                     if (zoneID >= 0)
 552                     {
 553                         const faceZone& fZone = faceZones[zoneID];
 554
 555                         zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)];
 556                     }
 557
 558                     // Get current connectivity
 559                     label own = faceOwner[faceI];
 560                     label nei = -1;
 561                     label patchID = -1;
 562
 563                     if (mesh_.isInternalFace(faceI))
 564                     {
 565                         nei = faceNeighbour[faceI];
 566                     }
 567                     else
 568                     {
 569                         patchID = boundaryMesh.whichPatch(faceI);
 570                     }
 571
 572                     meshMod.modifyFace
 573                     (
 574                         newFaces[faceI],            // face
 575                         faceI,                      // faceI to change
 576                         own,                        // owner
 577                         nei,                        // neighbour
 578                         false,                      // flipFaceFlux
 579                         patchID,                    // patch
 580                         zoneID,
 581                         zoneFlip
 582                     );
 583                     meshChanged = true;
 584                 }
 585             }
 586         }
 587     }
 588
 589     return meshChanged;
 590 }
 591
 592
 593 void Foam::edgeCollapser::updateMesh(const mapPolyMesh& map)
 594 {
 595     pointRegion_.setSize(mesh_.nPoints());
 596     pointRegion_ = -1;
 597     // Reset count, do not remove underlying storage
 598     pointRegionMaster_.clear();
 599     freeRegions_.clear();
 600 }
 601
 602
 603 // ************************************************************************* //