49 centralWeight_(cWeight),
50 #ifdef SPHERICAL_GEOMETRY
66 Info <<
"Contructing quadraticFitSnGradData" <<
endl;
70 if (centralWeight_ < 1 - SMALL)
72 FatalErrorIn(
"quadraticFitSnGradData::quadraticFitSnGradData")
73 <<
"centralWeight requested = " << centralWeight_
74 <<
" should not be less than one"
89 "quadraticFitSnGradPolySize",
112 snGradPolySize[faci] = calcFit(stencilPoints[faci], faci);
117 snGradPolySize.
write();
118 Info<<
"quadraticFitSnGradData::quadraticFitSnGradData() :"
119 <<
"Finished constructing polynomialFit data"
127 void Foam::quadraticFitSnGradData::findFaceDirs
136 idir = mesh.
Sf()[faci];
139 #ifndef SPHERICAL_GEOMETRY
157 const face&
f = mesh.
faces()[faci];
162 kdir = mesh.
Cf()[faci];
168 kdir -= (idir & kdir)*idir;
170 scalar magk =
mag(kdir);
174 FatalErrorIn(
"findFaceDirs") <<
" calculated kdir = zero"
187 Foam::label Foam::quadraticFitSnGradData::calcFit
189 const List<point>& C,
196 findFaceDirs(idir, jdir, kdir,
mesh(), faci);
199 wts[0] = centralWeight_;
200 wts[1] = centralWeight_;
208 for(label ip = 0; ip < C.size(); ip++)
212 scalar px = (p - p0)&idir;
213 scalar py = (p - p0)&jdir;
214 #ifdef SPHERICAL_GEOMETRY
215 scalar pz =
mag(p) -
mag(p0);
217 scalar pz = (p - p0)&kdir;
228 B[ip][is++] = wts[0]*wts[ip];
229 B[ip][is++] = wts[0]*wts[ip]*px;
230 B[ip][is++] = wts[ip]*
sqr(px);
234 B[ip][is++] = wts[ip]*py;
235 B[ip][is++] = wts[ip]*px*py;
236 B[ip][is++] = wts[ip]*
sqr(py);
240 B[ip][is++] = wts[ip]*pz;
241 B[ip][is++] = wts[ip]*px*pz;
243 B[ip][is++] = wts[ip]*
sqr(pz);
248 label stencilSize = C.size();
249 fit_[faci].setSize(stencilSize);
255 bool goodFit =
false;
256 for(
int iIt = 0; iIt < 10 && !goodFit; iIt++)
260 scalar fit0 = wts[0]*wts[0]*svd.VSinvUt()[1][0]/scale;
261 scalar fit1 = wts[0]*wts[1]*svd.VSinvUt()[1][1]/scale;
265 &&
mag(fit0 + deltaCoeff) < 0.5*deltaCoeff
266 &&
mag(fit1 - deltaCoeff) < 0.5*deltaCoeff;
270 fit_[faci][0] = fit0;
271 fit_[faci][1] = fit1;
272 for(label i = 2; i < stencilSize; i++)
274 fit_[faci][i] = wts[0]*wts[i]*svd.VSinvUt()[1][i]/scale;
277 nSVDzeros = svd.nZeros();
284 for(label i = 0; i < B.n(); i++)
290 for(label j = 0; j < B.m(); j++)
301 fit_[faci][0] += deltaCoeff;
302 fit_[faci][1] -= deltaCoeff;
306 Pout<<
"quadratifFitSnGradData could not fit face " << faci
307 <<
" fit_[faci][0] = " << fit_[faci][0]
308 <<
" fit_[faci][1] = " << fit_[faci][1]
309 <<
" deltaCoeff = " << deltaCoeff <<
endl;
313 return minSize_ - nSVDzeros;