cellQuality.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) 2011-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 "cellQuality.H"
27 #include "unitConversion.H"
28 #include "SubField.H"
29 
30 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
31 
32 Foam::cellQuality::cellQuality(const polyMesh& mesh)
33 :
34  mesh_(mesh)
35 {}
36 
37 
38 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
39 
40 Foam::tmp<Foam::scalarField> Foam::cellQuality::nonOrthogonality() const
41 {
42  tmp<scalarField> tresult
43  (
44  new scalarField
45  (
46  mesh_.nCells(), 0.0
47  )
48  );
49 
50  scalarField& result = tresult.ref();
51 
52  scalarField sumArea(mesh_.nCells(), 0.0);
53 
54  const vectorField& centres = mesh_.cellCentres();
55  const vectorField& areas = mesh_.faceAreas();
56 
57  const labelList& own = mesh_.faceOwner();
58  const labelList& nei = mesh_.faceNeighbour();
59 
60  forAll(nei, facei)
61  {
62  vector d = centres[nei[facei]] - centres[own[facei]];
63  vector s = areas[facei];
64  scalar magS = mag(s);
65 
66  scalar cosDDotS =
67  radToDeg(Foam::acos(min(1.0, (d & s)/(mag(d)*magS + vSmall))));
68 
69  result[own[facei]] = max(cosDDotS, result[own[facei]]);
70 
71  result[nei[facei]] = max(cosDDotS, result[nei[facei]]);
72  }
73 
74  forAll(mesh_.boundaryMesh(), patchi)
75  {
76  const labelUList& faceCells =
77  mesh_.boundaryMesh()[patchi].faceCells();
78 
79  const vectorField::subField faceCentres =
80  mesh_.boundaryMesh()[patchi].faceCentres();
81 
82  const vectorField::subField faceAreas =
83  mesh_.boundaryMesh()[patchi].faceAreas();
84 
85  forAll(faceCentres, facei)
86  {
87  vector d = faceCentres[facei] - centres[faceCells[facei]];
88  vector s = faceAreas[facei];
89  scalar magS = mag(s);
90 
91  scalar cosDDotS =
92  radToDeg(Foam::acos(min(1.0, (d & s)/(mag(d)*magS + vSmall))));
93 
94  result[faceCells[facei]] = max(cosDDotS, result[faceCells[facei]]);
95  }
96  }
97 
98  return tresult;
99 }
100 
101 
102 Foam::tmp<Foam::scalarField> Foam::cellQuality::skewness() const
103 {
104  tmp<scalarField> tresult
105  (
106  new scalarField
107  (
108  mesh_.nCells(), 0.0
109  )
110  );
111  scalarField& result = tresult.ref();
112 
113  scalarField sumArea(mesh_.nCells(), 0.0);
114 
115  const vectorField& cellCtrs = mesh_.cellCentres();
116  const vectorField& faceCtrs = mesh_.faceCentres();
117  const vectorField& areas = mesh_.faceAreas();
118 
119  const labelList& own = mesh_.faceOwner();
120  const labelList& nei = mesh_.faceNeighbour();
121 
122  forAll(nei, facei)
123  {
124  scalar dOwn = mag
125  (
126  (faceCtrs[facei] - cellCtrs[own[facei]]) & areas[facei]
127  )/mag(areas[facei]);
128 
129  scalar dNei = mag
130  (
131  (cellCtrs[nei[facei]] - faceCtrs[facei]) & areas[facei]
132  )/mag(areas[facei]);
133 
134  point faceIntersection =
135  cellCtrs[own[facei]]
136  + (dOwn/(dOwn+dNei))*(cellCtrs[nei[facei]] - cellCtrs[own[facei]]);
137 
138  scalar skewness =
139  mag(faceCtrs[facei] - faceIntersection)
140  /(mag(cellCtrs[nei[facei]] - cellCtrs[own[facei]]) + vSmall);
141 
142  result[own[facei]] = max(skewness, result[own[facei]]);
143 
144  result[nei[facei]] = max(skewness, result[nei[facei]]);
145  }
146 
147  forAll(mesh_.boundaryMesh(), patchi)
148  {
149  const labelUList& faceCells =
150  mesh_.boundaryMesh()[patchi].faceCells();
151 
152  const vectorField::subField faceCentres =
153  mesh_.boundaryMesh()[patchi].faceCentres();
154 
155  const vectorField::subField faceAreas =
156  mesh_.boundaryMesh()[patchi].faceAreas();
157 
158  forAll(faceCentres, facei)
159  {
160  vector n = faceAreas[facei]/mag(faceAreas[facei]);
161 
162  point faceIntersection =
163  cellCtrs[faceCells[facei]]
164  + ((faceCentres[facei] - cellCtrs[faceCells[facei]])&n)*n;
165 
166  scalar skewness =
167  mag(faceCentres[facei] - faceIntersection)
168  /(
169  mag(faceCentres[facei] - cellCtrs[faceCells[facei]])
170  + vSmall
171  );
172 
173  result[faceCells[facei]] = max(skewness, result[faceCells[facei]]);
174  }
175  }
176 
177  return tresult;
178 }
179 
180 
181 Foam::tmp<Foam::scalarField> Foam::cellQuality::faceNonOrthogonality() const
182 {
183  tmp<scalarField> tresult
184  (
185  new scalarField
186  (
187  mesh_.nFaces(), 0.0
188  )
189  );
190  scalarField& result = tresult.ref();
191 
192 
193  const vectorField& centres = mesh_.cellCentres();
194  const vectorField& areas = mesh_.faceAreas();
195 
196  const labelList& own = mesh_.faceOwner();
197  const labelList& nei = mesh_.faceNeighbour();
198 
199  forAll(nei, facei)
200  {
201  vector d = centres[nei[facei]] - centres[own[facei]];
202  vector s = areas[facei];
203  scalar magS = mag(s);
204 
205  scalar cosDDotS =
206  radToDeg(Foam::acos(min(1.0, (d & s)/(mag(d)*magS + vSmall))));
207 
208  result[facei] = cosDDotS;
209  }
210 
211  label globalFacei = mesh_.nInternalFaces();
212 
213  forAll(mesh_.boundaryMesh(), patchi)
214  {
215  const labelUList& faceCells =
216  mesh_.boundaryMesh()[patchi].faceCells();
217 
218  const vectorField::subField faceCentres =
219  mesh_.boundaryMesh()[patchi].faceCentres();
220 
221  const vectorField::subField faceAreas =
222  mesh_.boundaryMesh()[patchi].faceAreas();
223 
224  forAll(faceCentres, facei)
225  {
226  vector d = faceCentres[facei] - centres[faceCells[facei]];
227  vector s = faceAreas[facei];
228  scalar magS = mag(s);
229 
230  scalar cosDDotS =
231  radToDeg(Foam::acos(min(1.0, (d & s)/(mag(d)*magS + vSmall))));
232 
233  result[globalFacei++] = cosDDotS;
234  }
235  }
236 
237  return tresult;
238 }
239 
240 
241 Foam::tmp<Foam::scalarField> Foam::cellQuality::faceSkewness() const
242 {
243  tmp<scalarField> tresult
244  (
245  new scalarField
246  (
247  mesh_.nFaces(), 0.0
248  )
249  );
250  scalarField& result = tresult.ref();
251 
252 
253  const vectorField& cellCtrs = mesh_.cellCentres();
254  const vectorField& faceCtrs = mesh_.faceCentres();
255  const vectorField& areas = mesh_.faceAreas();
256 
257  const labelList& own = mesh_.faceOwner();
258  const labelList& nei = mesh_.faceNeighbour();
259 
260  forAll(nei, facei)
261  {
262  scalar dOwn = mag
263  (
264  (faceCtrs[facei] - cellCtrs[own[facei]]) & areas[facei]
265  )/mag(areas[facei]);
266 
267  scalar dNei = mag
268  (
269  (cellCtrs[nei[facei]] - faceCtrs[facei]) & areas[facei]
270  )/mag(areas[facei]);
271 
272  point faceIntersection =
273  cellCtrs[own[facei]]
274  + (dOwn/(dOwn+dNei))*(cellCtrs[nei[facei]] - cellCtrs[own[facei]]);
275 
276  result[facei] =
277  mag(faceCtrs[facei] - faceIntersection)
278  /(mag(cellCtrs[nei[facei]] - cellCtrs[own[facei]]) + vSmall);
279  }
280 
281 
282  label globalFacei = mesh_.nInternalFaces();
283 
284  forAll(mesh_.boundaryMesh(), patchi)
285  {
286  const labelUList& faceCells =
287  mesh_.boundaryMesh()[patchi].faceCells();
288 
289  const vectorField::subField faceCentres =
290  mesh_.boundaryMesh()[patchi].faceCentres();
291 
292  const vectorField::subField faceAreas =
293  mesh_.boundaryMesh()[patchi].faceAreas();
294 
295  forAll(faceCentres, facei)
296  {
297  vector n = faceAreas[facei]/mag(faceAreas[facei]);
298 
299  point faceIntersection =
300  cellCtrs[faceCells[facei]]
301  + ((faceCentres[facei] - cellCtrs[faceCells[facei]])&n)*n;
302 
303  result[globalFacei++] =
304  mag(faceCentres[facei] - faceIntersection)
305  /(
306  mag(faceCentres[facei] - cellCtrs[faceCells[facei]])
307  + vSmall
308  );
309  }
310  }
311 
312  return tresult;
313 }
314 
315 
316 // ************************************************************************* //
Foam::polyMesh::boundaryMesh
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:453
Foam::acos
dimensionedScalar acos(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:281
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
Foam::max
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
Definition: fvMeshStitchersMoving.C:102
Foam::cellQuality::nonOrthogonality
tmp< scalarField > nonOrthogonality() const
Return cell non-orthogonality.
Definition: cellQuality.C:40
Foam::radToDeg
scalar radToDeg(const scalar rad)
Conversion from radians to degrees.
Definition: unitConversion.H:51
Foam::primitiveMesh::nInternalFaces
label nInternalFaces() const
Definition: primitiveMeshI.H:82
Foam::polyMesh::faceNeighbour
virtual const labelList & faceNeighbour() const
Return face neighbour.
Definition: polyMesh.C:1243
Foam::primitiveMesh::nFaces
label nFaces() const
Definition: primitiveMeshI.H:88
Foam::tmp::ref
T & ref() const
Return non-const reference or generate a fatal error.
Definition: tmpI.H:181
unitConversion.H
Unit conversion functions.
Foam::primitiveMesh::nCells
label nCells() const
Definition: primitiveMeshI.H:94
Foam::cellQuality::faceSkewness
tmp< scalarField > faceSkewness() const
Return face skewness.
Definition: cellQuality.C:241
Foam::SubField
Pre-declare related SubField type.
Definition: Field.H:60
mesh
fvMesh & mesh
Definition: createMeshesPostProcess.H:9
s
gmvFile<< "tracers "<< particles.size()<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().x()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().y()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().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))
Foam::cellQuality::faceNonOrthogonality
tmp< scalarField > faceNonOrthogonality() const
Return face non-orthogonality.
Definition: cellQuality.C:181
Foam::cellQuality::cellQuality
cellQuality(const polyMesh &mesh)
Construct from mesh.
Definition: cellQuality.C:32
Foam::polyMesh::faceOwner
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1237
Foam::min
layerAndWeight min(const layerAndWeight &a, const layerAndWeight &b)
Definition: fvMeshStitchersMoving.C:108
Foam::primitiveMesh::cellCentres
const vectorField & cellCentres() const
Definition: primitiveMeshCellCentresAndVols.C:160
Foam::UList
A 1D vector of objects of type <T>, where the size of the vector is known and can be used for subscri...
Definition: HashTable.H:60
Foam::primitiveMesh::faceCentres
const vectorField & faceCentres() const
Definition: primitiveMeshFaceCentresAndAreas.C:162
patchi
label patchi
Definition: getPatchFieldScalar.H:1
Foam::primitiveMesh::faceAreas
const vectorField & faceAreas() const
Definition: primitiveMeshFaceCentresAndAreas.C:173
Foam::mag
dimensioned< scalar > mag(const dimensioned< Type > &)
n
label n
Definition: TABSMDCalcMethod2.H:31
Foam::polyMesh
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:76
Foam::tmp
A class for managing temporary objects.
Definition: PtrList.H:53
Foam::cellQuality::skewness
tmp< scalarField > skewness() const
Return cell skewness.
Definition: cellQuality.C:102