src/mesh/cfMesh/meshLibrary/utilities/meshes/polyMeshGenAddressing/polyMeshGenAddressingEdges.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 \*---------------------------------------------------------------------------*/
  25
  26 #include "polyMeshGenAddressing.H"
  27 #include "DynList.H"
  28 #include "helperFunctions.H"
  29 #include "demandDrivenData.H"
  30
  31 # ifdef USE_OMP
  32 #include <omp.h>
  33 # endif
  34
  35 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  36
  37 namespace Foam
  38 {
  39
  40 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  41
  42 void polyMeshGenAddressing::calcEdges() const
  43 {
  44     if( edgesPtr_ )
  45     {
  46         FatalErrorIn("polyMeshGenAddressing::calcEdges() const")
  47             << "edges already calculated"
  48             << abort(FatalError);
  49     }
  50     else
  51     {
  52         const faceListPMG& faces = mesh_.faces();
  53         const VRWGraph& pointFaces = this->pointFaces();
  54
  55         edgesPtr_ = new edgeList();
  56         label nEdges(0);
  57         labelList nfe(faces.size());
  58
  59         # ifdef USE_OMP
  60         const label nThreads = 3 * omp_get_num_procs();
  61         labelList nEdgesForThread(nThreads);
  62         # else
  63         labelList nEdgesForThread(1);
  64         # endif
  65
  66         # ifdef USE_OMP
  67         # pragma omp parallel num_threads(nThreads) if( faces.size() > 1000 )
  68         # endif
  69         {
  70             LongList<edge> edgesHelper;
  71
  72             # ifdef USE_OMP
  73             # pragma omp for schedule(static)
  74             # endif
  75             forAll(faces, faceI)
  76             {
  77                 const face& f = faces[faceI];
  78
  79                 forAll(f, pI)
  80                 {
  81                     const edge fe = f.faceEdge(pI);
  82                     const label s = fe.start();
  83                     const label e = fe.end();
  84
  85                     DynList<label> edgeFaces;
  86
  87                     bool store(true);
  88
  89                     //- find all faces attached to this edge
  90                     //- store the edge in case the face faceI is the face
  91                     //- with the smallest label
  92                     forAllRow(pointFaces, s, pfI)
  93                     {
  94                         const label ofI = pointFaces(s, pfI);
  95                         const face& of = faces[ofI];
  96
  97                         if( of.which(e) < 0 )
  98                             continue;
  99                         if( ofI < faceI )
 100                         {
 101                             store = false;
 102                             break;
 103                         }
 104
 105                         edgeFaces.append(ofI);
 106                     }
 107
 108                     if( store )
 109                         edgesHelper.append(fe);
 110                 }
 111             }
 112
 113             //- this enables other threads to see the number of edges
 114             //- generated by each thread
 115             # ifdef USE_OMP
 116             nEdgesForThread[omp_get_thread_num()] = edgesHelper.size();
 117             # else
 118             nEdgesForThread[0] = edgesHelper.size();
 119             # endif
 120
 121             # ifdef USE_OMP
 122             # pragma omp critical
 123             # endif
 124             nEdges += edgesHelper.size();
 125
 126             # ifdef USE_OMP
 127             # pragma omp barrier
 128
 129             # pragma omp master
 130             # endif
 131             edgesPtr_->setSize(nEdges);
 132
 133             # ifdef USE_OMP
 134             # pragma omp barrier
 135             # endif
 136
 137             //- find the starting position of the edges generated by this thread
 138             //- in the global list of edges
 139             label localStart(0);
 140             # ifdef USE_OMP
 141             const label threadI = omp_get_thread_num();
 142             # else
 143             const label threadI = 0;
 144             # endif
 145             for(label i=0;i<threadI;++i)
 146                 localStart += nEdgesForThread[i];
 147
 148             //- store edges into the global list
 149             forAll(edgesHelper, i)
 150                 edgesPtr_->operator[](localStart+i) = edgesHelper[i];
 151         }
 152     }
 153 }
 154
 155 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 156
 157 const edgeList& polyMeshGenAddressing::edges() const
 158 {
 159     if( !edgesPtr_ )
 160     {
 161         # ifdef USE_OMP
 162         if( omp_in_parallel() )
 163             FatalErrorIn
 164             (
 165                 "const edgeList& polyMeshGenAddressing::edges() const"
 166             ) << "Calculating addressing inside a parallel region."
 167                 << " This is not thread safe" << exit(FatalError);
 168         # endif
 169
 170         calcEdges();
 171     }
 172
 173     return *edgesPtr_;
 174 }
 175
 176 void polyMeshGenAddressing::clearOutEdges()
 177 {
 178     deleteDemandDrivenData(edgesPtr_);
 179 }
 180
 181
 182 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 183
 184 } // End namespace Foam
 185
 186 // ************************************************************************* //