src/meshTools/cellQuality/cellQuality.C

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