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