Go to the documentation of this file.00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
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
00060 const surfaceScalarField signF = pos(flux);
00061
00062
00063 GeometricField
00064 <
00065 typename outerProduct<vector, Type>::type, fvPatchField, volMesh
00066 > gradVf = fvc::grad(vf);
00067
00068
00069 forAll (owner, faceI)
00070 {
00071
00072 const vector& curFaceCentre = faceCentres[faceI];
00073
00074
00075 vector ownD = curFaceCentre - cellCentres[owner[faceI]];
00076
00077
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
00087 vector neiD = curFaceCentre - cellCentres[neighbour[faceI]];
00088
00089
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
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 }
00158
00159
00160
00161 }
00162
00163
00164
00165
00166
00167