src/finiteVolume/interpolation/surfaceInterpolation/schemes/skewCorrected/skewCorrectionVectors.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
   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 \*---------------------------------------------------------------------------*/
  26
  27 #include "skewCorrectionVectors.H"
  28 #include "surfaceFields.H"
  29 #include "volFields.H"
  30
  31 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  32
  33 namespace Foam
  34 {
  35     defineTypeNameAndDebug(skewCorrectionVectors, 0);
  36 }
  37
  38
  39 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
  40
  41 Foam::skewCorrectionVectors::skewCorrectionVectors(const fvMesh& mesh)
  42 :
  43     MeshObject<fvMesh, skewCorrectionVectors>(mesh),
  44     skew_(true),
  45     skewCorrectionVectors_(NULL)
  46 {}
  47
  48
  49 Foam::skewCorrectionVectors::~skewCorrectionVectors()
  50 {
  51     deleteDemandDrivenData(skewCorrectionVectors_);
  52 }
  53
  54
  55 void Foam::skewCorrectionVectors::makeSkewCorrectionVectors() const
  56 {
  57     if (debug)
  58     {
  59         Info<< "surfaceInterpolation::makeSkewCorrectionVectors() : "
  60             << "Constructing skew correction vectors"
  61             << endl;
  62     }
  63
  64     skewCorrectionVectors_ = new surfaceVectorField
  65     (
  66         IOobject
  67         (
  68             "skewCorrectionVectors",
  69             mesh_.pointsInstance(),
  70             mesh_,
  71             IOobject::NO_READ,
  72             IOobject::NO_WRITE,
  73             false
  74         ),
  75         mesh_,
  76         dimless
  77     );
  78     surfaceVectorField& SkewCorrVecs = *skewCorrectionVectors_;
  79
  80     // Set local references to mesh data
  81     const volVectorField& C = mesh_.C();
  82     const surfaceVectorField& Cf = mesh_.Cf();
  83     const surfaceVectorField& Sf = mesh_.Sf();
  84
  85     const unallocLabelList& owner = mesh_.owner();
  86
  87     // Build the d-vectors
  88     surfaceVectorField d = Sf/(mesh_.magSf()*mesh_.deltaCoeffs());
  89
  90     if (!mesh_.orthogonal())
  91     {
  92         d -= mesh_.correctionVectors()/mesh_.deltaCoeffs();
  93     }
  94
  95     forAll(owner, faceI)
  96     {
  97         vector Cpf = Cf[faceI] - C[owner[faceI]];
  98
  99         SkewCorrVecs[faceI] =
 100             Cpf - ((Sf[faceI] & Cpf)/(Sf[faceI] & d[faceI]))*d[faceI];
 101     }
 102
 103
 104     forAll(SkewCorrVecs.boundaryField(), patchI)
 105     {
 106         fvsPatchVectorField& patchSkewCorrVecs =
 107             SkewCorrVecs.boundaryField()[patchI];
 108
 109         if (!patchSkewCorrVecs.coupled())
 110         {
 111             patchSkewCorrVecs = vector::zero;
 112         }
 113         else
 114         {
 115             const fvPatch& p = patchSkewCorrVecs.patch();
 116             const unallocLabelList& faceCells = p.faceCells();
 117             const vectorField& patchFaceCentres = Cf.boundaryField()[patchI];
 118             const vectorField& patchSf = Sf.boundaryField()[patchI];
 119             const vectorField& patchD = d.boundaryField()[patchI];
 120
 121             forAll(p, patchFaceI)
 122             {
 123                 vector Cpf =
 124                     patchFaceCentres[patchFaceI] - C[faceCells[patchFaceI]];
 125
 126                 patchSkewCorrVecs[patchFaceI] =
 127                     Cpf
 128                   - (
 129                         (patchSf[patchFaceI] & Cpf)/
 130                         (patchSf[patchFaceI] & patchD[patchFaceI])
 131                     )*patchD[patchFaceI];
 132             }
 133         }
 134     }
 135
 136     scalar skewCoeff = 0.0;
 137
 138     if (Sf.internalField().size())
 139     {
 140         skewCoeff = max(mag(SkewCorrVecs)/mag(d)).value();
 141     }
 142
 143     if (debug)
 144     {
 145         Info<< "surfaceInterpolation::makeSkewCorrectionVectors() : "
 146             << "skew coefficient = " << skewCoeff << endl;
 147     }
 148
 149     //skewCoeff = 0.0;
 150
 151     if (skewCoeff < 1e-5)
 152     {
 153         skew_ = false;
 154         deleteDemandDrivenData(skewCorrectionVectors_);
 155     }
 156     else
 157     {
 158         skew_ = true;
 159     }
 160
 161     if (debug)
 162     {
 163         Info<< "surfaceInterpolation::makeSkewCorrectionVectors() : "
 164             << "Finished constructing skew correction vectors"
 165             << endl;
 166     }
 167 }
 168
 169
 170 bool Foam::skewCorrectionVectors::skew() const
 171 {
 172     if (skew_ == true && !skewCorrectionVectors_)
 173     {
 174         makeSkewCorrectionVectors();
 175     }
 176
 177     return skew_;
 178 }
 179
 180
 181 const Foam::surfaceVectorField& Foam::skewCorrectionVectors::operator()() const
 182 {
 183     if (!skew())
 184     {
 185         FatalErrorIn("skewCorrectionVectors::operator()()")
 186             << "Cannot return correctionVectors; mesh is not skewed"
 187             << abort(FatalError);
 188     }
 189
 190     return *skewCorrectionVectors_;
 191 }
 192
 193
 194 // Do what is neccessary if the mesh has moved
 195 bool Foam::skewCorrectionVectors::movePoints()
 196 {
 197     skew_ = true;
 198     deleteDemandDrivenData(skewCorrectionVectors_);
 199
 200     return true;
 201 }
 202
 203
 204 // ************************************************************************* //