src/mesh/cfMesh/utilities/meshes/polyMeshGenAddressing/polyMeshGenAddressingUpdateGeometry.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | cfMesh: A library for mesh generation
   4    \\    /   O peration     |
   5     \\  /    A nd           | Author: Franjo Juretic (franjo.juretic@c-fields.com)
   6      \\/     M anipulation  | Copyright (C) Creative Fields, Ltd.
   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 \*---------------------------------------------------------------------------*/
  25
  26 #include "polyMeshGenAddressing.H"
  27 #include "demandDrivenData.H"
  28
  29 # ifdef USE_OMP
  30 #include <omp.h>
  31 #endif
  32
  33 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  34
  35 namespace Foam
  36 {
  37
  38 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  39
  40 void polyMeshGenAddressing::updateGeometry
  41 (
  42     const boolList& changedFace
  43 )
  44 {
  45     const pointFieldPMG& p = mesh_.points();
  46     const faceListPMG& faces = mesh_.faces();
  47
  48     //- update face centres and face areas
  49     if( faceCentresPtr_ && faceAreasPtr_ )
  50     {
  51         vectorField& fCtrs = *faceCentresPtr_;
  52         vectorField& fAreas = *faceAreasPtr_;
  53
  54         # ifdef USE_OMP
  55         # pragma omp parallel for if( faces.size() > 100 ) \
  56         schedule(dynamic, 10)
  57         # endif
  58         forAll(faces, faceI)
  59             if( changedFace[faceI] )
  60             {
  61                 const face& f = faces[faceI];
  62                 const label nPoints = f.size();
  63
  64                 // If the face is a triangle, do a direct calculation for
  65                 // efficiency and to avoid round-off error-related problems
  66                 if (nPoints == 3)
  67                 {
  68                     fCtrs[faceI] = (1.0/3.0)*(p[f[0]] + p[f[1]] + p[f[2]]);
  69                     fAreas[faceI] =
  70                         0.5*((p[f[1]] - p[f[0]])^(p[f[2]] - p[f[0]]));
  71                 }
  72                 else
  73                 {
  74                     vector sumN = vector::zero;
  75                     scalar sumA = 0.0;
  76                     vector sumAc = vector::zero;
  77
  78                     point fCentre = p[f[0]];
  79                     for(label pI=1;pI<nPoints;++pI)
  80                     {
  81                         fCentre += p[f[pI]];
  82                     }
  83
  84                     fCentre /= nPoints;
  85
  86                     for(label pI=0;pI<nPoints;++pI)
  87                     {
  88                         const point& nextPoint = p[f.nextLabel(pI)];
  89
  90                         vector c = p[f[pI]] + nextPoint + fCentre;
  91                         vector n = (nextPoint - p[f[pI]])^(fCentre - p[f[pI]]);
  92                         scalar a = mag(n);
  93
  94                         sumN += n;
  95                         sumA += a;
  96                         sumAc += a*c;
  97                     }
  98
  99                     fCtrs[faceI] = (1.0/3.0)*sumAc/(sumA + VSMALL);
 100                     fAreas[faceI] = 0.5*sumN;
 101                 }
 102             }
 103     }
 104
 105     //- update cell centres and cell volumes
 106     if( cellCentresPtr_ && cellVolumesPtr_ && faceCentresPtr_ && faceAreasPtr_ )
 107     {
 108         const vectorField& fCtrs = *faceCentresPtr_;
 109         const vectorField& fAreas = *faceAreasPtr_;
 110         vectorField& cellCtrs = *cellCentresPtr_;
 111         scalarField& cellVols = *cellVolumesPtr_;
 112
 113         const labelList& own = mesh_.owner();
 114         const cellListPMG& cells = mesh_.cells();
 115
 116         # ifdef USE_OMP
 117         # pragma omp parallel for if( cells.size() > 100 ) \
 118         schedule(dynamic, 10)
 119         # endif
 120         forAll(cells, cellI)
 121         {
 122             const cell& c = cells[cellI];
 123
 124             bool update(false);
 125             forAll(c, fI)
 126                 if( changedFace[c[fI]] )
 127                 {
 128                     update = true;
 129                     break;
 130                 }
 131
 132             if( update )
 133             {
 134                 cellCtrs[cellI] = vector::zero;
 135                 cellVols[cellI] = 0.0;
 136
 137                 //- estimate position of cell centre
 138                 vector cEst(vector::zero);
 139                 forAll(c, fI)
 140                     cEst += fCtrs[c[fI]];
 141                 cEst /= c.size();
 142
 143                 forAll(c, fI)
 144                     if( own[c[fI]] == cellI )
 145                     {
 146                         // Calculate 3*face-pyramid volume
 147                         const scalar pyr3Vol =
 148                             max
 149                             (
 150                                 fAreas[c[fI]] &
 151                                 (
 152                                     fCtrs[c[fI]] -
 153                                     cEst
 154                                 ),
 155                                 VSMALL
 156                             );
 157
 158                         // Calculate face-pyramid centre
 159                         const vector pc =
 160                             (3.0/4.0)*fCtrs[c[fI]] + (1.0/4.0)*cEst;
 161
 162                         // Accumulate volume-weighted face-pyramid centre
 163                         cellCtrs[cellI] += pyr3Vol*pc;
 164
 165                         // Accumulate face-pyramid volume
 166                         cellVols[cellI] += pyr3Vol;
 167                     }
 168                     else
 169                     {
 170                         // Calculate 3*face-pyramid volume
 171                         const scalar pyr3Vol =
 172                             max
 173                             (
 174                                 fAreas[c[fI]] &
 175                                 (
 176                                     cEst - fCtrs[c[fI]]
 177                                 ),
 178                                 VSMALL
 179                             );
 180
 181                         // Calculate face-pyramid centre
 182                         const vector pc =
 183                             (3.0/4.0)*fCtrs[c[fI]] + (1.0/4.0)*cEst;
 184
 185                         // Accumulate volume-weighted face-pyramid centre
 186                         cellCtrs[cellI] += pyr3Vol*pc;
 187
 188                         // Accumulate face-pyramid volume
 189                         cellVols[cellI] += pyr3Vol;
 190                     }
 191
 192                 cellCtrs[cellI] /= cellVols[cellI];
 193                 cellVols[cellI] /= 3.0;
 194             }
 195         }
 196     }
 197 }
 198
 199 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 200
 201 } // End namespace Foam
 202
 203 // ************************************************************************* //