src/finiteVolume/fvMesh/fvPatches/constraint/cyclic/cyclicFvPatch.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 \*---------------------------------------------------------------------------*/
  25
  26 #include "cyclicFvPatch.H"
  27 #include "addToRunTimeSelectionTable.H"
  28 #include "fvMesh.H"
  29
  30 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  31
  32 namespace Foam
  33 {
  34
  35 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  36
  37 defineTypeNameAndDebug(cyclicFvPatch, 0);
  38 addToRunTimeSelectionTable(fvPatch, cyclicFvPatch, polyPatch);
  39
  40
  41 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  42
  43 // Make patch weighting factors
  44 void cyclicFvPatch::makeWeights(scalarField& w) const
  45 {
  46     const scalarField& magFa = magSf();
  47
  48     // Note: mag in the dot-product.
  49     // For all valid meshes, the non-orthogonality will be less than
  50     // 90 deg and the dot-product will be positive.  For invalid
  51     // meshes (d & s <= 0), this will stabilise the calculation
  52     // but the result will be poor.  HJ, 24/Aug/2011
  53     scalarField deltas = mag(nf() & fvPatch::delta());
  54     label sizeby2 = deltas.size()/2;
  55
  56     scalar maxMatchError = 0;
  57     label errorFace = -1;
  58
  59     for (label facei = 0; facei < sizeby2; facei++)
  60     {
  61         scalar avFa = (magFa[facei] + magFa[facei + sizeby2])/2.0;
  62
  63         if
  64         (
  65             mag(magFa[facei] - magFa[facei + sizeby2])/avFa
  66           > polyPatch::matchTol_()
  67         )
  68         {
  69             // Found error.  Look for largest matching error
  70             maxMatchError =
  71                 Foam::max
  72                 (
  73                     maxMatchError,
  74                     mag(magFa[facei] - magFa[facei + sizeby2])/avFa
  75                 );
  76
  77             errorFace = facei;
  78         }
  79
  80         scalar di = deltas[facei];
  81         scalar dni = deltas[facei + sizeby2];
  82
  83         w[facei] = dni/(di + dni);
  84         w[facei + sizeby2] = 1 - w[facei];
  85     }
  86
  87     // Check for error in matching
  88     if (maxMatchError > polyPatch::matchTol_())
  89     {
  90         scalar avFa = (magFa[errorFace] + magFa[errorFace + sizeby2])/2.0;
  91
  92         FatalErrorIn("cyclicFvPatch::makeWeights(scalarField& w) const")
  93             << "face " << errorFace << " and " << errorFace + sizeby2
  94             <<  " areas do not match by "
  95             << 100*mag(magFa[errorFace] - magFa[errorFace + sizeby2])/avFa
  96             << "% -- possible face ordering problem." << nl
  97             << "Cyclic area match tolerance = "
  98             << polyPatch::matchTol_() << " patch: " << name()
  99             << abort(FatalError);
 100     }
 101 }
 102
 103
 104 // Make patch face - neighbour cell distances
 105 void cyclicFvPatch::makeDeltaCoeffs(scalarField& dc) const
 106 {
 107     vectorField d = delta();
 108     vectorField n = nf();
 109     label sizeby2 = d.size()/2;
 110
 111     for (label facei = 0; facei < sizeby2; facei++)
 112     {
 113         // Stabilised form for bad meshes.  HJ, 24/Aug/2011
 114         dc[facei] = 1.0/max(n[facei] & d[facei], 0.05*mag(d[facei]));
 115         dc[facei + sizeby2] = dc[facei];
 116     }
 117 }
 118
 119
 120 // Return delta (P to N) vectors across coupled patch
 121 tmp<vectorField> cyclicFvPatch::delta() const
 122 {
 123     vectorField patchD = fvPatch::delta();
 124     label sizeby2 = patchD.size()/2;
 125
 126     tmp<vectorField> tpdv(new vectorField(patchD.size()));
 127     vectorField& pdv = tpdv();
 128
 129     // To the transformation if necessary
 130     if (parallel())
 131     {
 132         for (label facei = 0; facei < sizeby2; facei++)
 133         {
 134             vector ddi = patchD[facei];
 135             vector dni = patchD[facei + sizeby2];
 136
 137             pdv[facei] = ddi - dni;
 138             pdv[facei + sizeby2] = -pdv[facei];
 139         }
 140     }
 141     else
 142     {
 143         for (label facei = 0; facei < sizeby2; facei++)
 144         {
 145             vector ddi = patchD[facei];
 146             vector dni = patchD[facei + sizeby2];
 147
 148             pdv[facei] = ddi - transform(forwardT()[0], dni);
 149             pdv[facei + sizeby2] = -transform(reverseT()[0], pdv[facei]);
 150         }
 151     }
 152
 153     return tpdv;
 154 }
 155
 156
 157 tmp<labelField> cyclicFvPatch::interfaceInternalField
 158 (
 159     const unallocLabelList& internalData
 160 ) const
 161 {
 162     return patchInternalField(internalData);
 163 }
 164
 165
 166 tmp<labelField> cyclicFvPatch::transfer
 167 (
 168     const Pstream::commsTypes,
 169     const unallocLabelList& interfaceData
 170 ) const
 171 {
 172     tmp<labelField> tpnf(new labelField(this->size()));
 173     labelField& pnf = tpnf();
 174
 175     label sizeby2 = this->size()/2;
 176
 177     for (label facei=0; facei<sizeby2; facei++)
 178     {
 179         pnf[facei] = interfaceData[facei + sizeby2];
 180         pnf[facei + sizeby2] = interfaceData[facei];
 181     }
 182
 183     return tpnf;
 184 }
 185
 186
 187 tmp<labelField> cyclicFvPatch::internalFieldTransfer
 188 (
 189     const Pstream::commsTypes commsType,
 190     const unallocLabelList& iF
 191 ) const
 192 {
 193     const unallocLabelList& faceCells = this->patch().faceCells();
 194
 195     tmp<labelField> tpnf(new labelField(this->size()));
 196     labelField& pnf = tpnf();
 197
 198     label sizeby2 = this->size()/2;
 199
 200     for (label facei=0; facei<sizeby2; facei++)
 201     {
 202         pnf[facei] = iF[faceCells[facei + sizeby2]];
 203         pnf[facei + sizeby2] = iF[faceCells[facei]];
 204     }
 205
 206     return tpnf;
 207 }
 208
 209
 210 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 211
 212 } // End namespace Foam
 213
 214 // ************************************************************************* //