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-------------------------------------------------------------------------------
10License
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
34Foam::cellQuality::cellQuality(const polyMesh& mesh)
35:
36 mesh_(mesh)
37{}
38
39
40// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
41
42Foam::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
104Foam::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
183Foam::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
243Foam::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// ************************************************************************* //
n
label n
Definition: TABSMDCalcMethod2.H:31
Foam::SubField
SubField is a Field obtained as a section of another Field, without its own allocation....
Definition: SubField.H:62
Foam::cellQuality
Class calculates cell quality measures.
Definition: cellQuality.H:52
Foam::cellQuality::nonOrthogonality
tmp< scalarField > nonOrthogonality() const
Return cell non-orthogonality.
Definition: cellQuality.C:42
Foam::cellQuality::skewness
tmp< scalarField > skewness() const
Return cell skewness.
Definition: cellQuality.C:104
Foam::cellQuality::faceSkewness
tmp< scalarField > faceSkewness() const
Return face skewness.
Definition: cellQuality.C:243
Foam::cellQuality::faceNonOrthogonality
tmp< scalarField > faceNonOrthogonality() const
Return face non-orthogonality.
Definition: cellQuality.C:183
Foam::faceCells
Smooth ATC in cells next to a set of patches supplied by type.
Definition: faceCells.H:59
Foam::polyMesh
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:81
Foam::tmp
A class for managing temporary objects.
Definition: tmp.H:65
Foam::tmp::ref
T & ref() const
Definition: tmpI.H:227
mesh
dynamicFvMesh & mesh
Definition: createDynamicFvMesh.H:6
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))
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
Foam::mag
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
Foam::radToDeg
constexpr scalar radToDeg() noexcept
Multiplication factor for radians to degrees conversion.
Definition: unitConversion.H:66
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::Zero
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
Foam::acos
dimensionedScalar acos(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:268
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333
unitConversion.H
Unit conversion functions.