FitData.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2023 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  scalarRectangularMatrix invB(svd.VSinvUt());
212 
213  scalar maxCoeff = 0;
214  label maxCoeffi = 0;
215 
216  for (label i=0; i<stencilSize; i++)
217  {
218  coeffsi[i] = wts[0]*wts[i]*invB(0, i);
219  if (mag(coeffsi[i]) > maxCoeff)
220  {
221  maxCoeff = mag(coeffsi[i]);
222  maxCoeffi = i;
223  }
224  }
225 
226  if (linearCorrection_)
227  {
228  goodFit =
229  (mag(coeffsi[0] - wLin) < linearLimitFactor_*wLin)
230  && (mag(coeffsi[1] - (1 - wLin)) < linearLimitFactor_*(1 - wLin))
231  && maxCoeffi <= 1;
232  }
233  else
234  {
235  // Upwind: weight on face is 1.
236  goodFit =
237  (mag(coeffsi[0] - 1.0) < linearLimitFactor_*1.0)
238  && maxCoeffi <= 1;
239  }
240 
241  // if (goodFit && iIt > 0)
242  // {
243  // Info<< "FitData<Polynomial>::calcFit"
244  // << "(const List<point>& C, const label facei" << nl
245  // << "Can now fit face " << facei << " iteration " << iIt
246  // << " with sum of weights " << sum(coeffsi) << nl
247  // << " Weights " << coeffsi << nl
248  // << " Linear weights " << wLin << " " << 1 - wLin << nl
249  // << " sing vals " << svd.S() << endl;
250  // }
251 
252  if (!goodFit) // (not good fit so increase weight in the centre and
253  // weight for constant and linear terms)
254  {
255  // if (iIt == 7)
256  // {
257  // WarningInFunction
258  // << "Cannot fit face " << facei << " iteration " << iIt
259  // << " with sum of weights " << sum(coeffsi) << nl
260  // << " Weights " << coeffsi << nl
261  // << " Linear weights " << wLin << " " << 1 - wLin << nl
262  // << " sing vals " << svd.S() << endl;
263  // }
264 
265  wts[0] *= 10;
266  if (linearCorrection_)
267  {
268  wts[1] *= 10;
269  }
270 
271  for (label j = 0; j < B.n(); j++)
272  {
273  B(0, j) *= 10;
274  B(1, j) *= 10;
275  }
276 
277  for (label i = 0; i < B.m(); i++)
278  {
279  B(i, 0) *= 10;
280  B(i, 1) *= 10;
281  }
282  }
283  }
284 
285  if (goodFit)
286  {
287  if (linearCorrection_)
288  {
289  // Remove the uncorrected linear coefficients
290  coeffsi[0] -= wLin;
291  coeffsi[1] -= 1 - wLin;
292  }
293  else
294  {
295  // Remove the uncorrected upwind coefficients
296  coeffsi[0] -= 1.0;
297  }
298  }
299  else
300  {
301  // if (debug)
302  // {
304  << "Could not fit face " << facei
305  << " Weights = " << coeffsi
306  << ", reverting to linear." << nl
307  << " Linear weights " << wLin << " " << 1 - wLin << endl;
308  // }
309 
310  coeffsi = 0;
311  }
312 }
313 
314 
315 template<class FitDataType, class ExtendedStencil, class Polynomial>
317 {
318  calcFit();
319  return true;
320 }
321 
322 // ************************************************************************* //
static const Foam::dimensionedScalar B("B", Foam::dimless, 18.678)
static const Foam::dimensionedScalar C("C", Foam::dimTemperature, 234.5)
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
Graphite solid properties.
Definition: C.H:51
Templated abstract base-class for demand-driven mesh objects used to automate their allocation to the...
bool movePoints()
Recalculate weights (but not stencil) when the mesh moves.
Definition: FitData.C:316
scalar linearLimitFactor() const
Factor the fit is allowed to deviate from the base scheme.
Definition: FitData.H:131
FitData(const fvMesh &mesh, const ExtendedStencil &stencil, const bool linearCorrection, const scalar linearLimitFactor, const scalar centralWeight)
Construct from components.
Definition: FitData.C:35
virtual void calcFit()=0
Calculate the fit for all the faces.
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:281
Polynomial templated on size (order):
Definition: Polynomial.H:84
const Cmpt & z() const
Definition: VectorI.H:87
const Cmpt & y() const
Definition: VectorI.H:81
const Cmpt & x() const
Definition: VectorI.H:75
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:99
label nGeometricD() const
Return the number of valid geometric dimensions in the mesh.
Definition: polyMesh.C:1000
virtual const faceList & faces() const
Return raw faces.
Definition: polyMesh.C:1326
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:1313
const Vector< label > & geometricD() const
Return the vector of geometric directions in mesh.
Definition: polyMesh.C:989
const vectorField & faceCentres() const
const vectorField & faceAreas() const
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:334
#define WarningInFunction
Report a warning using Foam::Warning.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
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
RectangularMatrix< scalar > scalarRectangularMatrix
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:257
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
dimensioned< scalar > mag(const dimensioned< Type > &)
dimensionSet cmptMag(const dimensionSet &)
Definition: dimensionSet.C:296
void cmptMax(FieldField< Field, typename FieldField< Field, Type >::cmptType > &cf, const FieldField< Field, Type > &f)
error FatalError
static const char nl
Definition: Ostream.H:266
labelList f(nPoints)
volScalarField & p
Foam::surfaceFields.