FreeFOAM The Cross-Platform CFD Toolkit
leastSquaresVectors.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 "leastSquaresVectors.H"
28 #include <finiteVolume/volFields.H>
29 
30 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
31 
32 namespace Foam
33 {
34  defineTypeNameAndDebug(leastSquaresVectors, 0);
35 }
36 
37 
38 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
39 
41 :
43  pVectorsPtr_(NULL),
44  nVectorsPtr_(NULL)
45 {}
46 
47 
48 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
49 
51 {
52  deleteDemandDrivenData(pVectorsPtr_);
53  deleteDemandDrivenData(nVectorsPtr_);
54 }
55 
56 
57 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
58 
59 void Foam::leastSquaresVectors::makeLeastSquaresVectors() const
60 {
61  if (debug)
62  {
63  Info<< "leastSquaresVectors::makeLeastSquaresVectors() :"
64  << "Constructing least square gradient vectors"
65  << endl;
66  }
67 
68  const fvMesh& mesh = mesh_;
69 
70  pVectorsPtr_ = new surfaceVectorField
71  (
72  IOobject
73  (
74  "LeastSquaresP",
75  mesh_.pointsInstance(),
76  mesh_,
79  false
80  ),
81  mesh_,
83  );
84  surfaceVectorField& lsP = *pVectorsPtr_;
85 
86  nVectorsPtr_ = new surfaceVectorField
87  (
88  IOobject
89  (
90  "LeastSquaresN",
91  mesh_.pointsInstance(),
92  mesh_,
95  false
96  ),
97  mesh_,
99  );
100  surfaceVectorField& lsN = *nVectorsPtr_;
101 
102  // Set local references to mesh data
103  const unallocLabelList& owner = mesh_.owner();
104  const unallocLabelList& neighbour = mesh_.neighbour();
105 
106  const volVectorField& C = mesh.C();
107  const surfaceScalarField& w = mesh.weights();
108  const surfaceScalarField& magSf = mesh.magSf();
109 
110 
111  // Set up temporary storage for the dd tensor (before inversion)
112  symmTensorField dd(mesh_.nCells(), symmTensor::zero);
113 
114  forAll(owner, facei)
115  {
116  label own = owner[facei];
117  label nei = neighbour[facei];
118 
119  vector d = C[nei] - C[own];
120  symmTensor wdd = (magSf[facei]/magSqr(d))*sqr(d);
121 
122  dd[own] += (1 - w[facei])*wdd;
123  dd[nei] += w[facei]*wdd;
124  }
125 
126 
127  forAll(lsP.boundaryField(), patchi)
128  {
129  const fvsPatchScalarField& pw = w.boundaryField()[patchi];
130  const fvsPatchScalarField& pMagSf = magSf.boundaryField()[patchi];
131 
132  const fvPatch& p = pw.patch();
133  const unallocLabelList& faceCells = p.patch().faceCells();
134 
135  // Build the d-vectors
136  vectorField pd =
137  mesh.Sf().boundaryField()[patchi]
138  /(
139  mesh.magSf().boundaryField()[patchi]
140  *mesh.deltaCoeffs().boundaryField()[patchi]
141  );
142 
143  if (!mesh.orthogonal())
144  {
145  pd -= mesh.correctionVectors().boundaryField()[patchi]
146  /mesh.deltaCoeffs().boundaryField()[patchi];
147  }
148 
149  if (p.coupled())
150  {
151  forAll(pd, patchFacei)
152  {
153  const vector& d = pd[patchFacei];
154 
155  dd[faceCells[patchFacei]] +=
156  ((1 - pw[patchFacei])*pMagSf[patchFacei]/magSqr(d))*sqr(d);
157  }
158  }
159  else
160  {
161  forAll(pd, patchFacei)
162  {
163  const vector& d = pd[patchFacei];
164 
165  dd[faceCells[patchFacei]] +=
166  (pMagSf[patchFacei]/magSqr(d))*sqr(d);
167  }
168  }
169  }
170 
171 
172  // Invert the dd tensor
173  symmTensorField invDd = inv(dd);
174 
175 
176  // Revisit all faces and calculate the lsP and lsN vectors
177  forAll(owner, facei)
178  {
179  label own = owner[facei];
180  label nei = neighbour[facei];
181 
182  vector d = C[nei] - C[own];
183  scalar magSfByMagSqrd = magSf[facei]/magSqr(d);
184 
185  lsP[facei] = (1 - w[facei])*magSfByMagSqrd*(invDd[own] & d);
186  lsN[facei] = -w[facei]*magSfByMagSqrd*(invDd[nei] & d);
187  }
188 
189  forAll(lsP.boundaryField(), patchi)
190  {
191  fvsPatchVectorField& patchLsP = lsP.boundaryField()[patchi];
192 
193  const fvsPatchScalarField& pw = w.boundaryField()[patchi];
194  const fvsPatchScalarField& pMagSf = magSf.boundaryField()[patchi];
195 
196  const fvPatch& p = pw.patch();
197  const unallocLabelList& faceCells = p.faceCells();
198 
199  // Build the d-vectors
200  vectorField pd =
201  mesh.Sf().boundaryField()[patchi]
202  /(
203  mesh.magSf().boundaryField()[patchi]
204  *mesh.deltaCoeffs().boundaryField()[patchi]
205  );
206 
207  if (!mesh.orthogonal())
208  {
209  pd -= mesh.correctionVectors().boundaryField()[patchi]
210  /mesh.deltaCoeffs().boundaryField()[patchi];
211  }
212 
213 
214  if (p.coupled())
215  {
216  forAll(pd, patchFacei)
217  {
218  const vector& d = pd[patchFacei];
219 
220  patchLsP[patchFacei] =
221  ((1 - pw[patchFacei])*pMagSf[patchFacei]/magSqr(d))
222  *(invDd[faceCells[patchFacei]] & d);
223  }
224  }
225  else
226  {
227  forAll(pd, patchFacei)
228  {
229  const vector& d = pd[patchFacei];
230 
231  patchLsP[patchFacei] =
232  pMagSf[patchFacei]*(1.0/magSqr(d))
233  *(invDd[faceCells[patchFacei]] & d);
234  }
235  }
236  }
237 
238 
239  // For 3D meshes check the determinant of the dd tensor and switch to
240  // Gauss if it is less than 3
241  /* Currently the det(dd[celli]) criterion is incorrect: dd is weighted by Sf
242  if (mesh.nGeometricD() == 3)
243  {
244  label nBadCells = 0;
245 
246  const cellList& cells = mesh.cells();
247  const scalarField& V = mesh.V();
248  const surfaceVectorField& Sf = mesh.Sf();
249  const surfaceScalarField& w = mesh.weights();
250 
251  forAll (dd, celli)
252  {
253  if (det(dd[celli]) < 3)
254  {
255  nBadCells++;
256 
257  const cell& c = cells[celli];
258 
259  forAll(c, cellFacei)
260  {
261  label facei = c[cellFacei];
262 
263  if (mesh.isInternalFace(facei))
264  {
265  scalar wf = max(min(w[facei], 0.8), 0.2);
266 
267  if (celli == owner[facei])
268  {
269  lsP[facei] = (1 - wf)*Sf[facei]/V[celli];
270  }
271  else
272  {
273  lsN[facei] = -wf*Sf[facei]/V[celli];
274  }
275  }
276  else
277  {
278  label patchi = mesh.boundaryMesh().whichPatch(facei);
279 
280  if (mesh.boundary()[patchi].size())
281  {
282  label patchFacei =
283  facei - mesh.boundaryMesh()[patchi].start();
284 
285  if (mesh.boundary()[patchi].coupled())
286  {
287  scalar wf = max
288  (
289  min
290  (
291  w.boundaryField()[patchi][patchFacei],
292  0.8
293  ),
294  0.2
295  );
296 
297  lsP.boundaryField()[patchi][patchFacei] =
298  (1 - wf)
299  *Sf.boundaryField()[patchi][patchFacei]
300  /V[celli];
301  }
302  else
303  {
304  lsP.boundaryField()[patchi][patchFacei] =
305  Sf.boundaryField()[patchi][patchFacei]
306  /V[celli];
307  }
308  }
309  }
310  }
311  }
312  }
313 
314  if (debug)
315  {
316  InfoIn("leastSquaresVectors::makeLeastSquaresVectors()")
317  << "number of bad cells switched to Gauss = " << nBadCells
318  << endl;
319  }
320  }
321  */
322 
323  if (debug)
324  {
325  Info<< "leastSquaresVectors::makeLeastSquaresVectors() :"
326  << "Finished constructing least square gradient vectors"
327  << endl;
328  }
329 }
330 
331 
333 {
334  if (!pVectorsPtr_)
335  {
336  makeLeastSquaresVectors();
337  }
338 
339  return *pVectorsPtr_;
340 }
341 
342 
344 {
345  if (!nVectorsPtr_)
346  {
347  makeLeastSquaresVectors();
348  }
349 
350  return *nVectorsPtr_;
351 }
352 
353 
355 {
356  deleteDemandDrivenData(pVectorsPtr_);
357  deleteDemandDrivenData(nVectorsPtr_);
358 
359  return true;
360 }
361 
362 
363 // ************************ vim: set sw=4 sts=4 et: ************************ //