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