src/tetFiniteElement/tetFemMatrix/tetFemMatrixSolve.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 basic solvers.
  26
  27 \*---------------------------------------------------------------------------*/
  28
  29 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  30
  31 namespace Foam
  32 {
  33
  34 // * * * * * * * * * * * * * * * * Solvers * * * * * * * * * * * * * * * * * //
  35
  36 template<class Type>
  37 lduSolverPerformance tetFemMatrix<Type>::solve
  38 (
  39      const dictionary& solverControls
  40 )
  41 {
  42     if (debug)
  43     {
  44         Info<< "tetFemMatrix<Type>::solve(const dictionary&) : "
  45                "solving tetFemMatrix<Type>"
  46             << endl;
  47     }
  48
  49     // Check the matrix
  50     if (debug > 1)
  51     {
  52         this->check();
  53     }
  54
  55     lduSolverPerformance solverPerfVec
  56     (
  57         "tetFemMatrix<Type>::solve",
  58         psi_.name()
  59     );
  60
  61     // Add boundary source for gradient-type conditions
  62     addBoundarySourceDiag();
  63
  64     // Store the boundary coefficients for insertion of boundary conditions
  65     storeBoundaryCoeffs();
  66
  67     typename Type::labelType validComponents
  68     (
  69         pow
  70         (
  71             psi_.mesh()().solutionD(),
  72             pTraits<typename powProduct<Vector<label>, Type::rank>::type>::zero
  73         )
  74     );
  75
  76     // Make a copy of interfaces: no longer a reference
  77     // HJ, 20/Nov/2007
  78     lduInterfaceFieldPtrsList interfaces = psi_.boundaryField().interfaces();
  79
  80     for (direction cmpt = 0; cmpt < Type::nComponents; cmpt++)
  81     {
  82         if (validComponents[cmpt] == -1) continue;
  83
  84         scalarField psiCmpt = psi_.internalField().component(cmpt);
  85         scalarField sourceCmpt = source_.component(cmpt);
  86
  87         // Set component boundary conditions
  88         setComponentBoundaryConditions(cmpt, psiCmpt, sourceCmpt);
  89
  90         // Add the coupling coefficients
  91         addCouplingCoeffs();
  92
  93         addCouplingSource(sourceCmpt);
  94
  95         // Prepare for coupled interface update
  96         FieldField<Field, scalar> coupledBouCoeffs
  97         (
  98             psi_.boundaryField().size()
  99         );
 100
 101         FieldField<Field, scalar> coupledIntCoeffs
 102         (
 103             psi_.boundaryField().size()
 104         );
 105
 106         forAll(psi_.boundaryField(), patchI)
 107         {
 108             const tetPolyPatchField<Type>& ptf = psi_.boundaryField()[patchI];
 109
 110             coupledBouCoeffs.set
 111             (
 112                 patchI,
 113                 ptf.cutBouCoeffs(*this)
 114             );
 115
 116             coupledIntCoeffs.set
 117             (
 118                 patchI,
 119                 ptf.cutIntCoeffs(*this)
 120             );
 121         }
 122
 123         eliminateCouplingCoeffs();
 124
 125         scalarField res(psi_.size(), 0);
 126         lduMatrix::residual
 127         (
 128             res,
 129             psiCmpt,
 130             sourceCmpt,
 131             coupledBouCoeffs,
 132             interfaces,
 133             cmpt
 134         );
 135
 136         lduSolverPerformance solverPerf = lduSolver::New
 137         (
 138             psi_.name() + pTraits<Type>::componentNames[cmpt],
 139             *this,
 140             coupledBouCoeffs,
 141             coupledIntCoeffs,
 142             interfaces,
 143             solverControls
 144         )->solve(psiCmpt, sourceCmpt, cmpt);
 145
 146         solverPerf.print();
 147
 148         if
 149         (
 150             solverPerf.initialResidual() > solverPerfVec.initialResidual()
 151          && !solverPerf.singular()
 152         )
 153         {
 154             solverPerfVec = solverPerf;
 155         }
 156
 157         psi_.internalField().replace(cmpt, psiCmpt);
 158
 159         reconstructMatrix();
 160     }
 161
 162     if (debug)
 163     {
 164         Info<< "tetFemMatrix<Type>::solve : correcting boundary conditions"
 165             << endl;
 166     }
 167
 168     psi_.correctBoundaryConditions();
 169
 170     return solverPerfVec;
 171 }
 172
 173
 174 template<class Type>
 175 lduSolverPerformance tetFemMatrix<Type>::solve()
 176 {
 177     return solve(psi_.mesh().solutionDict().solver(psi_.name()));
 178 }
 179
 180
 181 // Return the matrix residual
 182 template<class Type>
 183 tmp<Field<Type> > tetFemMatrix<Type>::residual()
 184 {
 185     tmp<Field<Type> > tres(psi_.size());
 186
 187     // Store the boundary coefficients for insertion of boundary conditions
 188     storeBoundaryCoeffs();
 189
 190     // Make a copy of interfaces: no longer a reference
 191     // HJ, 20/Nov/2007
 192     lduInterfaceFieldPtrsList interfaces = psi_.boundaryField().interfaces();
 193
 194     // Loop over field components
 195     for (direction cmpt = 0; cmpt < Type::nComponents; cmpt++)
 196     {
 197         scalarField PsiInternalCmpt = psi_.internalField().component(cmpt);
 198         scalarField sourceCmpt = source_.component(cmpt);
 199
 200         setComponentBoundaryConditions(cmpt, PsiInternalCmpt, sourceCmpt);
 201
 202         // Add the coupling coefficients
 203         addCouplingCoeffs();
 204         addCouplingSource(sourceCmpt);
 205
 206         // Prepare for coupled interface update
 207         FieldField<Field, scalar> coupledBouCoeffs(psi_.boundaryField().size());
 208
 209         forAll(psi_.boundaryField(), patchI)
 210         {
 211             const tetPolyPatchField<Type>& ptf = psi_.boundaryField()[patchI];
 212             coupledBouCoeffs.set
 213             (
 214                 patchI,
 215                 new scalarField(ptf.cutBouCoeffs(*this))
 216             );
 217         }
 218
 219         eliminateCouplingCoeffs();
 220
 221         tres().replace
 222         (
 223             cmpt,
 224             lduMatrix::residual
 225             (
 226                 psi_.internalField().component(cmpt),
 227                 sourceCmpt,
 228                 coupledBouCoeffs,
 229                 interfaces,
 230                 cmpt
 231             )
 232         );
 233
 234         reconstructMatrix();
 235     }
 236
 237     return tres;
 238 }
 239
 240
 241 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 242
 243 } // End namespace Foam
 244
 245 // ************************************************************************* //