src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheckMotion.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2009 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 the
  13     Free Software Foundation; either version 2 of the License, or (at your
  14     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, write to the Free Software Foundation,
  23     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
  24
  25 Description
  26     Given a set of points, find out if the mesh resulting from point motion will
  27     be valid without actually changing the mesh.
  28
  29 \*---------------------------------------------------------------------------*/
  30
  31 #include "primitiveMesh.H"
  32 #include "pyramidPointFaceRef.H"
  33 #include "cell.H"
  34 #include "mathematicalConstants.H"
  35
  36 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  37
  38
  39 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  40
  41 bool Foam::primitiveMesh::checkMeshMotion
  42 (
  43     const pointField& newPoints,
  44     const bool report
  45 ) const
  46 {
  47     if (debug || report)
  48     {
  49         Pout<< "bool primitiveMesh::checkMeshMotion("
  50             << "const pointField& newPoints, const bool report) const: "
  51             << "checking mesh motion" << endl;
  52     }
  53
  54     bool error = false;
  55
  56     const faceList& f = faces();
  57
  58     const labelList& own = faceOwner();
  59     const labelList& nei = faceNeighbour();
  60
  61     vectorField fCtrs(nFaces());
  62     vectorField fAreas(nFaces());
  63
  64     makeFaceCentresAndAreas(newPoints, fCtrs, fAreas);
  65
  66     // Check cell volumes and calculate new cell centres
  67     vectorField cellCtrs(nCells());
  68     scalarField cellVols(nCells());
  69
  70     makeCellCentresAndVols(fCtrs, fAreas, cellCtrs, cellVols);
  71
  72     scalar minVolume = GREAT;
  73     label nNegVols = 0;
  74
  75     forAll (cellVols, cellI)
  76     {
  77         if (cellVols[cellI] < VSMALL)
  78         {
  79             if (debug || report)
  80             {
  81                 Pout<< "Zero or negative cell volume detected for cell "
  82                     << cellI << ".  Volume = " << cellVols[cellI] << endl;
  83             }
  84
  85             nNegVols++;
  86         }
  87
  88         minVolume = min(minVolume, cellVols[cellI]);
  89     }
  90
  91     if (nNegVols > 0)
  92     {
  93         error = true;
  94
  95         Pout<< "Zero or negative cell volume in mesh motion in " << nNegVols
  96             << " cells.  Min volume: " << minVolume << endl;
  97     }
  98     else
  99     {
 100         if (debug || report)
 101         {
 102             Pout<< "Min volume = " << minVolume
 103                 << ".  Total volume = " << sum(cellVols)
 104                 << ".  Cell volumes OK." << endl;
 105         }
 106     }
 107
 108     // Check face areas
 109
 110     scalar minArea = GREAT;
 111     label nNegAreas = 0;
 112     label nPyrErrors = 0;
 113     label nDotProductErrors = 0;
 114
 115     forAll (f, faceI)
 116     {
 117         const scalar a = Foam::mag(fAreas[faceI]);
 118
 119         if (a < VSMALL)
 120         {
 121             if (debug || report)
 122             {
 123                 if (isInternalFace(faceI))
 124                 {
 125                     Pout<< "Zero or negative face area detected for "
 126                         << "internal face "<< faceI << " between cells "
 127                         << own[faceI] << " and " << nei[faceI]
 128                         << ".  Face area magnitude = " << a << endl;
 129                 }
 130                 else
 131                 {
 132                     Pout<< "Zero or negative face area detected for "
 133                         << "boundary face " << faceI << " next to cell "
 134                         << own[faceI] << ".  Face area magnitude = "
 135                         << a << endl;
 136                 }
 137             }
 138
 139             nNegAreas++;
 140         }
 141
 142         minArea = min(minArea, a);
 143
 144         // Create the owner pyramid - it will have negative volume
 145         scalar pyrVol =
 146             pyramidPointFaceRef(f[faceI], cellCtrs[own[faceI]]).mag(newPoints);
 147
 148         if (pyrVol > SMALL)
 149         {
 150             if (debug || report)
 151             {
 152                 Pout<< "Negative pyramid volume: " << -pyrVol
 153                     << " for face " << faceI << " " << f[faceI]
 154                     << "  and owner cell: " << own[faceI] << endl
 155                     << "Owner cell vertex labels: "
 156                     << cells()[own[faceI]].labels(f)
 157                     << endl;
 158             }
 159
 160             nPyrErrors++;
 161         }
 162
 163         if (isInternalFace(faceI))
 164         {
 165             // Create the neighbour pyramid - it will have positive volume
 166             scalar pyrVol =
 167                 pyramidPointFaceRef
 168                 (
 169                     f[faceI],
 170                     cellCtrs[nei[faceI]]
 171                 ).mag(newPoints);
 172
 173             if (pyrVol < -SMALL)
 174             {
 175                 if (debug || report)
 176                 {
 177                     Pout<< "Negative pyramid volume: " << pyrVol
 178                         << " for face " << faceI << " " << f[faceI]
 179                         << "  and neighbour cell: " << nei[faceI] << nl
 180                         << "Neighbour cell vertex labels: "
 181                         << cells()[nei[faceI]].labels(f)
 182                         << endl;
 183                 }
 184
 185                 nPyrErrors++;
 186             }
 187
 188             const vector d = cellCtrs[nei[faceI]] - cellCtrs[own[faceI]];
 189             const vector& s = fAreas[faceI];
 190             scalar dDotS = (d & s)/(mag(d)*mag(s) + VSMALL);
 191
 192             // Only write full message the first time
 193             if (dDotS < SMALL && nDotProductErrors == 0)
 194             {
 195                 // Non-orthogonality greater than 90 deg
 196                 WarningIn
 197                 (
 198                     "primitiveMesh::checkMeshMotion"
 199                     "(const pointField& newPoints, const bool report) const"
 200                 )   << "Severe non-orthogonality in mesh motion for face "
 201                     << faceI
 202                     << " between cells " << own[faceI] << " and " << nei[faceI]
 203                     << ": Angle = " << ::acos(dDotS)/mathematicalConstant::pi*180.0
 204                     << " deg." << endl;
 205
 206                 nDotProductErrors++;
 207             }
 208         }
 209     }
 210
 211     if (nNegAreas > 0)
 212     {
 213         error = true;
 214
 215         WarningIn
 216         (
 217             "primitiveMesh::checkMeshMotion"
 218             "(const pointField& newPoints, const bool report) const"
 219         )   << "Zero or negative face area in mesh motion in " << nNegAreas
 220             << " faces.  Min area: " << minArea << endl;
 221     }
 222     else
 223     {
 224         if (debug || report)
 225         {
 226             Pout<< "Min area = " << minArea
 227                 << ".  Face areas OK." << endl;
 228         }
 229     }
 230
 231     if (nPyrErrors > 0)
 232     {
 233         Pout<< "Detected " << nPyrErrors
 234             << " negative pyramid volume in mesh motion" << endl;
 235
 236         error = true;
 237     }
 238     else
 239     {
 240         if (debug || report)
 241         {
 242             Pout<< "Pyramid volumes OK." << endl;
 243         }
 244     }
 245
 246     if (nDotProductErrors > 0)
 247     {
 248         Pout<< "Detected " << nDotProductErrors
 249             << " in non-orthogonality in mesh motion." << endl;
 250
 251         error = true;
 252     }
 253     else
 254     {
 255         if (debug || report)
 256         {
 257             Pout<< "Non-orthogonality check OK." << endl;
 258         }
 259     }
 260
 261     if (!error && (debug || report))
 262     {
 263         Pout << "Mesh motion check OK." << endl;
 264     }
 265
 266     return error;
 267 }
 268
 269
 270 // ************************************************************************* //