src/OpenFOAM/meshes/primitiveMesh/primitiveMeshFaceCentresAndAreas.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright held by original author
   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     Calulate the face centres and areas.
  27
  28     Calculate the centre by breaking the face into triangles using the face
  29     centre and area-weighted averaging their centres.  This method copes with
  30     small face-concavity.
  31
  32 \*---------------------------------------------------------------------------*/
  33
  34 #include "primitiveMesh.H"
  35
  36
  37 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  38
  39 void Foam::primitiveMesh::calcFaceCentresAndAreas() const
  40 {
  41     if (debug)
  42     {
  43         Pout<< "primitiveMesh::calcFaceCentresAndAreas() : "
  44             << "Calculating face centres and face areas"
  45             << endl;
  46     }
  47
  48     // It is an error to attempt to recalculate faceCentres
  49     // if the pointer is already set
  50     if (faceCentresPtr_ || faceAreasPtr_)
  51     {
  52         FatalErrorIn("primitiveMesh::calcFaceCentresAndAreas() const")
  53             << "Face centres or face areas already calculated"
  54             << abort(FatalError);
  55     }
  56
  57     faceCentresPtr_ = new vectorField(nFaces());
  58     vectorField& fCtrs = *faceCentresPtr_;
  59
  60     faceAreasPtr_ = new vectorField(nFaces());
  61     vectorField& fAreas = *faceAreasPtr_;
  62
  63     makeFaceCentresAndAreas(points(), fCtrs, fAreas);
  64
  65     if (debug)
  66     {
  67         Pout<< "primitiveMesh::calcFaceCentresAndAreas() : "
  68             << "Finished calculating face centres and face areas"
  69             << endl;
  70     }
  71 }
  72
  73
  74 void Foam::primitiveMesh::makeFaceCentresAndAreas
  75 (
  76     const pointField& p,
  77     vectorField& fCtrs,
  78     vectorField& fAreas
  79 ) const
  80 {
  81     const faceList& fs = faces();
  82
  83     forAll (fs, facei)
  84     {
  85         const labelList& f = fs[facei];
  86         label nPoints = f.size();
  87
  88         // If the face is a triangle, do a direct calculation for efficiency
  89         // and to avoid round-off error-related problems
  90         if (nPoints == 3)
  91         {
  92             fCtrs[facei] = (1.0/3.0)*(p[f[0]] + p[f[1]] + p[f[2]]);
  93             fAreas[facei] = 0.5*((p[f[1]] - p[f[0]])^(p[f[2]] - p[f[0]]));
  94         }
  95         else
  96         {
  97             vector sumN = vector::zero;
  98             scalar sumA = 0.0;
  99             vector sumAc = vector::zero;
 100
 101             point fCentre = p[f[0]];
 102             for (label pi = 1; pi < nPoints; pi++)
 103             {
 104                 fCentre += p[f[pi]];
 105             }
 106
 107             fCentre /= nPoints;
 108
 109             for (label pi = 0; pi < nPoints; pi++)
 110             {
 111                 const point& nextPoint = p[f[(pi + 1) % nPoints]];
 112
 113                 vector c = p[f[pi]] + nextPoint + fCentre;
 114                 vector n = (nextPoint - p[f[pi]])^(fCentre - p[f[pi]]);
 115                 scalar a = mag(n);
 116
 117                 sumN += n;
 118                 sumA += a;
 119                 sumAc += a*c;
 120             }
 121
 122             fCtrs[facei] = (1.0/3.0)*sumAc/(sumA + VSMALL);
 123             fAreas[facei] = 0.5*sumN;
 124         }
 125     }
 126 }
 127
 128
 129 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 130
 131 const Foam::vectorField& Foam::primitiveMesh::faceCentres() const
 132 {
 133     if (!faceCentresPtr_)
 134     {
 135         calcFaceCentresAndAreas();
 136     }
 137
 138     return *faceCentresPtr_;
 139 }
 140
 141
 142 const Foam::vectorField& Foam::primitiveMesh::faceAreas() const
 143 {
 144     if (!faceAreasPtr_)
 145     {
 146         calcFaceCentresAndAreas();
 147     }
 148
 149     return *faceAreasPtr_;
 150 }
 151
 152
 153 // ************************************************************************* //