FreeFOAM The Cross-Platform CFD Toolkit
Hosted by SourceForge:
Get FreeFOAM at SourceForge.net.
            Fast, secure and Free Open Source software downloads

resErrorDiv.C

Go to the documentation of this file.
00001 /*---------------------------------------------------------------------------*\
00002   =========                 |
00003   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
00004    \\    /   O peration     |
00005     \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
00006      \\/     M anipulation  |
00007 -------------------------------------------------------------------------------
00008 License
00009     This file is part of OpenFOAM.
00010 
00011     OpenFOAM is free software: you can redistribute it and/or modify it
00012     under the terms of the GNU General Public License as published by
00013     the Free Software Foundation, either version 3 of the License, or
00014     (at your option) any later version.
00015 
00016     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
00017     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
00019     for more details.
00020 
00021     You should have received a copy of the GNU General Public License
00022     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
00023 
00024 \*---------------------------------------------------------------------------*/
00025 
00026 #include "resErrorDiv.H"
00027 #include <finiteVolume/fvc.H>
00028 
00029 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
00030 
00031 namespace Foam
00032 {
00033 
00034 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
00035 
00036 namespace resError
00037 {
00038 
00039 template<class Type>
00040 tmp<errorEstimate<Type> >
00041 div
00042 (
00043     const surfaceScalarField& flux,
00044     const GeometricField<Type, fvPatchField, volMesh>& vf
00045 )
00046 {
00047     const fvMesh& mesh = vf.mesh();
00048 
00049     const scalarField& vols = mesh.V();
00050     const surfaceVectorField& faceCentres = mesh.Cf();
00051     const volVectorField& cellCentres = mesh.C();
00052     const fvPatchList& patches = mesh.boundary();
00053     const unallocLabelList& owner = mesh.owner();
00054     const unallocLabelList& neighbour = mesh.neighbour();
00055 
00056     Field<Type> res(vols.size(), pTraits<Type>::zero);
00057     scalarField aNorm(vols.size(), 0.0);
00058 
00059     // Get sign of flux
00060     const surfaceScalarField signF = pos(flux);
00061 
00062     // Calculate gradient of the solution
00063     GeometricField
00064     <
00065         typename outerProduct<vector, Type>::type, fvPatchField, volMesh
00066     > gradVf = fvc::grad(vf);
00067 
00068     // Internal faces
00069     forAll (owner, faceI)
00070     {
00071         // Calculate the centre of the face
00072         const vector& curFaceCentre = faceCentres[faceI];
00073 
00074         // Owner
00075         vector ownD = curFaceCentre - cellCentres[owner[faceI]];
00076 
00077         // Subtract convection
00078         res[owner[faceI]] -=
00079         (
00080             vf[owner[faceI]]
00081           + (ownD & gradVf[owner[faceI]])
00082         )*flux[faceI];
00083 
00084         aNorm[owner[faceI]] += signF[faceI]*flux[faceI];
00085 
00086         // Neighbour
00087         vector neiD = curFaceCentre - cellCentres[neighbour[faceI]];
00088 
00089         // Subtract convection
00090         res[neighbour[faceI]] +=
00091         (
00092             vf[neighbour[faceI]]
00093           + (neiD & gradVf[neighbour[faceI]])
00094         )*flux[faceI];
00095 
00096         aNorm[neighbour[faceI]] -= (1.0 - signF[faceI])*flux[faceI];
00097     }
00098 
00099     forAll (patches, patchI)
00100     {
00101         const vectorField& patchFaceCentres =
00102             faceCentres.boundaryField()[patchI];
00103 
00104         const scalarField& patchFlux = flux.boundaryField()[patchI];
00105         const scalarField& patchSignFlux = signF.boundaryField()[patchI];
00106 
00107         const labelList& fCells = patches[patchI].faceCells();
00108 
00109         forAll (fCells, faceI)
00110         {
00111             vector d =
00112                 patchFaceCentres[faceI] - cellCentres[fCells[faceI]];
00113 
00114             // Subtract convection
00115             res[fCells[faceI]] -=
00116             (
00117                 vf[fCells[faceI]]
00118               + (d & gradVf[fCells[faceI]])
00119             )*patchFlux[faceI];
00120 
00121             aNorm[fCells[faceI]] += patchSignFlux[faceI]*patchFlux[faceI];
00122         }
00123     }
00124 
00125     res /= vols;
00126     aNorm /= vols;
00127 
00128     return tmp<errorEstimate<Type> >
00129     (
00130         new errorEstimate<Type>
00131         (
00132             vf,
00133             flux.dimensions()*vf.dimensions(),
00134             res,
00135             aNorm
00136         )
00137     );
00138 }
00139 
00140 
00141 template<class Type>
00142 tmp<errorEstimate<Type> >
00143 div
00144 (
00145     const tmp<surfaceScalarField>& tflux,
00146     const GeometricField<Type, fvPatchField, volMesh>& vf
00147 )
00148 {
00149     tmp<errorEstimate<Type> > Div(resError::div(tflux(), vf));
00150     tflux.clear();
00151     return Div;
00152 }
00153 
00154 
00155 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
00156 
00157 } // End namespace resError
00158 
00159 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
00160 
00161 } // End namespace Foam
00162 
00163 // ************************************************************************* //
00164 
00165 
00166 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
00167 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines