src/tetDecompositionFiniteElement/tetFemMethod/tetFem.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 Description
  26     Class of static functions to calculate implicit finite element derivatives
  27     returning a matrix.
  28
  29 \*---------------------------------------------------------------------------*/
  30
  31 #include "tetCell.H"
  32 #include "tetCellList.H"
  33 #include "tetPointRef.H"
  34 #include "SquareMatrix.H"
  35
  36 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  37
  38 namespace Foam
  39 {
  40
  41 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  42
  43 template<class Type>
  44 tmp<tetFemMatrix<Type> > tetFem::laplacian
  45 (
  46     GeometricField<Type, tetPolyPatchField, tetPointMesh>& vf
  47 )
  48 {
  49     elementScalarField Gamma
  50     (
  51         IOobject
  52         (
  53             "gamma",
  54             vf.time().constant(),
  55             vf.db(),
  56             IOobject::NO_READ
  57         ),
  58         vf.mesh(),
  59         dimensionedScalar("1", dimless, 1.0)
  60     );
  61
  62     return tetFem::laplacian(Gamma, vf);
  63 }
  64
  65
  66 template<class Type>
  67 tmp<tetFemMatrix<Type> > tetFem::laplacian
  68 (
  69     const elementScalarField& gamma,
  70     GeometricField<Type, tetPolyPatchField, tetPointMesh>& vf
  71 )
  72 {
  73     const tetPolyMesh& mesh = vf.mesh();
  74
  75     tmp<tetFemMatrix<Type> > tfem
  76     (
  77         new tetFemMatrix<Type>
  78         (
  79             vf,
  80             gamma.dimensions()*vf.dimensions()/dimLength/dimLength
  81         )
  82     );
  83     tetFemMatrix<Type>& fem = tfem();
  84
  85     // Get reference to upper and diagonal
  86     scalarField& u = fem.upper();
  87     scalarField& d = fem.diag();
  88
  89     // In order to avoid excessive memory allocation and copying,
  90     // matrix coefficients and addressing are retrieved using a buffer.
  91     //
  92
  93     // Get reference to ldu addressing
  94     const lduAddressing& lduAddr = fem.lduAddr();
  95     const unallocLabelList& ownerStart = lduAddr.ownerStartAddr();
  96     const unallocLabelList& neighbour = lduAddr.upperAddr();
  97
  98     // Create local-to-global and global-to-local addressing arrays
  99     labelList localToGlobalBuffer(mesh.maxNPointsForCell());
 100     labelList globalToLocalBuffer(lduAddr.size(), -1);
 101
 102     SquareMatrix<scalar> denseMatrix
 103     (
 104         mesh.maxNPointsForCell(),
 105         0
 106     );
 107
 108     for (label cellI = 0; cellI < mesh.nCells(); cellI++)
 109     {
 110         scalar curGamma = gamma[cellI];
 111
 112         label nCellPoints =
 113             mesh.addressing
 114             (
 115                 cellI,
 116                 localToGlobalBuffer,
 117                 globalToLocalBuffer
 118             );
 119
 120         mesh.gradNiDotGradNj(cellI, denseMatrix, globalToLocalBuffer);
 121
 122         // Insertion and clearing of the dense matrix is done together!
 123
 124         for (label localI = 0; localI < nCellPoints; localI++)
 125         {
 126             label globalI = localToGlobalBuffer[localI];
 127
 128             // Insert the diagonal
 129             d[globalI] += curGamma*denseMatrix[localI][localI];
 130
 131             // Zero out the coefficient after insertion
 132             denseMatrix[localI][localI] = 0;
 133
 134             // Insert the off-diagonal.
 135             // Here, look for the known global neighbours of the point
 136             // using ldu owner start and grab their local index using
 137             // globalToLocalAddressing.  Then insert at the correct
 138             // place in the global matrix
 139             label startLabel = ownerStart[globalI];
 140             label endLabel = ownerStart[globalI + 1];
 141
 142             for
 143             (
 144                 label globalJEdge = startLabel;
 145                 globalJEdge < endLabel;
 146                 globalJEdge++
 147             )
 148             {
 149                 label localJ = globalToLocalBuffer[neighbour[globalJEdge]];
 150
 151                 // Not all neighbour of globalI are in the dense
 152                 // matrix; only if the global-to-local coefficient is not -1
 153                 if (localJ > -1)
 154                 {
 155                     // It is unknown if the localI or localJ is lower:
 156                     // that depends on how the local points have been
 157                     // selected.  Therefore, the matrix (which is
 158                     // symmetric!) needs to be "folded"
 159                     u[globalJEdge] +=
 160                         curGamma*
 161                         denseMatrix[min(localI, localJ)][max(localI, localJ)];
 162
 163                     // Zero out the coefficient after insertion
 164                     denseMatrix[min(localI, localJ)][max(localI, localJ)] = 0;
 165                 }
 166             }
 167         }
 168
 169         // Clear addressing for element
 170         mesh.clearAddressing
 171         (
 172             cellI,
 173             nCellPoints,
 174             localToGlobalBuffer,
 175             globalToLocalBuffer
 176         );
 177     }
 178
 179     return tfem;
 180 }
 181
 182
 183 template<class Type>
 184 tmp<tetFemMatrix<Type> > tetFem::smoother
 185 (
 186     GeometricField<Type, tetPolyPatchField, tetPointMesh>& vf
 187 )
 188 {
 189     tmp<tetFemMatrix<Type> > tfem
 190     (
 191         new tetFemMatrix<Type>
 192         (
 193             vf,
 194             vf.dimensions()
 195         )
 196     );
 197     tetFemMatrix<Type>& fem = tfem();
 198
 199     fem.upper() = 1.0;
 200
 201     fem.negSumDiag();
 202
 203     return tfem;
 204 }
 205
 206
 207 template<class Type>
 208 tmp<tetFemMatrix<Type> > tetFem::laplacian
 209 (
 210     const dimensionedScalar& gamma,
 211     GeometricField<Type, tetPolyPatchField, tetPointMesh>& vf
 212 )
 213 {
 214     elementScalarField Gamma
 215     (
 216         IOobject
 217         (
 218             "gamma",
 219             vf.instance(),
 220             vf.db(),
 221             IOobject::NO_READ
 222         ),
 223         vf.mesh(),
 224         gamma
 225     );
 226
 227     return tetFem::laplacian(Gamma, vf);
 228 }
 229
 230
 231 template<class Type>
 232 tmp<tetFemMatrix<Type> > tetFem::laplacianTranspose
 233 (
 234     GeometricField<Type, tetPolyPatchField, tetPointMesh>& vf
 235 )
 236 {
 237     elementScalarField Gamma
 238     (
 239         IOobject
 240         (
 241             "gamma",
 242             vf.time().constant(),
 243             vf.db(),
 244             IOobject::NO_READ
 245         ),
 246         vf.mesh(),
 247         dimensionedScalar("1", dimless, 1.0)
 248     );
 249
 250     return tetFem::laplacianTranspose(Gamma, vf);
 251 }
 252
 253
 254 template<class Type>
 255 tmp<tetFemMatrix<Type> > tetFem::laplacianTranspose
 256 (
 257     const elementScalarField& gamma,
 258     GeometricField<Type, tetPolyPatchField, tetPointMesh>& vf
 259 )
 260 {
 261     const tetPolyMesh& mesh = vf.mesh();
 262
 263     tmp<tetFemMatrix<Type> > tfem
 264     (
 265         new tetFemMatrix<Type>
 266         (
 267             vf,
 268             gamma.dimensions()*vf.dimensions()/dimLength/dimLength
 269         )
 270     );
 271     tetFemMatrix<Type>& fem = tfem();
 272
 273     // Get reference to internal field
 274     const Field<Type>& psi = fem.psi().internalField();
 275
 276     // Get reference to upper and diagonal
 277     scalarField& u = fem.upper();
 278     scalarField& d = fem.diag();
 279     Field<Type>& source = fem.source();
 280
 281     // In order to avoid excessive memory allocation and copying,
 282     // matrix coefficients and addressing are retrieved using a buffer.
 283     //
 284
 285     // Get reference to ldu addressing
 286     const lduAddressing& lduAddr = fem.lduAddr();
 287     const unallocLabelList& ownerStart = lduAddr.ownerStartAddr();
 288     const unallocLabelList& neighbour = lduAddr.upperAddr();
 289
 290     // Create local-to-global and global-to-local addressing arrays
 291     labelList localToGlobalBuffer(mesh.maxNPointsForCell());
 292     labelList globalToLocalBuffer(lduAddr.size(), -1);
 293
 294     SquareMatrix<tensor> denseMatrix
 295     (
 296         mesh.maxNPointsForCell(),
 297         tensor::zero
 298     );
 299
 300     for (label cellI = 0; cellI < mesh.nCells(); cellI++)
 301     {
 302         scalar curGamma = gamma[cellI];
 303
 304         label nCellPoints =
 305             mesh.addressing
 306             (
 307                 cellI,
 308                 localToGlobalBuffer,
 309                 globalToLocalBuffer
 310             );
 311
 312         mesh.gradNiGradNj(cellI, denseMatrix, globalToLocalBuffer);
 313
 314         // Insertion and clearing of the dense matrix is done together!
 315
 316         for (label localI = 0; localI < nCellPoints; localI++)
 317         {
 318             label globalI = localToGlobalBuffer[localI];
 319
 320             tensor& curDiagCoeff = denseMatrix[localI][localI];
 321
 322             // Insert the diagonal
 323             d[globalI] += curGamma*tr(curDiagCoeff);
 324
 325             source[globalI] +=
 326                 psi[globalI] &
 327                 (
 328                     tr(curDiagCoeff)*I
 329                   - curDiagCoeff
 330                 )*curGamma;
 331
 332             // Zero out the coefficient after insertion
 333             curDiagCoeff = tensor::zero;
 334
 335             // Insert the off-diagonal.
 336             // Here, look for the known global neighbours of the point
 337             // using ldu owner start and grab their local index using
 338             // globalToLocalAddressing.  Then insert at the correct
 339             // place in the global matrix
 340             label startLabel = ownerStart[globalI];
 341             label endLabel = ownerStart[globalI + 1];
 342
 343             for
 344             (
 345                 label globalJEdge = startLabel;
 346                 globalJEdge < endLabel;
 347                 globalJEdge++
 348             )
 349             {
 350                 label globalJ = neighbour[globalJEdge];
 351                 label localJ = globalToLocalBuffer[globalJ];
 352
 353                 // Not all neighbour of globalI are in the dense
 354                 // matrix; only if the global-to-local coefficient is not -1
 355                 if (localJ > -1)
 356                 {
 357                     // It is unknown if the localI or localJ is lower:
 358                     // that depends on how the local points have been
 359                     // selected.  Therefore, the matrix (which is
 360                     // symmetric!) needs to be "folded"
 361                     tensor& curOffDiagCoeff =
 362                         denseMatrix[min(localI, localJ)][max(localI, localJ)];
 363
 364                     u[globalJEdge] += curGamma*tr(curOffDiagCoeff);
 365
 366                     source[globalI] +=
 367                         curGamma*
 368                         (
 369                             tr(curOffDiagCoeff)*I
 370                           - curOffDiagCoeff
 371                         ).T() & psi[globalJ];
 372
 373                     source[globalJ] +=
 374                         curGamma*
 375                         (
 376                             tr(curOffDiagCoeff)*I
 377                           - curOffDiagCoeff
 378                         ) & psi[globalI];
 379
 380                     curOffDiagCoeff = tensor::zero;
 381                 }
 382             }
 383         }
 384
 385         // Clear addressing for element
 386         mesh.clearAddressing
 387         (
 388             cellI,
 389             nCellPoints,
 390             localToGlobalBuffer,
 391             globalToLocalBuffer
 392         );
 393     }
 394
 395     return tfem;
 396 }
 397
 398
 399 template<class Type>
 400 tmp<tetFemMatrix<Type> > tetFem::laplacianTranspose
 401 (
 402     const dimensionedScalar& gamma,
 403     GeometricField<Type, tetPolyPatchField, tetPointMesh>& vf
 404 )
 405 {
 406     elementScalarField Gamma
 407     (
 408         IOobject
 409         (
 410             gamma.name(),
 411             vf.instance(),
 412             vf.db(),
 413             IOobject::NO_READ
 414         ),
 415         vf.mesh(),
 416         gamma
 417     );
 418
 419     return tetFem::laplacianTranspose(Gamma, vf);
 420 }
 421
 422
 423 template<class Type>
 424 tmp<tetFemMatrix<Type> > tetFem::laplacianTrace
 425 (
 426     const elementScalarField& gamma,
 427     GeometricField<Type, tetPolyPatchField, tetPointMesh>& vf
 428 )
 429 {
 430     const tetPolyMesh& mesh = vf.mesh();
 431
 432     tmp<tetFemMatrix<Type> > tfem
 433     (
 434         new tetFemMatrix<Type>
 435         (
 436             vf,
 437             gamma.dimensions()*vf.dimensions()/dimLength/dimLength
 438         )
 439     );
 440     tetFemMatrix<Type>& fem = tfem();
 441
 442     // Get reference to internal field
 443     const Field<Type>& psi = fem.psi().internalField();
 444
 445     // Get reference to upper, diagonal and source
 446     scalarField& u = fem.upper();
 447     scalarField& d = fem.diag();
 448     Field<Type>& source = fem.source();
 449
 450     // In order to avoid excessive memory allocation and copying,
 451     // matrix coefficients and addressing are retrieved using a buffer.
 452     // HJ, date deleted
 453
 454     // Get reference to ldu addressing
 455     const lduAddressing& lduAddr = fem.lduAddr();
 456     const unallocLabelList& ownerStart = lduAddr.ownerStartAddr();
 457     const unallocLabelList& neighbour = lduAddr.upperAddr();
 458
 459     // Create local-to-global and global-to-local addressing arrays
 460     labelList localToGlobalBuffer(mesh.maxNPointsForCell());
 461     labelList globalToLocalBuffer(lduAddr.size(), -1);
 462
 463     SquareMatrix<tensor> denseMatrix
 464     (
 465         mesh.maxNPointsForCell(),
 466         tensor::zero
 467     );
 468
 469     for (label cellI = 0; cellI < mesh.nCells(); cellI++)
 470     {
 471         scalar curGamma = gamma[cellI];
 472
 473         label nCellPoints =
 474             mesh.addressing
 475             (
 476                 cellI,
 477                 localToGlobalBuffer,
 478                 globalToLocalBuffer
 479             );
 480
 481         mesh.gradNiGradNj(cellI, denseMatrix, globalToLocalBuffer);
 482
 483         // Insertion and clearing of the dense matrix is done together!
 484
 485         for (label localI = 0; localI < nCellPoints; localI++)
 486         {
 487             label globalI = localToGlobalBuffer[localI];
 488
 489             tensor& curDiagCoeff = denseMatrix[localI][localI];
 490
 491             d[globalI] += curGamma*tr(curDiagCoeff);
 492
 493             source[globalI] +=
 494                 curGamma*
 495                 ((tr(curDiagCoeff)*I - curDiagCoeff) & psi[globalI]);
 496
 497             // Zero out the coefficient after insertion
 498             curDiagCoeff = tensor::zero;
 499
 500             // Insert the off-diagonal.
 501             // Here, look for the known global neighbours of the point
 502             // using ldu owner start and grab their local index using
 503             // globalToLocalAddressing.  Then insert at the correct
 504             // place in the global matrix
 505             label startLabel = ownerStart[globalI];
 506             label endLabel = ownerStart[globalI + 1];
 507
 508             for
 509             (
 510                 label globalJEdge = startLabel;
 511                 globalJEdge < endLabel;
 512                 globalJEdge++
 513             )
 514             {
 515                 label globalJ = neighbour[globalJEdge];
 516                 label localJ = globalToLocalBuffer[globalJ];
 517
 518                 // Not all neighbour of globalI are in the dense
 519                 // matrix; only if the global-to-local coefficient is not -1
 520                 if (localJ > -1)
 521                 {
 522                     // It is unknown if the localI or localJ is lower:
 523                     // that depends on how the local points have been
 524                     // selected.  Therefore, the matrix (which is
 525                     // symmetric!) needs to be "folded"
 526                     tensor& curOffDiagCoeff =
 527                         denseMatrix[min(localI, localJ)][max(localI, localJ)];
 528
 529                     u[globalJEdge] += curGamma*tr(curOffDiagCoeff);
 530
 531                     source[globalI] +=
 532                         curGamma*
 533                         (
 534                             (tr(curOffDiagCoeff)*I - curOffDiagCoeff)
 535                           & psi[globalJ]
 536                         );
 537
 538                     source[globalJ] +=
 539                         curGamma*
 540                         (
 541                             (tr(curOffDiagCoeff)*I - curOffDiagCoeff).T()
 542                           & psi[globalI]
 543                         );
 544
 545                     curOffDiagCoeff = tensor::zero;
 546                 }
 547             }
 548         }
 549
 550         // Clear addressing for element
 551         mesh.clearAddressing
 552         (
 553             cellI,
 554             nCellPoints,
 555             localToGlobalBuffer,
 556             globalToLocalBuffer
 557         );
 558     }
 559
 560     return tfem;
 561 }
 562
 563
 564 template<class Type>
 565 tmp<tetFemMatrix<Type> > tetFem::laplacianTrace
 566 (
 567     const dimensionedScalar& gamma,
 568     GeometricField<Type, tetPolyPatchField, tetPointMesh>& vf
 569 )
 570 {
 571     elementScalarField Gamma
 572     (
 573         IOobject
 574         (
 575             gamma.name(),
 576             vf.instance(),
 577             vf.db(),
 578             IOobject::NO_READ
 579         ),
 580         vf.mesh(),
 581         gamma
 582     );
 583
 584     return tetFem::laplacianTrace(Gamma, vf);
 585 }
 586
 587
 588 template<class Type>
 589 tmp<tetFemMatrix<Type> > tetFem::laplacianTrace
 590 (
 591     GeometricField<Type, tetPolyPatchField, tetPointMesh>& vf
 592 )
 593 {
 594     elementScalarField Gamma
 595     (
 596         IOobject
 597         (
 598             "1",
 599             vf.instance(),
 600             vf.db(),
 601             IOobject::NO_READ
 602         ),
 603         vf.mesh(),
 604         dimensionedScalar("1", dimless, 1.0)
 605     );
 606
 607     return tetFem::laplacianTrace(Gamma, vf);
 608 }
 609
 610
 611 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 612
 613 } // End namespace Foam
 614
 615 // ************************************************************************* //