src/mesh/cfMesh/utilities/smoothers/geometry/meshOptimizer/meshOptimizerOptimizePointParallel.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 "meshOptimizer.H"
  30 #include "polyMeshGenAddressing.H"
  31
  32 #include "refLabelledPoint.H"
  33 #include "labelledPointScalar.H"
  34 #include "helperFunctionsPar.H"
  35
  36 #include <map>
  37
  38 //#define DEBUGSmooth
  39
  40 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  41
  42 namespace Foam
  43 {
  44
  45 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  46
  47 void meshOptimizer::laplaceSmoother::laplacianParallel
  48 (
  49     const labelLongList& procPoints,
  50     const bool smoothOnlySurfaceNodes
  51 )
  52 {
  53     if( !Pstream::parRun() )
  54         return;
  55
  56     if( returnReduce(procPoints.size(), sumOp<label>()) == 0 )
  57         return;
  58
  59     const polyMeshGenAddressing& addressing = mesh_.addressingData();
  60     const VRWGraph& pPoints = mesh_.addressingData().pointPoints();
  61     pointFieldPMG& points = mesh_.points();
  62
  63     //- exchange data between processors
  64     const labelLongList& globalPointLabel = addressing.globalPointLabel();
  65     const VRWGraph& pointAtProcs = addressing.pointAtProcs();
  66     const Map<label>& globalToLocal = addressing.globalToLocalPointAddressing();
  67
  68     //- create storage for the data
  69     std::map<label, labelledPoint> localData;
  70
  71     //- create data which shall be exchanged with other processors
  72     std::map<label, LongList<refLabelledPoint> > exchangeData;
  73     forAll(procPoints, pI)
  74     {
  75         const label pointI = procPoints[pI];
  76
  77         localData.insert(std::make_pair(pointI, labelledPoint(0,vector::zero)));
  78         labelledPoint& lpd = localData[pointI];
  79
  80         forAllRow(pPoints, pointI, ppI)
  81         {
  82             const label nei = pPoints(pointI, ppI);
  83
  84             if( smoothOnlySurfaceNodes && (vertexLocation_[nei] & INSIDE) )
  85                 continue;
  86
  87             if( pointAtProcs.sizeOfRow(nei) != 0 )
  88             {
  89                 label pMin(Pstream::nProcs());
  90                 forAllRow(pointAtProcs, nei, procI)
  91                 {
  92                     const label procJ = pointAtProcs(nei, procI);
  93                     if( (procJ < pMin) && pointAtProcs.contains(pointI, procJ) )
  94                         pMin = procJ;
  95                 }
  96
  97                 if( pMin != Pstream::myProcNo() )
  98                     continue;
  99             }
 100
 101             ++lpd.pointLabel();
 102             lpd.coordinates() += points[nei];
 103         }
 104
 105         forAllRow(pointAtProcs, pointI, procI)
 106         {
 107             const label neiProc = pointAtProcs(pointI, procI);
 108             if( neiProc == Pstream::myProcNo() )
 109                 continue;
 110
 111             if( exchangeData.find(neiProc) == exchangeData.end() )
 112             {
 113                 exchangeData.insert
 114                 (
 115                     std::make_pair(neiProc, LongList<refLabelledPoint>())
 116                 );
 117             }
 118
 119             //- add data to the list which will be sent to other processor
 120             LongList<refLabelledPoint>& dts = exchangeData[neiProc];
 121             dts.append(refLabelledPoint(globalPointLabel[pointI], lpd));
 122         }
 123     }
 124
 125     //- exchange data with other processors
 126     LongList<refLabelledPoint> receivedData;
 127     help::exchangeMap(exchangeData, receivedData);
 128
 129     forAll(receivedData, i)
 130     {
 131         const refLabelledPoint& lp = receivedData[i];
 132         const label pointI = globalToLocal[lp.objectLabel()];
 133
 134         labelledPoint& lpd = localData[pointI];
 135
 136         lpd.pointLabel() += lp.lPoint().pointLabel();
 137         lpd.coordinates() += lp.lPoint().coordinates();
 138     }
 139
 140     //- create new positions of nodes
 141     forAll(procPoints, pI)
 142     {
 143         const label pointI = procPoints[pI];
 144
 145         const labelledPoint& lp = localData[pointI];
 146
 147         if( lp.pointLabel() != 0 )
 148         {
 149             const point newP = lp.coordinates() / lp.pointLabel();
 150
 151             points[pointI] = newP;
 152         }
 153     }
 154
 155     # ifdef DEBUGSmooth
 156     //- check
 157     std::map<label, LongList<labelledPoint> > check;
 158     forAll(procPoints, pI)
 159     {
 160         const label pointI = procPoints[pI];
 161
 162         forAllRow(pointAtProcs, pointI, i)
 163         {
 164             const label procI = pointAtProcs(pointI, i);
 165             if( procI == Pstream::myProcNo() )
 166                 continue;
 167             if( check.find(procI) == check.end() )
 168             {
 169                 check.insert(std::make_pair(procI, LongList<labelledPoint>()));
 170             }
 171
 172             LongList<labelledPoint>& data = check[procI];
 173             data.append(labelledPoint(globalPointLabel[pointI],points[pointI]));
 174         }
 175     }
 176
 177     LongList<labelledPoint> data;
 178     help::exchangeMap(check, data);
 179
 180     forAll(data, i)
 181     {
 182         const label pointI = globalToLocal[data[i].pointLabel()];
 183
 184         if( mag(points[pointI] - data[i].coordinates()) > SMALL )
 185             Pout << "Crap " << globalPointLabel[pointI] << " coordinates "
 186                 << points[pointI] << " point there " << data[i] << endl;
 187     }
 188     # endif
 189 }
 190
 191 void meshOptimizer::laplaceSmoother::laplacianPCParallel
 192 (
 193     const labelLongList& procPoints
 194 )
 195 {
 196     if( !Pstream::parRun() )
 197         return;
 198
 199     if( returnReduce(procPoints.size(), sumOp<label>()) == 0 )
 200         return;
 201
 202     const polyMeshGenAddressing& addressing = mesh_.addressingData();
 203     const VRWGraph& pCells = mesh_.addressingData().pointCells();
 204     const vectorField& centres = mesh_.addressingData().cellCentres();
 205     pointFieldPMG& points = mesh_.points();
 206
 207     //- exchange data between processors
 208     const labelLongList& globalPointLabel = addressing.globalPointLabel();
 209     const VRWGraph& pointAtProcs = addressing.pointAtProcs();
 210     const Map<label>& globalToLocal = addressing.globalToLocalPointAddressing();
 211
 212     //- create storage for the data
 213     std::map<label, labelledPoint> localData;
 214
 215     //- create data which shall be exchanged with other processors
 216     std::map<label, LongList<refLabelledPoint> > exchangeData;
 217     forAll(procPoints, pI)
 218     {
 219         const label pointI = procPoints[pI];
 220
 221         if( vertexLocation_[pointI] & LOCKED )
 222             continue;
 223
 224         localData.insert(std::make_pair(pointI, labelledPoint(0,vector::zero)));
 225         labelledPoint& lpd = localData[pointI];
 226
 227         forAllRow(pCells, pointI, pcI)
 228         {
 229             const label cellI = pCells(pointI, pcI);
 230
 231             ++lpd.pointLabel();
 232             lpd.coordinates() += centres[cellI];
 233         }
 234
 235         forAllRow(pointAtProcs, pointI, procI)
 236         {
 237             const label neiProc = pointAtProcs(pointI, procI);
 238             if( neiProc == Pstream::myProcNo() )
 239                 continue;
 240
 241             if( exchangeData.find(neiProc) == exchangeData.end() )
 242             {
 243                 exchangeData.insert
 244                 (
 245                     std::make_pair(neiProc, LongList<refLabelledPoint>())
 246                 );
 247             }
 248
 249             //- add data to the list which will be sent to other processor
 250             LongList<refLabelledPoint>& dts = exchangeData[neiProc];
 251             dts.append(refLabelledPoint(globalPointLabel[pointI], lpd));
 252         }
 253     }
 254
 255     //- exchange data with other processors
 256     LongList<refLabelledPoint> receivedData;
 257     help::exchangeMap(exchangeData, receivedData);
 258
 259     forAll(receivedData, i)
 260     {
 261         const refLabelledPoint& lp = receivedData[i];
 262         const label pointI = globalToLocal[lp.objectLabel()];
 263
 264         labelledPoint& lpd = localData[pointI];
 265
 266         lpd.pointLabel() += lp.lPoint().pointLabel();
 267         lpd.coordinates() += lp.lPoint().coordinates();
 268     }
 269
 270     //- create new positions of nodes
 271     forAll(procPoints, pI)
 272     {
 273         const label pointI = procPoints[pI];
 274
 275         const labelledPoint& lp = localData[pointI];
 276
 277         if( lp.pointLabel() != 0 )
 278         {
 279             const point newP = lp.coordinates() / lp.pointLabel();
 280
 281             points[pointI] = newP;
 282         }
 283     }
 284
 285     # ifdef DEBUGSmooth
 286     //- check
 287     std::map<label, LongList<labelledPoint> > check;
 288     forAll(procPoints, pI)
 289     {
 290         const label pointI = procPoints[pI];
 291
 292         forAllRow(pointAtProcs, pointI, i)
 293         {
 294             const label procI = pointAtProcs(pointI, i);
 295             if( procI == Pstream::myProcNo() )
 296                 continue;
 297             if( check.find(procI) == check.end() )
 298             {
 299                 check.insert(std::make_pair(procI, LongList<labelledPoint>()));
 300             }
 301
 302             LongList<labelledPoint>& data = check[procI];
 303             data.append(labelledPoint(globalPointLabel[pointI],points[pointI]));
 304         }
 305     }
 306
 307     LongList<labelledPoint> data;
 308     help::exchangeMap(check, data);
 309
 310     forAll(data, i)
 311     {
 312         const label pointI = globalToLocal[data[i].pointLabel()];
 313
 314         if( mag(points[pointI] - data[i].coordinates()) > SMALL )
 315             Pout << "Crap " << globalPointLabel[pointI] << " coordinates "
 316                 << points[pointI] << " point there " << data[i] << endl;
 317     }
 318     # endif
 319 }
 320
 321 void meshOptimizer::laplaceSmoother::laplacianWPCParallel
 322 (
 323     const labelLongList& procPoints
 324 )
 325 {
 326     if( !Pstream::parRun() )
 327         return;
 328
 329     if( returnReduce(procPoints.size(), sumOp<label>()) == 0 )
 330         return;
 331
 332     const polyMeshGenAddressing& addressing = mesh_.addressingData();
 333     const VRWGraph& pCells = mesh_.addressingData().pointCells();
 334     const vectorField& centres = mesh_.addressingData().cellCentres();
 335     const scalarField& volumes = mesh_.addressingData().cellVolumes();
 336     pointFieldPMG& points = mesh_.points();
 337
 338     //- exchange data between processors
 339     const labelLongList& globalPointLabel = addressing.globalPointLabel();
 340     const VRWGraph& pointAtProcs = addressing.pointAtProcs();
 341     const Map<label>& globalToLocal = addressing.globalToLocalPointAddressing();
 342
 343     //- create storage for the data
 344     std::map<label, labelledPointScalar> localData;
 345
 346     //- create data which shall be exchanged with other processors
 347     std::map<label, LongList<labelledPointScalar> > exchangeData;
 348     forAll(procPoints, pI)
 349     {
 350         const label pointI = procPoints[pI];
 351
 352         if( vertexLocation_[pointI] & LOCKED )
 353             continue;
 354
 355         localData.insert
 356         (
 357             std::make_pair
 358             (
 359                 pointI,
 360                 labelledPointScalar(globalPointLabel[pointI], vector::zero, 0.)
 361             )
 362         );
 363         labelledPointScalar& lps = localData[pointI];
 364
 365         forAllRow(pCells, pointI, pcI)
 366         {
 367             const label cellI = pCells(pointI, pcI);
 368
 369             const scalar w = Foam::max(volumes[cellI], VSMALL);
 370             lps.coordinates() += w * centres[cellI];
 371             lps.scalarValue() += w;
 372         }
 373
 374         forAllRow(pointAtProcs, pointI, procI)
 375         {
 376             const label neiProc = pointAtProcs(pointI, procI);
 377             if( neiProc == Pstream::myProcNo() )
 378                 continue;
 379
 380             if( exchangeData.find(neiProc) == exchangeData.end() )
 381             {
 382                 exchangeData.insert
 383                 (
 384                     std::make_pair(neiProc, LongList<labelledPointScalar>())
 385                 );
 386             }
 387
 388             //- add data to the list which will be sent to other processor
 389             LongList<labelledPointScalar>& dts = exchangeData[neiProc];
 390             dts.append(lps);
 391         }
 392     }
 393
 394     //- exchange data with other processors
 395     LongList<labelledPointScalar> receivedData;
 396     help::exchangeMap(exchangeData, receivedData);
 397
 398     forAll(receivedData, i)
 399     {
 400         const labelledPointScalar& lps = receivedData[i];
 401         const label pointI = globalToLocal[lps.pointLabel()];
 402
 403         labelledPointScalar& lp = localData[pointI];
 404
 405         lp.scalarValue() += lps.scalarValue();
 406         lp.coordinates() += lps.coordinates();
 407     }
 408
 409     //- create new positions of nodes
 410     forAll(procPoints, pI)
 411     {
 412         const label pointI = procPoints[pI];
 413
 414         if( vertexLocation_[pointI] & LOCKED )
 415             continue;
 416
 417         const labelledPointScalar& lp = localData[pointI];
 418
 419         if( lp.pointLabel() != 0 )
 420         {
 421             const point newP = lp.coordinates() / lp.scalarValue();
 422
 423             points[pointI] = newP;
 424         }
 425     }
 426
 427     # ifdef DEBUGSmooth
 428     //- check
 429     std::map<label, LongList<labelledPoint> > check;
 430     forAll(procPoints, pI)
 431     {
 432         const label pointI = procPoints[pI];
 433
 434         forAllRow(pointAtProcs, pointI, i)
 435         {
 436             const label procI = pointAtProcs(pointI, i);
 437             if( procI == Pstream::myProcNo() )
 438                 continue;
 439             if( check.find(procI) == check.end() )
 440             {
 441                 check.insert(std::make_pair(procI, LongList<labelledPoint>()));
 442             }
 443
 444             LongList<labelledPoint>& data = check[procI];
 445             data.append(labelledPoint(globalPointLabel[pointI],points[pointI]));
 446         }
 447     }
 448
 449     LongList<labelledPoint> data;
 450     help::exchangeMap(check, data);
 451
 452     forAll(data, i)
 453     {
 454         const label pointI = globalToLocal[data[i].pointLabel()];
 455
 456         if( mag(points[pointI] - data[i].coordinates()) > SMALL )
 457             Pout << "Crap " << globalPointLabel[pointI] << " coordinates "
 458                 << points[pointI] << " point there " << data[i] << endl;
 459     }
 460     # endif
 461 }
 462
 463 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 464
 465 } // End namespace Foam
 466
 467 // ************************************************************************* //