cellQuality.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 -------------------------------------------------------------------------------
10 License
11  This file is part of OpenFOAM.
12 
13  OpenFOAM is free software: you can redistribute it and/or modify it
14  under the terms of the GNU General Public License as published by
15  the Free Software Foundation, either version 3 of the License, or
16  (at your option) any later version.
17 
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 
23  You should have received a copy of the GNU General Public License
24  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
25 
26 \*---------------------------------------------------------------------------*/
27 
28 #include "cellQuality.H"
29 #include "unitConversion.H"
30 #include "SubField.H"
31 
32 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
33 
34 Foam::cellQuality::cellQuality(const polyMesh& mesh)
35 :
36  mesh_(mesh)
37 {}
38 
39 
40 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
41 
42 Foam::tmp<Foam::scalarField> Foam::cellQuality::nonOrthogonality() const
43 {
44  tmp<scalarField> tresult
45  (
46  new scalarField
47  (
48  mesh_.nCells(), 0.0
49  )
50  );
51 
52  scalarField& result = tresult.ref();
53 
54  scalarField sumArea(mesh_.nCells(), Zero);
55 
56  const vectorField& centres = mesh_.cellCentres();
57  const vectorField& areas = mesh_.faceAreas();
58 
59  const labelList& own = mesh_.faceOwner();
60  const labelList& nei = mesh_.faceNeighbour();
61 
62  forAll(nei, facei)
63  {
64  vector d = centres[nei[facei]] - centres[own[facei]];
65  vector s = areas[facei];
66  scalar magS = mag(s);
67 
68  scalar cosDDotS =
69  radToDeg(Foam::acos(min(1.0, (d & s)/(mag(d)*magS + VSMALL))));
70 
71  result[own[facei]] = max(cosDDotS, result[own[facei]]);
72 
73  result[nei[facei]] = max(cosDDotS, result[nei[facei]]);
74  }
75 
76  forAll(mesh_.boundaryMesh(), patchi)
77  {
78  const labelUList& faceCells =
79  mesh_.boundaryMesh()[patchi].faceCells();
80 
81  const vectorField::subField faceCentres =
82  mesh_.boundaryMesh()[patchi].faceCentres();
83 
84  const vectorField::subField faceAreas =
85  mesh_.boundaryMesh()[patchi].faceAreas();
86 
87  forAll(faceCentres, facei)
88  {
89  vector d = faceCentres[facei] - centres[faceCells[facei]];
90  vector s = faceAreas[facei];
91  scalar magS = mag(s);
92 
93  scalar cosDDotS =
94  radToDeg(Foam::acos(min(1.0, (d & s)/(mag(d)*magS + VSMALL))));
95 
96  result[faceCells[facei]] = max(cosDDotS, result[faceCells[facei]]);
97  }
98  }
99 
100  return tresult;
101 }
102 
103 
104 Foam::tmp<Foam::scalarField> Foam::cellQuality::skewness() const
105 {
106  tmp<scalarField> tresult
107  (
108  new scalarField
109  (
110  mesh_.nCells(), 0.0
111  )
112  );
113  scalarField& result = tresult.ref();
114 
115  scalarField sumArea(mesh_.nCells(), Zero);
116 
117  const vectorField& cellCtrs = mesh_.cellCentres();
118  const vectorField& faceCtrs = mesh_.faceCentres();
119  const vectorField& areas = mesh_.faceAreas();
120 
121  const labelList& own = mesh_.faceOwner();
122  const labelList& nei = mesh_.faceNeighbour();
123 
124  forAll(nei, facei)
125  {
126  scalar dOwn = mag
127  (
128  (faceCtrs[facei] - cellCtrs[own[facei]]) & areas[facei]
129  )/mag(areas[facei]);
130 
131  scalar dNei = mag
132  (
133  (cellCtrs[nei[facei]] - faceCtrs[facei]) & areas[facei]
134  )/mag(areas[facei]);
135 
136  point faceIntersection =
137  cellCtrs[own[facei]]
138  + (dOwn/(dOwn+dNei))*(cellCtrs[nei[facei]] - cellCtrs[own[facei]]);
139 
140  scalar skewness =
141  mag(faceCtrs[facei] - faceIntersection)
142  /(mag(cellCtrs[nei[facei]] - cellCtrs[own[facei]]) + VSMALL);
143 
144  result[own[facei]] = max(skewness, result[own[facei]]);
145 
146  result[nei[facei]] = max(skewness, result[nei[facei]]);
147  }
148 
149  forAll(mesh_.boundaryMesh(), patchi)
150  {
151  const labelUList& faceCells =
152  mesh_.boundaryMesh()[patchi].faceCells();
153 
154  const vectorField::subField faceCentres =
155  mesh_.boundaryMesh()[patchi].faceCentres();
156 
157  const vectorField::subField faceAreas =
158  mesh_.boundaryMesh()[patchi].faceAreas();
159 
160  forAll(faceCentres, facei)
161  {
162  vector n = faceAreas[facei]/mag(faceAreas[facei]);
163 
164  point faceIntersection =
165  cellCtrs[faceCells[facei]]
166  + ((faceCentres[facei] - cellCtrs[faceCells[facei]])&n)*n;
167 
168  scalar skewness =
169  mag(faceCentres[facei] - faceIntersection)
170  /(
171  mag(faceCentres[facei] - cellCtrs[faceCells[facei]])
172  + VSMALL
173  );
174 
175  result[faceCells[facei]] = max(skewness, result[faceCells[facei]]);
176  }
177  }
178 
179  return tresult;
180 }
181 
182 
183 Foam::tmp<Foam::scalarField> Foam::cellQuality::faceNonOrthogonality() const
184 {
185  tmp<scalarField> tresult
186  (
187  new scalarField
188  (
189  mesh_.nFaces(), 0.0
190  )
191  );
192  scalarField& result = tresult.ref();
193 
194 
195  const vectorField& centres = mesh_.cellCentres();
196  const vectorField& areas = mesh_.faceAreas();
197 
198  const labelList& own = mesh_.faceOwner();
199  const labelList& nei = mesh_.faceNeighbour();
200 
201  forAll(nei, facei)
202  {
203  vector d = centres[nei[facei]] - centres[own[facei]];
204  vector s = areas[facei];
205  scalar magS = mag(s);
206 
207  scalar cosDDotS =
208  radToDeg(Foam::acos(min(1.0, (d & s)/(mag(d)*magS + VSMALL))));
209 
210  result[facei] = cosDDotS;
211  }
212 
213  label globalFacei = mesh_.nInternalFaces();
214 
215  forAll(mesh_.boundaryMesh(), patchi)
216  {
217  const labelUList& faceCells =
218  mesh_.boundaryMesh()[patchi].faceCells();
219 
220  const vectorField::subField faceCentres =
221  mesh_.boundaryMesh()[patchi].faceCentres();
222 
223  const vectorField::subField faceAreas =
224  mesh_.boundaryMesh()[patchi].faceAreas();
225 
226  forAll(faceCentres, facei)
227  {
228  vector d = faceCentres[facei] - centres[faceCells[facei]];
229  vector s = faceAreas[facei];
230  scalar magS = mag(s);
231 
232  scalar cosDDotS =
233  radToDeg(Foam::acos(min(1.0, (d & s)/(mag(d)*magS + VSMALL))));
234 
235  result[globalFacei++] = cosDDotS;
236  }
237  }
238 
239  return tresult;
240 }
241 
242 
243 Foam::tmp<Foam::scalarField> Foam::cellQuality::faceSkewness() const
244 {
245  tmp<scalarField> tresult
246  (
247  new scalarField
248  (
249  mesh_.nFaces(), 0.0
250  )
251  );
252  scalarField& result = tresult.ref();
253 
254 
255  const vectorField& cellCtrs = mesh_.cellCentres();
256  const vectorField& faceCtrs = mesh_.faceCentres();
257  const vectorField& areas = mesh_.faceAreas();
258 
259  const labelList& own = mesh_.faceOwner();
260  const labelList& nei = mesh_.faceNeighbour();
261 
262  forAll(nei, facei)
263  {
264  scalar dOwn = mag
265  (
266  (faceCtrs[facei] - cellCtrs[own[facei]]) & areas[facei]
267  )/mag(areas[facei]);
268 
269  scalar dNei = mag
270  (
271  (cellCtrs[nei[facei]] - faceCtrs[facei]) & areas[facei]
272  )/mag(areas[facei]);
273 
274  point faceIntersection =
275  cellCtrs[own[facei]]
276  + (dOwn/(dOwn+dNei))*(cellCtrs[nei[facei]] - cellCtrs[own[facei]]);
277 
278  result[facei] =
279  mag(faceCtrs[facei] - faceIntersection)
280  /(mag(cellCtrs[nei[facei]] - cellCtrs[own[facei]]) + VSMALL);
281  }
282 
283 
284  label globalFacei = mesh_.nInternalFaces();
285 
286  forAll(mesh_.boundaryMesh(), patchi)
287  {
288  const labelUList& faceCells =
289  mesh_.boundaryMesh()[patchi].faceCells();
290 
291  const vectorField::subField faceCentres =
292  mesh_.boundaryMesh()[patchi].faceCentres();
293 
294  const vectorField::subField faceAreas =
295  mesh_.boundaryMesh()[patchi].faceAreas();
296 
297  forAll(faceCentres, facei)
298  {
299  vector n = faceAreas[facei]/mag(faceAreas[facei]);
300 
301  point faceIntersection =
302  cellCtrs[faceCells[facei]]
303  + ((faceCentres[facei] - cellCtrs[faceCells[facei]])&n)*n;
304 
305  result[globalFacei++] =
306  mag(faceCentres[facei] - faceIntersection)
307  /(
308  mag(faceCentres[facei] - cellCtrs[faceCells[facei]])
309  + VSMALL
310  );
311  }
312  }
313 
314  return tresult;
315 }
316 
317 
318 // ************************************************************************* //
Foam::cellQuality::faceNonOrthogonality
tmp< scalarField > faceNonOrthogonality() const
Return face non-orthogonality.
Definition: cellQuality.C:183
SubField.H
s
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
Definition: gmvOutputSpray.H:25
Foam::tmp
A class for managing temporary objects.
Definition: PtrList.H:61
Foam::Zero
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
Foam::cellQuality::nonOrthogonality
tmp< scalarField > nonOrthogonality() const
Return cell non-orthogonality.
Definition: cellQuality.C:42
unitConversion.H
Unit conversion functions.
Foam::cellQuality::skewness
tmp< scalarField > skewness() const
Return cell skewness.
Definition: cellQuality.C:104
Foam::min
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33
Foam::polyMesh
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:77
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:296
Foam::SubField
SubField is a Field obtained as a section of another Field, without its own allocation....
Definition: Field.H:64
n
label n
Definition: TABSMDCalcMethod2.H:31
Foam::tmp::ref
T & ref() const
Definition: tmpI.H:227
Foam::Field< scalar >
Foam::radToDeg
constexpr scalar radToDeg(const scalar rad) noexcept
Conversion from radians to degrees.
Definition: unitConversion.H:54
Foam::max
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:47
mesh
dynamicFvMesh & mesh
Definition: createDynamicFvMesh.H:6
Foam::cellQuality::faceSkewness
tmp< scalarField > faceSkewness() const
Return face skewness.
Definition: cellQuality.C:243
cellQuality.H
Foam::Vector< scalar >
Foam::List< label >
Foam::acos
dimensionedScalar acos(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:268
Foam::mag
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
Foam::UList< label >
Foam::faceCells
Smooth ATC in cells next to a set of patches supplied by type.
Definition: faceCells.H:56