src/tetFiniteElement/tetFemMatrix/tetFemMatrixCheck.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 check for diagonal dominance.
  26     Especially useful for debugging complex interfaces
  27
  28 \*---------------------------------------------------------------------------*/
  29
  30 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  31
  32 namespace Foam
  33 {
  34
  35 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
  36
  37 template<class Type>
  38 void tetFemMatrix<Type>::check()
  39 {
  40     if (debug)
  41     {
  42         Info<< "tetFemMatrix<Type>::check : checking tetFemMatrix<Type>"
  43             << endl;
  44     }
  45
  46     Pout << "First check diagonal dominance" << endl;
  47
  48     // Get constant matrix access
  49     const tetFemMatrix<Type>& constMatrix = *this;
  50
  51     const scalarField& lower = constMatrix.lower();
  52     const scalarField& upper = constMatrix.upper();
  53     const scalarField& dd = constMatrix.diag();
  54
  55     // Calculate local matrix off-diag sum
  56     scalarField dummySource(lduAddr().size(), 0.0);
  57
  58     const scalarField oldLower = lower;
  59     const scalarField oldUpper = upper;
  60     const scalarField oldDiag = dd;
  61
  62     // Get matrix addressing
  63     const unallocLabelList& L = lduAddr().lowerAddr();
  64     const unallocLabelList& U = lduAddr().upperAddr();
  65
  66     {
  67         scalarField matrixSumOffDiag(lduAddr().size(), 0.0);
  68         scalarField dummyInternal(lduAddr().size(), 1.0);
  69
  70         forAll (L, face)
  71         {
  72             matrixSumOffDiag[L[face]] += oldLower[face];
  73             matrixSumOffDiag[U[face]] += oldUpper[face];
  74         }
  75
  76         Pout<< "void tetFemMatrix<Type>::check() : "
  77             << "Raw matrix difference: "
  78             << sum(mag(matrixSumOffDiag + oldDiag)) << endl;
  79     }
  80
  81     // Add the coupling coefficients
  82     addCouplingCoeffs();
  83
  84     addCouplingSource(dummySource);
  85
  86     // Make a copy of interfaces: no longer a reference
  87     // HJ, 20/Nov/2007
  88     lduInterfaceFieldPtrsList interfaces = psi_.boundaryField().interfaces();
  89
  90     // Prepare for coupled interface update
  91     FieldField<Field, scalar> coupledBouCoeffs(psi_.boundaryField().size());
  92     FieldField<Field, scalar> coupledIntCoeffs(psi_.boundaryField().size());
  93
  94     forAll(psi_.boundaryField(), patchI)
  95     {
  96         const tetPolyPatchField<Type>& ptf = psi_.boundaryField()[patchI];
  97
  98         coupledBouCoeffs.set(patchI, ptf.cutBouCoeffs(*this));
  99
 100         coupledIntCoeffs.set(patchI, ptf.cutIntCoeffs(*this));
 101     }
 102
 103     eliminateCouplingCoeffs();
 104
 105     Pout << "Second check diagonal dominance" << endl;
 106     // Calculate local matrix off-diag sum
 107     {
 108         scalarField matrixSumOffDiag(lduAddr().size(), 0.0);
 109         scalarField dummyInternal(lduAddr().size(), 1.0);
 110
 111         forAll (L, face)
 112         {
 113             matrixSumOffDiag[L[face]] += lower[face];
 114             matrixSumOffDiag[U[face]] += upper[face];
 115         }
 116
 117         // Initialise coupling
 118         forAll (interfaces, interfaceI)
 119         {
 120             // Check if the pointer is set
 121             if (interfaces.set(interfaceI))
 122             {
 123                 interfaces[interfaceI].initInterfaceMatrixUpdate
 124                 (
 125                     dummyInternal,
 126                     matrixSumOffDiag,
 127                     *this,
 128                     coupledBouCoeffs[interfaceI],
 129                     0,
 130                     static_cast<Pstream::commsTypes>
 131                     (
 132                         Pstream::defaultCommsType()
 133                     ),
 134                     false                       // Do not switch to lhs
 135                 );
 136             }
 137         }
 138
 139         // Update coupled interface
 140         forAll (interfaces, interfaceI)
 141         {
 142             if (interfaces.set(interfaceI))
 143             {
 144                 interfaces[interfaceI].updateInterfaceMatrix
 145                 (
 146                     dummyInternal,
 147                     matrixSumOffDiag,
 148                     *this,
 149                     coupledBouCoeffs[interfaceI],
 150                     0,
 151                     static_cast<Pstream::commsTypes>(Pstream::defaultCommsType()),
 152                     false                       // Do not switch to lhs
 153                 );
 154             }
 155         }
 156
 157         Pout<< "void tetFemMatrix<Type>::check() : "
 158             << "parallel matrix difference: "
 159             << sum(mag(matrixSumOffDiag + dd)) << endl;
 160     }
 161
 162     // restore the matrix in the original state
 163     if (hasUpper())
 164     {
 165         this->upper() = oldUpper;
 166     }
 167
 168     if (hasLower())
 169     {
 170         this->lower() = oldLower;
 171     }
 172
 173     if (hasDiag())
 174     {
 175         this->diag() = oldDiag;
 176     }
 177 }
 178
 179
 180 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 181
 182 } // End namespace Foam
 183
 184 // ************************************************************************* //