FitData.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) 2011-2016 OpenFOAM Foundation
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 "FitData.H"
27 #include "surfaceFields.H"
28 #include "volFields.H"
29 #include "SVD.H"
30 
31 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
32 
33 template<class Form, class ExtendedStencil, class Polynomial>
35 (
36  const fvMesh& mesh,
37  const ExtendedStencil& stencil,
38  const bool linearCorrection,
39  const scalar linearLimitFactor,
40  const scalar centralWeight
41 )
42 :
44  stencil_(stencil),
45  linearCorrection_(linearCorrection),
46  linearLimitFactor_(linearLimitFactor),
47  centralWeight_(centralWeight),
48  #ifdef SPHERICAL_GEOMETRY
49  dim_(2),
50  #else
51  dim_(mesh.nGeometricD()),
52  #endif
53  minSize_(Polynomial::nTerms(dim_))
54 {
55  // Check input
56  if (linearLimitFactor <= SMALL || linearLimitFactor > 3)
57  {
59  << "linearLimitFactor requested = " << linearLimitFactor
60  << " should be between zero and 3"
61  << exit(FatalError);
62  }
63 }
64 
65 
66 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
67 
68 template<class FitDataType, class ExtendedStencil, class Polynomial>
70 (
71  vector& idir, // value changed in return
72  vector& jdir, // value changed in return
73  vector& kdir, // value changed in return
74  const label facei
75 )
76 {
77  const fvMesh& mesh = this->mesh();
78 
79  idir = mesh.faceAreas()[facei];
80  idir /= mag(idir);
81 
82  #ifndef SPHERICAL_GEOMETRY
83  if (mesh.nGeometricD() <= 2) // find the normal direction
84  {
85  if (mesh.geometricD()[0] == -1)
86  {
87  kdir = vector(1, 0, 0);
88  }
89  else if (mesh.geometricD()[1] == -1)
90  {
91  kdir = vector(0, 1, 0);
92  }
93  else
94  {
95  kdir = vector(0, 0, 1);
96  }
97  }
98  else // 3D so find a direction in the plane of the face
99  {
100  const face& f = mesh.faces()[facei];
101  kdir = mesh.points()[f[0]] - mesh.faceCentres()[facei];
102  }
103  #else
104  // Spherical geometry so kdir is the radial direction
105  kdir = mesh.faceCentres()[facei];
106  #endif
107 
108  if (mesh.nGeometricD() == 3)
109  {
110  // Remove the idir component from kdir and normalise
111  kdir -= (idir & kdir)*idir;
112 
113  scalar magk = mag(kdir);
114 
115  if (magk < SMALL)
116  {
118  << exit(FatalError);
119  }
120  else
121  {
122  kdir /= magk;
123  }
124  }
125 
126  jdir = kdir ^ idir;
127 }
128 
129 
130 template<class FitDataType, class ExtendedStencil, class Polynomial>
132 (
133  scalarList& coeffsi,
134  const List<point>& C,
135  const scalar wLin,
136  const label facei
137 )
138 {
139  vector idir(1,0,0);
140  vector jdir(0,1,0);
141  vector kdir(0,0,1);
142  findFaceDirs(idir, jdir, kdir, facei);
143 
144  // Setup the point weights
145  scalarList wts(C.size(), scalar(1));
146  wts[0] = centralWeight_;
147  if (linearCorrection_)
148  {
149  wts[1] = centralWeight_;
150  }
151 
152  // Reference point
153  point p0 = this->mesh().faceCentres()[facei];
154 
155  // Info<< "Face " << facei << " at " << p0 << " stencil points at:\n"
156  // << C - p0 << endl;
157 
158  // p0 -> p vector in the face-local coordinate system
159  vector d;
160 
161  // Local coordinate scaling
162  scalar scale = 1;
163 
164  // Matrix of the polynomial components
165  scalarRectangularMatrix B(C.size(), minSize_, scalar(0));
166 
167  forAll(C, ip)
168  {
169  const point& p = C[ip];
170 
171  d.x() = (p - p0)&idir;
172  d.y() = (p - p0)&jdir;
173  #ifndef SPHERICAL_GEOMETRY
174  d.z() = (p - p0)&kdir;
175  #else
176  d.z() = mag(p) - mag(p0);
177  #endif
178 
179  if (ip == 0)
180  {
181  scale = cmptMax(cmptMag((d)));
182  }
183 
184  // Scale the radius vector
185  d /= scale;
186 
187  Polynomial::addCoeffs
188  (
189  B[ip],
190  d,
191  wts[ip],
192  dim_
193  );
194  }
195 
196  // Additional weighting for constant and linear terms
197  for (label i = 0; i < B.m(); i++)
198  {
199  B(i, 0) *= wts[0];
200  B(i, 1) *= wts[0];
201  }
202 
203  // Set the fit
204  label stencilSize = C.size();
205  coeffsi.setSize(stencilSize);
206 
207  bool goodFit = false;
208  for (int iIt = 0; iIt < 8 && !goodFit; iIt++)
209  {
210  SVD svd(B, SMALL);
211 
212  scalar maxCoeff = 0;
213  label maxCoeffi = 0;
214 
215  for (label i=0; i<stencilSize; i++)
216  {
217  coeffsi[i] = wts[0]*wts[i]*svd.VSinvUt()(0, i);
218  if (mag(coeffsi[i]) > maxCoeff)
219  {
220  maxCoeff = mag(coeffsi[i]);
221  maxCoeffi = i;
222  }
223  }
224 
225  if (linearCorrection_)
226  {
227  goodFit =
228  (mag(coeffsi[0] - wLin) < linearLimitFactor_*wLin)
229  && (mag(coeffsi[1] - (1 - wLin)) < linearLimitFactor_*(1 - wLin))
230  && maxCoeffi <= 1;
231  }
232  else
233  {
234  // Upwind: weight on face is 1.
235  goodFit =
236  (mag(coeffsi[0] - 1.0) < linearLimitFactor_*1.0)
237  && maxCoeffi <= 1;
238  }
239 
240  // if (goodFit && iIt > 0)
241  // {
242  // Info<< "FitData<Polynomial>::calcFit"
243  // << "(const List<point>& C, const label facei" << nl
244  // << "Can now fit face " << facei << " iteration " << iIt
245  // << " with sum of weights " << sum(coeffsi) << nl
246  // << " Weights " << coeffsi << nl
247  // << " Linear weights " << wLin << " " << 1 - wLin << nl
248  // << " sing vals " << svd.S() << endl;
249  // }
250 
251  if (!goodFit) // (not good fit so increase weight in the centre and
252  // weight for constant and linear terms)
253  {
254  // if (iIt == 7)
255  // {
256  // WarningInFunction
257  // << "Cannot fit face " << facei << " iteration " << iIt
258  // << " with sum of weights " << sum(coeffsi) << nl
259  // << " Weights " << coeffsi << nl
260  // << " Linear weights " << wLin << " " << 1 - wLin << nl
261  // << " sing vals " << svd.S() << endl;
262  // }
263 
264  wts[0] *= 10;
265  if (linearCorrection_)
266  {
267  wts[1] *= 10;
268  }
269 
270  for (label j = 0; j < B.n(); j++)
271  {
272  B(0, j) *= 10;
273  B(1, j) *= 10;
274  }
275 
276  for (label i = 0; i < B.m(); i++)
277  {
278  B(i, 0) *= 10;
279  B(i, 1) *= 10;
280  }
281  }
282  }
283 
284  if (goodFit)
285  {
286  if (linearCorrection_)
287  {
288  // Remove the uncorrected linear coefficients
289  coeffsi[0] -= wLin;
290  coeffsi[1] -= 1 - wLin;
291  }
292  else
293  {
294  // Remove the uncorrected upwind coefficients
295  coeffsi[0] -= 1.0;
296  }
297  }
298  else
299  {
300  // if (debug)
301  // {
303  << "Could not fit face " << facei
304  << " Weights = " << coeffsi
305  << ", reverting to linear." << nl
306  << " Linear weights " << wLin << " " << 1 - wLin << endl;
307  // }
308 
309  coeffsi = 0;
310  }
311 }
312 
313 
314 template<class FitDataType, class ExtendedStencil, class Polynomial>
316 {
317  calcFit();
318  return true;
319 }
320 
321 // ************************************************************************* //
Foam::surfaceFields.
const Cmpt & z() const
Definition: VectorI.H:87
void cmptMax(FieldField< Field, typename FieldField< Field, Type >::cmptType > &cf, const FieldField< Field, Type > &f)
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
virtual void calcFit()=0
Calculate the fit for all the faces.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
const Cmpt & x() const
Definition: VectorI.H:75
error FatalError
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
const vectorField & faceAreas() const
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:76
const Vector< label > & geometricD() const
Return the vector of geometric directions in mesh.
Definition: polyMesh.C:784
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:253
const scalarRectangularMatrix & VSinvUt() const
Return VSinvUt (the pseudo inverse)
Definition: SVDI.H:52
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
const vectorField & faceCentres() const
label nGeometricD() const
Return the number of valid geometric dimensions in the mesh.
Definition: polyMesh.C:795
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:979
Templated abstract base-class for optional mesh objects used to automate their allocation to the mesh...
Definition: MeshObject.H:85
dynamicFvMesh & mesh
const Cmpt & y() const
Definition: VectorI.H:81
Singular value decomposition of a rectangular matrix.
Definition: SVD.H:51
void cmptMag(FieldField< Field, Type > &cf, const FieldField< Field, Type > &f)
static const char nl
Definition: Ostream.H:262
labelList f(nPoints)
void findFaceDirs(vector &idir, vector &jdir, vector &kdir, const label faci)
Find the normal direction (i) and j and k directions for face faci.
Definition: FitData.C:70
void setSize(const label)
Reset size of List.
Definition: List.C:295
#define WarningInFunction
Report a warning using Foam::Warning.
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:78
FitData(const fvMesh &mesh, const ExtendedStencil &stencil, const bool linearCorrection, const scalar linearLimitFactor, const scalar centralWeight)
Construct from components.
Definition: FitData.C:35
dimensioned< scalar > mag(const dimensioned< Type > &)
volScalarField & p
virtual const faceList & faces() const
Return raw faces.
Definition: polyMesh.C:1004
bool movePoints()
Recalculate weights (but not stencil) when the mesh moves.
Definition: FitData.C:315