src/finiteVolume/interpolation/surfaceInterpolation/schemes/pointLinear/pointLinear.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2004-2011 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
  13     the Free Software Foundation, either version 3 of the License, or
  14     (at your 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, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "pointLinear.H"
  27 #include "fvMesh.H"
  28 #include "volPointInterpolation.H"
  29
  30 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  31
  32 template<class Type>
  33 Foam::tmp<Foam::GeometricField<Type, Foam::fvsPatchField, Foam::surfaceMesh> >
  34 Foam::pointLinear<Type>::
  35 correction
  36 (
  37     const GeometricField<Type, fvPatchField, volMesh>& vf
  38 ) const
  39 {
  40     const fvMesh& mesh = this->mesh();
  41
  42     GeometricField<Type, pointPatchField, pointMesh> pvf
  43     (
  44         volPointInterpolation::New(mesh).interpolate(vf)
  45     );
  46
  47     tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > tsfCorr =
  48         linearInterpolate(vf);
  49
  50     Field<Type>& sfCorr = tsfCorr().internalField();
  51
  52     const pointField& points = mesh.points();
  53     const pointField& C = mesh.C().internalField();
  54     const faceList& faces = mesh.faces();
  55     const scalarField& w = mesh.weights().internalField();
  56     const labelList& owner = mesh.owner();
  57     const labelList& neighbour = mesh.neighbour();
  58
  59     forAll(sfCorr, facei)
  60     {
  61         point pi =
  62             w[owner[facei]]*C[owner[facei]]
  63           + (1.0 - w[owner[facei]])*C[neighbour[facei]];
  64
  65         scalar at = triangle<point, const point&>
  66         (
  67             pi,
  68             points[faces[facei][0]],
  69             points[faces[facei][faces[facei].size()-1]]
  70         ).mag();
  71
  72         scalar sumAt = at;
  73         Type sumPsip = at*(1.0/3.0)*
  74         (
  75             sfCorr[facei]
  76           + pvf[faces[facei][0]]
  77           + pvf[faces[facei][faces[facei].size()-1]]
  78         );
  79
  80         for (label pointi=1; pointi<faces[facei].size(); pointi++)
  81         {
  82             at = triangle<point, const point&>
  83             (
  84                 pi,
  85                 points[faces[facei][pointi]],
  86                 points[faces[facei][pointi-1]]
  87             ).mag();
  88
  89             sumAt += at;
  90             sumPsip += at*(1.0/3.0)*
  91             (
  92                 sfCorr[facei]
  93               + pvf[faces[facei][pointi]]
  94               + pvf[faces[facei][pointi-1]]
  95             );
  96
  97         }
  98
  99         sfCorr[facei] = sumPsip/sumAt - sfCorr[facei];
 100     }
 101
 102     tsfCorr().boundaryField() = pTraits<Type>::zero;
 103
 104     return tsfCorr;
 105 }
 106
 107
 108 namespace Foam
 109 {
 110     makeSurfaceInterpolationScheme(pointLinear);
 111 }
 112
 113 // ************************************************************************* //