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