src/mesh/cfMesh/utilities/meshes/triSurf/triSurfAddressing.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 Description
  25
  26 \*---------------------------------------------------------------------------*/
  27
  28 #include "triSurfAddressing.H"
  29 #include "VRWGraphSMPModifier.H"
  30 #include "demandDrivenData.H"
  31
  32 #include <set>
  33
  34 # ifdef USE_OMP
  35 #include <omp.h>
  36 # endif
  37
  38 namespace Foam
  39 {
  40
  41 void triSurfAddressing::calculatePointFacets() const
  42 {
  43     pointFacetsPtr_ = new VRWGraph();
  44
  45     VRWGraphSMPModifier(*pointFacetsPtr_).reverseAddressing(facets_);
  46 }
  47
  48 void triSurfAddressing::calculateEdges() const
  49 {
  50     edgesPtr_ = new edgeLongList();
  51
  52     const VRWGraph& pFacets = pointFacets();
  53
  54     # ifdef USE_OMP
  55     label nThreads = 3 * omp_get_num_procs();
  56     if( pFacets.size() < 1000 )
  57         nThreads = 1;
  58     # else
  59     const label nThreads(1);
  60     # endif
  61
  62     labelList nEdgesForThread(nThreads);
  63
  64     label edgeI(0);
  65
  66     # ifdef USE_OMP
  67     # pragma omp parallel num_threads(nThreads)
  68     # endif
  69     {
  70         edgeLongList edgesHelper;
  71
  72         # ifdef USE_OMP
  73         # pragma omp for schedule(static)
  74         # endif
  75         forAll(pFacets, pI)
  76         {
  77             std::set<std::pair<label, label> > edgesAtPoint;
  78
  79             forAllRow(pFacets, pI, pfI)
  80             {
  81                 const label triI = pFacets(pI, pfI);
  82                 const labelledTri& tri = facets_[triI];
  83
  84                 label pos(-1);
  85                 forAll(tri, i)
  86                 {
  87                     if( tri[i] == pI )
  88                     {
  89                         pos = i;
  90
  91
  92                         if( tri[(pos+1)%3] >= pI )
  93                             edgesAtPoint.insert
  94                             (
  95                                 std::make_pair(pI, tri[(pos+1)%3])
  96                             );
  97                         if( tri[(pos+2)%3] >= pI )
  98                             edgesAtPoint.insert
  99                             (
 100                                 std::make_pair(pI, tri[(pos+2)%3])
 101                             );
 102                     }
 103                 }
 104
 105
 106             }
 107
 108             std::set<std::pair<label, label> >::const_iterator it;
 109             for(it=edgesAtPoint.begin();it!=edgesAtPoint.end();++it)
 110                 edgesHelper.append(edge(it->first, it->second));
 111         }
 112
 113         //- this enables other threads to see the number of edges
 114         //- generated by each thread
 115         # ifdef USE_OMP
 116         const label threadI = omp_get_thread_num();
 117         # else
 118         const label threadI(0);
 119         # endif
 120         nEdgesForThread[threadI] = edgesHelper.size();
 121
 122
 123         # ifdef USE_OMP
 124         # pragma omp critical
 125         # endif
 126         edgeI += edgesHelper.size();
 127
 128         # ifdef USE_OMP
 129         # pragma omp barrier
 130
 131         # pragma omp master
 132         # endif
 133         edgesPtr_->setSize(edgeI);
 134
 135         # ifdef USE_OMP
 136         # pragma omp barrier
 137         # endif
 138
 139         //- find the starting position of the edges generated by this thread
 140         //- in the global list of edges
 141         label localStart(0);
 142         for(label i=0;i<threadI;++i)
 143             localStart += nEdgesForThread[i];
 144
 145         //- store edges into the global list
 146         forAll(edgesHelper, i)
 147             edgesPtr_->operator[](localStart+i) = edgesHelper[i];
 148     }
 149 }
 150
 151 void triSurfAddressing::calculateFacetEdges() const
 152 {
 153     const edgeLongList& edges = this->edges();
 154     const VRWGraph& pointFaces = this->pointFacets();
 155
 156     facetEdgesPtr_ = new VRWGraph(facets_.size(), 3, -1);
 157     VRWGraph& faceEdges = *facetEdgesPtr_;
 158
 159     # ifdef USE_OMP
 160     # pragma omp parallel for schedule(dynamic, 100)
 161     # endif
 162     forAll(edges, edgeI)
 163     {
 164         const edge ee = edges[edgeI];
 165         const label pI = ee.start();
 166
 167         forAllRow(pointFaces, pI, pfI)
 168         {
 169             const label fI = pointFaces(pI, pfI);
 170
 171             const labelledTri& tri = facets_[fI];
 172             forAll(tri, eI)
 173             {
 174                 const edge e(tri[eI], tri[(eI+1)%3]);
 175
 176                 if( e == ee )
 177                 {
 178                     faceEdges(fI, eI) = edgeI;
 179                 }
 180             }
 181         }
 182     }
 183
 184     # ifdef DEBUGTriSurfAddressing
 185     forAll(faceEdges, triI)
 186     {
 187         forAllRow(faceEdges, triI, feI)
 188         {
 189             if( facets_[triI][feI] < 0 || facets_[triI][feI] >= points_.size() )
 190                 FatalErrorIn
 191                 (
 192                     "void triSurfAddressing::calculateFacetEdges() const"
 193                 ) << "Invalid entry in triangle " << triI
 194                   << " " << facets_[triI] << abort(FatalError);
 195
 196             const label edgeI = faceEdges(triI, feI);
 197
 198             if( edgeI < 0 || edgeI >= edges.size() )
 199                 FatalErrorIn
 200                 (
 201                     "void triSurfAddressing::calculateFacetEdges() const"
 202                 ) << "Invalid entry in face " << triI << " "
 203                      << facets_[triI] << " edges "
 204                      << faceEdges[triI] << abort(FatalError);
 205         }
 206     }
 207     # endif
 208 }
 209
 210 void triSurfAddressing::calculateEdgeFacets() const
 211 {
 212     const edgeLongList& edges = this->edges();
 213     const VRWGraph& faceEdges = this->facetEdges();
 214
 215     edgeFacetsPtr_ = new VRWGraph(edges.size());
 216
 217     VRWGraphSMPModifier(*edgeFacetsPtr_).reverseAddressing(faceEdges);
 218 }
 219
 220 void triSurfAddressing::calculatePointEdges() const
 221 {
 222     const edgeLongList& edges = this->edges();
 223
 224     pointEdgesPtr_ = new VRWGraph(points_.size());
 225
 226     pointEdgesPtr_->reverseAddressing(edges);
 227 }
 228
 229 void triSurfAddressing::calculateFacetFacetsEdges() const
 230 {
 231     facetFacetsEdgesPtr_ = new VRWGraph();
 232
 233     const VRWGraph& facetEdges = this->facetEdges();
 234     const VRWGraph& edgeFacets = this->edgeFacets();
 235
 236     facetFacetsEdgesPtr_->setSize(facets_.size());
 237
 238     forAll(facetEdges, facetI)
 239     {
 240         labelHashSet fLabels;
 241
 242         forAllRow(facetEdges, facetI, feI)
 243         {
 244             const label edgeI = facetEdges(facetI, feI);
 245
 246             forAllRow(edgeFacets, edgeI, efI)
 247                 fLabels.insert(edgeFacets(edgeI, efI));
 248         }
 249
 250         facetFacetsEdgesPtr_->setRowSize(facetI, fLabels.size());
 251
 252         label counter(0);
 253         forAllConstIter(labelHashSet, fLabels, iter)
 254             facetFacetsEdgesPtr_->operator()(facetI, counter++) = iter.key();
 255     }
 256 }
 257
 258 void triSurfAddressing::calculatePointNormals() const
 259 {
 260     const VRWGraph& pFacets = pointFacets();
 261     const vectorField& fNormals = facetNormals();
 262
 263     pointNormalsPtr_ = new vectorField(pFacets.size());
 264
 265     const label size = pFacets.size();
 266     # ifdef USE_OMP
 267     # pragma omp parallel for if( size > 1000 )
 268     # endif
 269     for(label pI=0;pI<size;++pI)
 270     {
 271         vector normal(vector::zero);
 272
 273         forAllRow(pFacets, pI, pfI)
 274             normal += fNormals[pFacets(pI, pfI)];
 275
 276         const scalar d = mag(normal);
 277         if( d > VSMALL )
 278         {
 279             normal /= d;
 280         }
 281         else
 282         {
 283             normal = vector::zero;
 284         }
 285
 286         (*pointNormalsPtr_)[pI] = normal;
 287     }
 288 }
 289
 290 void triSurfAddressing::calculateFacetNormals() const
 291 {
 292     facetNormalsPtr_ = new vectorField(facets_.size());
 293
 294     # ifdef USE_OMP
 295     # pragma omp parallel for schedule(dynamic, 40)
 296     # endif
 297     forAll(facets_, fI)
 298     {
 299         vector v = facets_[fI].normal(points_);
 300         v /= (mag(v) + VSMALL);
 301         (*facetNormalsPtr_)[fI] = v;
 302     }
 303 }
 304
 305 void triSurfAddressing::calculateFacetCentres() const
 306 {
 307     facetCentresPtr_ = new vectorField(facets_.size());
 308
 309     # ifdef USE_OMP
 310     # pragma omp parallel for schedule(dynamic, 40)
 311     # endif
 312     forAll(facets_, fI)
 313         (*facetCentresPtr_)[fI] = facets_[fI].centre(points_);
 314 }
 315
 316 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 317
 318 triSurfAddressing::triSurfAddressing
 319 (
 320     const pointField& points,
 321     const LongList<labelledTri>& triangles)
 322 :
 323     points_(points),
 324     facets_(triangles),
 325     pointFacetsPtr_(NULL),
 326     edgesPtr_(NULL),
 327     facetEdgesPtr_(NULL),
 328     edgeFacetsPtr_(NULL),
 329     pointEdgesPtr_(NULL),
 330     facetFacetsEdgesPtr_(NULL),
 331     pointNormalsPtr_(NULL),
 332     facetNormalsPtr_(NULL),
 333     facetCentresPtr_(NULL)
 334
 335 {}
 336
 337 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 338
 339 triSurfAddressing::~triSurfAddressing()
 340 {
 341     clearAddressing();
 342     clearGeometry();
 343 }
 344
 345 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 346
 347 void triSurfAddressing::clearAddressing()
 348 {
 349     deleteDemandDrivenData(pointFacetsPtr_);
 350     deleteDemandDrivenData(edgesPtr_);
 351     deleteDemandDrivenData(facetEdgesPtr_);
 352     deleteDemandDrivenData(edgeFacetsPtr_);
 353     deleteDemandDrivenData(pointEdgesPtr_);
 354     deleteDemandDrivenData(facetFacetsEdgesPtr_);
 355 }
 356
 357 void triSurfAddressing::clearGeometry()
 358 {
 359     deleteDemandDrivenData(pointNormalsPtr_);
 360     deleteDemandDrivenData(facetNormalsPtr_);
 361     deleteDemandDrivenData(facetCentresPtr_);
 362 }
 363
 364 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 365
 366 } // End namespace Foam
 367
 368 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //