All Classes Namespaces Files Functions Variables Typedefs Friends Macros Groups
extendedFaceStencilCalculateWeights.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 | www.openfoam.com
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8  Copyright (C) 2011-2016 OpenFOAM Foundation
9  Copyright (C) 2019 OpenCFD Ltd.
10  Copyright (C) 2016-2019 ISP RAS (www.ispras.ru) UniCFD Group (www.unicfd.ru)
11 -------------------------------------------------------------------------------
12 License
13  This file is part of QGDsolver library, based on OpenFOAM+.
14  OpenFOAM is free software: you can redistribute it and/or modify it
15  under the terms of the GNU General Public License as published by
16  the Free Software Foundation, either version 3 of the License, or
17  (at your option) any later version.
18  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
19  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
20  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
21  for more details.
22  You should have received a copy of the GNU General Public License
23  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
24 Class
25  Foam::fvsc::leastSquares::extendedFaceStencilCalculateWeights
26 Description
27  Base methods for calculating weights
28 \*---------------------------------------------------------------------------*/
29 
30 
31 
32 #include "leastSquaresBase.H"
33 #include "polyMesh.H"
34 #include "fvMesh.H"
35 #include "word.H"
36 #include "IOstream.H"
37 #include "Ostream.H"
38 #include <HashTable.H>
39 
40 //- Compute weights for least squares scheme for gradient calculation.
41 void Foam::fvsc::leastSquaresBase::calculateWeights()
42 {
43  //Pout << "Start calculateWeights()" << endl;
44 
45  const faceList& faces = cMesh_.faces();
46  GdfAll_.resize(faces.size());
47  wf2All_.resize(faces.size());
48  label cellDim = 3;
49 
50  //create list of underdetermined cells:
51  //cellSet udCells(mesh_, "udCells", mesh_.nCells()/100);
52  //mesh_.checkCellDeterminant(true, &udCells);
53  label nDegFaces = 0;
54  scalar minDet = GREAT;
55  scalar detG = 0.0;
56 
57  //for internal faces
58  forAll(faces, facei)
59  {
60  if (cMesh_.isInternalFace(facei))
61  {
62  List<vector> df(neighbourCells_[facei].size());
63  scalarList wf2(neighbourCells_[facei].size());
64  symmTensor G(0);
65 
66  vector Cf = cMesh_.faceCentres()[facei];
67 
68  forAll(neighbourCells_[facei], i)
69  {
70  df[i] = cMesh_.cellCentres()[neighbourCells_[facei][i]] - Cf;
71  wf2[i] = 1/magSqr(df[i]);
72  symmTensor addToG(0);
73  addToG = sqr(df[i]);
74  addToG = addToG * wf2[i];
75  G += addToG;
76  }
77 
78  symmTensor G0(0);
79  cellDim = 3;
80 
81  //correct G tensor for zero directions
82  {
83  if (mag(G.xx()) < SMALL)
84  {
85  G0.xx() = 1;
86  cellDim--;
87  }
88  if (mag(G.yy()) < SMALL)
89  {
90  G0.yy() = 1;
91  cellDim--;
92  }
93  if (mag(G.zz()) < SMALL)
94  {
95  G0.zz() = 1;
96  cellDim--;
97  }
98 
99  if (cellDim != cMesh_.nGeometricD())
100  {
101  WarningInFunction
102  << "face " << facei << " with center "
103  << cMesh_.faceCentres()[facei] << nl
104  << " connected to cells with dimensions " << cellDim
105  << " less then geometric " << cMesh_.nGeometricD()
106  << nl << endl;
107  Pout << "Degenerate face: " << facei << endl;
108  }
109  }
110 
111  /*
112  if(mesh_.nGeometricD()==1)
113  {
114  symmTensor G01(1, 0, 0, 1, 0, 1);
115  symmTensor G02(sqr((Vector<label>::one + mesh_.geometricD())/2));
116  G0 = G01 - G02;
117  }
118  else
119  {
120  if(mesh_.nGeometricD()==2)
121  {
122  G0 = sqr((Vector<label>::one - mesh_.geometricD())/2);
123  }
124  };
125  */
126 
127  G = G + G0;
128  detG = det(G);
129  if (detG < minDet)
130  {
131  minDet = detG;
132  }
133 
134  if (detG < 1)
135  {
136  nDegFaces++;
137  internalDegFaces_.append(facei);
138  }
139  else
140  {
141  G = inv(G);
142  G = G - G0;
143  }
144 
145  forAll(df, i)
146  {
147  df[i] = G&df[i];
148  }
149 
150  GdfAll_[facei] = df;
151  wf2All_[facei] = wf2;
152  }
153  }
154 
155  if (nDegFaces > 0)
156  {
157  Pout << "Min determinant : " << minDet << endl;
158  Pout << "Total # of deg. faces: " << nDegFaces << endl;
159  }
160 
161  //Info << "End for not parallel" << endl;
162 
163  if (Pstream::parRun())
164  {
165  List3<vector> ownCellCenters(nProcPatches_);
166  List3<vector> neiCellCenters(nProcPatches_);
167  List3<vector> corCellCenters(nProcPatches_);
168  //List<label> nOwnCells(nProcPatches_, 0);
169 
170  label cellId = -1;
171  label nCorCells = -1;
172 
173  //Step 1. Set cell centers at my patches
174  //Cell centers are stored only per processor patch
175  forAll(ownCellCenters, iProcPatch)
176  {
177  ownCellCenters[iProcPatch].resize(myProcPatchCells_[iProcPatch].size());
178  neiCellCenters[iProcPatch].resize(myProcPatchCells_[iProcPatch].size());
179  corCellCenters[iProcPatch].resize(myProcPatchCells_[iProcPatch].size());
180  forAll(ownCellCenters[iProcPatch], iFace)
181  {
182  ownCellCenters[iProcPatch][iFace].resize
183  (
184  myProcPatchCells_[iProcPatch][iFace].size()
185  );
186 
187  nCorCells = corEnd_[iProcPatch][iFace] - corStart_[iProcPatch][iFace] + 1;
188 
189  //Pout << "corEnd = " << corEnd_[iProcPatch][iFace] << " nCorCell = " << nCorCells << endl;
190  corCellCenters[iProcPatch][iFace].resize
191  (
192  nCorCells
193  );
194 
195  forAll(ownCellCenters[iProcPatch][iFace], iCell)
196  {
197  cellId = myProcPatchCells_[iProcPatch][iFace][iCell];
198  ownCellCenters[iProcPatch][iFace][iCell] = cMesh_.C()[cellId];
199  }
200  //nOwnCells[iProcPatch] += ownCellCenters[iProcPatch][iFace].size();
201  }
202  }
203 
204  // Step 2. Loop over all neighboring processors and send/receive cell centers
205  {
206  PstreamBuffers pBuffers(Pstream::commsTypes::nonBlocking);
207 
208  forAll(neigProcs_, iProcPair)
209  {
210  if (procPairs_[iProcPair] > -1) //send cell centers for patch-neighbouring processes
211  {
212  label procId = neigProcs_[iProcPair];
213  UOPstream oProcStr(procId, pBuffers);
214  oProcStr << ownCellCenters[iProcPair];
215  }
216  }
217 
218  pBuffers.finishedSends();
219 
220  forAll(neigProcs_, iProcPair)
221  {
222  if (procPairs_[iProcPair] > -1) //recieve cell centers for patch-neighbouring processes
223  {
224  label procId = neigProcs_[iProcPair];
225  UIPstream iProcStr(procId, pBuffers);
226  iProcStr >> neiCellCenters[iProcPair];
227  }
228  }
229 
230  //Pout << "neiCellCenters = " << neiCellCenters << endl;
231  }
232 
233  // Step 3. Loop over all corner neigbouring processors and send/receive cell centers
234  {
235  PstreamBuffers pBuffers(Pstream::commsTypes::nonBlocking);
236 
237  // Send
238  forAll(neigProcs_, iProcPair)
239  {
240  if (procPairs_[iProcPair] < 0)
241  {
242  label procId = neigProcs_[iProcPair];
243  UOPstream oProcStr(procId, pBuffers);
244  label id = corProcIds_[procId];
245 
246  List<vector> locCc (corCellIds_[id].size());
247 
248  forAll(locCc, iCell)
249  {
250  label cellId = corCellIds_[id][iCell];
251  locCc[iCell] = cMesh_.C()[cellId];
252  }
253  oProcStr << locCc;
254  //Pout << "Sending " << locCc << " to " << procId << endl;
255  }
256  }
257 
258  // Recieve
259  pBuffers.finishedSends();
260  label iCorProc = 0;
261  forAll(neigProcs_, iProcPair)
262  {
263  if (procPairs_[iProcPair] < 0)
264  {
265  label procId = neigProcs_[iProcPair];
266  UIPstream iProcStr(procId, pBuffers);
267 
268  List<vector> corCc (iProcStr);
269 
270  //Pout << "Received from " << procId << " cell centers " << corCc << endl;
271 
272  const List<Triple<label> > & addr = corAddr_[iCorProc];
273  label patchNo = -1;
274  label faceNo = -1;
275  label cellNo = -1;
276 
277  forAll(corCc, iCell)
278  {
279  patchNo = addr[iCell][0];
280  faceNo = addr[iCell][1];
281  cellNo = addr[iCell][2];
282 
283  corCellCenters[patchNo][faceNo][cellNo] = corCc[iCell];
284  }
285 
286  iCorProc++;
287  }
288  }
289 
290  //Pout << "corCellCenters = " << corCellCenters << endl;
291  }
292 
293 
294  // Step 4. Calculate weights
295  //Pout << "neiCellCenters = " << neiCellCenters << endl;
296  //Pout << "corCellCenters = " << corCellCenters << endl;
297  procDegFaces_.resize(nProcPatches_);
298  forAll(myProcPatchCells_,iProcPatch)
299  {
300  if (procPairs_[iProcPatch] > -1)
301  {
302  const label iProcPatchId = procPairs_[iProcPatch];
303  const fvPatch& fvp = cMesh_.boundary()[iProcPatchId];
304 
305  label nFaceCells = 0;
306 
307  forAll(ownCellCenters[iProcPatch], facei)
308  {
309  nFaceCells = corEnd_[iProcPatch][facei] + 1;
310 
311  List<vector> df (nFaceCells, vector::zero);
312  List<scalar> wf2(nFaceCells, 0.0);
313 
314  symmTensor G(0);
315 
316  forAll(ownCellCenters[iProcPatch][facei], i)
317  {
318  df[i] = ownCellCenters[iProcPatch][facei][i] - fvp.Cf()[facei];
319  wf2[i] = 1/magSqr(df[i]);
320  symmTensor addToG(0);
321  addToG = sqr(df[i]);
322  addToG = addToG * wf2[i];
323  G += addToG;
324  }
325 
326  label k = neiStart_[iProcPatch][facei];
327  forAll(neiCellCenters[iProcPatch][facei], i)
328  {
329  df[k] = neiCellCenters[iProcPatch][facei][i] - fvp.Cf()[facei];
330  wf2[k] = 1/magSqr(df[k]);
331  symmTensor addToG(0);
332  addToG = sqr(df[k]);
333  addToG = addToG * wf2[k];
334  G += addToG;
335  k++;
336  }
337 
338  label l = corStart_[iProcPatch][facei];
339  forAll(corCellCenters[iProcPatch][facei], i)
340  {
341  df[l] = corCellCenters[iProcPatch][facei][i] - fvp.Cf()[facei];
342  wf2[l] = 1/magSqr(df[l]);
343  symmTensor addToG(0);
344  addToG = sqr(df[l]);
345  addToG = addToG * wf2[l];
346  G += addToG;
347  l++;
348  }
349 
350  symmTensor G0(0);
351  cellDim = 3;
352  /*
353  if(mesh_.nGeometricD()==1)
354  {
355  symmTensor G01(1, 0, 0, 1, 0, 1);
356  symmTensor G02(sqr((Vector<label>::one + mesh_.geometricD())/2));
357  G0 = G01 - G02;
358  }
359  else
360  {
361  if(mesh_.nGeometricD()==2)
362  {
363  G0 = sqr((Vector<label>::one - mesh_.geometricD())/2);
364  }
365  };
366  */
367 
368  //correct G tensor for zero directions
369  {
370  if (mag(G.xx()) < SMALL)
371  {
372  G0.xx() = 1;
373  cellDim--;
374  }
375  if (mag(G.yy()) < SMALL)
376  {
377  G0.yy() = 1;
378  cellDim--;
379  }
380  if (mag(G.zz()) < SMALL)
381  {
382  G0.zz() = 1;
383  cellDim--;
384  }
385  }
386 
387  if (cellDim != cMesh_.nGeometricD())
388  {
389  WarningInFunction
390  << "face " << facei << " with center "
391  << fvp.Cf()[facei] << nl
392  << " connected to cells with dimensions " << cellDim
393  << " less then geometric " << cMesh_.nGeometricD()
394  << nl << endl;
395  }
396 
397  G = G + G0;
398 
399  detG = det(G);
400  if (detG < 1)
401  {
402  procDegFaces_[iProcPatch].append(facei);
403  }
404  else
405  {
406  G = inv(G);
407  G = G - G0;
408  }
409 
410  forAll(df, i)
411  {
412  df[i] = G&df[i];
413  }
414 
415  procGdf_[iProcPatch][facei] = df;
416  procWf2_[iProcPatch][facei] = wf2;
417  }
418  }
419  }
420  }
421 
422  //Pout << "End calculateWeights()" << endl;
423 };
424 
Foam::fvsc::leastSquaresBase::procDegFaces_
List< DynamicList< label > > procDegFaces_
Definition: leastSquaresBase.H:107
Foam::fvsc::leastSquaresBase::calculateWeights
void calculateWeights()
Compute weights for least squares scheme for gradient calculation.
Definition: extendedFaceStencilCalculateWeights.C:36
Foam::fvsc::leastSquaresBase::corProcIds_
labelHashTable< label > corProcIds_
Definition: leastSquaresBase.H:105
Foam::faceList
List< face > faceList
Definition: GaussVolPointBase3D.H:40
Foam::fvsc::leastSquaresBase::corAddr_
List2< Triple< label > > corAddr_
Definition: leastSquaresBase.H:103
forAll
forAll(Y, i)
Definition: QGDYEqn.H:36
Foam::fvsc::leastSquaresBase::internalDegFaces_
DynamicList< label > internalDegFaces_
Definition: leastSquaresBase.H:76

Generated by   doxygen 1.8.5