src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2010-2011 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
  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 \*---------------------------------------------------------------------------*/
  25
  26 #include "polyMeshTetDecomposition.H"
  27
  28 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  29
  30 // Note: the use of this tolerance is ad-hoc, there may be extreme
  31 // cases where the resulting tetrahedra still have particle tracking
  32 // problems, or tets with lower quality may track OK.
  33 const Foam::scalar Foam::polyMeshTetDecomposition::minTetQuality = sqr(SMALL);
  34
  35
  36 // * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
  37
  38 Foam::label Foam::polyMeshTetDecomposition::findSharedBasePoint
  39 (
  40     const polyMesh& mesh,
  41     label fI,
  42     const point& nCc,
  43     scalar tol,
  44     bool report
  45 )
  46 {
  47     const faceList& pFaces = mesh.faces();
  48     const pointField& pPts = mesh.points();
  49     const vectorField& pC = mesh.cellCentres();
  50     const labelList& pOwner = mesh.faceOwner();
  51
  52     const face& f = pFaces[fI];
  53
  54     label oCI = pOwner[fI];
  55
  56     const point& oCc = pC[oCI];
  57
  58     List<scalar> tetQualities(2, 0.0);
  59
  60     forAll(f, faceBasePtI)
  61     {
  62         scalar thisBaseMinTetQuality = VGREAT;
  63
  64         const point& tetBasePt = pPts[f[faceBasePtI]];
  65
  66         for (label tetPtI = 1; tetPtI < f.size() - 1; tetPtI++)
  67         {
  68             label facePtI = (tetPtI + faceBasePtI) % f.size();
  69             label otherFacePtI = f.fcIndex(facePtI);
  70
  71             {
  72                 // owner cell tet
  73                 label ptAI = f[facePtI];
  74                 label ptBI = f[otherFacePtI];
  75
  76                 const point& pA = pPts[ptAI];
  77                 const point& pB = pPts[ptBI];
  78
  79                 tetPointRef tet(oCc, tetBasePt, pA, pB);
  80
  81                 tetQualities[0] = tet.quality();
  82             }
  83
  84             {
  85                 // neighbour cell tet
  86                 label ptAI = f[otherFacePtI];
  87                 label ptBI = f[facePtI];
  88
  89                 const point& pA = pPts[ptAI];
  90                 const point& pB = pPts[ptBI];
  91
  92                 tetPointRef tet(nCc, tetBasePt, pA, pB);
  93
  94                 tetQualities[1] = tet.quality();
  95             }
  96
  97             if (min(tetQualities) < thisBaseMinTetQuality)
  98             {
  99                 thisBaseMinTetQuality = min(tetQualities);
 100             }
 101         }
 102
 103         if (thisBaseMinTetQuality > tol)
 104         {
 105             return faceBasePtI;
 106         }
 107     }
 108
 109     // If a base point hasn't triggered a return by now, then there is
 110     // none that can produce a good decomposition
 111     return -1;
 112 }
 113
 114
 115 Foam::label Foam::polyMeshTetDecomposition::findSharedBasePoint
 116 (
 117     const polyMesh& mesh,
 118     label fI,
 119     scalar tol,
 120     bool report
 121 )
 122 {
 123     return findSharedBasePoint
 124     (
 125         mesh,
 126         fI,
 127         mesh.cellCentres()[mesh.faceNeighbour()[fI]],
 128         tol,
 129         report
 130     );
 131 }
 132
 133
 134 Foam::label Foam::polyMeshTetDecomposition::findBasePoint
 135 (
 136     const polyMesh& mesh,
 137     label fI,
 138     scalar tol,
 139     bool report
 140 )
 141 {
 142     const faceList& pFaces = mesh.faces();
 143     const pointField& pPts = mesh.points();
 144     const vectorField& pC = mesh.cellCentres();
 145     const labelList& pOwner = mesh.faceOwner();
 146
 147     const face& f = pFaces[fI];
 148
 149     label cI = pOwner[fI];
 150
 151     bool own = (pOwner[fI] == cI);
 152
 153     const point& cC = pC[cI];
 154
 155     forAll(f, faceBasePtI)
 156     {
 157         scalar thisBaseMinTetQuality = VGREAT;
 158
 159         const point& tetBasePt = pPts[f[faceBasePtI]];
 160
 161         for (label tetPtI = 1; tetPtI < f.size() - 1; tetPtI++)
 162         {
 163             label facePtI = (tetPtI + faceBasePtI) % f.size();
 164             label otherFacePtI = f.fcIndex(facePtI);
 165
 166             label ptAI = -1;
 167             label ptBI = -1;
 168
 169             if (own)
 170             {
 171                 ptAI = f[facePtI];
 172                 ptBI = f[otherFacePtI];
 173             }
 174             else
 175             {
 176                 ptAI = f[otherFacePtI];
 177                 ptBI = f[facePtI];
 178             }
 179
 180             const point& pA = pPts[ptAI];
 181             const point& pB = pPts[ptBI];
 182
 183             tetPointRef tet(cC, tetBasePt, pA, pB);
 184
 185             scalar tetQuality = tet.quality();
 186
 187             if (tetQuality < thisBaseMinTetQuality)
 188             {
 189                 thisBaseMinTetQuality = tetQuality;
 190             }
 191         }
 192
 193         if (thisBaseMinTetQuality > tol)
 194         {
 195             return faceBasePtI;
 196         }
 197     }
 198
 199     // If a base point hasn't triggered a return by now, then there is
 200     // non that can produce a good decomposition
 201     return -1;
 202 }
 203
 204
 205 Foam::labelList Foam::polyMeshTetDecomposition::findFaceBasePts
 206 (
 207     const polyMesh& mesh,
 208     scalar tol,
 209     bool report
 210 )
 211 {
 212     const labelList& pOwner = mesh.faceOwner();
 213
 214     // Find a suitable base point for each face, considering both
 215     // cells for interface faces or those on coupled patches
 216
 217     labelList tetBasePtIs(mesh.nFaces(), -1);
 218
 219     label nInternalFaces = mesh.nInternalFaces();
 220
 221     for (label fI = 0; fI < nInternalFaces; fI++)
 222     {
 223         tetBasePtIs[fI] = findSharedBasePoint(mesh, fI, tol, report);
 224     }
 225
 226     pointField neighbourCellCentres(mesh.nFaces() - nInternalFaces);
 227
 228     for(label faceI = nInternalFaces; faceI < mesh.nFaces(); faceI++)
 229     {
 230         neighbourCellCentres[faceI - nInternalFaces] =
 231             mesh.cellCentres()[pOwner[faceI]];
 232     }
 233
 234     syncTools::swapBoundaryFacePositions(mesh, neighbourCellCentres);
 235
 236     const polyBoundaryMesh& patches = mesh.boundaryMesh();
 237
 238     SubList<label> boundaryFaceTetBasePtIs
 239     (
 240         tetBasePtIs,
 241         mesh.nFaces() - nInternalFaces,
 242         nInternalFaces
 243     );
 244
 245     for
 246     (
 247         label fI = nInternalFaces, bFI = 0;
 248         fI < mesh.nFaces();
 249         fI++, bFI++
 250     )
 251     {
 252         label patchI =
 253             mesh.boundaryMesh().patchID()[bFI];
 254
 255         if (patches[patchI].coupled())
 256         {
 257             const coupledPolyPatch& pp =
 258                 refCast<const coupledPolyPatch>(patches[patchI]);
 259
 260             if (pp.owner())
 261             {
 262                 boundaryFaceTetBasePtIs[bFI] = findSharedBasePoint
 263                 (
 264                     mesh,
 265                     fI,
 266                     neighbourCellCentres[bFI],
 267                     tol,
 268                     report
 269                 );
 270             }
 271             else
 272             {
 273                 // Assign -2, to distinguish from a failed base point
 274                 // find, which returns -1.
 275                 boundaryFaceTetBasePtIs[bFI] = -2;
 276             }
 277         }
 278         else
 279         {
 280             boundaryFaceTetBasePtIs[bFI] = findBasePoint
 281             (
 282                 mesh,
 283                 fI,
 284                 tol,
 285                 report
 286             );
 287         }
 288     }
 289
 290     // maxEqOp will replace the -2 values on the neighbour patches
 291     // with the result from the owner base point find.
 292
 293     syncTools::syncBoundaryFaceList
 294     (
 295         mesh,
 296         boundaryFaceTetBasePtIs,
 297         maxEqOp<label>()
 298     );
 299
 300     for
 301     (
 302         label fI = nInternalFaces, bFI = 0;
 303         fI < mesh.nFaces();
 304         fI++, bFI++
 305     )
 306     {
 307         label& bFTetBasePtI = boundaryFaceTetBasePtIs[bFI];
 308
 309         if (bFTetBasePtI == -2)
 310         {
 311             FatalErrorIn
 312             (
 313                 "Foam::labelList"
 314                 "Foam::polyMeshTetDecomposition::findFaceBasePts"
 315                 "("
 316                     "const polyMesh& mesh, "
 317                     "scalar tol, "
 318                     "bool report"
 319                 ")"
 320             )
 321                 << "Coupled face base point exchange failure for face "
 322                 << fI
 323                 << abort(FatalError);
 324         }
 325
 326         if (bFTetBasePtI < 1)
 327         {
 328             // If the base point is -1, it should be left as such to
 329             // indicate a problem, if it is 0, then no action is required.
 330
 331             continue;
 332         }
 333
 334         label patchI = mesh.boundaryMesh().patchID()[bFI];
 335
 336         if (patches[patchI].coupled())
 337         {
 338             const coupledPolyPatch& pp =
 339                 refCast<const coupledPolyPatch>(patches[patchI]);
 340
 341             // Calculated base points on coupled faces are those of
 342             // the owner patch face. They need to be reindexed to for
 343             // the non-owner face, which has the opposite order.
 344
 345             // So, for fPtI_o != 0, fPtI_n = f.size() - fPtI_o
 346
 347             // i.e.:
 348
 349             // owner coupledPolyPatch face
 350             // face    (a b c d e f)
 351             // fPtI     0 1 2 3 4 5
 352             //            +
 353             //              #
 354
 355             // neighbour coupledPolyPatch face
 356             // face    (a f e d c b)
 357             // fPtI     0 1 2 3 4 5
 358             //                    +
 359             //                  #
 360             // +: 6 - 1 = 5
 361             // #: 6 - 2 = 4
 362
 363             if (!pp.owner())
 364             {
 365                 bFTetBasePtI = mesh.faces()[fI].size() - bFTetBasePtI;
 366             }
 367         }
 368     }
 369
 370     return tetBasePtIs;
 371 }
 372
 373
 374 bool Foam::polyMeshTetDecomposition::checkFaceTets
 375 (
 376     const polyMesh& mesh,
 377     scalar tol,
 378     const bool report,
 379     labelHashSet* setPtr
 380 )
 381 {
 382     const labelList& own = mesh.faceOwner();
 383     const labelList& nei = mesh.faceNeighbour();
 384     const polyBoundaryMesh& patches = mesh.boundaryMesh();
 385
 386     const vectorField& cc = mesh.cellCentres();
 387     const vectorField& fc = mesh.faceCentres();
 388
 389     // Calculate coupled cell centre
 390     pointField neiCc(mesh.nFaces() - mesh.nInternalFaces());
 391
 392     for (label faceI = mesh.nInternalFaces(); faceI < mesh.nFaces(); faceI++)
 393     {
 394         neiCc[faceI - mesh.nInternalFaces()] = cc[own[faceI]];
 395     }
 396
 397     syncTools::swapBoundaryFacePositions(mesh, neiCc);
 398
 399     const faceList& fcs = mesh.faces();
 400
 401     const pointField& p = mesh.points();
 402
 403     label nErrorTets = 0;
 404
 405     forAll(fcs, faceI)
 406     {
 407         const face& f = fcs[faceI];
 408
 409         forAll(f, fPtI)
 410         {
 411             scalar tetQual = tetPointRef
 412             (
 413                 p[f[fPtI]],
 414                 p[f.nextLabel(fPtI)],
 415                 fc[faceI],
 416                 cc[own[faceI]]
 417             ).quality();
 418
 419             if (tetQual > -tol)
 420             {
 421                 if (setPtr)
 422                 {
 423                     setPtr->insert(faceI);
 424                 }
 425
 426                 nErrorTets++;
 427                 break;              // no need to check other tets
 428             }
 429         }
 430
 431         if (mesh.isInternalFace(faceI))
 432         {
 433             // Create the neighbour tet - it will have positive volume
 434             const face& f = fcs[faceI];
 435
 436             forAll(f, fPtI)
 437             {
 438                 scalar tetQual = tetPointRef
 439                 (
 440                     p[f[fPtI]],
 441                     p[f.nextLabel(fPtI)],
 442                     fc[faceI],
 443                     cc[nei[faceI]]
 444                 ).quality();
 445
 446                 if (tetQual < tol)
 447                 {
 448                     if (setPtr)
 449                     {
 450                         setPtr->insert(faceI);
 451                     }
 452
 453                     nErrorTets++;
 454                     break;
 455                 }
 456             }
 457
 458             if (findSharedBasePoint(mesh, faceI, tol, report) == -1)
 459             {
 460                 if (setPtr)
 461                 {
 462                     setPtr->insert(faceI);
 463                 }
 464
 465                 nErrorTets++;
 466             }
 467         }
 468         else
 469         {
 470             label patchI = patches.patchID()[faceI - mesh.nInternalFaces()];
 471
 472             if (patches[patchI].coupled())
 473             {
 474                 if
 475                 (
 476                     findSharedBasePoint
 477                     (
 478                         mesh,
 479                         faceI,
 480                         neiCc[faceI - mesh.nInternalFaces()],
 481                         tol,
 482                         report
 483                     ) == -1
 484                 )
 485                 {
 486                     if (setPtr)
 487                     {
 488                         setPtr->insert(faceI);
 489                     }
 490
 491                     nErrorTets++;
 492                 }
 493             }
 494             else
 495             {
 496                 if (findBasePoint(mesh, faceI, tol, report) == -1)
 497                 {
 498                     if (setPtr)
 499                     {
 500                         setPtr->insert(faceI);
 501                     }
 502
 503                     nErrorTets++;
 504                 }
 505             }
 506         }
 507     }
 508
 509     reduce(nErrorTets, sumOp<label>());
 510
 511     if (nErrorTets > 0)
 512     {
 513         if (report)
 514         {
 515             Info<< " ***Error in face tets: "
 516                 << nErrorTets << " faces with low quality or negative volume "
 517                 << "decomposition tets." << endl;
 518         }
 519
 520         return true;
 521     }
 522     else
 523     {
 524         if (report)
 525         {
 526             Info<< "    Face tets OK." << endl;
 527         }
 528
 529         return false;
 530     }
 531 }
 532
 533
 534 Foam::List<Foam::tetIndices> Foam::polyMeshTetDecomposition::faceTetIndices
 535 (
 536     const polyMesh& mesh,
 537     label fI,
 538     label cI
 539 )
 540 {
 541     const faceList& pFaces = mesh.faces();
 542     const labelList& pOwner = mesh.faceOwner();
 543
 544     const face& f = pFaces[fI];
 545
 546     label nTets = f.size() - 2;
 547
 548     List<tetIndices> faceTets(nTets);
 549
 550     bool own = (pOwner[fI] == cI);
 551
 552     label tetBasePtI = mesh.tetBasePtIs()[fI];
 553
 554     if (tetBasePtI == -1)
 555     {
 556         WarningIn
 557         (
 558             "Foam::List<Foam::FixedList<Foam::label, 4> >"
 559             "Foam::Cloud<ParticleType>::"
 560             "faceTetIndices(label fI, label cI) const"
 561         )
 562             << "No base point for face " << fI << ", " << f
 563             << ", produces a valid tet decomposition."
 564             << endl;
 565
 566         tetBasePtI = 0;
 567     }
 568
 569     for (label tetPtI = 1; tetPtI < f.size() - 1; tetPtI++)
 570     {
 571         tetIndices& faceTetIs = faceTets[tetPtI - 1];
 572
 573         label facePtI = (tetPtI + tetBasePtI) % f.size();
 574         label otherFacePtI = f.fcIndex(facePtI);
 575
 576         faceTetIs.cell() = cI;
 577
 578         faceTetIs.face() = fI;
 579
 580         faceTetIs.faceBasePt() = tetBasePtI;
 581
 582         if (own)
 583         {
 584             faceTetIs.facePtA() = facePtI;
 585             faceTetIs.facePtB() = otherFacePtI;
 586         }
 587         else
 588         {
 589             faceTetIs.facePtA() = otherFacePtI;
 590             faceTetIs.facePtB() = facePtI;
 591         }
 592
 593         faceTetIs.tetPt() = tetPtI;
 594     }
 595
 596     return faceTets;
 597 }
 598
 599
 600 Foam::List<Foam::tetIndices> Foam::polyMeshTetDecomposition::cellTetIndices
 601 (
 602     const polyMesh& mesh,
 603     label cI
 604 )
 605 {
 606     const faceList& pFaces = mesh.faces();
 607     const cellList& pCells = mesh.cells();
 608
 609     const cell& thisCell = pCells[cI];
 610
 611     label nTets = 0;
 612
 613     forAll(thisCell, cFI)
 614     {
 615         nTets += pFaces[thisCell[cFI]].size() - 2;
 616     }
 617
 618     DynamicList<tetIndices> cellTets(nTets);
 619
 620     forAll(thisCell, cFI)
 621     {
 622         label fI = thisCell[cFI];
 623
 624         cellTets.append(faceTetIndices(mesh, fI, cI));
 625     }
 626
 627     return cellTets;
 628 }
 629
 630
 631 // ************************************************************************* //