src/triSurface/meshTriangulation/meshTriangulation.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
   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 "meshTriangulation.H"
  27 #include "polyMesh.H"
  28 #include "faceTriangulation.H"
  29
  30 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  31
  32 bool Foam::meshTriangulation::isInternalFace
  33 (
  34     const primitiveMesh& mesh,
  35     const boolList& includedCell,
  36     const label faceI
  37 )
  38 {
  39     if (mesh.isInternalFace(faceI))
  40     {
  41         label own = mesh.faceOwner()[faceI];
  42         label nei = mesh.faceNeighbour()[faceI];
  43
  44         if (includedCell[own] && includedCell[nei])
  45         {
  46             // Neighbouring cell will get included in subset
  47             // as well so face is internal.
  48             return true;
  49         }
  50         else
  51         {
  52             return false;
  53         }
  54     }
  55     else
  56     {
  57         return false;
  58     }
  59 }
  60
  61
  62 void Foam::meshTriangulation::getFaces
  63 (
  64     const primitiveMesh& mesh,
  65     const boolList& includedCell,
  66     boolList& faceIsCut,
  67     label& nFaces,
  68     label& nInternalFaces
  69 )
  70 {
  71     // All faces to be triangulated.
  72     faceIsCut.setSize(mesh.nFaces());
  73     faceIsCut = false;
  74
  75     nFaces = 0;
  76     nInternalFaces = 0;
  77
  78     forAll(includedCell, cellI)
  79     {
  80         // Include faces of cut cells only.
  81         if (includedCell[cellI])
  82         {
  83             const labelList& cFaces = mesh.cells()[cellI];
  84
  85             forAll(cFaces, i)
  86             {
  87                 label faceI = cFaces[i];
  88
  89                 if (!faceIsCut[faceI])
  90                 {
  91                     // First visit of face.
  92                     nFaces++;
  93                     faceIsCut[faceI] = true;
  94
  95                     // See if would become internal or external face
  96                     if (isInternalFace(mesh, includedCell, faceI))
  97                     {
  98                         nInternalFaces++;
  99                     }
 100                 }
 101             }
 102         }
 103     }
 104
 105     Pout<< "Subset consists of " << nFaces << " faces out of " << mesh.nFaces()
 106         << " of which " << nInternalFaces << " are internal" << endl;
 107 }
 108
 109
 110 void Foam::meshTriangulation::insertTriangles
 111 (
 112     const triFaceList& faceTris,
 113     const label faceI,
 114     const label regionI,
 115     const bool reverse,
 116
 117     List<labelledTri>& triangles,
 118     label& triI
 119 )
 120 {
 121     // Copy triangles. Optionally reverse them
 122     forAll(faceTris, i)
 123     {
 124         const triFace& f = faceTris[i];
 125
 126         labelledTri& tri = triangles[triI];
 127
 128         if (reverse)
 129         {
 130             tri[0] = f[0];
 131             tri[2] = f[1];
 132             tri[1] = f[2];
 133         }
 134         else
 135         {
 136             tri[0] = f[0];
 137             tri[1] = f[1];
 138             tri[2] = f[2];
 139         }
 140
 141         tri.region() = regionI;
 142
 143         faceMap_[triI] = faceI;
 144
 145         triI++;
 146     }
 147 }
 148
 149
 150 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 151
 152 // Null constructor
 153 Foam::meshTriangulation::meshTriangulation()
 154 :
 155     triSurface(),
 156     nInternalFaces_(0),
 157     faceMap_()
 158 {}
 159
 160
 161 // Construct from faces of cells
 162 Foam::meshTriangulation::meshTriangulation
 163 (
 164     const polyMesh& mesh,
 165     const label internalFacesPatch,
 166     const boolList& includedCell,
 167     const bool faceCentreDecomposition
 168 )
 169 :
 170     triSurface(),
 171     nInternalFaces_(0),
 172     faceMap_()
 173 {
 174     const faceList& faces = mesh.faces();
 175     const pointField& points = mesh.points();
 176     const polyBoundaryMesh& patches = mesh.boundaryMesh();
 177
 178     // All faces to be triangulated.
 179     boolList faceIsCut;
 180     label nFaces, nInternalFaces;
 181
 182     getFaces
 183     (
 184         mesh,
 185         includedCell,
 186         faceIsCut,
 187         nFaces,
 188         nInternalFaces
 189     );
 190
 191
 192     // Find upper limit for number of triangles
 193     // (can be less if triangulation fails)
 194     label nTotTri = 0;
 195
 196     if (faceCentreDecomposition)
 197     {
 198         forAll(faceIsCut, faceI)
 199         {
 200             if (faceIsCut[faceI])
 201             {
 202                 nTotTri += faces[faceI].size();
 203             }
 204         }
 205     }
 206     else
 207     {
 208         forAll(faceIsCut, faceI)
 209         {
 210             if (faceIsCut[faceI])
 211             {
 212                 nTotTri += faces[faceI].nTriangles(points);
 213             }
 214         }
 215     }
 216     Pout<< "nTotTri : " << nTotTri << endl;
 217
 218
 219     // Storage for new and old points (only for faceCentre decomposition;
 220     // for triangulation uses only existing points)
 221     pointField newPoints;
 222
 223     if (faceCentreDecomposition)
 224     {
 225         newPoints.setSize(mesh.nPoints() + faces.size());
 226         forAll(mesh.points(), pointI)
 227         {
 228             newPoints[pointI] = mesh.points()[pointI];
 229         }
 230         // Face centres
 231         forAll(faces, faceI)
 232         {
 233             newPoints[mesh.nPoints() + faceI] = mesh.faceCentres()[faceI];
 234         }
 235     }
 236
 237     // Storage for all triangles
 238     List<labelledTri> triangles(nTotTri);
 239     faceMap_.setSize(nTotTri);
 240     label triI = 0;
 241
 242
 243     if (faceCentreDecomposition)
 244     {
 245         // Decomposition around face centre
 246         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 247
 248         // Triangulate internal faces
 249         forAll(faceIsCut, faceI)
 250         {
 251             if (faceIsCut[faceI] && isInternalFace(mesh, includedCell, faceI))
 252             {
 253                 // Face was internal to the mesh and will be 'internal' to
 254                 // the surface.
 255
 256                 // Triangulate face
 257                 const face& f = faces[faceI];
 258
 259                 forAll(f, fp)
 260                 {
 261                     faceMap_[triI] = faceI;
 262
 263                     triangles[triI++] =
 264                         labelledTri
 265                         (
 266                             f[fp],
 267                             f.nextLabel(fp),
 268                             mesh.nPoints() + faceI,     // face centre
 269                             internalFacesPatch
 270                         );
 271                 }
 272             }
 273         }
 274         nInternalFaces_ = triI;
 275
 276
 277         // Triangulate external faces
 278         forAll(faceIsCut, faceI)
 279         {
 280             if (faceIsCut[faceI] && !isInternalFace(mesh, includedCell, faceI))
 281             {
 282                 // Face will become outside of the surface.
 283
 284                 label patchI = -1;
 285                 bool reverse = false;
 286
 287                 if (mesh.isInternalFace(faceI))
 288                 {
 289                     patchI = internalFacesPatch;
 290
 291                     // Check orientation. Check which side of the face gets
 292                     // included (note: only one side is).
 293                     if (includedCell[mesh.faceOwner()[faceI]])
 294                     {
 295                         reverse = false;
 296                     }
 297                     else
 298                     {
 299                         reverse = true;
 300                     }
 301                 }
 302                 else
 303                 {
 304                     // Face was already outside so orientation ok.
 305
 306                     patchI = patches.whichPatch(faceI);
 307
 308                     reverse = false;
 309                 }
 310
 311
 312                 // Triangulate face
 313                 const face& f = faces[faceI];
 314
 315                 if (reverse)
 316                 {
 317                     forAll(f, fp)
 318                     {
 319                         faceMap_[triI] = faceI;
 320
 321                         triangles[triI++] =
 322                             labelledTri
 323                             (
 324                                 f.nextLabel(fp),
 325                                 f[fp],
 326                                 mesh.nPoints() + faceI,     // face centre
 327                                 patchI
 328                             );
 329                     }
 330                 }
 331                 else
 332                 {
 333                     forAll(f, fp)
 334                     {
 335                         faceMap_[triI] = faceI;
 336
 337                         triangles[triI++] =
 338                             labelledTri
 339                             (
 340                                 f[fp],
 341                                 f.nextLabel(fp),
 342                                 mesh.nPoints() + faceI,     // face centre
 343                                 patchI
 344                             );
 345                     }
 346                 }
 347             }
 348         }
 349     }
 350     else
 351     {
 352         // Triangulation using existing vertices
 353         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 354
 355         // Triangulate internal faces
 356         forAll(faceIsCut, faceI)
 357         {
 358             if (faceIsCut[faceI] && isInternalFace(mesh, includedCell, faceI))
 359             {
 360                 // Face was internal to the mesh and will be 'internal' to
 361                 // the surface.
 362
 363                 // Triangulate face. Fall back to naive triangulation if failed.
 364                 faceTriangulation faceTris(points, faces[faceI], true);
 365
 366                 if (faceTris.empty())
 367                 {
 368                     WarningIn("meshTriangulation::meshTriangulation")
 369                         << "Could not find triangulation for face " << faceI
 370                         << " vertices " << faces[faceI] << " coords "
 371                         << IndirectList<point>(points, faces[faceI])() << endl;
 372                 }
 373                 else
 374                 {
 375                     // Copy triangles. Make them internalFacesPatch
 376                     insertTriangles
 377                     (
 378                         faceTris,
 379                         faceI,
 380                         internalFacesPatch,
 381                         false,                  // no reverse
 382
 383                         triangles,
 384                         triI
 385                     );
 386                 }
 387             }
 388         }
 389         nInternalFaces_ = triI;
 390
 391
 392         // Triangulate external faces
 393         forAll(faceIsCut, faceI)
 394         {
 395             if (faceIsCut[faceI] && !isInternalFace(mesh, includedCell, faceI))
 396             {
 397                 // Face will become outside of the surface.
 398
 399                 label patchI = -1;
 400                 bool reverse = false;
 401
 402                 if (mesh.isInternalFace(faceI))
 403                 {
 404                     patchI = internalFacesPatch;
 405
 406                     // Check orientation. Check which side of the face gets
 407                     // included (note: only one side is).
 408                     if (includedCell[mesh.faceOwner()[faceI]])
 409                     {
 410                         reverse = false;
 411                     }
 412                     else
 413                     {
 414                         reverse = true;
 415                     }
 416                 }
 417                 else
 418                 {
 419                     // Face was already outside so orientation ok.
 420
 421                     patchI = patches.whichPatch(faceI);
 422
 423                     reverse = false;
 424                 }
 425
 426                 // Triangulate face
 427                 faceTriangulation faceTris(points, faces[faceI], true);
 428
 429                 if (faceTris.empty())
 430                 {
 431                     WarningIn("meshTriangulation::meshTriangulation")
 432                         << "Could not find triangulation for face " << faceI
 433                         << " vertices " << faces[faceI] << " coords "
 434                         << IndirectList<point>(points, faces[faceI])() << endl;
 435                 }
 436                 else
 437                 {
 438                     // Copy triangles. Optionally reverse them
 439                     insertTriangles
 440                     (
 441                         faceTris,
 442                         faceI,
 443                         patchI,
 444                         reverse,    // whether to reverse
 445
 446                         triangles,
 447                         triI
 448                     );
 449                 }
 450             }
 451         }
 452     }
 453
 454     // Shrink if nessecary (because of invalid triangulations)
 455     triangles.setSize(triI);
 456     faceMap_.setSize(triI);
 457
 458     Pout<< "nInternalFaces_:" << nInternalFaces_ << endl;
 459     Pout<< "triangles:" << triangles.size() << endl;
 460
 461
 462     geometricSurfacePatchList surfPatches(patches.size());
 463
 464     forAll(patches, patchI)
 465     {
 466         surfPatches[patchI] =
 467             geometricSurfacePatch
 468             (
 469                 patches[patchI].physicalType(),
 470                 patches[patchI].name(),
 471                 patchI
 472             );
 473     }
 474
 475     // Create globally numbered tri surface
 476     if (faceCentreDecomposition)
 477     {
 478         // Use newPoints (mesh points + face centres)
 479         triSurface globalSurf(triangles, surfPatches, newPoints);
 480
 481         // Create locally numbered tri surface
 482         triSurface::operator=
 483         (
 484             triSurface
 485             (
 486                 globalSurf.localFaces(),
 487                 surfPatches,
 488                 globalSurf.localPoints()
 489             )
 490         );
 491     }
 492     else
 493     {
 494         // Use mesh points
 495         triSurface globalSurf(triangles, surfPatches, mesh.points());
 496
 497         // Create locally numbered tri surface
 498         triSurface::operator=
 499         (
 500             triSurface
 501             (
 502                 globalSurf.localFaces(),
 503                 surfPatches,
 504                 globalSurf.localPoints()
 505             )
 506         );
 507     }
 508 }
 509
 510
 511 // ************************************************************************* //