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 "FitData.H"
00027 #include <finiteVolume/surfaceFields.H>
00028 #include <finiteVolume/volFields.H>
00029 #include <OpenFOAM/SVD.H>
00030
00031
00032
00033 template<class Form, class ExtendedStencil, class Polynomial>
00034 Foam::FitData<Form, ExtendedStencil, Polynomial>::FitData
00035 (
00036 const fvMesh& mesh,
00037 const ExtendedStencil& stencil,
00038 const bool linearCorrection,
00039 const scalar linearLimitFactor,
00040 const scalar centralWeight
00041 )
00042 :
00043 MeshObject<fvMesh, Form>(mesh),
00044 stencil_(stencil),
00045 linearCorrection_(linearCorrection),
00046 linearLimitFactor_(linearLimitFactor),
00047 centralWeight_(centralWeight),
00048 # ifdef SPHERICAL_GEOMETRY
00049 dim_(2),
00050 # else
00051 dim_(mesh.nGeometricD()),
00052 # endif
00053 minSize_(Polynomial::nTerms(dim_))
00054 {
00055
00056 if (linearLimitFactor <= SMALL || linearLimitFactor > 3)
00057 {
00058 FatalErrorIn("FitData<Polynomial>::FitData(..)")
00059 << "linearLimitFactor requested = " << linearLimitFactor
00060 << " should be between zero and 3"
00061 << exit(FatalError);
00062 }
00063 }
00064
00065
00066
00067
00068 template<class FitDataType, class ExtendedStencil, class Polynomial>
00069 void Foam::FitData<FitDataType, ExtendedStencil, Polynomial>::findFaceDirs
00070 (
00071 vector& idir,
00072 vector& jdir,
00073 vector& kdir,
00074 const label facei
00075 )
00076 {
00077 const fvMesh& mesh = this->mesh();
00078
00079 idir = mesh.faceAreas()[facei];
00080 idir /= mag(idir);
00081
00082 # ifndef SPHERICAL_GEOMETRY
00083 if (mesh.nGeometricD() <= 2)
00084 {
00085 if (mesh.geometricD()[0] == -1)
00086 {
00087 kdir = vector(1, 0, 0);
00088 }
00089 else if (mesh.geometricD()[1] == -1)
00090 {
00091 kdir = vector(0, 1, 0);
00092 }
00093 else
00094 {
00095 kdir = vector(0, 0, 1);
00096 }
00097 }
00098 else
00099 {
00100 const face& f = mesh.faces()[facei];
00101 kdir = mesh.points()[f[0]] - mesh.faceCentres()[facei];
00102 }
00103 # else
00104
00105 kdir = mesh.faceCentres()[facei];
00106 # endif
00107
00108 if (mesh.nGeometricD() == 3)
00109 {
00110
00111 kdir -= (idir & kdir)*idir;
00112
00113 scalar magk = mag(kdir);
00114
00115 if (magk < SMALL)
00116 {
00117 FatalErrorIn("findFaceDirs(..)") << " calculated kdir = zero"
00118 << exit(FatalError);
00119 }
00120 else
00121 {
00122 kdir /= magk;
00123 }
00124 }
00125
00126 jdir = kdir ^ idir;
00127 }
00128
00129
00130 template<class FitDataType, class ExtendedStencil, class Polynomial>
00131 void Foam::FitData<FitDataType, ExtendedStencil, Polynomial>::calcFit
00132 (
00133 scalarList& coeffsi,
00134 const List<point>& C,
00135 const scalar wLin,
00136 const label facei
00137 )
00138 {
00139 vector idir(1,0,0);
00140 vector jdir(0,1,0);
00141 vector kdir(0,0,1);
00142 findFaceDirs(idir, jdir, kdir, facei);
00143
00144
00145 scalarList wts(C.size(), scalar(1));
00146 wts[0] = centralWeight_;
00147 if (linearCorrection_)
00148 {
00149 wts[1] = centralWeight_;
00150 }
00151
00152
00153 point p0 = this->mesh().faceCentres()[facei];
00154
00155
00156
00157
00158
00159 vector d;
00160
00161
00162 scalar scale = 1;
00163
00164
00165 scalarRectangularMatrix B(C.size(), minSize_, scalar(0));
00166
00167 for(label ip = 0; ip < C.size(); ip++)
00168 {
00169 const point& p = C[ip];
00170
00171 d.x() = (p - p0)&idir;
00172 d.y() = (p - p0)&jdir;
00173 # ifndef SPHERICAL_GEOMETRY
00174 d.z() = (p - p0)&kdir;
00175 # else
00176 d.z() = mag(p) - mag(p0);
00177 # endif
00178
00179 if (ip == 0)
00180 {
00181 scale = cmptMax(cmptMag((d)));
00182 }
00183
00184
00185 d /= scale;
00186
00187 Polynomial::addCoeffs
00188 (
00189 B[ip],
00190 d,
00191 wts[ip],
00192 dim_
00193 );
00194 }
00195
00196
00197 for(label i = 0; i < B.n(); i++)
00198 {
00199 B[i][0] *= wts[0];
00200 B[i][1] *= wts[0];
00201 }
00202
00203
00204 label stencilSize = C.size();
00205 coeffsi.setSize(stencilSize);
00206
00207 bool goodFit = false;
00208 for(int iIt = 0; iIt < 8 && !goodFit; iIt++)
00209 {
00210 SVD svd(B, SMALL);
00211
00212 scalar maxCoeff = 0;
00213 label maxCoeffi = 0;
00214
00215 for(label i=0; i<stencilSize; i++)
00216 {
00217 coeffsi[i] = wts[0]*wts[i]*svd.VSinvUt()[0][i];
00218 if (mag(coeffsi[i]) > maxCoeff)
00219 {
00220 maxCoeff = mag(coeffsi[i]);
00221 maxCoeffi = i;
00222 }
00223 }
00224
00225 if (linearCorrection_)
00226 {
00227 goodFit =
00228 (mag(coeffsi[0] - wLin) < linearLimitFactor_*wLin)
00229 && (mag(coeffsi[1] - (1 - wLin)) < linearLimitFactor_*(1 - wLin))
00230 && maxCoeffi <= 1;
00231 }
00232 else
00233 {
00234
00235 goodFit =
00236 (mag(coeffsi[0] - 1.0) < linearLimitFactor_*1.0)
00237 && maxCoeffi <= 1;
00238 }
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250
00251 if (!goodFit)
00252
00253 {
00254
00255
00256
00257
00258
00259
00260
00261
00262
00263
00264
00265
00266
00267 wts[0] *= 10;
00268 if (linearCorrection_)
00269 {
00270 wts[1] *= 10;
00271 }
00272
00273 for(label j = 0; j < B.m(); j++)
00274 {
00275 B[0][j] *= 10;
00276 B[1][j] *= 10;
00277 }
00278
00279 for(label i = 0; i < B.n(); i++)
00280 {
00281 B[i][0] *= 10;
00282 B[i][1] *= 10;
00283 }
00284 }
00285 }
00286
00287 if (goodFit)
00288 {
00289 if (linearCorrection_)
00290 {
00291
00292 coeffsi[0] -= wLin;
00293 coeffsi[1] -= 1 - wLin;
00294 }
00295 else
00296 {
00297
00298 coeffsi[0] -= 1.0;
00299 }
00300 }
00301 else
00302 {
00303
00304
00305 WarningIn
00306 (
00307 "FitData<Polynomial>::calcFit(..)"
00308 ) << "Could not fit face " << facei
00309 << " Weights = " << coeffsi
00310 << ", reverting to linear." << nl
00311 << " Linear weights " << wLin << " " << 1 - wLin << endl;
00312
00313
00314 coeffsi = 0;
00315 }
00316 }
00317
00318
00319 template<class FitDataType, class ExtendedStencil, class Polynomial>
00320 bool Foam::FitData<FitDataType, ExtendedStencil, Polynomial>::movePoints()
00321 {
00322 calcFit();
00323 return true;
00324 }
00325
00326