src/mesh/cfMesh/utilities/surfaceTools/meshSurfaceCheckInvertedVertices/meshSurfaceCheckInvertedVertices.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 "meshSurfaceCheckInvertedVertices.H"
  29 #include "meshSurfacePartitioner.H"
  30 #include "boolList.H"
  31 #include "demandDrivenData.H"
  32 #include "refLabelledPoint.H"
  33 #include "helperFunctions.H"
  34
  35 #include <map>
  36
  37 # ifdef USE_OMP
  38 #include <omp.h>
  39 # endif
  40
  41 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  42
  43 namespace Foam
  44 {
  45
  46 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  47
  48 void meshSurfaceCheckInvertedVertices::checkVertices()
  49 {
  50     const labelList& facePatch = surfacePartitioner_.boundaryFacePatches();
  51     const meshSurfaceEngine& mse = surfacePartitioner_.surfaceEngine();
  52     const pointFieldPMG& points = mse.points();
  53     const labelList& bp = mse.bp();
  54     const VRWGraph& pointFaces = mse.pointFaces();
  55     const VRWGraph& pointInFaces = mse.pointInFaces();
  56     const faceList::subList& bFaces = mse.boundaryFaces();
  57     const vectorField& pNormals = mse.pointNormals();
  58     const vectorField& fCentres = mse.faceCentres();
  59     const vectorField& fNormals = mse.faceNormals();
  60
  61     const labelHashSet& corners = surfacePartitioner_.corners();
  62     const labelHashSet& edgePoints = surfacePartitioner_.edgePoints();
  63
  64     typedef std::map<label, vector> ltvMap;
  65     typedef std::map<label, ltvMap> lltvMap;
  66     lltvMap pointPatchNormal;
  67
  68     forAllConstIter(labelHashSet, corners, it)
  69     {
  70         const label bpI = it.key();
  71
  72         if( activePointsPtr_ && !activePointsPtr_->operator[](bpI))
  73             continue;
  74
  75         ltvMap& patchNormal = pointPatchNormal[bpI];
  76
  77         forAllRow(pointFaces, bpI, pfI)
  78         {
  79             const label bfI = pointFaces(bpI, pfI);
  80             const label patchI = facePatch[bfI];
  81
  82             if( patchNormal.find(patchI) == patchNormal.end() )
  83             {
  84                 patchNormal[patchI] = fNormals[bfI];
  85             }
  86             else
  87             {
  88                 patchNormal[patchI] += fNormals[bfI];
  89             }
  90         }
  91     }
  92
  93     forAllConstIter(labelHashSet, edgePoints, it)
  94     {
  95         const label bpI = it.key();
  96
  97         if( activePointsPtr_ && !activePointsPtr_->operator[](bpI))
  98             continue;
  99
 100         ltvMap& patchNormal = pointPatchNormal[bpI];
 101
 102         forAllRow(pointFaces, bpI, pfI)
 103         {
 104             const label bfI = pointFaces(bpI, pfI);
 105             const label patchI = facePatch[bfI];
 106
 107             if( patchNormal.find(patchI) == patchNormal.end() )
 108             {
 109                 patchNormal[patchI] = fNormals[bfI];
 110             }
 111             else
 112             {
 113                 patchNormal[patchI] += fNormals[bfI];
 114             }
 115         }
 116     }
 117
 118     if( Pstream::parRun() )
 119     {
 120         const Map<label>& globalToLocal = mse.globalToLocalBndPointAddressing();
 121         const DynList<label>& neiProcs = mse.bpNeiProcs();
 122         const VRWGraph& bpAtProcs = mse.bpAtProcs();
 123
 124         std::map<label, LongList<refLabelledPoint> > exchangeData;
 125         forAll(neiProcs, i)
 126             exchangeData[neiProcs[i]].clear();
 127
 128         forAllConstIter(Map<label>, globalToLocal, it)
 129         {
 130             const label bpI = it();
 131
 132             if( pointPatchNormal.find(bpI) != pointPatchNormal.end() )
 133             {
 134                 const ltvMap& patchNormal = pointPatchNormal[bpI];
 135
 136                 forAllRow(bpAtProcs, bpI, i)
 137                 {
 138                     const label neiProc = bpAtProcs(bpI, i);
 139
 140                     if( neiProc == Pstream::myProcNo() )
 141                         continue;
 142
 143                     forAllConstIter(ltvMap, patchNormal, pIt)
 144                         exchangeData[neiProc].append
 145                         (
 146                             refLabelledPoint
 147                             (
 148                                 it.key(),
 149                                 labelledPoint(pIt->first, pIt->second)
 150                             )
 151                         );
 152                 }
 153             }
 154         }
 155
 156         LongList<refLabelledPoint> receivedData;
 157         help::exchangeMap(exchangeData, receivedData);
 158
 159         forAll(receivedData, i)
 160         {
 161             const refLabelledPoint& rlp = receivedData[i];
 162             const label bpI = globalToLocal[rlp.objectLabel()];
 163
 164             ltvMap& patchNormal = pointPatchNormal[bpI];
 165
 166             const labelledPoint& lp = rlp.lPoint();
 167             patchNormal[lp.pointLabel()] += lp.coordinates();
 168         }
 169     }
 170
 171     forAllIter(lltvMap, pointPatchNormal, it)
 172     {
 173         ltvMap& patchNormal = pointPatchNormal[it->first];
 174
 175         forAllIter(ltvMap, patchNormal, pIt)
 176         {
 177             const scalar magv = mag(pIt->second) + VSMALL;
 178
 179             pIt->second /= magv;
 180         }
 181     }
 182
 183     invertedVertices_.clear();
 184
 185     # ifdef USE_OMP
 186     # pragma omp parallel for if( pointFaces.size() > 100 ) \
 187     schedule(dynamic, 20)
 188     # endif
 189     forAll(pointFaces, bpI)
 190     {
 191         if( activePointsPtr_ && !activePointsPtr_->operator[](bpI) )
 192             continue;
 193
 194         forAllRow(pointFaces, bpI, pfI)
 195         {
 196             const label pI = pointInFaces(bpI, pfI);
 197             const label bfI = pointFaces(bpI, pfI);
 198
 199             vector pNormal = pNormals[bpI];
 200
 201             if( pointPatchNormal.find(bpI) != pointPatchNormal.end() )
 202                 pNormal = pointPatchNormal[bpI][facePatch[bfI]];
 203
 204             const face& bf = bFaces[bfI];
 205
 206             //- chech the first triangle (with the next node)
 207             triangle<point, point> triNext
 208             (
 209                 points[bf[pI]],
 210                 points[bf.nextLabel(pI)],
 211                 fCentres[bfI]
 212             );
 213
 214             vector nNext = triNext.normal();
 215             scalar mNext = mag(nNext);
 216
 217             //- face has zero area
 218             if( mNext < VSMALL )
 219             {
 220                 # ifdef USE_OMP
 221                 # pragma omp critical
 222                 # endif
 223                 invertedVertices_.insert(bf[pI]);
 224
 225                 break;
 226             }
 227             else
 228             {
 229                 nNext /= mNext;
 230             }
 231
 232             //- collocated points
 233             if( magSqr(triNext.a() - triNext.b()) < VSMALL )
 234             {
 235                 # ifdef USE_OMP
 236                 # pragma omp critical
 237                 # endif
 238                 invertedVertices_.insert(bf[pI]);
 239
 240                 break;
 241             }
 242             if( magSqr(triNext.c() - triNext.a()) < VSMALL )
 243             {
 244                 # ifdef USE_OMP
 245                 # pragma omp critical
 246                 # endif
 247                 invertedVertices_.insert(bf[pI]);
 248
 249                 break;
 250             }
 251
 252             //- normal vector is not visible
 253             if( (nNext & pNormal) < 0.0 )
 254             {
 255                 # ifdef USE_OMP
 256                 # pragma omp critical
 257                 # endif
 258                 invertedVertices_.insert(bf[pI]);
 259
 260                 break;
 261             }
 262
 263             //- check the second triangle (with previous node)
 264             triangle<point, point> triPrev
 265             (
 266                 points[bf[pI]],
 267                 fCentres[bfI],
 268                 points[bf.prevLabel(pI)]
 269             );
 270
 271             vector nPrev = triPrev.normal();
 272             scalar mPrev = mag(nPrev);
 273
 274             //- face has zero area
 275             if( mPrev < VSMALL )
 276             {
 277                 # ifdef USE_OMP
 278                 # pragma omp critical
 279                 # endif
 280                 invertedVertices_.insert(bf[pI]);
 281
 282                 break;
 283             }
 284             else
 285             {
 286                 nPrev /= mPrev;
 287             }
 288
 289             //- collocated points
 290             if( magSqr(triPrev.a() - triPrev.b()) < VSMALL )
 291             {
 292                 # ifdef USE_OMP
 293                 # pragma omp critical
 294                 # endif
 295                 invertedVertices_.insert(bf[pI]);
 296
 297                 break;
 298             }
 299             if( magSqr(triPrev.c() - triPrev.a()) < VSMALL )
 300             {
 301                 # ifdef USE_OMP
 302                 # pragma omp critical
 303                 # endif
 304                 invertedVertices_.insert(bf[pI]);
 305
 306                 break;
 307             }
 308
 309             //- normal vector is not visible
 310             if( (nPrev & pNormal) < 0.0 )
 311             {
 312                 # ifdef USE_OMP
 313                 # pragma omp critical
 314                 # endif
 315                 invertedVertices_.insert(bf[pI]);
 316
 317                 break;
 318             }
 319
 320             //- check whether the normals of both triangles
 321             //- point in the same direction
 322             if( (nNext & nPrev) < 0.0 )
 323             {
 324                 # ifdef USE_OMP
 325                 # pragma omp critical
 326                 # endif
 327                 invertedVertices_.insert(bf[pI]);
 328
 329                 break;
 330             }
 331         }
 332     }
 333
 334     //- check if there exist concave faces
 335     # ifdef USE_OMP
 336     # pragma omp parallel for schedule(dynamic, 50)
 337     # endif
 338     forAll(bFaces, bfI)
 339     {
 340         const face& bf = bFaces[bfI];
 341
 342         DynList<bool> OkPoints;
 343         if( !help::isFaceConvexAndOk(bf, points, OkPoints) )
 344         {
 345             forAll(OkPoints, pI)
 346             {
 347                 if( activePointsPtr_ && !(*activePointsPtr_)[bp[bf[pI]]] )
 348                     continue;
 349
 350                 if( !OkPoints[pI] )
 351                 {
 352                     # ifdef USE_OMP
 353                     # pragma omp critical
 354                     # endif
 355                     {
 356                         invertedVertices_.insert(bf[pI]);
 357                     }
 358                 }
 359             }
 360         }
 361     }
 362
 363     if( Pstream::parRun() )
 364     {
 365         //- exchange global labels of inverted points
 366         const labelList& bPoints = mse.boundaryPoints();
 367         const Map<label>& globalToLocal =
 368             mse.globalToLocalBndPointAddressing();
 369         const VRWGraph& bpAtProcs = mse.bpAtProcs();
 370         const DynList<label>& neiProcs = mse.bpNeiProcs();
 371
 372         std::map<label, labelLongList> shareData;
 373         forAll(neiProcs, i)
 374             shareData.insert(std::make_pair(neiProcs[i], labelLongList()));
 375
 376         forAllConstIter(Map<label>, globalToLocal, iter)
 377         {
 378             const label bpI = iter();
 379
 380             if( !invertedVertices_.found(bPoints[bpI]) )
 381                 continue;
 382
 383             forAllRow(bpAtProcs, bpI, procI)
 384             {
 385                 const label neiProc = bpAtProcs(bpI, procI);
 386
 387                 if( neiProc == Pstream::myProcNo() )
 388                     continue;
 389
 390                 shareData[neiProc].append(iter.key());
 391             }
 392         }
 393
 394         //- exchange data with other processors
 395         labelLongList receivedData;
 396         help::exchangeMap(shareData, receivedData);
 397
 398         forAll(receivedData, i)
 399         {
 400             const label bpI = globalToLocal[receivedData[i]];
 401             invertedVertices_.insert(bPoints[bpI]);
 402         }
 403     }
 404 }
 405
 406 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 407
 408 meshSurfaceCheckInvertedVertices::meshSurfaceCheckInvertedVertices
 409 (
 410     const meshSurfacePartitioner& mpart
 411 )
 412 :
 413     surfacePartitioner_(mpart),
 414     activePointsPtr_(NULL),
 415     invertedVertices_()
 416 {
 417     checkVertices();
 418 }
 419
 420 meshSurfaceCheckInvertedVertices::meshSurfaceCheckInvertedVertices
 421 (
 422     const meshSurfacePartitioner& mpart,
 423     const boolList& activePoints
 424 )
 425 :
 426     surfacePartitioner_(mpart),
 427     activePointsPtr_(&activePoints),
 428     invertedVertices_()
 429 {
 430     checkVertices();
 431 }
 432
 433 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 434
 435 meshSurfaceCheckInvertedVertices::~meshSurfaceCheckInvertedVertices()
 436 {}
 437
 438 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 439
 440 } // End namespace Foam
 441
 442 // ************************************************************************* //