src/tetFiniteElement/tetFemMatrix/tetFemScalarMatrix.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      Tet Finite Element scalar matrix member functions and operators
  26
  27 \*---------------------------------------------------------------------------*/
  28
  29 #include "objectRegistry.H"
  30 #include "foamTime.H"
  31 #include "tetFemScalarMatrix.H"
  32 #include "tetPointFields.H"
  33 #include "tetPolyPatchFields.H"
  34
  35 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  36
  37 namespace Foam
  38 {
  39
  40 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  41
  42 template<>
  43 lduSolverPerformance tetFemMatrix<scalar>::solve
  44 (
  45     const dictionary& solverControls
  46 )
  47 {
  48     if (debug)
  49     {
  50         Info<< "tetFemMatrix<scalar>::solve(const dictionary&) : "
  51             << "solving tetFemMatrix<scalar>"
  52             << endl;
  53     }
  54
  55     // Add boundary source for gradient-type conditions
  56     addBoundarySourceDiag();
  57
  58     // Store the boundary coefficients for insertion of boundary conditions
  59     storeBoundaryCoeffs();
  60
  61     // Set component boundary conditions
  62     scalarField sourceCpy = source_;
  63
  64     // Store the boundary coefficients for insertion of boundary conditions
  65     setComponentBoundaryConditions(0, psi_, sourceCpy);
  66
  67     // Add the coupling coefficients
  68     addCouplingCoeffs();
  69     addCouplingSource(sourceCpy);
  70
  71     // prepare for coupled interface update
  72     FieldField<Field, scalar> coupledBouCoeffs(psi_.boundaryField().size());
  73     FieldField<Field, scalar> coupledIntCoeffs(psi_.boundaryField().size());
  74
  75     forAll(psi_.boundaryField(), patchI)
  76     {
  77         const tetPolyPatchScalarField& ptf = psi_.boundaryField()[patchI];
  78
  79         coupledBouCoeffs.set
  80         (
  81             patchI,
  82             new scalarField(ptf.cutBouCoeffs(*this))
  83         );
  84
  85         coupledIntCoeffs.set
  86         (
  87             patchI,
  88             new scalarField(ptf.cutIntCoeffs(*this))
  89         );
  90     }
  91
  92     eliminateCouplingCoeffs();
  93
  94     lduSolverPerformance solverPerf
  95     (
  96         "tetFemMatrix<scalar>::solve",
  97         psi_.name()
  98     );
  99
 100     // Make a copy of interfaces: no longer a reference
 101     // HJ, 20/Nov/2007
 102     lduInterfaceFieldPtrsList interfaces = psi_.boundaryField().interfaces();
 103
 104     solverPerf = lduSolver::New
 105     (
 106         psi_.name(),
 107         *this,
 108         coupledBouCoeffs,
 109         coupledIntCoeffs,
 110         interfaces,
 111         solverControls
 112     )->solve(psi_.internalField(), sourceCpy);
 113
 114     solverPerf.print();
 115
 116     reconstructMatrix();
 117
 118     if (debug)
 119     {
 120         Info<< "tetFemMatrix<scalar>::solve(const scalar, const scalar) : "
 121             << "correcting boundary conditions"
 122             << endl;
 123     }
 124
 125     psi_.correctBoundaryConditions();
 126
 127     return solverPerf;
 128 }
 129
 130
 131 // Return the matrix residual
 132 template<>
 133 tmp<scalarField> tetFemMatrix<scalar>::residual()
 134 {
 135     // Store the boundary coefficients for insertion of boundary conditions
 136     storeBoundaryCoeffs();
 137
 138     // Set component boundary conditions
 139     scalarField sourceCpy = source_;
 140
 141     setComponentBoundaryConditions(0, psi_, sourceCpy);
 142
 143     // Add the coupling coefficients
 144     addCouplingCoeffs();
 145     addCouplingSource(sourceCpy);
 146
 147     // Prepare for coupled interface update
 148     FieldField<Field, scalar> coupledBouCoeffs(psi_.boundaryField().size());
 149
 150     forAll(psi_.boundaryField(), patchI)
 151     {
 152         const tetPolyPatchScalarField& ptf = psi_.boundaryField()[patchI];
 153         coupledBouCoeffs.set
 154         (
 155             patchI,
 156             new scalarField(ptf.cutBouCoeffs(*this))
 157         );
 158     }
 159
 160     eliminateCouplingCoeffs();
 161
 162     // Make a copy of interfaces: no longer a reference
 163     // HJ, 20/Nov/2007
 164     lduInterfaceFieldPtrsList interfaces = psi_.boundaryField().interfaces();
 165
 166     tmp<scalarField> tres
 167     (
 168         lduMatrix::residual
 169         (
 170             psi_.internalField(),
 171             sourceCpy,
 172             coupledBouCoeffs,
 173             interfaces,
 174             0
 175         )
 176     );
 177
 178     reconstructMatrix();
 179
 180     return tres;
 181 }
 182
 183
 184 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 185
 186 } // End namespace Foam
 187
 188 // ************************************************************************* //