src/tetDecompositionFiniteElement/tetPolyMeshCellDecomp/calcTetPolyMeshCellDecompGeometry.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 "tetPolyMeshCellDecomp.H"
  28 #include "tetPointRef.H"
  29 #include "SquareMatrix.H"
  30
  31 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  32
  33 namespace Foam
  34 {
  35
  36 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  37
  38 label tetPolyMeshCellDecomp::nPoints() const
  39 {
  40     if (nPoints_ < 0)
  41     {
  42         nPoints_ = mesh_.nPoints() + mesh_.nCells();
  43     }
  44
  45     return nPoints_;
  46 }
  47
  48
  49 label tetPolyMeshCellDecomp::nEdges() const
  50 {
  51     if (nEdges_ < 0)
  52     {
  53         nEdges_ = mesh_.nEdges();
  54
  55         // add edges going across the face
  56         const faceList& f  = mesh_.faces();
  57
  58         forAll (f, fI)
  59         {
  60             nEdges_ += f[fI].size() - 3;
  61         }
  62
  63         // add point-to-cell edges
  64         const labelListList& pc = mesh_.pointCells();
  65
  66         forAll (pc, pointI)
  67         {
  68             nEdges_ += pc[pointI].size();
  69         }
  70     }
  71
  72     return nEdges_;
  73 }
  74
  75
  76 label tetPolyMeshCellDecomp::nTets() const
  77 {
  78     if (nTets_ < 0)
  79     {
  80         // Count the cells
  81         nTets_ = 0;
  82
  83         const cellList& polyCells = mesh_.cells();
  84
  85         forAll (polyCells, cellI)
  86         {
  87             nTets_ += nTetsForCell(cellI);
  88         }
  89     }
  90
  91     return nTets_;
  92 }
  93
  94
  95 tmp<pointField> tetPolyMeshCellDecomp::points() const
  96 {
  97     tmp<pointField> ttetPoints(new pointField(nPoints()));
  98     pointField& tetPoints = ttetPoints();
  99
 100     const pointField& points = mesh_.points();
 101
 102     const pointField& cellCentres = mesh_.cellCentres();
 103
 104     label tetPointI = 0;
 105
 106     forAll (points, pointI)
 107     {
 108         tetPoints[tetPointI] = points[pointI];
 109         tetPointI++;
 110     }
 111
 112     forAll (cellCentres, cellI)
 113     {
 114         tetPoints[tetPointI] = cellCentres[cellI];
 115         tetPointI++;
 116     }
 117
 118     return ttetPoints;
 119 }
 120
 121
 122 cellShapeList tetPolyMeshCellDecomp::tetCells() const
 123 {
 124     cellShapeList t(nTets());
 125
 126     const cellList& polyCells = mesh_.cells();
 127
 128     label nTetCells = 0;
 129
 130     forAll (polyCells, cellI)
 131     {
 132         tetCellList cellTets = tets(cellI);
 133
 134         forAll (cellTets, tetI)
 135         {
 136             t[nTetCells] = cellTets[tetI].tetCellShape();
 137             nTetCells++;
 138         }
 139     }
 140
 141     return t;
 142 }
 143
 144
 145 // Return tetrahedral decomposition for cell
 146 tetCellList tetPolyMeshCellDecomp::tets(const label cellID) const
 147 {
 148     const unallocFaceList& meshFaces = mesh_.faces();
 149
 150     const unallocLabelList& owner = mesh_.faceOwner();
 151
 152     // Initialise the size of the return
 153     tetCellList t(nTetsForCell(cellID));
 154
 155     label nTetras = 0;
 156
 157     const labelList& cellFaces = mesh_.cells()[cellID];
 158
 159     forAll (cellFaces, faceI)
 160     {
 161         const label curFaceID = cellFaces[faceI];
 162
 163         face f = meshFaces[curFaceID];
 164
 165         // Take care of owner/neighbour face reversal
 166         if (cellID == owner[curFaceID])
 167         {
 168             f = f.reverseFace();
 169         }
 170
 171         // Grab the point zero of the face and create triangles by grabbing the
 172         // labels around the perimeter
 173         label pointZeroID = f[0];
 174
 175         for (label i = 1; i < f.size() - 1; i++)
 176         {
 177             // Face is assumed to be inward-pointing
 178             t[nTetras] =
 179                 tetCell
 180                 (
 181                     pointZeroID,
 182                     f[i],
 183                     f[i + 1],
 184                     cellID + cellOffset_
 185                 );
 186
 187             nTetras++;
 188         }
 189     }
 190
 191     return t;
 192 }
 193
 194
 195 void tetPolyMeshCellDecomp::gradNiDotGradNj
 196 (
 197     const label cellID,
 198     SquareMatrix<scalar>& denseMatrix,
 199     const labelList& globalToLocalBuffer
 200 ) const
 201 {
 202     const unallocFaceList& meshFaces = mesh_.faces();
 203     const unallocLabelList& owner = mesh_.faceOwner();
 204     const labelList& cellFaces = mesh_.cells()[cellID];
 205
 206     const pointField& points = mesh_.points();
 207     const pointField& cellCentres = mesh_.cellCentres();
 208
 209     scalarField tetBuffer(tetPointRef::nEdges);
 210
 211     // For optimisation, point allocation is taken out of the loop
 212     label vertexI, edgeI, localI, localJ;
 213     edge curEdge;
 214
 215     forAll (cellFaces, faceI)
 216     {
 217         const label curFaceID = cellFaces[faceI];
 218
 219         face f = meshFaces[curFaceID];
 220
 221         if (cellID == owner[curFaceID])
 222         {
 223             f = f.reverseFace();
 224         }
 225
 226         label pointZeroID = f[0];
 227
 228         for (label i = 1; i < f.size() - 1; i++)
 229         {
 230             tetCell curTet
 231             (
 232                 pointZeroID,
 233                 f[i],
 234                 f[i + 1],
 235                 cellOffset() + cellID
 236             );
 237
 238             tetPointRef tpr
 239             (
 240                 points[pointZeroID],
 241                 points[f[i]],
 242                 points[f[i + 1]],
 243                 cellCentres[cellID]
 244             );
 245
 246             // Calculate the diagonal and insert into local matrix
 247             tpr.gradNiSquared(tetBuffer);
 248
 249             for
 250             (
 251                 vertexI = 0;
 252                 vertexI < tetPointRef::nVertices;
 253                 vertexI++
 254             )
 255             {
 256                 localI = globalToLocalBuffer[curTet[vertexI]];
 257
 258                 denseMatrix[localI][localI] += tetBuffer[vertexI];
 259             }
 260
 261             // Note: the coefficient is symmetric, i.e. equal in the upper and
 262             // lower triangle
 263             tpr.gradNiDotGradNj(tetBuffer);
 264
 265             for (edgeI = 0; edgeI < tetPointRef::nEdges; edgeI++)
 266             {
 267                 edge curEdge = curTet.tetEdge(edgeI);
 268
 269                 localI = globalToLocalBuffer[curEdge.start()];
 270                 localJ = globalToLocalBuffer[curEdge.end()];
 271
 272                 // Insert the coefficient into upper and lower triangle
 273                 // HJ, 4/Dec/2006
 274                 denseMatrix[localI][localJ] += tetBuffer[edgeI];
 275                 denseMatrix[localJ][localI] += tetBuffer[edgeI];
 276             }
 277         }
 278     }
 279 }
 280
 281
 282 void tetPolyMeshCellDecomp::gradNiGradNj
 283 (
 284     const label cellID,
 285     SquareMatrix<tensor>& denseMatrix,
 286     const labelList& globalToLocalBuffer
 287 ) const
 288 {
 289     const unallocFaceList& meshFaces = mesh_.faces();
 290     const unallocLabelList& owner = mesh_.faceOwner();
 291     const labelList& cellFaces = mesh_.cells()[cellID];
 292
 293     const pointField& points = mesh_.points();
 294     const pointField& cellCentres = mesh_.cellCentres();
 295
 296     tensorField tetBuffer(tetPointRef::nEdges);
 297
 298     // For optimisation, point allocation is taken out of the loop
 299     label vertexI, edgeI, localI, localJ, minIJ, maxIJ;
 300     edge curEdge;
 301
 302     forAll (cellFaces, faceI)
 303     {
 304         const label curFaceID = cellFaces[faceI];
 305
 306         face f = meshFaces[curFaceID];
 307
 308         if (cellID == owner[curFaceID])
 309         {
 310             f = f.reverseFace();
 311         }
 312
 313         label pointZeroID = f[0];
 314
 315         for (label i = 1; i < f.size() - 1; i++)
 316         {
 317             tetCell curTet
 318             (
 319                 pointZeroID,
 320                 f[i],
 321                 f[i + 1],
 322                 cellOffset() + cellID
 323             );
 324
 325             tetPointRef tpr
 326             (
 327                 points[pointZeroID],
 328                 points[f[i]],
 329                 points[f[i + 1]],
 330                 cellCentres[cellID]
 331             );
 332
 333             // Calculate the diagonal and insert into local matrix
 334             tpr.gradNiGradNi(tetBuffer);
 335
 336             for
 337             (
 338                 vertexI = 0;
 339                 vertexI < tetPointRef::nVertices;
 340                 vertexI++
 341             )
 342             {
 343                 localI = globalToLocalBuffer[curTet[vertexI]];
 344
 345                 denseMatrix[localI][localI] += tetBuffer[vertexI];
 346             }
 347
 348             tpr.gradNiGradNj(tetBuffer);
 349
 350             for (edgeI = 0; edgeI < tetPointRef::nEdges; edgeI++)
 351             {
 352                 curEdge = curTet.tetEdge(edgeI);
 353
 354                 localI = globalToLocalBuffer[curEdge.start()];
 355                 localJ = globalToLocalBuffer[curEdge.end()];
 356
 357                 // Accumulate the off-diagonal; folding the matrix.
 358                 // It is unknown if the localI or localJ is lower:
 359                 // that depends on how the local points have been
 360                 // selected.  Therefore, the matrix (which is
 361                 // symmetric!) needs to be "folded"
 362
 363                 minIJ = min(localI, localJ);
 364                 maxIJ = max(localI, localJ);
 365
 366                 // Take care of the directionality of the edge
 367                 if (curEdge.start() < curEdge.end())
 368                 {
 369                     denseMatrix[minIJ][maxIJ] += tetBuffer[edgeI];
 370                 }
 371                 else
 372                 {
 373                     denseMatrix[minIJ][maxIJ] += tetBuffer[edgeI].T();
 374                 }
 375             }
 376         }
 377     }
 378 }
 379
 380
 381 // Fill buffer with the volume integral distributed into vertices
 382 void tetPolyMeshCellDecomp::volIntegral
 383 (
 384     const label cellID,
 385     scalarField& buffer,
 386     const labelList& globalToLocalBuffer
 387 ) const
 388 {
 389     // The addressing and volume distribution has been done tet-by-tet
 390     // and the combination is done later for the whole cell.
 391
 392     const unallocFaceList& meshFaces = mesh_.faces();
 393     const unallocLabelList& owner = mesh_.faceOwner();
 394     const labelList& cellFaces = mesh_.cells()[cellID];
 395
 396     const pointField& points = mesh_.points();
 397     const pointField& cellCentres = mesh_.cellCentres();
 398
 399     label vertexI, localI;
 400
 401     forAll (cellFaces, faceI)
 402     {
 403         const label curFaceID = cellFaces[faceI];
 404
 405         face f = meshFaces[curFaceID];
 406
 407         // Take care of owner/neighbour face reversal
 408         if (cellID == owner[curFaceID])
 409         {
 410             f = f.reverseFace();
 411         }
 412
 413         // Grab the point zero of the face and create triangles by grabbing the
 414         // labels around the perimeter
 415         label pointZeroID = f[0];
 416
 417         for (label i = 1; i < f.size() - 1; i++)
 418         {
 419             // Face is assumed to be inward-pointing
 420             tetCell curTet
 421             (
 422                 pointZeroID,
 423                 f[i],
 424                 f[i + 1],
 425                 cellOffset() + cellID
 426             );
 427
 428             scalar quarterVolume =
 429                 0.25*
 430                 tetPointRef
 431                 (
 432                     points[pointZeroID],
 433                     points[f[i]],
 434                     points[f[i + 1]],
 435                     cellCentres[cellID]
 436                 ).mag();
 437
 438             for
 439             (
 440                 vertexI = 0;
 441                 vertexI < tetPointRef::nVertices;
 442                 vertexI++
 443             )
 444             {
 445                 localI = globalToLocalBuffer[curTet[vertexI]];
 446
 447                 // Accumulate the diagonal
 448                 buffer[localI] += quarterVolume;
 449             }
 450         }
 451     }
 452 }
 453
 454
 455 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 456
 457 } // End namespace Foam
 458
 459 // ************************************************************************* //