src/mesh/cfMesh/meshLibrary/utilities/surfaceTools/meshSurfaceEngine/meshSurfaceEngineModifier.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 "meshSurfaceEngineModifier.H"
  29 #include "polyMeshGenModifier.H"
  30 #include "demandDrivenData.H"
  31
  32 #include "labelledPoint.H"
  33 #include "helperFunctionsPar.H"
  34
  35 // #define DEBUGSearch
  36
  37 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  38
  39 namespace Foam
  40 {
  41
  42 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  43
  44 meshSurfaceEngineModifier::meshSurfaceEngineModifier
  45 (
  46     meshSurfaceEngine& surfaceEngine
  47 )
  48 :
  49     surfaceEngine_(surfaceEngine)
  50 {}
  51
  52 meshSurfaceEngineModifier::meshSurfaceEngineModifier
  53 (
  54     const meshSurfaceEngine& surfaceEngine
  55 )
  56 :
  57     surfaceEngine_(const_cast<meshSurfaceEngine&>(surfaceEngine))
  58 {}
  59
  60 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
  61
  62 meshSurfaceEngineModifier::~meshSurfaceEngineModifier()
  63 {}
  64
  65 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  66
  67 void meshSurfaceEngineModifier::moveBoundaryVertexNoUpdate
  68 (
  69     const label bpI,
  70     const point& newP
  71 )
  72 {
  73     surfaceEngine_.mesh_.points()[surfaceEngine_.boundaryPoints()[bpI]] = newP;
  74 }
  75
  76 void meshSurfaceEngineModifier::moveBoundaryVertex
  77 (
  78     const label bpI,
  79     const point& newP
  80 )
  81 {
  82     const labelList& bPoints = surfaceEngine_.boundaryPoints();
  83     pointFieldPMG& points = surfaceEngine_.mesh_.points();
  84     points[bPoints[bpI]] = newP;
  85
  86     if( surfaceEngine_.faceCentresPtr_ )
  87     {
  88         vectorField& faceCentres = *surfaceEngine_.faceCentresPtr_;
  89         const VRWGraph& pFaces = surfaceEngine_.pointFaces();
  90         const faceList::subList& bFaces = surfaceEngine_.boundaryFaces();
  91
  92         forAllRow(pFaces, bpI, pfI)
  93         {
  94             const label bfI = pFaces(bpI, pfI);
  95
  96             faceCentres[bfI] = bFaces[bfI].centre(points);
  97         }
  98     }
  99
 100     if( surfaceEngine_.faceNormalsPtr_ )
 101     {
 102         vectorField& faceNormals = *surfaceEngine_.faceNormalsPtr_;
 103         const VRWGraph& pFaces = surfaceEngine_.pointFaces();
 104         const faceList::subList& bFaces = surfaceEngine_.boundaryFaces();
 105
 106         forAllRow(pFaces, bpI, pfI)
 107         {
 108             const label bfI = pFaces(bpI, pfI);
 109
 110             faceNormals[bfI] = bFaces[bfI].normal(points);
 111         }
 112     }
 113
 114     if( surfaceEngine_.pointNormalsPtr_ )
 115     {
 116         const vectorField& faceNormals = *surfaceEngine_.faceNormalsPtr_;
 117         const VRWGraph& pFaces = surfaceEngine_.pointFaces();
 118         const VRWGraph& pPoints = surfaceEngine_.pointPoints();
 119
 120         vectorField& pn = *surfaceEngine_.pointNormalsPtr_;
 121         vector n(vector::zero);
 122         forAllRow(pFaces, bpI, pfI)
 123             n += faceNormals[pFaces(bpI, pfI)];
 124
 125         const scalar l = mag(n);
 126         if( l > VSMALL )
 127         {
 128             n /= l;
 129         }
 130         else
 131         {
 132             n = vector::zero;
 133         }
 134
 135         pn[bpI] = n;
 136
 137         //- change normal of vertices connected to bpI
 138         forAllRow(pPoints, bpI, ppI)
 139         {
 140             const label bpJ = pPoints(bpI, ppI);
 141             n = vector::zero;
 142             forAllRow(pFaces, bpJ, pfI)
 143                 n += faceNormals[pFaces(bpJ, pfI)];
 144
 145             const scalar d = mag(n);
 146             if( d > VSMALL )
 147             {
 148                 n /= d;
 149             }
 150             else
 151             {
 152                 n = vector::zero;
 153             }
 154
 155             pn[bpJ] = n;
 156         }
 157     }
 158 }
 159
 160 void meshSurfaceEngineModifier::syncVerticesAtParallelBoundaries()
 161 {
 162     if( !Pstream::parRun() )
 163         return;
 164
 165     const Map<label>& globalToLocal =
 166         surfaceEngine_.globalToLocalBndPointAddressing();
 167     labelLongList syncNodes;
 168     forAllConstIter(Map<label>, globalToLocal, it)
 169         syncNodes.append(it());
 170
 171     syncVerticesAtParallelBoundaries(syncNodes);
 172 }
 173
 174 void meshSurfaceEngineModifier::syncVerticesAtParallelBoundaries
 175 (
 176     const labelLongList& syncNodes
 177 )
 178 {
 179     if( !Pstream::parRun() )
 180         return;
 181
 182     const VRWGraph& bpAtProcs = surfaceEngine_.bpAtProcs();
 183     const labelList& globalLabel =
 184         surfaceEngine_.globalBoundaryPointLabel();
 185     const Map<label>& globalToLocal =
 186         surfaceEngine_.globalToLocalBndPointAddressing();
 187     const DynList<label>& neiProcs = surfaceEngine_.bpNeiProcs();
 188     const labelList& bPoints = surfaceEngine_.boundaryPoints();
 189     const pointFieldPMG& points = surfaceEngine_.mesh().points();
 190
 191     std::map<label, LongList<labelledPoint> > exchangeData;
 192     forAll(neiProcs, i)
 193         exchangeData.insert
 194         (
 195             std::make_pair(neiProcs[i], LongList<labelledPoint>())
 196         );
 197
 198     //- construct the map
 199     forAll(syncNodes, snI)
 200     {
 201         const label bpI = syncNodes[snI];
 202
 203         if( bpAtProcs.sizeOfRow(bpI) == 0 )
 204             continue;
 205
 206         point p = points[bPoints[bpI]] / bpAtProcs.sizeOfRow(bpI);
 207         moveBoundaryVertexNoUpdate(bpI, p);
 208
 209         forAllRow(bpAtProcs, bpI, i)
 210         {
 211             const label neiProc = bpAtProcs(bpI, i);
 212             if( neiProc == Pstream::myProcNo() )
 213                 continue;
 214
 215             exchangeData[neiProc].append(labelledPoint(globalLabel[bpI], p));
 216         }
 217     }
 218
 219     //- exchange the data with other processors
 220     LongList<labelledPoint> receivedData;
 221     help::exchangeMap(exchangeData, receivedData);
 222
 223     //- adjust the coordinates
 224     forAll(receivedData, i)
 225     {
 226         const labelledPoint& lp = receivedData[i];
 227         const label bpI = globalToLocal[lp.pointLabel()];
 228         const point newP = points[bPoints[bpI]] + lp.coordinates();
 229         moveBoundaryVertexNoUpdate(bpI, newP);
 230     }
 231 }
 232
 233 void meshSurfaceEngineModifier::updateGeometry
 234 (
 235     const labelLongList& updateBndNodes
 236 )
 237 {
 238     const pointFieldPMG& points = surfaceEngine_.points();
 239     const faceList::subList& bFaces = surfaceEngine_.boundaryFaces();
 240     const VRWGraph& pFaces = surfaceEngine_.pointFaces();
 241     const labelList& bp = surfaceEngine_.bp();
 242
 243     boolList updateFaces(bFaces.size(), false);
 244     # ifdef USE_OMP
 245     # pragma omp parallel for if( updateBndNodes.size() > 1000 )
 246     # endif
 247     forAll(updateBndNodes, i)
 248     {
 249         const label bpI = updateBndNodes[i];
 250         forAllRow(pFaces, bpI, j)
 251             updateFaces[pFaces(bpI, j)] = true;
 252     }
 253
 254     if( surfaceEngine_.faceCentresPtr_ )
 255     {
 256         vectorField& faceCentres = *surfaceEngine_.faceCentresPtr_;
 257
 258         # ifdef USE_OMP
 259         # pragma omp parallel for if( updateFaces.size() > 1000 ) \
 260         schedule(dynamic, 100)
 261         # endif
 262         forAll(updateFaces, bfI)
 263         {
 264             if( updateFaces[bfI] )
 265                 faceCentres[bfI] = bFaces[bfI].centre(points);
 266         }
 267     }
 268
 269     if( surfaceEngine_.faceNormalsPtr_ )
 270     {
 271         vectorField& faceNormals = *surfaceEngine_.faceNormalsPtr_;
 272
 273         # ifdef USE_OMP
 274         # pragma omp parallel for if( updateFaces.size() > 1000 ) \
 275         schedule(dynamic, 100)
 276         # endif
 277         forAll(updateFaces, bfI)
 278         {
 279             if( updateFaces[bfI] )
 280                 faceNormals[bfI] = bFaces[bfI].normal(points);
 281         }
 282     }
 283
 284     if( surfaceEngine_.pointNormalsPtr_ )
 285     {
 286         const vectorField& faceNormals = surfaceEngine_.faceNormals();
 287
 288         boolList updateBndPoint(pFaces.size(), false);
 289         # ifdef USE_OMP
 290         # pragma omp parallel for schedule(dynamic, 50)
 291         # endif
 292         forAll(updateBndNodes, i)
 293         {
 294             const label bpI = updateBndNodes[i];
 295
 296             forAllRow(pFaces, bpI, pfI)
 297             {
 298                 const face& bf = bFaces[pFaces(bpI, pfI)];
 299
 300                 forAll(bf, pI)
 301                     updateBndPoint[bp[bf[pI]]] = true;
 302             }
 303         }
 304
 305         vectorField& pn = *surfaceEngine_.pointNormalsPtr_;
 306         # ifdef USE_OMP
 307         # pragma omp parallel for schedule(dynamic, 100)
 308         # endif
 309         forAll(updateBndPoint, bpI)
 310         {
 311             if( !updateBndPoint[bpI] )
 312                 continue;
 313
 314             vector n(vector::zero);
 315             forAllRow(pFaces, bpI, pfI)
 316                 n += faceNormals[pFaces(bpI, pfI)];
 317
 318             const scalar l = mag(n);
 319             if( l > VSMALL )
 320             {
 321                 n /= l;
 322             }
 323             else
 324             {
 325                 n = vector::zero;
 326             }
 327
 328             pn[bpI] = n;
 329         }
 330
 331         if( Pstream::parRun() )
 332         {
 333             //- update point normals at inter-processor boundaries
 334             const Map<label>& globalToLocal =
 335                 surfaceEngine_.globalToLocalBndPointAddressing();
 336             const VRWGraph& bpAtProcs = surfaceEngine_.bpAtProcs();
 337             const DynList<label>& neiProcs = surfaceEngine_.bpNeiProcs();
 338
 339             //- make sure that the points ar updated on all processors
 340             std::map<label, labelLongList> exchangeNodeLabels;
 341             forAll(neiProcs, i)
 342                 exchangeNodeLabels[neiProcs[i]].clear();
 343
 344             forAllConstIter(Map<label>, globalToLocal, it)
 345             {
 346                 const label bpI = it();
 347
 348                 if( updateBndPoint[bpI] )
 349                 {
 350                     forAllRow(bpAtProcs, bpI, i)
 351                     {
 352                         const label neiProc = bpAtProcs(bpI, i);
 353
 354                         if( neiProc == Pstream::myProcNo() )
 355                             continue;
 356
 357                         exchangeNodeLabels[neiProc].append(it.key());
 358                     }
 359                 }
 360             }
 361
 362             labelLongList receivedNodes;
 363             help::exchangeMap(exchangeNodeLabels, receivedNodes);
 364
 365             forAll(receivedNodes, i)
 366                 updateBndPoint[globalToLocal[receivedNodes[i]]] = true;
 367
 368
 369             //- start updating point normals
 370             std::map<label, LongList<labelledPoint> > exchangeData;
 371             forAll(neiProcs, i)
 372                 exchangeData[neiProcs[i]].clear();
 373
 374             //- prepare data for sending
 375             forAllConstIter(Map<label>, globalToLocal, iter)
 376             {
 377                 const label bpI = iter();
 378
 379                 if( !updateBndPoint[bpI] )
 380                     continue;
 381
 382                 vector& n = pn[bpI];
 383                 n = vector::zero;
 384
 385                 forAllRow(pFaces, bpI, pfI)
 386                     n += faceNormals[pFaces(bpI, pfI)];
 387
 388                 forAllRow(bpAtProcs, bpI, procI)
 389                 {
 390                     const label neiProc = bpAtProcs(bpI, procI);
 391                     if( neiProc == Pstream::myProcNo() )
 392                         continue;
 393
 394                     exchangeData[neiProc].append(labelledPoint(iter.key(), n));
 395                 }
 396             }
 397
 398             //- exchange data with other procs
 399             LongList<labelledPoint> receivedData;
 400             help::exchangeMap(exchangeData, receivedData);
 401
 402             forAll(receivedData, i)
 403             {
 404                 const label bpI = globalToLocal[receivedData[i].pointLabel()];
 405                 pn[bpI] += receivedData[i].coordinates();
 406             }
 407
 408             //- normalize vectors
 409             forAllConstIter(Map<label>, globalToLocal, it)
 410             {
 411                 const label bpI = it();
 412
 413                 if( !updateBndPoint[bpI] )
 414                     continue;
 415
 416                 vector normal = pn[bpI];
 417                 const scalar d = mag(normal);
 418                 if( d > VSMALL )
 419                 {
 420                     normal /= d;
 421                 }
 422                 else
 423                 {
 424                     normal = vector::zero;
 425                 }
 426
 427                 pn[bpI] = normal;
 428             }
 429         }
 430     }
 431 }
 432
 433 void meshSurfaceEngineModifier::updateGeometry()
 434 {
 435     labelLongList updateBndNodes(surfaceEngine_.boundaryPoints().size());
 436
 437     # ifdef USE_OMP
 438     # pragma omp parallel for if( updateBndNodes.size() > 10000 )
 439     # endif
 440     forAll(updateBndNodes, bpI)
 441         updateBndNodes[bpI] = bpI;
 442
 443     updateGeometry(updateBndNodes);
 444 }
 445
 446 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 447
 448 } // End namespace Foam
 449
 450 // ************************************************************************* //