CentredFitSnGradData.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) 2013 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 "CentredFitSnGradData.H"
27 #include "surfaceFields.H"
28 #include "SVD.H"
30 
31 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
32 
33 template<class Polynomial>
35 (
36  const fvMesh& mesh,
37  const extendedCentredCellToFaceStencil& stencil,
38  const scalar linearLimitFactor,
39  const scalar centralWeight
40 )
41 :
42  FitData
43  <
44  CentredFitSnGradData<Polynomial>,
45  extendedCentredCellToFaceStencil,
47  >
48  (
49  mesh, stencil, true, linearLimitFactor, centralWeight
50  ),
51  coeffs_(mesh.nFaces())
52 {
53  if (debug)
54  {
55  Info<< "Contructing CentredFitSnGradData<Polynomial>" << endl;
56  }
57 
58  calcFit();
59 
60  if (debug)
61  {
62  Info<< "CentredFitSnGradData<Polynomial>::CentredFitSnGradData() :"
63  << "Finished constructing polynomialFit data"
64  << endl;
65  }
66 }
67 
68 
69 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
70 
71 template<class Polynomial>
73 (
74  scalarList& coeffsi,
75  const List<point>& C,
76  const scalar wLin,
77  const scalar deltaCoeff,
78  const label facei
79 )
80 {
81  vector idir(1,0,0);
82  vector jdir(0,1,0);
83  vector kdir(0,0,1);
84  this->findFaceDirs(idir, jdir, kdir, facei);
85 
86  // Setup the point weights
87  scalarList wts(C.size(), scalar(1));
88  wts[0] = this->centralWeight();
89  wts[1] = this->centralWeight();
90 
91  // Reference point
92  point p0 = this->mesh().faceCentres()[facei];
93 
94  // p0 -> p vector in the face-local coordinate system
95  vector d;
96 
97  // Local coordinate scaling
98  scalar scale = 1;
99 
100  // Matrix of the polynomial components
101  scalarRectangularMatrix B(C.size(), this->minSize(), scalar(0));
102 
103  forAll(C, ip)
104  {
105  const point& p = C[ip];
106  const vector p0p = p - p0;
107 
108  d.x() = p0p & idir;
109  d.y() = p0p & jdir;
110  d.z() = p0p & kdir;
111 
112  if (ip == 0)
113  {
114  scale = cmptMax(cmptMag((d)));
115  }
116 
117  // Scale the radius vector
118  d /= scale;
119 
120  Polynomial::addCoeffs(B[ip], d, wts[ip], this->dim());
121  }
122 
123  // Additional weighting for constant and linear terms
124  for (label i = 0; i < B.n(); i++)
125  {
126  B[i][0] *= wts[0];
127  B[i][1] *= wts[0];
128  }
129 
130  // Set the fit
131  label stencilSize = C.size();
132  coeffsi.setSize(stencilSize);
133 
134  bool goodFit = false;
135  for (int iIt = 0; iIt < 8 && !goodFit; iIt++)
136  {
137  SVD svd(B, SMALL);
138 
139  for (label i=0; i<stencilSize; i++)
140  {
141  coeffsi[i] = wts[1]*wts[i]*svd.VSinvUt()[1][i]/scale;
142  }
143 
144  goodFit =
145  (
146  mag(wts[0]*wts[0]*svd.VSinvUt()[0][0] - wLin)
147  < this->linearLimitFactor()*wLin)
148  && (mag(wts[0]*wts[1]*svd.VSinvUt()[0][1] - (1 - wLin)
149  ) < this->linearLimitFactor()*(1 - wLin))
150  && coeffsi[0] < 0 && coeffsi[1] > 0
151  && mag(coeffsi[0] + deltaCoeff) < 0.5*deltaCoeff
152  && mag(coeffsi[1] - deltaCoeff) < 0.5*deltaCoeff;
153 
154  if (!goodFit)
155  {
156  // (not good fit so increase weight in the centre and weight
157  // for constant and linear terms)
158 
159  WarningIn
160  (
161  "CentredFitSnGradData<Polynomial>::calcFit"
162  "(const List<point>& C, const label facei"
163  ) << "Cannot fit face " << facei << " iteration " << iIt
164  << " with sum of weights " << sum(coeffsi) << nl
165  << " Weights " << coeffsi << nl
166  << " Linear weights " << wLin << " " << 1 - wLin << nl
167  << " deltaCoeff " << deltaCoeff << nl
168  << " sing vals " << svd.S() << nl
169  << "Components of goodFit:\n"
170  << " wts[0]*wts[0]*svd.VSinvUt()[0][0] = "
171  << wts[0]*wts[0]*svd.VSinvUt()[0][0] << nl
172  << " wts[0]*wts[1]*svd.VSinvUt()[0][1] = "
173  << wts[0]*wts[1]*svd.VSinvUt()[0][1]
174  << " dim = " << this->dim() << endl;
175 
176  wts[0] *= 10;
177  wts[1] *= 10;
178 
179  for (label j = 0; j < B.m(); j++)
180  {
181  B[0][j] *= 10;
182  B[1][j] *= 10;
183  }
184 
185  for (label i = 0; i < B.n(); i++)
186  {
187  B[i][0] *= 10;
188  B[i][1] *= 10;
189  }
190  }
191  }
192 
193  if (goodFit)
194  {
195  // Remove the uncorrected coefficients
196  coeffsi[0] += deltaCoeff;
197  coeffsi[1] -= deltaCoeff;
198  }
199  else
200  {
201  WarningIn
202  (
203  "CentredFitSnGradData<Polynomial>::calcFit(..)"
204  ) << "Could not fit face " << facei
205  << " Coefficients = " << coeffsi
206  << ", reverting to uncorrected." << endl;
207 
208  coeffsi = 0;
209  }
210 }
211 
212 
213 template<class Polynomial>
215 {
216  const fvMesh& mesh = this->mesh();
217 
218  // Get the cell/face centres in stencil order.
219  // Centred face stencils no good for triangles or tets.
220  // Need bigger stencils
221  List<List<point> > stencilPoints(mesh.nFaces());
222  this->stencil().collectData(mesh.C(), stencilPoints);
223 
224  // find the fit coefficients for every face in the mesh
225 
226  const surfaceScalarField& w = mesh.surfaceInterpolation::weights();
227  const surfaceScalarField& dC = mesh.nonOrthDeltaCoeffs();
228 
229  for (label facei = 0; facei < mesh.nInternalFaces(); facei++)
230  {
231  calcFit
232  (
233  coeffs_[facei],
234  stencilPoints[facei],
235  w[facei],
236  dC[facei],
237  facei
238  );
239  }
240 
241  const surfaceScalarField::GeometricBoundaryField& bw = w.boundaryField();
242  const surfaceScalarField::GeometricBoundaryField& bdC = dC.boundaryField();
243 
244  forAll(bw, patchi)
245  {
246  const fvsPatchScalarField& pw = bw[patchi];
247  const fvsPatchScalarField& pdC = bdC[patchi];
248 
249  if (pw.coupled())
250  {
251  label facei = pw.patch().start();
252 
253  forAll(pw, i)
254  {
255  calcFit
256  (
257  coeffs_[facei],
258  stencilPoints[facei],
259  pw[i],
260  pdC[i],
261  facei
262  );
263  facei++;
264  }
265  }
266  }
267 }
268 
269 
270 // ************************************************************************* //
Data for the upwinded and centred polynomial fit interpolation schemes. The linearCorrection_ determi...
Definition: FitData.H:54
virtual bool coupled() const
Return true if this patch field is coupled.
Polynomial templated on size (order):
Definition: Polynomial.H:63
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:78
Foam::surfaceFields.
label nFaces() const
dimensioned< scalar > mag(const dimensioned< Type > &)
GeometricBoundaryField & boundaryField()
Return reference to GeometricBoundaryField.
const fvPatch & patch() const
Return patch.
void cmptMax(FieldField< Field, typename FieldField< Field, Type >::cmptType > &cf, const FieldField< Field, Type > &f)
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
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:76
messageStream Info
dynamicFvMesh & mesh
static const char nl
Definition: Ostream.H:260
void setSize(const label)
Reset size of List.
Definition: List.C:318
const Cmpt & y() const
Definition: VectorI.H:71
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
#define WarningIn(functionName)
Report a warning using Foam::Warning.
const volVectorField & C() const
Return cell centres as volVectorField.
volScalarField & p
Definition: createFields.H:51
Singular value decomposition of a rectangular matrix.
Definition: SVD.H:52
#define forAll(list, i)
Definition: UList.H:421
label patchi
const Cmpt & x() const
Definition: VectorI.H:65
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
const Cmpt & z() const
Definition: VectorI.H:77
const surfaceScalarField & nonOrthDeltaCoeffs() const
Return reference to non-orthogonal cell-centre difference.
label nInternalFaces() const
const scalarDiagonalMatrix & S() const
Return the singular values.
Definition: SVDI.H:47
CentredFitSnGradData(const fvMesh &mesh, const extendedCentredCellToFaceStencil &stencil, const scalar linearLimitFactor, const scalar centralWeight)
Construct from components.
const scalarRectangularMatrix & VSinvUt() const
Return VSinvUt (the pseudo inverse)
Definition: SVDI.H:52
void cmptMag(FieldField< Field, Type > &cf, const FieldField< Field, Type > &f)
label start() const
Return start label of this patch in the polyMesh face list.
Definition: fvPatch.H:155
An abstract base class with a fat-interface to all derived classes covering all possible ways in whic...
Definition: fvsPatchField.H:65