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