src/tetFiniteElement/tetPolyMesh/calcTetPolyMeshGeometry.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend 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 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "tetPolyMesh.H"
  27 #include "tetPointRef.H"
  28 #include "tetPointRef.H"
  29 #include "SquareMatrix.H"
  30
  31 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  32
  33 namespace Foam
  34 {
  35
  36 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  37
  38 label tetPolyMesh::nPoints() const
  39 {
  40     if (nPoints_ < 0)
  41     {
  42         nPoints_ =
  43             mesh_.nPoints()
  44           + mesh_.nFaces()
  45           + mesh_.nCells();
  46     }
  47
  48     return nPoints_;
  49 }
  50
  51
  52 label tetPolyMesh::nEdges() const
  53 {
  54     if (nEdges_ < 0)
  55     {
  56         nEdges_ = mesh_.nEdges();
  57
  58         const labelListList& pf = mesh_.pointFaces();
  59
  60         const labelListList& pc = mesh_.pointCells();
  61
  62         forAll (pf, pointI)
  63         {
  64             nEdges_ += pf[pointI].size();
  65
  66             nEdges_ += pc[pointI].size();
  67         }
  68
  69         const cellList& cellFaces = mesh_.cells();
  70
  71         forAll (cellFaces, cellI)
  72         {
  73             nEdges_ += cellFaces[cellI].size();
  74         }
  75     }
  76
  77     return nEdges_;
  78 }
  79
  80
  81 label tetPolyMesh::nTets() const
  82 {
  83     if (nTets_ < 0)
  84     {
  85         // count the cells
  86         nTets_ = 0;
  87
  88         const cellList& polyCells = mesh_.cells();
  89
  90         forAll (polyCells, cellI)
  91         {
  92             nTets_ += nTetsForCell(cellI);
  93         }
  94     }
  95
  96     return nTets_;
  97 }
  98
  99
 100 // Return number of tetrahedra in decomposition for cell
 101 label tetPolyMesh::nTetsForCell(const label cellID) const
 102 {
 103     const unallocFaceList& f = mesh_.faces();
 104
 105     label nTetrasForCell = 0;
 106
 107     const labelList& cellFaces = mesh_.cells()[cellID];
 108
 109     forAll (cellFaces, faceI)
 110     {
 111         nTetrasForCell += f[cellFaces[faceI]].nEdges();
 112     }
 113
 114     return nTetrasForCell;
 115 }
 116
 117
 118 tmp<pointField> tetPolyMesh::points() const
 119 {
 120     tmp<pointField> ttetPoints(new pointField(nPoints()));
 121     pointField& tetPoints = ttetPoints();
 122
 123     const pointField& points = mesh_.points();
 124
 125     const pointField& faceCentres = mesh_.faceCentres();
 126
 127     const pointField& cellCentres = mesh_.cellCentres();
 128
 129     label tetPointI = 0;
 130
 131     forAll (points, pointI)
 132     {
 133         tetPoints[tetPointI] = points[pointI];
 134         tetPointI++;
 135     }
 136
 137     forAll (faceCentres, faceI)
 138     {
 139         tetPoints[tetPointI] = faceCentres[faceI];
 140         tetPointI++;
 141     }
 142
 143     forAll (cellCentres, cellI)
 144     {
 145         tetPoints[tetPointI] = cellCentres[cellI];
 146         tetPointI++;
 147     }
 148
 149     return ttetPoints;
 150 }
 151
 152
 153 cellShapeList tetPolyMesh::tetCells() const
 154 {
 155     cellShapeList t(nTets());
 156
 157     const cellList& polyCells = mesh_.cells();
 158
 159     label nTetCells = 0;
 160
 161     forAll (polyCells, cellI)
 162     {
 163         tetCellList cellTets = tets(cellI);
 164
 165         forAll (cellTets, tetI)
 166         {
 167             t[nTetCells] = cellTets[tetI].tetCellShape();
 168             nTetCells++;
 169         }
 170     }
 171
 172     return t;
 173 }
 174
 175
 176 // Return tetrahedral decomposition for cell
 177 tetCellList tetPolyMesh::tets(const label cellID) const
 178 {
 179     const unallocFaceList& f = mesh_.faces();
 180
 181     const unallocLabelList& owner = mesh_.faceOwner();
 182
 183     // Initialise the size of the return
 184     tetCellList t(nTetsForCell(cellID));
 185
 186     label nTetras = 0;
 187
 188     const labelList& cellFaces = mesh_.cells()[cellID];
 189
 190     // A tet is created from the face edge, face centre and cell centre
 191     forAll (cellFaces, faceI)
 192     {
 193         const label curFace = cellFaces[faceI];
 194
 195         edgeList faceEdges;
 196
 197         // Take care of owner/neighbour face reversal
 198         if (cellID == owner[curFace])
 199         {
 200             faceEdges = f[curFace].reverseFace().edges();
 201         }
 202         else
 203         {
 204             faceEdges = f[curFace].edges();
 205         }
 206
 207         forAll (faceEdges, edgeI)
 208         {
 209             // Face is assumed to be inward-pointing
 210             t[nTetras] =
 211                 tetCell
 212                 (
 213                     faceEdges[edgeI].start(),
 214                     faceEdges[edgeI].end(),
 215                     curFace + faceOffset_,
 216                     cellID + cellOffset_
 217                 );
 218
 219             nTetras++;
 220         }
 221     }
 222
 223     return t;
 224 }
 225
 226
 227 void tetPolyMesh::gradNiDotGradNj
 228 (
 229     const label cellID,
 230     SquareMatrix<scalar>& denseMatrix,
 231     const labelList& globalToLocalBuffer
 232 ) const
 233 {
 234     const unallocFaceList& meshFaces = mesh_.faces();
 235     const unallocLabelList& owner = mesh_.faceOwner();
 236     const labelList& cellFaces = mesh_.cells()[cellID];
 237
 238     const pointField& points = mesh_.points();
 239     const pointField& faceCentres = mesh_.faceCentres();
 240     const pointField& cellCentres = mesh_.cellCentres();
 241
 242     scalarField tetBuffer(tetPointRef::nEdges);
 243
 244     // For optimisation, point allocation is taken out of the loop
 245     label vertexI, edgeI, localI, localJ;
 246     edge curEdge;
 247
 248     forAll (cellFaces, faceI)
 249     {
 250         const label curFaceID = cellFaces[faceI];
 251
 252         edgeList faceEdges;
 253
 254         if (cellID == owner[curFaceID])
 255         {
 256             faceEdges = meshFaces[curFaceID].reverseFace().edges();
 257         }
 258         else
 259         {
 260             faceEdges = meshFaces[curFaceID].edges();
 261         }
 262
 263         forAll (faceEdges, i)
 264         {
 265             tetCell curTet
 266             (
 267                 faceEdges[i].start(),
 268                 faceEdges[i].end(),
 269                 faceOffset() + curFaceID,
 270                 cellOffset() + cellID
 271             );
 272
 273             tetPointRef tpr
 274             (
 275                 points[faceEdges[i].start()],
 276                 points[faceEdges[i].end()],
 277                 faceCentres[curFaceID],
 278                 cellCentres[cellID]
 279             );
 280
 281             // Calculate the diagonal and insert into local matrix
 282             tpr.gradNiSquared(tetBuffer);
 283
 284             for
 285             (
 286                 vertexI = 0;
 287                 vertexI < tetPointRef::nVertices;
 288                 vertexI++
 289             )
 290             {
 291                 localI = globalToLocalBuffer[curTet[vertexI]];
 292
 293                 denseMatrix[localI][localI] += tetBuffer[vertexI];
 294             }
 295
 296             tpr.gradNiDotGradNj(tetBuffer);
 297
 298             for (edgeI = 0; edgeI < tetPointRef::nEdges; edgeI++)
 299             {
 300                 curEdge = curTet.tetEdge(edgeI);
 301
 302                 localI = globalToLocalBuffer[curEdge.start()];
 303                 localJ = globalToLocalBuffer[curEdge.end()];
 304
 305                 // Insert the coefficient into upper and lower triangle
 306                 // HJ, 4/Dec/2006
 307                 denseMatrix[localI][localJ] += tetBuffer[edgeI];
 308                 denseMatrix[localJ][localI] += tetBuffer[edgeI];
 309             }
 310         }
 311     }
 312 }
 313
 314
 315 void tetPolyMesh::gradNiGradNj
 316 (
 317     const label cellID,
 318     SquareMatrix<tensor>& denseMatrix,
 319     const labelList& globalToLocalBuffer
 320 ) const
 321 {
 322     const unallocFaceList& meshFaces = mesh_.faces();
 323     const unallocLabelList& owner = mesh_.faceOwner();
 324     const labelList& cellFaces = mesh_.cells()[cellID];
 325
 326     const pointField& points = mesh_.points();
 327     const pointField& faceCentres = mesh_.faceCentres();
 328     const pointField& cellCentres = mesh_.cellCentres();
 329
 330     tensorField tetBuffer(tetPointRef::nEdges);
 331
 332     // For optimisation, point allocation is taken out of the loop
 333     label vertexI, edgeI, localI, localJ;
 334     edge curEdge;
 335
 336     forAll (cellFaces, faceI)
 337     {
 338         const label curFaceID = cellFaces[faceI];
 339
 340         edgeList faceEdges;
 341
 342         if (cellID == owner[curFaceID])
 343         {
 344             faceEdges = meshFaces[curFaceID].reverseFace().edges();
 345         }
 346         else
 347         {
 348             faceEdges = meshFaces[curFaceID].edges();
 349         }
 350
 351         forAll (faceEdges, i)
 352         {
 353             tetCell curTet
 354             (
 355                 faceEdges[i].start(),
 356                 faceEdges[i].end(),
 357                 faceOffset() + curFaceID,
 358                 cellOffset() + cellID
 359             );
 360
 361             tetPointRef tpr
 362             (
 363                 points[faceEdges[i].start()],
 364                 points[faceEdges[i].end()],
 365                 faceCentres[curFaceID],
 366                 cellCentres[cellID]
 367             );
 368
 369             // Calculate the diagonal and insert into local matrix
 370             tpr.gradNiGradNi(tetBuffer);
 371
 372             for
 373             (
 374                 vertexI = 0;
 375                 vertexI < tetPointRef::nVertices;
 376                 vertexI++
 377             )
 378             {
 379                 localI = globalToLocalBuffer[curTet[vertexI]];
 380
 381                 denseMatrix[localI][localI] += tetBuffer[vertexI];
 382             }
 383
 384             tpr.gradNiGradNj(tetBuffer);
 385
 386             for (edgeI = 0; edgeI < tetPointRef::nEdges; edgeI++)
 387             {
 388                 curEdge = curTet.tetEdge(edgeI);
 389
 390                 localI = globalToLocalBuffer[curEdge.start()];
 391                 localJ = globalToLocalBuffer[curEdge.end()];
 392
 393                 // Add coefficient and its transpose to upper and lower
 394                 // triangle.  HJ, 4/Dec/2006
 395                 denseMatrix[localI][localJ] += tetBuffer[edgeI];
 396                 denseMatrix[localJ][localI] += tetBuffer[edgeI].T();
 397             }
 398         }
 399     }
 400 }
 401
 402
 403 // Fill buffer with the volume integral distributed into vertices
 404 void tetPolyMesh::volIntegral
 405 (
 406     const label cellID,
 407     scalarField& buffer,
 408     const labelList& globalToLocalBuffer
 409 ) const
 410 {
 411     // The addressing and volume distribution has been done tet-by-tet
 412     // and the combination is done later for the whole cell.
 413
 414     const unallocFaceList& meshFaces = mesh_.faces();
 415     const unallocLabelList& owner = mesh_.faceOwner();
 416     const labelList& cellFaces = mesh_.cells()[cellID];
 417
 418     const pointField& points = mesh_.points();
 419     const pointField& faceCentres = mesh_.faceCentres();
 420     const pointField& cellCentres = mesh_.cellCentres();
 421
 422     label vertexI, localI;
 423
 424     forAll (cellFaces, faceI)
 425     {
 426         const label curFaceID = cellFaces[faceI];
 427
 428         edgeList faceEdges;
 429
 430         // Take care of owner/neighbour face reversal
 431         if (cellID == owner[curFaceID])
 432         {
 433             faceEdges = meshFaces[curFaceID].reverseFace().edges();
 434         }
 435         else
 436         {
 437             faceEdges = meshFaces[curFaceID].edges();
 438         }
 439
 440         forAll (faceEdges, edgeI)
 441         {
 442             // Face is assumed to be inward-pointing
 443             tetCell curTet
 444             (
 445                 faceEdges[edgeI].start(),
 446                 faceEdges[edgeI].end(),
 447                 faceOffset() + curFaceID,
 448                 cellOffset() + cellID
 449             );
 450
 451             scalar quarterVolume =
 452                 0.25*
 453                 tetPointRef
 454                 (
 455                     points[faceEdges[edgeI].start()],
 456                     points[faceEdges[edgeI].end()],
 457                     faceCentres[curFaceID],
 458                     cellCentres[cellID]
 459                 ).mag();
 460
 461             for
 462             (
 463                 vertexI = 0;
 464                 vertexI < tetPointRef::nVertices;
 465                 vertexI++
 466             )
 467             {
 468                 localI = globalToLocalBuffer[curTet[vertexI]];
 469
 470                 // Accumulate the diagonal
 471                 buffer[localI] += quarterVolume;
 472             }
 473         }
 474     }
 475 }
 476
 477
 478 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 479
 480 } // End namespace Foam
 481
 482 // ************************************************************************* //