src/tetFiniteElement/tetFemMatrix/tetFemMatrixTools.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 Description
  25     Tetrahedral Finite Element matrix boundary treatment tools
  26
  27 \*---------------------------------------------------------------------------*/
  28
  29 #include "tetFemMatrices.H"
  30 #include "tetPointRef.H"
  31 #include "constraint.H"
  32
  33 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  34
  35 namespace Foam
  36 {
  37
  38 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  39
  40 // Add boundary source for gradient-type conditions
  41 template<class Type>
  42 void tetFemMatrix<Type>::addBoundarySourceDiag()
  43 {
  44     // Loop through all boundaries and fix the condition
  45     const FieldField<tetPolyPatchField, Type>& patches =
  46         psi().boundaryField();
  47
  48     forAll (patches, patchI)
  49     {
  50         patches[patchI].addBoundarySourceDiag(*this);
  51     }
  52 }
  53
  54
  55 template<class Type>
  56 void tetFemMatrix<Type>::storeBoundaryCoeffs() const
  57 {
  58     // Make a map of all possible equations and only fix the ones that are
  59     // requested. Once the list of requested fixes is complete, collapse the
  60     // list. Note the use of list of pointers for the first part of operation
  61     // and creation of collapsed list by copying the pointers into a pointer
  62     // list.   HJ, 28/Feb/2001
  63
  64     if (!boundaryConditionsSet_)
  65     {
  66         boundaryConditionsSet_ = true;
  67
  68         // Loop through all boundaries and fix the condition
  69         const FieldField<tetPolyPatchField, Type>& patches =
  70             psi().boundaryField();
  71
  72         forAll (patches, patchI)
  73         {
  74             patches[patchI].setBoundaryCondition(fixedEqns_);
  75         }
  76
  77         // Loop through all fixed equations and grab the matrix
  78         labelList toc = fixedEqns_.toc();
  79
  80         forAll (toc, eqnI)
  81         {
  82             fixedEqns_[toc[eqnI]].setMatrix(*this);
  83         }
  84     }
  85 }
  86
  87
  88 template<class Type>
  89 void tetFemMatrix<Type>::setComponentBoundaryConditions
  90 (
  91     const direction d,
  92     scalarField& psiCmpt,
  93     scalarField& sourceCmpt
  94 )
  95 {
  96     if (!boundaryConditionsSet_)
  97     {
  98         FatalErrorIn
  99         (
 100             "void tetFemMatrix<Type>::setComponentBoundaryConditions("
 101             "const direction& d, scalarField& psiCmpt, "
 102             "scalarField& sourceCmpt)"
 103         )   << "cannot reconstruct matrix: boundary conditions not set"
 104             << abort(FatalError);
 105     }
 106
 107     // Loop through all fixed equations and set boundary condition
 108     labelList toc = fixedEqns_.toc();
 109
 110     forAll (toc, eqnI)
 111     {
 112         fixedEqns_[toc[eqnI]].eliminateEquation(*this, d, sourceCmpt);
 113     }
 114
 115     forAll (toc, eqnI)
 116     {
 117         fixedEqns_[toc[eqnI]].setSourceDiag(*this, d, psiCmpt, sourceCmpt);
 118     }
 119 }
 120
 121
 122 template<class Type>
 123 void tetFemMatrix<Type>::reconstructMatrix()
 124 {
 125     if (!boundaryConditionsSet_)
 126     {
 127         FatalErrorIn
 128         (
 129             "void tetFemMatrix<Type>::reconstructMatrix()"
 130         )   << "cannot reconstruct matrix: boundary conditions not set"
 131             << abort(FatalError);
 132     }
 133
 134     // Loop through all fixed equations and reconstruct the matrix
 135     labelList toc = fixedEqns_.toc();
 136
 137     forAll (toc, eqnI)
 138     {
 139         fixedEqns_[toc[eqnI]].reconstructMatrix(*this);
 140     }
 141 }
 142
 143
 144 template<class Type>
 145 void tetFemMatrix<Type>::addCouplingCoeffs()
 146 {
 147     // Diagonal
 148
 149     if (hasDiag())
 150     {
 151         // Initialise diagonal transfer
 152         forAll (psi_.boundaryField(), patchI)
 153         {
 154             const tetPolyPatchField<Type>& ptf = psi_.boundaryField()[patchI];
 155
 156             if (ptf.coupled())
 157             {
 158                 ptf.initAddDiag(diag());
 159             }
 160         }
 161
 162         // Add the diagonal
 163         forAll (psi_.boundaryField(), patchI)
 164         {
 165             const tetPolyPatchField<Type>& ptf = psi_.boundaryField()[patchI];
 166
 167             if (ptf.coupled())
 168             {
 169                 ptf.addDiag(diag());
 170             }
 171         }
 172     }
 173
 174     // Upper triangle
 175
 176     if (hasUpper())
 177     {
 178         // Initialise upper transfer
 179         forAll (psi_.boundaryField(), patchI)
 180         {
 181             const tetPolyPatchField<Type>& ptf = psi_.boundaryField()[patchI];
 182
 183             if (ptf.coupled())
 184             {
 185                 ptf.initAddUpperLower(upper());
 186             }
 187         }
 188
 189         // Add the upper
 190         forAll (psi_.boundaryField(), patchI)
 191         {
 192             const tetPolyPatchField<Type>& ptf = psi_.boundaryField()[patchI];
 193
 194             if (ptf.coupled())
 195             {
 196                 ptf.addUpperLower(upper());
 197             }
 198         }
 199     }
 200
 201     // Lower triangle
 202
 203     if (hasLower())
 204     {
 205         // Initialise lower transfer
 206         forAll (psi_.boundaryField(), patchI)
 207         {
 208             const tetPolyPatchField<Type>& ptf = psi_.boundaryField()[patchI];
 209
 210             if (ptf.coupled())
 211             {
 212                 ptf.initAddUpperLower(lower());
 213             }
 214         }
 215
 216         // Add the lower
 217         forAll (psi_.boundaryField(), patchI)
 218         {
 219             const tetPolyPatchField<Type>& ptf = psi_.boundaryField()[patchI];
 220
 221             if (ptf.coupled())
 222             {
 223                 ptf.addUpperLower(lower());
 224             }
 225         }
 226     }
 227 }
 228
 229
 230 template<class Type>
 231 void tetFemMatrix<Type>::addCouplingSource(scalarField& sourceCmpt) const
 232 {
 233     // Initialise source transfer
 234     forAll (psi_.boundaryField(), patchI)
 235     {
 236         const tetPolyPatchField<Type>& ptf = psi_.boundaryField()[patchI];
 237
 238         if (ptf.coupled())
 239         {
 240             ptf.initAddSource(sourceCmpt);
 241         }
 242     }
 243
 244     // Add the source
 245     forAll (psi_.boundaryField(), patchI)
 246     {
 247         const tetPolyPatchField<Type>& ptf = psi_.boundaryField()[patchI];
 248
 249         if (ptf.coupled())
 250         {
 251             ptf.addSource(sourceCmpt);
 252         }
 253     }
 254 }
 255
 256
 257 template<class Type>
 258 void tetFemMatrix<Type>::eliminateCouplingCoeffs()
 259 {
 260     // Upper triangle
 261     if (hasUpper())
 262     {
 263         // Eliminate the upper
 264         forAll (psi_.boundaryField(), patchI)
 265         {
 266             const tetPolyPatchField<Type>& ptf = psi_.boundaryField()[patchI];
 267
 268             if (ptf.coupled())
 269             {
 270                 ptf.eliminateUpperLower(upper());
 271             }
 272         }
 273     }
 274
 275     // Lower triangle
 276
 277     if (hasLower())
 278     {
 279         // Eliminate the lower
 280         forAll (psi_.boundaryField(), patchI)
 281         {
 282             const tetPolyPatchField<Type>& ptf = psi_.boundaryField()[patchI];
 283
 284             if (ptf.coupled())
 285             {
 286                 ptf.eliminateUpperLower(lower());
 287             }
 288         }
 289     }
 290 }
 291
 292
 293 template<class Type>
 294 tmp<Field<Type> > tetFemMatrix<Type>::distributeSource
 295 (
 296     const Field<Type>& ef
 297 ) const
 298 {
 299     const tetPolyMesh& mesh = psi().mesh();
 300
 301     // Make a field over all points
 302     tmp<Field<Type> > tdistSu
 303     (
 304         new Field<Type>
 305         (
 306             mesh.nPoints(),
 307             pTraits<Type>::zero
 308         )
 309     );
 310     Field<Type>& distSu = tdistSu();
 311
 312     // Go through all the tets and add a quarter of the volume times ef
 313     // into each of the vertices. The cell integrals are prepared by the mesh
 314     labelList localToGlobalBuffer(mesh.maxNPointsForCell());
 315     labelList globalToLocalBuffer(lduAddr().size(), -1);
 316
 317     scalarField coeffsBuffer
 318     (
 319         mesh.maxNPointsForCell()*tetPointRef::nVertices
 320     );
 321
 322     for (label cellI = 0; cellI < mesh.nCells(); cellI++)
 323     {
 324         label nCellPoints =
 325             mesh.addressing
 326             (
 327                 cellI,
 328                 localToGlobalBuffer,
 329                 globalToLocalBuffer
 330             );
 331
 332         mesh.volIntegral(cellI, coeffsBuffer, globalToLocalBuffer);
 333
 334         const Type& curEf = ef[cellI];
 335
 336         for (label localI = 0; localI < nCellPoints; localI++)
 337         {
 338             label globalI = localToGlobalBuffer[localI];
 339
 340             // Insert the source
 341             distSu[globalI] += coeffsBuffer[localI]*curEf;
 342         }
 343
 344         // Clear addressing for element
 345         mesh.clearAddressing
 346         (
 347             cellI,
 348             nCellPoints,
 349             localToGlobalBuffer,
 350             globalToLocalBuffer
 351         );
 352     }
 353
 354     return tdistSu;
 355 }
 356
 357
 358 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 359
 360 } // End namespace Foam
 361
 362 // ************************************************************************* //