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