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 | Copyright (C) 2011-2013 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 
41 {
42  tmp<scalarField> tresult
43  (
44  new scalarField
45  (
46  mesh_.nCells(), 0.0
47  )
48  );
49 
50  scalarField& result = tresult();
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 
103 {
104  tmp<scalarField> tresult
105  (
106  new scalarField
107  (
108  mesh_.nCells(), 0.0
109  )
110  );
111  scalarField& result = tresult();
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 
182 {
183  tmp<scalarField> tresult
184  (
185  new scalarField
186  (
187  mesh_.nFaces(), 0.0
188  )
189  );
190  scalarField& result = tresult();
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 
242 {
243  tmp<scalarField> tresult
244  (
245  new scalarField
246  (
247  mesh_.nFaces(), 0.0
248  )
249  );
250  scalarField& result = tresult();
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 // ************************************************************************* //
const vectorField & faceAreas() const
Pre-declare related SubField type.
Definition: Field.H:61
label nFaces() const
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 ))
dimensioned< scalar > mag(const dimensioned< Type > &)
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
label nCells() const
const vectorField & cellCentres() const
dynamicFvMesh & mesh
tmp< scalarField > skewness() const
Return cell skewness.
Definition: cellQuality.C:102
label n
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define forAll(list, i)
Definition: UList.H:421
scalar radToDeg(const scalar rad)
Conversion from radians to degrees.
Unit conversion functions.
dimensionedScalar acos(const dimensionedScalar &ds)
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1073
tmp< scalarField > nonOrthogonality() const
Return cell non-orthogonality.
Definition: cellQuality.C:40
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:421
label nInternalFaces() const
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
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
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
virtual const labelList & faceNeighbour() const
Return face neighbour.
Definition: polyMesh.C:1079
tmp< scalarField > faceNonOrthogonality() const
Return face non-orthogonality.
Definition: cellQuality.C:181
A class for managing temporary objects.
Definition: PtrList.H:118
const vectorField & faceCentres() const
tmp< scalarField > faceSkewness() const
Return face skewness.
Definition: cellQuality.C:241