src/meshTools/cellQuality/cellQuality.C

   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"
  30 #include "mathematicalConstants.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
  43 Foam::tmp<Foam::scalarField> Foam::cellQuality::nonOrthogonality() const
  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)))
  71             *180.0/mathematicalConstant::pi;
  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)))
  97                *180.0/mathematicalConstant::pi;
  98
  99             result[faceCells[faceI]] = max(cosDDotS, result[faceCells[faceI]]);
 100         }
 101     }
 102
 103     return tresult;
 104 }
 105
 106
 107 Foam::tmp<Foam::scalarField> Foam::cellQuality::skewness() const
 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
 186 Foam::tmp<Foam::scalarField> Foam::cellQuality::faceNonOrthogonality() const
 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)))
 212             *180.0/mathematicalConstant::pi;
 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)))
 238                *180.0/mathematicalConstant::pi;
 239
 240             result[globalFaceI++] = cosDDotS;
 241         }
 242     }
 243
 244     return tresult;
 245 }
 246
 247
 248 Foam::tmp<Foam::scalarField> Foam::cellQuality::faceSkewness() const
 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 // ************************************************************************* //