FreeFOAM The Cross-Platform CFD Toolkit
extendedLeastSquaresVectors.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 
28 #include <finiteVolume/volFields.H>
29 
30 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
31 
32 namespace Foam
33 {
34  defineTypeNameAndDebug(extendedLeastSquaresVectors, 0);
35 }
36 
37 
38 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
39 
41 (
42  const fvMesh& mesh,
43  const scalar minDet
44 )
45 :
47  minDet_(minDet),
48  pVectorsPtr_(NULL),
49  nVectorsPtr_(NULL),
50  additionalCellsPtr_(NULL),
51  additionalVectorsPtr_(NULL)
52 {}
53 
54 
55 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
56 
58 {
59  deleteDemandDrivenData(pVectorsPtr_);
60  deleteDemandDrivenData(nVectorsPtr_);
61 
62  deleteDemandDrivenData(additionalCellsPtr_);
63  deleteDemandDrivenData(additionalVectorsPtr_);
64 }
65 
66 
67 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
68 
69 void Foam::extendedLeastSquaresVectors::makeLeastSquaresVectors() const
70 {
71  if (debug)
72  {
73  Info<< "extendedLeastSquaresVectors::makeLeastSquaresVectors() :"
74  << "Constructing least square gradient vectors"
75  << endl;
76  }
77 
78  const fvMesh& mesh = mesh_;
79 
80  pVectorsPtr_ = new surfaceVectorField
81  (
82  IOobject
83  (
84  "LeastSquaresP",
85  mesh_.pointsInstance(),
86  mesh_,
89  false
90  ),
91  mesh_,
93  );
94  surfaceVectorField& lsP = *pVectorsPtr_;
95 
96  nVectorsPtr_ = new surfaceVectorField
97  (
98  IOobject
99  (
100  "LeastSquaresN",
101  mesh_.pointsInstance(),
102  mesh_,
105  false
106  ),
107  mesh_,
109  );
110  surfaceVectorField& lsN = *nVectorsPtr_;
111 
112  // Set local references to mesh data
113  const unallocLabelList& owner = mesh_.owner();
114  const unallocLabelList& neighbour = mesh_.neighbour();
115 
116  // Build the d-vectors
117  surfaceVectorField d = mesh_.Sf()/(mesh_.magSf()*mesh_.deltaCoeffs());
118 
119  if (!mesh_.orthogonal())
120  {
121  d -= mesh_.correctionVectors()/mesh_.deltaCoeffs();
122  }
123 
124 
125  // Set up temporary storage for the dd tensor (before inversion)
126  symmTensorField dd(mesh_.nCells(), symmTensor::zero);
127 
128  forAll(owner, faceI)
129  {
130  symmTensor wdd = 1.0/magSqr(d[faceI])*sqr(d[faceI]);
131 
132  dd[owner[faceI]] += wdd;
133  dd[neighbour[faceI]] += wdd;
134  }
135 
136  // Visit the boundaries. Coupled boundaries are taken into account
137  // in the construction of d vectors.
138  forAll(d.boundaryField(), patchI)
139  {
140  const fvsPatchVectorField& pd = d.boundaryField()[patchI];
141 
142  const fvPatch& p = pd.patch();
143  const unallocLabelList& faceCells = p.faceCells();
144 
145  forAll(pd, patchFaceI)
146  {
147  dd[faceCells[patchFaceI]] +=
148  (1.0/magSqr(pd[patchFaceI]))*sqr(pd[patchFaceI]);
149  }
150  }
151 
152  scalarField detdd = det(dd);
153 
154  Info<< "max(detdd) = " << max(detdd) << endl;
155  Info<< "min(detdd) = " << min(detdd) << endl;
156  Info<< "average(detdd) = " << average(detdd) << endl;
157 
158  label nAddCells = 0;
159  label maxNaddCells = 4*detdd.size();
160  additionalCellsPtr_ = new List<labelPair>(maxNaddCells);
161  List<labelPair>& additionalCells_ = *additionalCellsPtr_;
162 
163  forAll (detdd, i)
164  {
165  label count = 0;
166 
167  while (++count < 100 && detdd[i] < minDet_)
168  {
169  if (nAddCells == maxNaddCells)
170  {
172  (
173  "extendedLeastSquaresVectors::"
174  "makeLeastSquaresVectors() const"
175  ) << "nAddCells exceeds maxNaddCells"
176  << exit(FatalError);
177  }
178 
179  labelList pointLabels = mesh.cells()[i].labels(mesh.faces());
180 
181  scalar maxDetddij = 0.0;
182 
183  label addCell = -1;
184 
185  forAll(pointLabels, j)
186  {
187  forAll(mesh.pointCells()[pointLabels[j]], k)
188  {
189  label cellj = mesh.pointCells()[pointLabels[j]][k];
190 
191  if (cellj != i)
192  {
193  if (cellj != -1)
194  {
195  vector dCij = (mesh.C()[cellj] - mesh.C()[i]);
196 
197  symmTensor ddij =
198  dd[i] + (1.0/magSqr(dCij))*sqr(dCij);
199 
200  scalar detddij = det(ddij);
201 
202  if (detddij > maxDetddij)
203  {
204  addCell = cellj;
205  maxDetddij = detddij;
206  }
207  }
208  }
209  }
210  }
211 
212  if (addCell != -1)
213  {
214  additionalCells_[nAddCells][0] = i;
215  additionalCells_[nAddCells++][1] = addCell;
216  vector dCij = mesh.C()[addCell] - mesh.C()[i];
217  dd[i] += (1.0/magSqr(dCij))*sqr(dCij);
218  detdd[i] = det(dd[i]);
219  }
220  }
221  }
222 
223  additionalCells_.setSize(nAddCells);
224 
225  Info<< "max(detdd) = " << max(detdd) << endl;
226  Info<< "min(detdd) = " << min(detdd) << endl;
227  Info<< "average(detdd) = " << average(detdd) << endl;
228  Info<< "nAddCells/nCells = " << scalar(nAddCells)/mesh.nCells() << endl;
229 
230  // Invert the dd tensor
231  symmTensorField invDd = inv(dd);
232 
233 
234  // Revisit all faces and calculate the lsP and lsN vectors
235  forAll(owner, faceI)
236  {
237  lsP[faceI] =
238  (1.0/magSqr(d[faceI]))*(invDd[owner[faceI]] & d[faceI]);
239 
240  lsN[faceI] =
241  ((-1.0)/magSqr(d[faceI]))*(invDd[neighbour[faceI]] & d[faceI]);
242  }
243 
244  forAll(lsP.boundaryField(), patchI)
245  {
246  const fvsPatchVectorField& pd = d.boundaryField()[patchI];
247 
248  fvsPatchVectorField& patchLsP = lsP.boundaryField()[patchI];
249 
250  const fvPatch& p = patchLsP.patch();
251  const unallocLabelList& faceCells = p.faceCells();
252 
253  forAll(p, patchFaceI)
254  {
255  patchLsP[patchFaceI] =
256  (1.0/magSqr(pd[patchFaceI]))
257  *(invDd[faceCells[patchFaceI]] & pd[patchFaceI]);
258  }
259  }
260 
261 
262  additionalVectorsPtr_ = new vectorField(additionalCells_.size());
263  vectorField& additionalVectors_ = *additionalVectorsPtr_;
264 
265  forAll(additionalCells_, i)
266  {
267  vector dCij =
268  mesh.C()[additionalCells_[i][1]] - mesh.C()[additionalCells_[i][0]];
269 
270  additionalVectors_[i] =
271  (1.0/magSqr(dCij))*(invDd[additionalCells_[i][0]] & dCij);
272  }
273 
274  if (debug)
275  {
276  Info<< "extendedLeastSquaresVectors::makeLeastSquaresVectors() :"
277  << "Finished constructing least square gradient vectors"
278  << endl;
279  }
280 }
281 
282 
285 {
286  if (!pVectorsPtr_)
287  {
288  makeLeastSquaresVectors();
289  }
290 
291  return *pVectorsPtr_;
292 }
293 
294 
297 {
298  if (!nVectorsPtr_)
299  {
300  makeLeastSquaresVectors();
301  }
302 
303  return *nVectorsPtr_;
304 }
305 
306 
309 {
310  if (!additionalCellsPtr_)
311  {
312  makeLeastSquaresVectors();
313  }
314 
315  return *additionalCellsPtr_;
316 }
317 
318 
319 const Foam::vectorField&
321 {
322  if (!additionalVectorsPtr_)
323  {
324  makeLeastSquaresVectors();
325  }
326 
327  return *additionalVectorsPtr_;
328 }
329 
330 
332 {
333  deleteDemandDrivenData(pVectorsPtr_);
334  deleteDemandDrivenData(nVectorsPtr_);
335 
336  deleteDemandDrivenData(additionalCellsPtr_);
337  deleteDemandDrivenData(additionalVectorsPtr_);
338 
339  return true;
340 }
341 
342 
343 // ************************ vim: set sw=4 sts=4 et: ************************ //