src/errorEstimation/errorEstimate/resErrorLaplacian.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2008 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 "resErrorLaplacian.H"
  28
  29 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  30
  31 namespace Foam
  32 {
  33
  34 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  35
  36 namespace resError
  37 {
  38
  39 template<class Type>
  40 tmp<errorEstimate<Type> >
  41 laplacian
  42 (
  43     const GeometricField<Type, fvPatchField, volMesh>& vf
  44 )
  45 {
  46     surfaceScalarField Gamma
  47     (
  48         IOobject
  49         (
  50             "gamma",
  51             vf.time().constant(),
  52             vf.db(),
  53             IOobject::NO_READ
  54         ),
  55         vf.mesh(),
  56         dimensionedScalar("1", dimless, 1.0)
  57     );
  58
  59     return resError::laplacian(Gamma, vf);
  60 }
  61
  62
  63 template<class Type>
  64 tmp<errorEstimate<Type> >
  65 laplacian
  66 (
  67     const dimensionedScalar& gamma,
  68     const GeometricField<Type, fvPatchField, volMesh>& vf
  69 )
  70 {
  71     surfaceScalarField Gamma
  72     (
  73         IOobject
  74         (
  75             gamma.name(),
  76             vf.time().timeName(),
  77             vf.db(),
  78             IOobject::NO_READ
  79         ),
  80         vf.mesh(),
  81         gamma
  82     );
  83
  84     return resError::laplacian(Gamma, vf);
  85 }
  86
  87
  88 template<class Type>
  89 tmp<errorEstimate<Type> >
  90 laplacian
  91 (
  92     const volScalarField& gamma,
  93     const GeometricField<Type, fvPatchField, volMesh>& vf
  94 )
  95 {
  96     return resError::laplacian(fvc::interpolate(gamma), vf);
  97 }
  98
  99
 100 template<class Type>
 101 tmp<errorEstimate<Type> >
 102 laplacian
 103 (
 104     const tmp<volScalarField>& tgamma,
 105     const GeometricField<Type, fvPatchField, volMesh>& vf
 106 )
 107 {
 108     tmp<errorEstimate<Type> > Laplacian(resError::laplacian(tgamma(), vf));
 109     tgamma.clear();
 110     return Laplacian;
 111 }
 112
 113
 114 template<class Type>
 115 tmp<errorEstimate<Type> >
 116 laplacian
 117 (
 118     const surfaceScalarField& gamma,
 119     const GeometricField<Type, fvPatchField, volMesh>& vf
 120 )
 121 {
 122     const fvMesh& mesh = vf.mesh();
 123
 124     const scalarField& vols = mesh.V();
 125     const surfaceVectorField& Sf = mesh.Sf();
 126     const surfaceScalarField magSf = mesh.magSf();
 127     const fvPatchList& patches = mesh.boundary();
 128     const unallocLabelList& owner = mesh.owner();
 129     const unallocLabelList& neighbour = mesh.neighbour();
 130
 131     const surfaceScalarField& delta =
 132         mesh.surfaceInterpolation::deltaCoeffs();
 133
 134     Field<Type> res(vols.size(), pTraits<Type>::zero);
 135     scalarField aNorm(vols.size(), 0.0);
 136
 137     // Calculate gradient of the solution
 138     GeometricField
 139     <
 140         typename outerProduct<vector, Type>::type, fvPatchField, volMesh
 141     > gradVf = fvc::grad(vf);
 142
 143     // Internal faces
 144     forAll (owner, faceI)
 145     {
 146         // Owner
 147
 148         // Subtract diffusion
 149         res[owner[faceI]] -=
 150             gamma[faceI]*(Sf[faceI] & gradVf[owner[faceI]]);
 151
 152         aNorm[owner[faceI]] += delta[faceI]*gamma[faceI]*magSf[faceI];
 153
 154         // Neighbour
 155
 156         // Subtract diffusion
 157         res[neighbour[faceI]] +=
 158             gamma[faceI]*(Sf[faceI] & gradVf[neighbour[faceI]]);
 159
 160         aNorm[neighbour[faceI]] += delta[faceI]*gamma[faceI]*magSf[faceI];
 161
 162     }
 163
 164     forAll (patches, patchI)
 165     {
 166         const vectorField& patchSf = Sf.boundaryField()[patchI];
 167         const scalarField& patchMagSf = magSf.boundaryField()[patchI];
 168         const scalarField& patchGamma = gamma.boundaryField()[patchI];
 169         const scalarField& patchDelta = delta.boundaryField()[patchI];
 170
 171         const labelList& fCells = patches[patchI].faceCells();
 172
 173         forAll (fCells, faceI)
 174         {
 175             // Subtract diffusion
 176             res[fCells[faceI]] -=
 177                 patchGamma[faceI]*
 178                 (
 179                     patchSf[faceI] & gradVf[fCells[faceI]]
 180                 );
 181
 182             aNorm[fCells[faceI]] +=
 183                 patchDelta[faceI]*patchGamma[faceI]*patchMagSf[faceI];
 184         }
 185     }
 186
 187     res /= vols;
 188     aNorm /= vols;
 189
 190     return tmp<errorEstimate<Type> >
 191     (
 192         new errorEstimate<Type>
 193         (
 194             vf,
 195             delta.dimensions()*gamma.dimensions()*magSf.dimensions()
 196             *vf.dimensions(),
 197             res,
 198             aNorm
 199         )
 200     );
 201 }
 202
 203 template<class Type>
 204 tmp<errorEstimate<Type> >
 205 laplacian
 206 (
 207     const tmp<surfaceScalarField>& tgamma,
 208     const GeometricField<Type, fvPatchField, volMesh>& vf
 209 )
 210 {
 211     tmp<errorEstimate<Type> > tresError(resError::laplacian(tgamma(), vf));
 212     tgamma.clear();
 213     return tresError;
 214 }
 215
 216
 217 template<class Type>
 218 tmp<errorEstimate<Type> >
 219 laplacian
 220 (
 221     const volTensorField& gamma,
 222     const GeometricField<Type, fvPatchField, volMesh>& vf
 223 )
 224 {
 225     const fvMesh& mesh = vf.mesh();
 226
 227     return resError::laplacian
 228     (
 229         (mesh.Sf() & fvc::interpolate(gamma) & mesh.Sf())
 230         /sqr(mesh.magSf()),
 231         vf
 232     );
 233 }
 234
 235 template<class Type>
 236 tmp<errorEstimate<Type> >
 237 laplacian
 238 (
 239     const tmp<volTensorField>& tgamma,
 240     const GeometricField<Type, fvPatchField, volMesh>& vf
 241 )
 242 {
 243     tmp<errorEstimate<Type> > Laplacian = resError::laplacian(tgamma(), vf);
 244     tgamma.clear();
 245     return Laplacian;
 246 }
 247
 248
 249 template<class Type>
 250 tmp<errorEstimate<Type> >
 251 laplacian
 252 (
 253     const surfaceTensorField& gamma,
 254     const GeometricField<Type, fvPatchField, volMesh>& vf
 255 )
 256 {
 257     const fvMesh& mesh = vf.mesh();
 258
 259     return resError::laplacian
 260     (
 261         (mesh.Sf() & gamma & mesh.Sf())/sqr(mesh.magSf()),
 262         vf
 263     );
 264 }
 265
 266 template<class Type>
 267 tmp<errorEstimate<Type> >
 268 laplacian
 269 (
 270     const tmp<surfaceTensorField>& tgamma,
 271     const GeometricField<Type, fvPatchField, volMesh>& vf
 272 )
 273 {
 274     tmp<errorEstimate<Type> > Laplacian = resError::laplacian(tgamma(), vf);
 275     tgamma.clear();
 276     return Laplacian;
 277 }
 278
 279
 280 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 281
 282 } // End namespace resError
 283
 284 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 285
 286 } // End namespace Foam
 287
 288 // ************************************************************************* //
 289
 290
 291 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 292