src/mesh/cfMesh/meshLibrary/utilities/meshes/polyMeshGenAddressing/polyMeshGenAddressingCentresAndAreas.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | cfMesh: A library for mesh generation
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright held by the original author
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of cfMesh.
  10
  11     cfMesh 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 3 of the License, or (at your
  14     option) any later version.
  15
  16     cfMesh 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 cfMesh.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 Description
  25     Calulate the face centres and areas.
  26
  27     Calculate the centre by breaking the face into triangles using the face
  28     centre and area-weighted averaging their centres.  This method copes with
  29     small face-concavity.
  30
  31 \*---------------------------------------------------------------------------*/
  32
  33 #include "polyMeshGenAddressing.H"
  34
  35 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  36
  37 namespace Foam
  38 {
  39
  40 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  41
  42 void polyMeshGenAddressing::calcFaceCentresAndAreas() const
  43 {
  44     if( faceCentresPtr_ || faceAreasPtr_ )
  45     {
  46         FatalErrorIn("polyMeshGenAddressing::calcFaceCentresAndAreas() const")
  47             << "Face centres or face areas already calculated"
  48             << abort(FatalError);
  49     }
  50
  51     const pointFieldPMG& points = mesh_.points();
  52     const faceListPMG& faces = mesh_.faces();
  53
  54     faceCentresPtr_ = new vectorField(faces.size());
  55     vectorField& fCtrs = *faceCentresPtr_;
  56
  57     faceAreasPtr_ = new vectorField(faces.size());
  58     vectorField& fAreas = *faceAreasPtr_;
  59
  60     makeFaceCentresAndAreas(points, fCtrs, fAreas);
  61 }
  62
  63 void polyMeshGenAddressing::makeFaceCentresAndAreas
  64 (
  65     const pointFieldPMG& p,
  66     vectorField& fCtrs,
  67     vectorField& fAreas
  68 ) const
  69 {
  70     const faceListPMG& fs = mesh_.faces();
  71     const label nFaces = fs.size();
  72
  73     # ifdef USE_OMP
  74     # pragma omp parallel for if( nFaces > 1000 )
  75     # endif
  76     for(label faceI=0;faceI<nFaces;++faceI)
  77     {
  78         const face& f = fs[faceI];
  79         label nPoints = f.size();
  80
  81         // If the face is a triangle, do a direct calculation for efficiency
  82         // and to avoid round-off error-related problems
  83         if (nPoints == 3)
  84         {
  85             fCtrs[faceI] = (1.0/3.0)*(p[f[0]] + p[f[1]] + p[f[2]]);
  86             fAreas[faceI] = 0.5*((p[f[1]] - p[f[0]])^(p[f[2]] - p[f[0]]));
  87         }
  88         else
  89         {
  90             vector sumN = vector::zero;
  91             scalar sumA = 0.0;
  92             vector sumAc = vector::zero;
  93
  94             point fCentre = p[f[0]];
  95             for(label pi=1;pi<nPoints;++pi)
  96             {
  97                 fCentre += p[f[pi]];
  98             }
  99
 100             fCentre /= nPoints;
 101
 102             for(label pi=0;pi<nPoints;++pi)
 103             {
 104                 const point& nextPoint = p[f.nextLabel(pi)];
 105
 106                 vector c = p[f[pi]] + nextPoint + fCentre;
 107                 vector n = (nextPoint - p[f[pi]])^(fCentre - p[f[pi]]);
 108                 scalar a = mag(n);
 109
 110                 sumN += n;
 111                 sumA += a;
 112                 sumAc += a*c;
 113             }
 114
 115             fCtrs[faceI] = (1.0/3.0)*sumAc/(sumA + VSMALL);
 116             fAreas[faceI] = 0.5*sumN;
 117         }
 118     }
 119 }
 120
 121 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 122
 123 const vectorField& polyMeshGenAddressing::faceCentres() const
 124 {
 125     if( !faceCentresPtr_ )
 126     {
 127         # ifdef USE_OMP
 128         if( omp_in_parallel() )
 129             FatalErrorIn
 130             (
 131                 "const vectorField& polyMeshGenAddressing::faceCentres() const"
 132             ) << "Calculating addressing inside a parallel region."
 133                 << " This is not thread safe" << exit(FatalError);
 134         # endif
 135
 136         calcFaceCentresAndAreas();
 137     }
 138
 139     return *faceCentresPtr_;
 140 }
 141
 142 const vectorField& polyMeshGenAddressing::faceAreas() const
 143 {
 144     if( !faceAreasPtr_ )
 145     {
 146         # ifdef USE_OMP
 147         if( omp_in_parallel() )
 148             FatalErrorIn
 149             (
 150                 "const vectorField& polyMeshGenAddressing::faceAreas() const"
 151             ) << "Calculating addressing inside a parallel region."
 152                 << " This is not thread safe" << exit(FatalError);
 153         # endif
 154
 155         calcFaceCentresAndAreas();
 156     }
 157
 158     return *faceAreasPtr_;
 159 }
 160
 161 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 162
 163 } // End namespace Foam
 164
 165 // ************************************************************************* //