src/triSurface/meshTriangulation/meshTriangulation.C

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