FreeFOAM The Cross-Platform CFD Toolkit
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) 1991-2010 OpenCFD Ltd.
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 Description
25  Class calculates cell quality measures.
26 
27 \*---------------------------------------------------------------------------*/
28 
29 #include "cellQuality.H"
31 
32 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
33 
34 // Construct from mesh
35 Foam::cellQuality::cellQuality(const polyMesh& mesh)
36 :
37  mesh_(mesh)
38 {}
39 
40 
41 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
42 
44 {
45  tmp<scalarField> tresult
46  (
47  new scalarField
48  (
49  mesh_.nCells(), 0.0
50  )
51  );
52 
53  scalarField& result = tresult();
54 
55  scalarField sumArea(mesh_.nCells(), 0.0);
56 
57  const vectorField& centres = mesh_.cellCentres();
58  const vectorField& areas = mesh_.faceAreas();
59 
60  const labelList& own = mesh_.faceOwner();
61  const labelList& nei = mesh_.faceNeighbour();
62 
63  forAll (nei, faceI)
64  {
65  vector d = centres[nei[faceI]] - centres[own[faceI]];
66  vector s = areas[faceI];
67  scalar magS = mag(s);
68 
69  scalar cosDDotS =
70  Foam::acos(min(1.0, (d & s)/(mag(d)*magS + VSMALL)))
72 
73  result[own[faceI]] = max(cosDDotS, result[own[faceI]]);
74 
75  result[nei[faceI]] = max(cosDDotS, result[nei[faceI]]);
76  }
77 
78  forAll (mesh_.boundaryMesh(), patchI)
79  {
80  const unallocLabelList& faceCells =
81  mesh_.boundaryMesh()[patchI].faceCells();
82 
83  const vectorField::subField faceCentres =
84  mesh_.boundaryMesh()[patchI].faceCentres();
85 
86  const vectorField::subField faceAreas =
87  mesh_.boundaryMesh()[patchI].faceAreas();
88 
89  forAll(faceCentres, faceI)
90  {
91  vector d = faceCentres[faceI] - centres[faceCells[faceI]];
92  vector s = faceAreas[faceI];
93  scalar magS = mag(s);
94 
95  scalar cosDDotS =
96  Foam::acos(min(1.0, (d & s)/(mag(d)*magS + VSMALL)))
98 
99  result[faceCells[faceI]] = max(cosDDotS, result[faceCells[faceI]]);
100  }
101  }
102 
103  return tresult;
104 }
105 
106 
108 {
109  tmp<scalarField> tresult
110  (
111  new scalarField
112  (
113  mesh_.nCells(), 0.0
114  )
115  );
116  scalarField& result = tresult();
117 
118  scalarField sumArea(mesh_.nCells(), 0.0);
119 
120  const vectorField& cellCtrs = mesh_.cellCentres();
121  const vectorField& faceCtrs = mesh_.faceCentres();
122  const vectorField& areas = mesh_.faceAreas();
123 
124  const labelList& own = mesh_.faceOwner();
125  const labelList& nei = mesh_.faceNeighbour();
126 
127  forAll (nei, faceI)
128  {
129  scalar dOwn = mag
130  (
131  (faceCtrs[faceI] - cellCtrs[own[faceI]]) & areas[faceI]
132  )/mag(areas[faceI]);
133 
134  scalar dNei = mag
135  (
136  (cellCtrs[nei[faceI]] - faceCtrs[faceI]) & areas[faceI]
137  )/mag(areas[faceI]);
138 
139  point faceIntersection =
140  cellCtrs[own[faceI]]
141  + (dOwn/(dOwn+dNei))*(cellCtrs[nei[faceI]] - cellCtrs[own[faceI]]);
142 
143  scalar skewness =
144  mag(faceCtrs[faceI] - faceIntersection)
145  /(mag(cellCtrs[nei[faceI]] - cellCtrs[own[faceI]]) + VSMALL);
146 
147  result[own[faceI]] = max(skewness, result[own[faceI]]);
148 
149  result[nei[faceI]] = max(skewness, result[nei[faceI]]);
150  }
151 
152  forAll (mesh_.boundaryMesh(), patchI)
153  {
154  const unallocLabelList& faceCells =
155  mesh_.boundaryMesh()[patchI].faceCells();
156 
157  const vectorField::subField faceCentres =
158  mesh_.boundaryMesh()[patchI].faceCentres();
159 
160  const vectorField::subField faceAreas =
161  mesh_.boundaryMesh()[patchI].faceAreas();
162 
163  forAll(faceCentres, faceI)
164  {
165  vector n = faceAreas[faceI]/mag(faceAreas[faceI]);
166 
167  point faceIntersection =
168  cellCtrs[faceCells[faceI]]
169  + ((faceCentres[faceI] - cellCtrs[faceCells[faceI]])&n)*n;
170 
171  scalar skewness =
172  mag(faceCentres[faceI] - faceIntersection)
173  /(
174  mag(faceCentres[faceI] - cellCtrs[faceCells[faceI]])
175  + VSMALL
176  );
177 
178  result[faceCells[faceI]] = max(skewness, result[faceCells[faceI]]);
179  }
180  }
181 
182  return tresult;
183 }
184 
185 
187 {
188  tmp<scalarField> tresult
189  (
190  new scalarField
191  (
192  mesh_.nFaces(), 0.0
193  )
194  );
195  scalarField& result = tresult();
196 
197 
198  const vectorField& centres = mesh_.cellCentres();
199  const vectorField& areas = mesh_.faceAreas();
200 
201  const labelList& own = mesh_.faceOwner();
202  const labelList& nei = mesh_.faceNeighbour();
203 
204  forAll (nei, faceI)
205  {
206  vector d = centres[nei[faceI]] - centres[own[faceI]];
207  vector s = areas[faceI];
208  scalar magS = mag(s);
209 
210  scalar cosDDotS =
211  Foam::acos(min(1.0, (d & s)/(mag(d)*magS + VSMALL)))
213 
214  result[faceI] = cosDDotS;
215  }
216 
217  label globalFaceI = mesh_.nInternalFaces();
218 
219  forAll (mesh_.boundaryMesh(), patchI)
220  {
221  const unallocLabelList& faceCells =
222  mesh_.boundaryMesh()[patchI].faceCells();
223 
224  const vectorField::subField faceCentres =
225  mesh_.boundaryMesh()[patchI].faceCentres();
226 
227  const vectorField::subField faceAreas =
228  mesh_.boundaryMesh()[patchI].faceAreas();
229 
230  forAll(faceCentres, faceI)
231  {
232  vector d = faceCentres[faceI] - centres[faceCells[faceI]];
233  vector s = faceAreas[faceI];
234  scalar magS = mag(s);
235 
236  scalar cosDDotS =
237  Foam::acos(min(1.0, (d & s)/(mag(d)*magS + VSMALL)))
239 
240  result[globalFaceI++] = cosDDotS;
241  }
242  }
243 
244  return tresult;
245 }
246 
247 
249 {
250  tmp<scalarField> tresult
251  (
252  new scalarField
253  (
254  mesh_.nFaces(), 0.0
255  )
256  );
257  scalarField& result = tresult();
258 
259 
260  const vectorField& cellCtrs = mesh_.cellCentres();
261  const vectorField& faceCtrs = mesh_.faceCentres();
262  const vectorField& areas = mesh_.faceAreas();
263 
264  const labelList& own = mesh_.faceOwner();
265  const labelList& nei = mesh_.faceNeighbour();
266 
267  forAll (nei, faceI)
268  {
269  scalar dOwn = mag
270  (
271  (faceCtrs[faceI] - cellCtrs[own[faceI]]) & areas[faceI]
272  )/mag(areas[faceI]);
273 
274  scalar dNei = mag
275  (
276  (cellCtrs[nei[faceI]] - faceCtrs[faceI]) & areas[faceI]
277  )/mag(areas[faceI]);
278 
279  point faceIntersection =
280  cellCtrs[own[faceI]]
281  + (dOwn/(dOwn+dNei))*(cellCtrs[nei[faceI]] - cellCtrs[own[faceI]]);
282 
283  result[faceI] =
284  mag(faceCtrs[faceI] - faceIntersection)
285  /(mag(cellCtrs[nei[faceI]] - cellCtrs[own[faceI]]) + VSMALL);
286  }
287 
288 
289  label globalFaceI = mesh_.nInternalFaces();
290 
291  forAll (mesh_.boundaryMesh(), patchI)
292  {
293  const unallocLabelList& faceCells =
294  mesh_.boundaryMesh()[patchI].faceCells();
295 
296  const vectorField::subField faceCentres =
297  mesh_.boundaryMesh()[patchI].faceCentres();
298 
299  const vectorField::subField faceAreas =
300  mesh_.boundaryMesh()[patchI].faceAreas();
301 
302  forAll(faceCentres, faceI)
303  {
304  vector n = faceAreas[faceI]/mag(faceAreas[faceI]);
305 
306  point faceIntersection =
307  cellCtrs[faceCells[faceI]]
308  + ((faceCentres[faceI] - cellCtrs[faceCells[faceI]])&n)*n;
309 
310  result[globalFaceI++] =
311  mag(faceCentres[faceI] - faceIntersection)
312  /(
313  mag(faceCentres[faceI] - cellCtrs[faceCells[faceI]])
314  + VSMALL
315  );
316  }
317  }
318 
319  return tresult;
320 }
321 
322 
323 // ************************ vim: set sw=4 sts=4 et: ************************ //