src/dynamicMesh/polyTopoChange/polyTopoChange/localPointRegion.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 "localPointRegion.H"
  28 #include "syncTools.H"
  29 #include "polyMesh.H"
  30 #include "mapPolyMesh.H"
  31 #include "globalIndex.H"
  32 #include "indirectPrimitivePatch.H"
  33
  34 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  35
  36 namespace Foam
  37 {
  38
  39 defineTypeNameAndDebug(localPointRegion, 0);
  40
  41 // Reduction class to get minimum value over face.
  42 class minEqOpFace
  43 {
  44 public:
  45
  46     void operator()(face& x, const face& y) const
  47     {
  48         if (x.size())
  49         {
  50             label j = 0;
  51             forAll(x, i)
  52             {
  53                 x[i] = min(x[i], y[j]);
  54
  55                 j = y.rcIndex(j);
  56             }
  57         }
  58     };
  59 };
  60
  61
  62 // Dummy transform for faces. Used in synchronisation
  63 void transformList
  64 (
  65     const tensorField& rotTensor,
  66     UList<face>& field
  67 )
  68 {};
  69
  70 }
  71
  72 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  73
  74 // Are two lists identical either in forward or in reverse order.
  75 bool Foam::localPointRegion::isDuplicate
  76 (
  77     const face& f0,
  78     const face& f1,
  79     const bool forward
  80 )
  81 {
  82     label fp1 = findIndex(f1, f0[0]);
  83
  84     if (fp1 == -1)
  85     {
  86         return false;
  87     }
  88
  89     forAll(f0, fp0)
  90     {
  91         if (f0[fp0] != f1[fp1])
  92         {
  93             return false;
  94         }
  95
  96         if (forward)
  97         {
  98             fp1 = f1.fcIndex(fp1);
  99         }
 100         else
 101         {
 102             fp1 = f1.rcIndex(fp1);
 103         }
 104     }
 105     return true;
 106 }
 107
 108
 109 // Count regions per point
 110 void Foam::localPointRegion::countPointRegions
 111 (
 112     const polyMesh& mesh,
 113     const boolList& candidatePoint,
 114     const Map<label>& candidateFace,
 115     faceList& minRegion
 116 )
 117 {
 118     // Almost all will have only one so only
 119     // populate Map if more than one.
 120     labelList minPointRegion(mesh.nPoints(), -1);
 121     // From global point to local (multi-region) point numbering
 122     meshPointMap_.resize(candidateFace.size()/100);
 123     // From local (multi-region) point to regions
 124     DynamicList<labelList> pointRegions(meshPointMap_.size());
 125
 126     // From faces with any duplicated point on it to local face
 127     meshFaceMap_.resize(meshPointMap_.size());
 128
 129     forAllConstIter(Map<label>, candidateFace, iter)
 130     {
 131         label faceI = iter.key();
 132
 133         if (!mesh.isInternalFace(faceI))
 134         {
 135             const face& f = mesh.faces()[faceI];
 136
 137             if (minRegion[faceI].empty())
 138             {
 139                 FatalErrorIn("localPointRegion::countPointRegions(..)")
 140                     << "Face from candidateFace without minRegion set." << endl
 141                     << "Face:" << faceI << " fc:" << mesh.faceCentres()[faceI]
 142                     << " verts:" << f << abort(FatalError);
 143             }
 144
 145             forAll(f, fp)
 146             {
 147                 label pointI = f[fp];
 148
 149                 // Even points which were not candidates for splitting might
 150                 // be on multiple baffles that are being split so check.
 151
 152                 if (candidatePoint[pointI])
 153                 {
 154                     label region = minRegion[faceI][fp];
 155
 156                     if (minPointRegion[pointI] == -1)
 157                     {
 158                         minPointRegion[pointI] = region;
 159                     }
 160                     else if (minPointRegion[pointI] != region)
 161                     {
 162                         // Multiple regions for this point. Add.
 163                         Map<label>::iterator iter = meshPointMap_.find(pointI);
 164                         if (iter != meshPointMap_.end())
 165                         {
 166                             labelList& regions = pointRegions[iter()];
 167                             if (findIndex(regions, region) == -1)
 168                             {
 169                                 label sz = regions.size();
 170                                 regions.setSize(sz+1);
 171                                 regions[sz] = region;
 172                             }
 173                         }
 174                         else
 175                         {
 176                             label localPointI = meshPointMap_.size();
 177                             meshPointMap_.insert(pointI, localPointI);
 178                             labelList regions(2);
 179                             regions[0] = minPointRegion[pointI];
 180                             regions[1] = region;
 181                             pointRegions.append(regions);
 182                         }
 183
 184                         label meshFaceMapI = meshFaceMap_.size();
 185                         meshFaceMap_.insert(faceI, meshFaceMapI);
 186                     }
 187                 }
 188             }
 189         }
 190     }
 191     minPointRegion.clear();
 192
 193     // Add internal faces that use any duplicated point. Can only have one
 194     // region!
 195     forAllConstIter(Map<label>, candidateFace, iter)
 196     {
 197         label faceI = iter.key();
 198
 199         if (mesh.isInternalFace(faceI))
 200         {
 201             const face& f = mesh.faces()[faceI];
 202
 203             forAll(f, fp)
 204             {
 205                 // Note: candidatePoint test not really necessary but
 206                 // speeds up rejection.
 207                 if (candidatePoint[f[fp]] && meshPointMap_.found(f[fp]))
 208                 {
 209                     label meshFaceMapI = meshFaceMap_.size();
 210                     meshFaceMap_.insert(faceI, meshFaceMapI);
 211                 }
 212             }
 213         }
 214     }
 215
 216
 217     // Transfer to member data
 218     pointRegions.shrink();
 219     pointRegions_.setSize(pointRegions.size());
 220     forAll(pointRegions, i)
 221     {
 222         pointRegions_[i].transfer(pointRegions[i]);
 223     }
 224
 225     // Compact minRegion
 226     faceRegions_.setSize(meshFaceMap_.size());
 227     forAllConstIter(Map<label>, meshFaceMap_, iter)
 228     {
 229         faceRegions_[iter()].labelList::transfer(minRegion[iter.key()]);
 230
 231         //// Print a bit
 232         //{
 233         //    label faceI = iter.key();
 234         //    const face& f = mesh.faces()[faceI];
 235         //    Pout<< "Face:" << faceI << " fc:" << mesh.faceCentres()[faceI]
 236         //        << " verts:" << f << endl;
 237         //    forAll(f, fp)
 238         //    {
 239         //        Pout<< "    " << f[fp] << " min:" << faceRegions_[iter()][fp]
 240         //            << endl;
 241         //    }
 242         //    Pout<< endl;
 243         //}
 244     }
 245
 246     // Compact region numbering
 247     // ? TBD.
 248 }
 249
 250
 251 void Foam::localPointRegion::calcPointRegions
 252 (
 253     const polyMesh& mesh,
 254     boolList& candidatePoint
 255 )
 256 {
 257     label nBnd = mesh.nFaces()-mesh.nInternalFaces();
 258     const labelList& faceOwner = mesh.faceOwner();
 259     const labelList& faceNeighbour = mesh.faceNeighbour();
 260
 261
 262     syncTools::syncPointList
 263     (
 264         mesh,
 265         candidatePoint,
 266         orEqOp<bool>(),
 267         false,              // nullValue
 268         false               // applySeparation
 269     );
 270
 271
 272     // Mark any face/boundaryFace/cell with a point on a candidate point.
 273     // - candidateFace does not necessary have to be a baffle!
 274     // - candidateFace is synchronised (since candidatePoint is)
 275     Map<label> candidateFace(2*nBnd);
 276     label candidateFaceI = 0;
 277
 278     Map<label> candidateCell(nBnd);
 279     label candidateCellI = 0;
 280
 281     forAll(mesh.faces(), faceI)
 282     {
 283         const face& f = mesh.faces()[faceI];
 284
 285         forAll(f, fp)
 286         {
 287             if (candidatePoint[f[fp]])
 288             {
 289                 // Mark face
 290                 if (candidateFace.insert(faceI, candidateFaceI))
 291                 {
 292                     candidateFaceI++;
 293                 }
 294
 295                 // Mark cells
 296                 if (candidateCell.insert(faceOwner[faceI], candidateCellI))
 297                 {
 298                     candidateCellI++;
 299                 }
 300
 301                 if (mesh.isInternalFace(faceI))
 302                 {
 303                     label nei = faceNeighbour[faceI];
 304                     if (candidateCell.insert(nei, candidateCellI))
 305                     {
 306                         candidateCellI++;
 307                     }
 308                 }
 309
 310                 break;
 311             }
 312         }
 313     }
 314
 315
 316
 317     // Get global indices for cells
 318     globalIndex globalCells(mesh.nCells());
 319
 320
 321     // Determine for every candidate face per point the minimum region
 322     // (global cell) it is connected to. (candidateFaces are the
 323     // only ones using a
 324     // candidate point so the only ones that can be affected)
 325     faceList minRegion(mesh.nFaces());
 326     forAllConstIter(Map<label>, candidateFace, iter)
 327     {
 328         label faceI = iter.key();
 329         const face& f = mesh.faces()[faceI];
 330
 331         if (mesh.isInternalFace(faceI))
 332         {
 333             label globOwn = globalCells.toGlobal(faceOwner[faceI]);
 334             label globNei = globalCells.toGlobal(faceNeighbour[faceI]);
 335             minRegion[faceI].setSize(f.size(), min(globOwn, globNei));
 336         }
 337         else
 338         {
 339             label globOwn = globalCells.toGlobal(faceOwner[faceI]);
 340             minRegion[faceI].setSize(f.size(), globOwn);
 341         }
 342     }
 343
 344     // Now minimize over all faces that are connected through internal
 345     // faces or through cells. This loop iterates over the max number of
 346     // cells connected to a point (=8 for hex mesh)
 347
 348     while (true)
 349     {
 350         // Transport minimum from face across cell
 351         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 352
 353         Map<label> minPointValue(100);
 354         label nChanged = 0;
 355         forAllConstIter(Map<label>, candidateCell, iter)
 356         {
 357             minPointValue.clear();
 358
 359             label cellI = iter.key();
 360             const cell& cFaces = mesh.cells()[cellI];
 361
 362             // Determine minimum per point
 363             forAll(cFaces, cFaceI)
 364             {
 365                 label faceI = cFaces[cFaceI];
 366
 367                 if (minRegion[faceI].size())
 368                 {
 369                     const face& f = mesh.faces()[faceI];
 370
 371                     forAll(f, fp)
 372                     {
 373                         label pointI = f[fp];
 374                         Map<label>::iterator iter = minPointValue.find(pointI);
 375
 376                         if (iter == minPointValue.end())
 377                         {
 378                             minPointValue.insert(pointI, minRegion[faceI][fp]);
 379                         }
 380                         else
 381                         {
 382                             label currentMin = iter();
 383                             iter() = min(currentMin, minRegion[faceI][fp]);
 384                         }
 385                     }
 386                 }
 387             }
 388
 389             // Set face minimum from point minimum
 390             forAll(cFaces, cFaceI)
 391             {
 392                 label faceI = cFaces[cFaceI];
 393
 394                 if (minRegion[faceI].size())
 395                 {
 396                     const face& f = mesh.faces()[faceI];
 397
 398                     forAll(f, fp)
 399                     {
 400                         label minVal = minPointValue[f[fp]];
 401
 402                         if (minVal != minRegion[faceI][fp])
 403                         {
 404                             minRegion[faceI][fp] = minVal;
 405                             nChanged++;
 406                         }
 407                     }
 408                 }
 409             }
 410         }
 411
 412         //Pout<< "nChanged:" << nChanged << endl;
 413
 414         if (returnReduce(nChanged, sumOp<label>()) == 0)
 415         {
 416             break;
 417         }
 418
 419
 420         // Transport minimum across coupled faces
 421         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 422
 423         syncTools::syncFaceList
 424         (
 425             mesh,
 426             minRegion,
 427             minEqOpFace(),
 428             false               // applySeparation
 429         );
 430     }
 431
 432
 433     // Count regions per point
 434     countPointRegions(mesh, candidatePoint, candidateFace, minRegion);
 435     minRegion.clear();
 436
 437
 438     //// Print points with multiple regions. These points need to be duplicated.
 439     //forAllConstIter(Map<label>, meshPointMap_, iter)
 440     //{
 441     //    Pout<< "point:" << iter.key()
 442     //        << " coord:" << mesh.points()[iter.key()]
 443     //        << " regions:" << pointRegions_[iter()] << endl;
 444     //}
 445 }
 446
 447
 448 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 449
 450 Foam::localPointRegion::localPointRegion(const polyMesh& mesh)
 451 :
 452     //nRegions_(0),
 453     meshPointMap_(0),
 454     pointRegions_(0),
 455     meshFaceMap_(0),
 456     faceRegions_(0)
 457 {
 458     const polyBoundaryMesh& patches = mesh.boundaryMesh();
 459
 460     // Get any point on the outside which is on a non-coupled boundary
 461     boolList candidatePoint(mesh.nPoints(), false);
 462
 463     forAll(patches, patchI)
 464     {
 465         if (!patches[patchI].coupled())
 466         {
 467             const polyPatch& pp = patches[patchI];
 468
 469             forAll(pp.meshPoints(), i)
 470             {
 471                 candidatePoint[pp.meshPoints()[i]] = true;
 472             }
 473         }
 474     }
 475
 476     calcPointRegions(mesh, candidatePoint);
 477 }
 478
 479
 480 Foam::localPointRegion::localPointRegion
 481 (
 482     const polyMesh& mesh,
 483     const labelList& candidatePoints
 484 )
 485 :
 486     //nRegions_(0),
 487     meshPointMap_(0),
 488     pointRegions_(0),
 489     meshFaceMap_(0),
 490     faceRegions_(0)
 491 {
 492     boolList candidatePoint(mesh.nPoints(), false);
 493
 494     forAll(candidatePoints, i)
 495     {
 496         candidatePoint[candidatePoints[i]] = true;
 497     }
 498
 499     calcPointRegions(mesh, candidatePoint);
 500 }
 501
 502
 503 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 504
 505 // Return a list (in allPatch indices) with either -1 or the face label
 506 // of the face that uses the same vertices.
 507 Foam::labelList Foam::localPointRegion::findDuplicateFaces
 508 (
 509     const primitiveMesh& mesh,
 510     const labelList& boundaryFaces
 511 )
 512 {
 513     // Addressing engine for all boundary faces.
 514     indirectPrimitivePatch allPatch
 515     (
 516         IndirectList<face>(mesh.faces(), boundaryFaces),
 517         mesh.points()
 518     );
 519
 520     labelList duplicateFace(allPatch.size(), -1);
 521     label nDuplicateFaces = 0;
 522
 523     // Find all duplicate faces.
 524     forAll(allPatch, bFaceI)
 525     {
 526         const face& f = allPatch.localFaces()[bFaceI];
 527
 528         // Get faces connected to f[0].
 529         // Check whether share all points with f.
 530         const labelList& pFaces = allPatch.pointFaces()[f[0]];
 531
 532         forAll(pFaces, i)
 533         {
 534             label otherFaceI = pFaces[i];
 535
 536             if (otherFaceI > bFaceI)
 537             {
 538                 const face& otherF = allPatch.localFaces()[otherFaceI];
 539
 540                 if (isDuplicate(f, otherF, true))
 541                 {
 542                     FatalErrorIn
 543                     (
 544                         "findDuplicateFaces(const primitiveMesh&"
 545                         ", const labelList&)"
 546                     )   << "Face:" << bFaceI + mesh.nInternalFaces()
 547                         << " has local points:" << f
 548                         << " which are in same order as face:"
 549                         << otherFaceI + mesh.nInternalFaces()
 550                         << " with local points:" << otherF
 551                         << abort(FatalError);
 552                 }
 553                 else if (isDuplicate(f, otherF, false))
 554                 {
 555                     label meshFace0 = bFaceI + mesh.nInternalFaces();
 556                     label meshFace1 = otherFaceI + mesh.nInternalFaces();
 557
 558                     if
 559                     (
 560                         duplicateFace[bFaceI] != -1
 561                      || duplicateFace[otherFaceI] != -1
 562                     )
 563                     {
 564                         FatalErrorIn
 565                         (
 566                             "findDuplicateFaces(const primitiveMesh&"
 567                             ", const labelList&)"
 568                         )   << "One of two duplicate faces already marked"
 569                             << " as duplicate." << nl
 570                             << "This means that three or more faces share"
 571                             << " the same points and this is illegal." << nl
 572                             << "Face:" << meshFace0
 573                             << " with local points:" << f
 574                             << " which are in same order as face:"
 575                             << meshFace1
 576                             << " with local points:" << otherF
 577                             << abort(FatalError);
 578                     }
 579
 580                     duplicateFace[bFaceI] = otherFaceI;
 581                     duplicateFace[otherFaceI] = bFaceI;
 582                     nDuplicateFaces++;
 583                 }
 584             }
 585         }
 586     }
 587
 588     return duplicateFace;
 589 }
 590
 591
 592 void Foam::localPointRegion::updateMesh(const mapPolyMesh& map)
 593 {
 594     {
 595         Map<label> newMap(meshFaceMap_.size());
 596
 597         forAllConstIter(Map<label>, meshFaceMap_, iter)
 598         {
 599             label newFaceI = map.reverseFaceMap()[iter.key()];
 600
 601             if (newFaceI >= 0)
 602             {
 603                 newMap.insert(newFaceI, iter());
 604             }
 605         }
 606         meshFaceMap_.transfer(newMap);
 607     }
 608     {
 609         Map<label> newMap(meshPointMap_.size());
 610
 611         forAllConstIter(Map<label>, meshPointMap_, iter)
 612         {
 613             label newPointI = map.reversePointMap()[iter.key()];
 614
 615             if (newPointI >= 0)
 616             {
 617                 newMap.insert(newPointI, iter());
 618             }
 619         }
 620
 621         meshPointMap_.transfer(newMap);
 622     }
 623 }
 624
 625
 626 // ************************************************************************* //