FreeFOAM The Cross-Platform CFD Toolkit
quadraticFitSnGradData.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 1991-2010 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
13  the Free Software Foundation, either version 3 of the License, or
14  (at your 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, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "quadraticFitSnGradData.H"
28 #include <finiteVolume/volFields.H>
29 #include <OpenFOAM/SVD.H>
30 #include <OpenFOAM/syncTools.H>
31 
32 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
33 
34 namespace Foam
35 {
37 }
38 
39 
40 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
41 
43 (
44  const fvMesh& mesh,
45  const scalar cWeight
46 )
47 :
49  centralWeight_(cWeight),
50  #ifdef SPHERICAL_GEOMETRY
51  dim_(2),
52  #else
53  dim_(mesh.nGeometricD()),
54  #endif
55  minSize_
56  (
57  dim_ == 1 ? 3 :
58  dim_ == 2 ? 6 :
59  dim_ == 3 ? 9 : 0
60  ),
61  stencil_(mesh),
62  fit_(mesh.nInternalFaces())
63 {
64  if (debug)
65  {
66  Info << "Contructing quadraticFitSnGradData" << endl;
67  }
68 
69  // check input
70  if (centralWeight_ < 1 - SMALL)
71  {
72  FatalErrorIn("quadraticFitSnGradData::quadraticFitSnGradData")
73  << "centralWeight requested = " << centralWeight_
74  << " should not be less than one"
75  << exit(FatalError);
76  }
77 
78  if (minSize_ == 0)
79  {
80  FatalErrorIn("quadraticFitSnGradData")
81  << " dimension must be 1,2 or 3, not" << dim_ << exit(FatalError);
82  }
83 
84  // store the polynomial size for each face to write out
85  surfaceScalarField snGradPolySize
86  (
87  IOobject
88  (
89  "quadraticFitSnGradPolySize",
90  "constant",
91  mesh,
92  IOobject::NO_READ,
93  IOobject::NO_WRITE
94  ),
95  mesh,
96  dimensionedScalar("quadraticFitSnGradPolySize", dimless, scalar(0))
97  );
98 
99  // Get the cell/face centres in stencil order.
100  // Centred face stencils no good for triangles of tets. Need bigger stencils
101  List<List<point> > stencilPoints(stencil_.stencil().size());
102  stencil_.collectData
103  (
104  mesh.C(),
105  stencilPoints
106  );
107 
108  // find the fit coefficients for every face in the mesh
109 
110  for(label faci = 0; faci < mesh.nInternalFaces(); faci++)
111  {
112  snGradPolySize[faci] = calcFit(stencilPoints[faci], faci);
113  }
114 
115  if (debug)
116  {
117  snGradPolySize.write();
118  Info<< "quadraticFitSnGradData::quadraticFitSnGradData() :"
119  << "Finished constructing polynomialFit data"
120  << endl;
121  }
122 }
123 
124 
125 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
126 
127 void Foam::quadraticFitSnGradData::findFaceDirs
128 (
129  vector& idir, // value changed in return
130  vector& jdir, // value changed in return
131  vector& kdir, // value changed in return
132  const fvMesh& mesh,
133  const label faci
134 )
135 {
136  idir = mesh.Sf()[faci];
137  idir /= mag(idir);
138 
139  #ifndef SPHERICAL_GEOMETRY
140  if (mesh.nGeometricD() <= 2) // find the normal direcion
141  {
142  if (mesh.geometricD()[0] == -1)
143  {
144  kdir = vector(1, 0, 0);
145  }
146  else if (mesh.geometricD()[1] == -1)
147  {
148  kdir = vector(0, 1, 0);
149  }
150  else
151  {
152  kdir = vector(0, 0, 1);
153  }
154  }
155  else // 3D so find a direction in the plane of the face
156  {
157  const face& f = mesh.faces()[faci];
158  kdir = mesh.points()[f[0]] - mesh.points()[f[1]];
159  }
160  #else
161  // Spherical geometry so kdir is the radial direction
162  kdir = mesh.Cf()[faci];
163  #endif
164 
165  if (mesh.nGeometricD() == 3)
166  {
167  // Remove the idir component from kdir and normalise
168  kdir -= (idir & kdir)*idir;
169 
170  scalar magk = mag(kdir);
171 
172  if (magk < SMALL)
173  {
174  FatalErrorIn("findFaceDirs") << " calculated kdir = zero"
175  << exit(FatalError);
176  }
177  else
178  {
179  kdir /= magk;
180  }
181  }
182 
183  jdir = kdir ^ idir;
184 }
185 
186 
187 Foam::label Foam::quadraticFitSnGradData::calcFit
188 (
189  const List<point>& C,
190  const label faci
191 )
192 {
193  vector idir(1,0,0);
194  vector jdir(0,1,0);
195  vector kdir(0,0,1);
196  findFaceDirs(idir, jdir, kdir, mesh(), faci);
197 
198  scalarList wts(C.size(), scalar(1));
199  wts[0] = centralWeight_;
200  wts[1] = centralWeight_;
201 
202  point p0 = mesh().faceCentres()[faci];
203  scalar scale = 0;
204 
205  // calculate the matrix of the polynomial components
206  scalarRectangularMatrix B(C.size(), minSize_, scalar(0));
207 
208  for(label ip = 0; ip < C.size(); ip++)
209  {
210  const point& p = C[ip];
211 
212  scalar px = (p - p0)&idir;
213  scalar py = (p - p0)&jdir;
214  #ifdef SPHERICAL_GEOMETRY
215  scalar pz = mag(p) - mag(p0);
216  #else
217  scalar pz = (p - p0)&kdir;
218  #endif
219 
220  if (ip == 0) scale = max(max(mag(px), mag(py)), mag(pz));
221 
222  px /= scale;
223  py /= scale;
224  pz /= scale;
225 
226  label is = 0;
227 
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);
231 
232  if (dim_ >= 2)
233  {
234  B[ip][is++] = wts[ip]*py;
235  B[ip][is++] = wts[ip]*px*py;
236  B[ip][is++] = wts[ip]*sqr(py);
237  }
238  if (dim_ == 3)
239  {
240  B[ip][is++] = wts[ip]*pz;
241  B[ip][is++] = wts[ip]*px*pz;
242  //B[ip][is++] = wts[ip]*py*pz;
243  B[ip][is++] = wts[ip]*sqr(pz);
244  }
245  }
246 
247  // Set the fit
248  label stencilSize = C.size();
249  fit_[faci].setSize(stencilSize);
250  scalarList singVals(minSize_);
251  label nSVDzeros = 0;
252 
253  const scalar& deltaCoeff = mesh().deltaCoeffs()[faci];
254 
255  bool goodFit = false;
256  for(int iIt = 0; iIt < 10 && !goodFit; iIt++)
257  {
258  SVD svd(B, SMALL);
259 
260  scalar fit0 = wts[0]*wts[0]*svd.VSinvUt()[1][0]/scale;
261  scalar fit1 = wts[0]*wts[1]*svd.VSinvUt()[1][1]/scale;
262 
263  goodFit =
264  fit0 < 0 && fit1 > 0
265  && mag(fit0 + deltaCoeff) < 0.5*deltaCoeff
266  && mag(fit1 - deltaCoeff) < 0.5*deltaCoeff;
267 
268  if (goodFit)
269  {
270  fit_[faci][0] = fit0;
271  fit_[faci][1] = fit1;
272  for(label i = 2; i < stencilSize; i++)
273  {
274  fit_[faci][i] = wts[0]*wts[i]*svd.VSinvUt()[1][i]/scale;
275  }
276  singVals = svd.S();
277  nSVDzeros = svd.nZeros();
278  }
279  else // (not good fit so increase weight in the centre and for linear)
280  {
281  wts[0] *= 10;
282  wts[1] *= 10;
283 
284  for(label i = 0; i < B.n(); i++)
285  {
286  B[i][0] *= 10;
287  B[i][1] *= 10;
288  }
289 
290  for(label j = 0; j < B.m(); j++)
291  {
292  B[0][j] *= 10;
293  B[1][j] *= 10;
294  }
295  }
296  }
297 
298  if (goodFit)
299  {
300  // remove the uncorrected snGradScheme coefficients
301  fit_[faci][0] += deltaCoeff;
302  fit_[faci][1] -= deltaCoeff;
303  }
304  else
305  {
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;
310  fit_[faci] = 0;
311  }
312 
313  return minSize_ - nSVDzeros;
314 }
315 
317 {
318  notImplemented("quadraticFitSnGradData::movePoints()");
319 
320  return true;
321 }
322 
323 
324 // ************************ vim: set sw=4 sts=4 et: ************************ //