src/mesh/cfMesh/utilities/smoothers/geometry/meshSurfaceOptimizer/meshSurfaceOptimizerOptimizePointParallel.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 "demandDrivenData.H"
  29 #include "meshSurfaceOptimizer.H"
  30 #include "meshSurfaceEngineModifier.H"
  31 #include "plane.H"
  32 #include "meshSurfaceMapper.H"
  33 #include "surfaceOptimizer.H"
  34 #include "refLabelledPoint.H"
  35 #include "helperFunctions.H"
  36
  37 //#define DEBUGSearch
  38
  39 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  40
  41 namespace Foam
  42 {
  43
  44 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  45
  46 void meshSurfaceOptimizer::nodeDisplacementLaplacianParallel
  47 (
  48     const labelLongList& nodesToSmooth,
  49     const bool transformIntoPlane
  50 )
  51 {
  52     if( returnReduce(nodesToSmooth.size(), sumOp<label>()) == 0 )
  53         return;
  54
  55     //- create storage for data
  56     std::map<label, labelledPoint> localData;
  57
  58     //- exchange data with other processors
  59     std::map<label, LongList<refLabelledPoint> > exchangeData;
  60
  61     const pointField& points = surfaceEngine_.points();
  62     const labelList& bPoints = surfaceEngine_.boundaryPoints();
  63     const VRWGraph& pPoints = surfaceEngine_.pointPoints();
  64     const vectorField& pNormals = surfaceEngine_.pointNormals();
  65     const labelList& globalPointLabel =
  66         surfaceEngine_.globalBoundaryPointLabel();
  67     const VRWGraph& bpAtProcs = surfaceEngine_.bpAtProcs();
  68     const Map<label>& globalToLocal =
  69         surfaceEngine_.globalToLocalBndPointAddressing();
  70
  71     //- perform smoothing
  72     forAll(nodesToSmooth, pI)
  73     {
  74         const label bpI = nodesToSmooth[pI];
  75
  76         if( vertexType_[bpI] & LOCKED )
  77             continue;
  78
  79         if( magSqr(pNormals[bpI]) < VSMALL )
  80             continue;
  81
  82         const plane pl(points[bPoints[bpI]], pNormals[bpI]);
  83
  84         //- project points onto the plane
  85         localData.insert(std::make_pair(bpI, labelledPoint(0, vector::zero)));
  86         labelledPoint& lpd = localData[bpI];
  87
  88         forAllRow(pPoints, bpI, ppI)
  89         {
  90             const label nei = pPoints(bpI, ppI);
  91
  92             if( bpAtProcs.sizeOfRow(nei) != 0 )
  93             {
  94                 label pMin(Pstream::nProcs());
  95                 forAllRow(bpAtProcs, nei, procI)
  96                 {
  97                     const label procJ = bpAtProcs(nei, procI);
  98                     if( (procJ < pMin) && bpAtProcs.contains(bpI, procJ) )
  99                         pMin = procJ;
 100                 }
 101
 102                 if( pMin != Pstream::myProcNo() )
 103                     continue;
 104             }
 105
 106             const point& p = points[bPoints[nei]];
 107             ++lpd.pointLabel();
 108             lpd.coordinates() += transformIntoPlane?pl.nearestPoint(p):p;
 109         }
 110
 111         forAllRow(bpAtProcs, bpI, procI)
 112         {
 113             const label neiProc = bpAtProcs(bpI, procI);
 114             if( neiProc == Pstream::myProcNo() )
 115                 continue;
 116
 117             if( exchangeData.find(neiProc) == exchangeData.end() )
 118             {
 119                 exchangeData.insert
 120                 (
 121                     std::make_pair(neiProc, LongList<refLabelledPoint>())
 122                 );
 123             }
 124
 125             //- add data to the list which will be sent to other processor
 126             LongList<refLabelledPoint>& dts = exchangeData[neiProc];
 127             dts.append(refLabelledPoint(globalPointLabel[bpI], lpd));
 128         }
 129     }
 130
 131     //- exchange data with other processors
 132     LongList<refLabelledPoint> receivedData;
 133     help::exchangeMap(exchangeData, receivedData);
 134
 135     forAll(receivedData, i)
 136     {
 137         const refLabelledPoint& lp = receivedData[i];
 138         const label bpI = globalToLocal[lp.objectLabel()];
 139
 140         labelledPoint& lpd = localData[bpI];
 141
 142         lpd.pointLabel() += lp.lPoint().pointLabel();
 143         lpd.coordinates() += lp.lPoint().coordinates();
 144     }
 145
 146     forAll(nodesToSmooth, pI)
 147     {
 148         const label bpI = nodesToSmooth[pI];
 149
 150         if( localData.find(bpI) == localData.end() )
 151             continue;
 152
 153         //- create new point position
 154         const labelledPoint& lp = localData[bpI];
 155         const point newP = lp.coordinates() / lp.pointLabel();
 156
 157         meshSurfaceEngineModifier surfaceModifier(surfaceEngine_);
 158         surfaceModifier.moveBoundaryVertexNoUpdate
 159         (
 160             bpI,
 161             newP
 162         );
 163     }
 164 }
 165
 166 void meshSurfaceOptimizer::nodeDisplacementLaplacianFCParallel
 167 (
 168     const labelLongList& nodesToSmooth,
 169     const bool transformIntoPlane
 170 )
 171 {
 172     if( returnReduce(nodesToSmooth.size(), sumOp<label>()) == 0 )
 173         return;
 174
 175     //- create storage for data
 176     std::map<label, labelledPoint> localData;
 177
 178     //- exchange data with other processors
 179     std::map<label, LongList<refLabelledPoint> > exchangeData;
 180
 181     const pointField& points = surfaceEngine_.points();
 182     const labelList& bPoints = surfaceEngine_.boundaryPoints();
 183     const VRWGraph& pFaces = surfaceEngine_.pointFaces();
 184     const vectorField& faceCentres = surfaceEngine_.faceCentres();
 185     const vectorField& pNormals = surfaceEngine_.pointNormals();
 186
 187     const labelList& globalPointLabel =
 188         surfaceEngine_.globalBoundaryPointLabel();
 189     const VRWGraph& bpAtProcs = surfaceEngine_.bpAtProcs();
 190     const Map<label>& globalToLocal =
 191         surfaceEngine_.globalToLocalBndPointAddressing();
 192
 193     //- perform smoothing
 194     forAll(nodesToSmooth, pI)
 195     {
 196         const label bpI = nodesToSmooth[pI];
 197
 198         if( vertexType_[bpI] & LOCKED )
 199             continue;
 200
 201         if( magSqr(pNormals[bpI]) < VSMALL )
 202             continue;
 203
 204         const plane pl(points[bPoints[bpI]], pNormals[bpI]);
 205
 206         //- project points onto the plane
 207         localData.insert(std::make_pair(bpI, labelledPoint(0, vector::zero)));
 208         labelledPoint& lpd = localData[bpI];
 209
 210         forAllRow(pFaces, bpI, pfI)
 211         {
 212             const point& p = faceCentres[pFaces(bpI, pfI)];
 213             ++lpd.pointLabel();
 214             lpd.coordinates() += transformIntoPlane?pl.nearestPoint(p):p;
 215         }
 216
 217         forAllRow(bpAtProcs, bpI, procI)
 218         {
 219             const label neiProc = bpAtProcs(bpI, procI);
 220             if( neiProc == Pstream::myProcNo() )
 221                 continue;
 222
 223             if( exchangeData.find(neiProc) == exchangeData.end() )
 224             {
 225                 exchangeData.insert
 226                 (
 227                     std::make_pair(neiProc, LongList<refLabelledPoint>())
 228                 );
 229             }
 230
 231             //- add data to the list which will be sent to other processor
 232             LongList<refLabelledPoint>& dts = exchangeData[neiProc];
 233             dts.append(refLabelledPoint(globalPointLabel[bpI], lpd));
 234         }
 235     }
 236
 237     //- exchange data with other processors
 238     LongList<refLabelledPoint> receivedData;
 239     help::exchangeMap(exchangeData, receivedData);
 240
 241     forAll(receivedData, i)
 242     {
 243         const refLabelledPoint& lp = receivedData[i];
 244         const label bpI = globalToLocal[lp.objectLabel()];
 245
 246         labelledPoint& lpd = localData[bpI];
 247
 248         lpd.pointLabel() += lp.lPoint().pointLabel();
 249         lpd.coordinates() += lp.lPoint().coordinates();
 250     }
 251
 252     pointField newPositions(nodesToSmooth.size());
 253
 254     # ifdef USE_OMP
 255     # pragma omp parallel for schedule(dynamic, 20)
 256     # endif
 257     forAll(nodesToSmooth, pI)
 258     {
 259         const label bpI = nodesToSmooth[pI];
 260
 261         if( localData.find(bpI) == localData.end() )
 262         {
 263             newPositions[pI] = points[bPoints[bpI]];
 264             continue;
 265         }
 266
 267         //- create new point position
 268         const labelledPoint& lp = localData[bpI];
 269         const point newP = lp.coordinates() / lp.pointLabel();
 270
 271         newPositions[pI] = newP;
 272     }
 273
 274     meshSurfaceEngineModifier surfaceModifier(surfaceEngine_);
 275     # ifdef USE_OMP
 276     # pragma omp parallel for schedule(dynamic, 20)
 277     # endif
 278     forAll(newPositions, pI)
 279     {
 280         surfaceModifier.moveBoundaryVertexNoUpdate
 281         (
 282             nodesToSmooth[pI],
 283             newPositions[pI]
 284         );
 285     }
 286 }
 287
 288 void meshSurfaceOptimizer::edgeNodeDisplacementParallel
 289 (
 290     const labelLongList& nodesToSmooth
 291 )
 292 {
 293     std::map<label, DynList<labelledPoint, 2> > mPts;
 294
 295     //- insert labels into the map
 296     const labelList& globalPointLabel =
 297         surfaceEngine_.globalBoundaryPointLabel();
 298
 299     forAll(nodesToSmooth, nI)
 300     {
 301         mPts.insert
 302         (
 303             std::make_pair
 304             (
 305                 globalPointLabel[nodesToSmooth[nI]],
 306                 DynList<labelledPoint, 2>()
 307             )
 308         );
 309     }
 310
 311     //- add local edge points
 312     const pointField& points = surfaceEngine_.points();
 313     const labelList& bPoints = surfaceEngine_.boundaryPoints();
 314     const edgeList& edges = surfaceEngine_.edges();
 315     const VRWGraph& bpEdges = surfaceEngine_.boundaryPointEdges();
 316     const  labelList& bp = surfaceEngine_.bp();
 317
 318     forAll(nodesToSmooth, nI)
 319     {
 320         const label bpI = nodesToSmooth[nI];
 321
 322         if( vertexType_[bpI] & LOCKED)
 323             continue;
 324
 325         DynList<labelledPoint, 2>& neiPoints = mPts[globalPointLabel[bpI]];
 326
 327         forAllRow(bpEdges, bpI, epI)
 328         {
 329             const edge& e = edges[bpEdges(bpI, epI)];
 330             const label pI = bp[e.otherVertex(bPoints[bpI])];
 331             if( vertexType_[pI] & (EDGE+CORNER) )
 332             {
 333                 neiPoints.append
 334                 (
 335                     labelledPoint(globalPointLabel[pI], points[bPoints[pI]])
 336                 );
 337             }
 338         }
 339     }
 340
 341     //- start preparing data which can be exchanged with other processors
 342     const DynList<label>& neiProcs = surfaceEngine_.bpNeiProcs();
 343     const VRWGraph& bpAtProcs = surfaceEngine_.bpAtProcs();
 344
 345     std::map<label, LongList<refLabelledPoint> > mProcs;
 346     forAll(neiProcs, procI)
 347     {
 348         const label neiProc = neiProcs[procI];
 349
 350         if( neiProc == Pstream::myProcNo() )
 351             continue;
 352
 353         mProcs.insert(std::make_pair(neiProc, LongList<refLabelledPoint>()));
 354     }
 355
 356     forAll(nodesToSmooth, nI)
 357     {
 358         const label bpI = nodesToSmooth[nI];
 359
 360         const DynList<labelledPoint, 2>& neiPoints =
 361             mPts[globalPointLabel[bpI]];
 362
 363         forAllRow(bpAtProcs, bpI, procI)
 364         {
 365             const label neiProc = bpAtProcs(bpI, procI);
 366
 367             if( neiProc == Pstream::myProcNo() )
 368                 continue;
 369
 370             forAll(neiPoints, npI)
 371             {
 372                 LongList<refLabelledPoint>& neiProcPts = mProcs[neiProc];
 373                 neiProcPts.append
 374                 (
 375                     refLabelledPoint(globalPointLabel[bpI], neiPoints[npI])
 376                 );
 377             }
 378         }
 379     }
 380
 381     //- exchange data with other processors
 382     LongList<refLabelledPoint> receivedData;
 383     help::exchangeMap(mProcs, receivedData);
 384
 385     forAll(receivedData, prI)
 386     {
 387         const refLabelledPoint& lp = receivedData[prI];
 388         DynList<labelledPoint, 2>& lPts = mPts[lp.objectLabel()];
 389         lPts.appendIfNotIn(receivedData[prI].lPoint());
 390     }
 391
 392     //- Finally, the data is ready to start smoothing
 393     meshSurfaceEngineModifier sm(surfaceEngine_);
 394
 395     pointField newPositions(nodesToSmooth.size());
 396     # ifdef USE_OMP
 397     # pragma omp parallel for schedule(dynamic, 20)
 398     # endif
 399     forAll(nodesToSmooth, pI)
 400     {
 401         const label bpI = nodesToSmooth[pI];
 402
 403         const DynList<labelledPoint, 2>& nPts = mPts[globalPointLabel[bpI]];
 404         point newP(vector::zero);
 405         forAll(nPts, ppI)
 406             newP += nPts[ppI].coordinates();
 407
 408         if( nPts.size() == 2 )
 409         {
 410             newP /= nPts.size();
 411             newPositions[pI] = newP;
 412         }
 413         else
 414         {
 415             newPositions[pI] = points[bPoints[bpI]];
 416         }
 417     }
 418
 419     # ifdef USE_OMP
 420     # pragma omp parallel for schedule(dynamic, 20)
 421     # endif
 422     forAll(newPositions, pI)
 423         sm.moveBoundaryVertexNoUpdate(nodesToSmooth[pI], newPositions[pI]);
 424 }
 425
 426 void meshSurfaceOptimizer::exchangeData
 427 (
 428     const labelLongList& nodesToSmooth,
 429     std::map<label, DynList<parTriFace> >& m
 430 ) const
 431 {
 432     /*
 433     const VRWGraph& bpAtProcs = surfaceEngine_.bpAtProcs();
 434     const DynList<label>& neiProcs = surfaceEngine_.bpNeiProcs();
 435
 436     std::map<label, LongList<parTriFace> > shareData;
 437     forAll(neiProcs, procI)
 438     {
 439         const label neiProc = neiProcs[procI];
 440
 441         if( neiProc == Pstream::myProcNo() )
 442             continue;
 443
 444         shareData.insert(std::make_pair(neiProc, LongList<parTriFace>()));
 445     }
 446
 447     //- create data which will be sent to other processors
 448     const pointField& points = surfaceEngine_.points();
 449     const labelList& globalPointLabel =
 450         surfaceEngine_.globalBoundaryPointLabel();
 451     const labelList& bPoints = surfaceEngine_.boundaryPoints();
 452     const LongList<triFace>& triangles = this->triangles();
 453     const VRWGraph& pTriangles = pointTriangles();
 454
 455     m.clear();
 456     std::map<label, DynList<parTriFace> >::iterator pIter;
 457     forAll(nodesToSmooth, nI)
 458     {
 459         const label bpI = nodesToSmooth[nI];
 460
 461         pIter = m.find(globalPointLabel[bpI]);
 462         if( pIter == m.end() )
 463         {
 464             m.insert
 465             (
 466                 std::make_pair(globalPointLabel[bpI], DynList<parTriFace>())
 467             );
 468
 469             pIter = m.find(globalPointLabel[bpI]);
 470         }
 471
 472         forAll(pTriangles[bpI], ptI)
 473         {
 474             const triFace& tri = triangles[pTriangles[bpI][ptI]];
 475
 476             parTriFace ptf
 477             (
 478                 globalPointLabel[tri[0]],
 479                 globalPointLabel[tri[1]],
 480                 globalPointLabel[tri[2]],
 481                 triangle<point, point>
 482                 (
 483                     points[bPoints[tri[0]]],
 484                     points[bPoints[tri[1]]],
 485                     points[bPoints[tri[2]]]
 486                 )
 487             );
 488
 489             pIter->second.append(ptf);
 490         }
 491
 492         forAllRow(bpAtProcs, bpI, procI)
 493         {
 494             const label neiProc = bpAtProcs(bpI, procI);
 495             if( neiProc == Pstream::myProcNo() )
 496                 continue;
 497
 498             LongList<parTriFace>& shData = shareData[neiProc];
 499
 500             const DynList<parTriFace>& localData = pIter->second;
 501
 502             forAll(localData, i)
 503                 shData.append(localData[i]);
 504         }
 505     }
 506
 507     //- send data to other processors
 508     LongList<parTriFace> receivedData;
 509     help::exchangeMap(shareData, receivedData);
 510
 511     forAll(receivedData, tI)
 512     {
 513         const label gpI = receivedData[tI].globalLabelOfPoint(0);
 514
 515         pIter = m.find(gpI);
 516         if( pIter == m.end() )
 517         {
 518             FatalErrorIn
 519             (
 520                 "void meshSurfaceOptimizer::exchangeData("
 521                 "const labelLongList& nodesToSmooth,"
 522                 "map<label, DynList<parTriFace> >& m"
 523                 ") const"
 524             ) << "Unknown point " << gpI << abort(FatalError);
 525         }
 526
 527         pIter->second.append(receivedData[tI]);
 528     }
 529     */
 530 }
 531
 532 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 533
 534 } // End namespace Foam
 535
 536 // ************************************************************************* //