src/mesh/cfMesh/meshLibrary/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         localData.insert(std::make_pair(pointI, labelledPoint(0,vector::zero)));
 222         labelledPoint& lpd = localData[pointI];
 223
 224         forAllRow(pCells, pointI, pcI)
 225         {
 226             const label cellI = pCells(pointI, pcI);
 227
 228             ++lpd.pointLabel();
 229             lpd.coordinates() += centres[cellI];
 230         }
 231
 232         forAllRow(pointAtProcs, pointI, procI)
 233         {
 234             const label neiProc = pointAtProcs(pointI, procI);
 235             if( neiProc == Pstream::myProcNo() )
 236                 continue;
 237
 238             if( exchangeData.find(neiProc) == exchangeData.end() )
 239             {
 240                 exchangeData.insert
 241                 (
 242                     std::make_pair(neiProc, LongList<refLabelledPoint>())
 243                 );
 244             }
 245
 246             //- add data to the list which will be sent to other processor
 247             LongList<refLabelledPoint>& dts = exchangeData[neiProc];
 248             dts.append(refLabelledPoint(globalPointLabel[pointI], lpd));
 249         }
 250     }
 251
 252     //- exchange data with other processors
 253     LongList<refLabelledPoint> receivedData;
 254     help::exchangeMap(exchangeData, receivedData);
 255
 256     forAll(receivedData, i)
 257     {
 258         const refLabelledPoint& lp = receivedData[i];
 259         const label pointI = globalToLocal[lp.objectLabel()];
 260
 261         labelledPoint& lpd = localData[pointI];
 262
 263         lpd.pointLabel() += lp.lPoint().pointLabel();
 264         lpd.coordinates() += lp.lPoint().coordinates();
 265     }
 266
 267     //- create new positions of nodes
 268     forAll(procPoints, pI)
 269     {
 270         const label pointI = procPoints[pI];
 271
 272         const labelledPoint& lp = localData[pointI];
 273
 274         if( lp.pointLabel() != 0 )
 275         {
 276             const point newP = lp.coordinates() / lp.pointLabel();
 277
 278             points[pointI] = newP;
 279         }
 280     }
 281
 282     # ifdef DEBUGSmooth
 283     //- check
 284     std::map<label, LongList<labelledPoint> > check;
 285     forAll(procPoints, pI)
 286     {
 287         const label pointI = procPoints[pI];
 288
 289         forAllRow(pointAtProcs, pointI, i)
 290         {
 291             const label procI = pointAtProcs(pointI, i);
 292             if( procI == Pstream::myProcNo() )
 293                 continue;
 294             if( check.find(procI) == check.end() )
 295             {
 296                 check.insert(std::make_pair(procI, LongList<labelledPoint>()));
 297             }
 298
 299             LongList<labelledPoint>& data = check[procI];
 300             data.append(labelledPoint(globalPointLabel[pointI],points[pointI]));
 301         }
 302     }
 303
 304     LongList<labelledPoint> data;
 305     help::exchangeMap(check, data);
 306
 307     forAll(data, i)
 308     {
 309         const label pointI = globalToLocal[data[i].pointLabel()];
 310
 311         if( mag(points[pointI] - data[i].coordinates()) > SMALL )
 312             Pout << "Crap " << globalPointLabel[pointI] << " coordinates "
 313                 << points[pointI] << " point there " << data[i] << endl;
 314     }
 315     # endif
 316 }
 317
 318 void meshOptimizer::laplaceSmoother::laplacianWPCParallel
 319 (
 320     const labelLongList& procPoints
 321 )
 322 {
 323     if( !Pstream::parRun() )
 324         return;
 325
 326     if( returnReduce(procPoints.size(), sumOp<label>()) == 0 )
 327         return;
 328
 329     const polyMeshGenAddressing& addressing = mesh_.addressingData();
 330     const VRWGraph& pCells = mesh_.addressingData().pointCells();
 331     const vectorField& centres = mesh_.addressingData().cellCentres();
 332     const scalarField& volumes = mesh_.addressingData().cellVolumes();
 333     pointFieldPMG& points = mesh_.points();
 334
 335     //- exchange data between processors
 336     const labelLongList& globalPointLabel = addressing.globalPointLabel();
 337     const VRWGraph& pointAtProcs = addressing.pointAtProcs();
 338     const Map<label>& globalToLocal = addressing.globalToLocalPointAddressing();
 339
 340     //- create storage for the data
 341     std::map<label, labelledPointScalar> localData;
 342
 343     //- create data which shall be exchanged with other processors
 344     std::map<label, LongList<labelledPointScalar> > exchangeData;
 345     forAll(procPoints, pI)
 346     {
 347         const label pointI = procPoints[pI];
 348
 349         localData.insert
 350         (
 351             std::make_pair
 352             (
 353                 pointI,
 354                 labelledPointScalar(globalPointLabel[pointI], vector::zero, 0.)
 355             )
 356         );
 357         labelledPointScalar& lps = localData[pointI];
 358
 359         forAllRow(pCells, pointI, pcI)
 360         {
 361             const label cellI = pCells(pointI, pcI);
 362
 363             const scalar w = Foam::max(volumes[cellI], VSMALL);
 364             lps.coordinates() += w * centres[cellI];
 365             lps.scalarValue() += w;
 366         }
 367
 368         forAllRow(pointAtProcs, pointI, procI)
 369         {
 370             const label neiProc = pointAtProcs(pointI, procI);
 371             if( neiProc == Pstream::myProcNo() )
 372                 continue;
 373
 374             if( exchangeData.find(neiProc) == exchangeData.end() )
 375             {
 376                 exchangeData.insert
 377                 (
 378                     std::make_pair(neiProc, LongList<labelledPointScalar>())
 379                 );
 380             }
 381
 382             //- add data to the list which will be sent to other processor
 383             LongList<labelledPointScalar>& dts = exchangeData[neiProc];
 384             dts.append(lps);
 385         }
 386     }
 387
 388     //- exchange data with other processors
 389     LongList<labelledPointScalar> receivedData;
 390     help::exchangeMap(exchangeData, receivedData);
 391
 392     forAll(receivedData, i)
 393     {
 394         const labelledPointScalar& lps = receivedData[i];
 395         const label pointI = globalToLocal[lps.pointLabel()];
 396
 397         labelledPointScalar& lp = localData[pointI];
 398
 399         lp.scalarValue() += lps.scalarValue();
 400         lp.coordinates() += lps.coordinates();
 401     }
 402
 403     //- create new positions of nodes
 404     forAll(procPoints, pI)
 405     {
 406         const label pointI = procPoints[pI];
 407
 408         const labelledPointScalar& lp = localData[pointI];
 409
 410         if( lp.pointLabel() != 0 )
 411         {
 412             const point newP = lp.coordinates() / lp.scalarValue();
 413
 414             points[pointI] = newP;
 415         }
 416     }
 417
 418     # ifdef DEBUGSmooth
 419     //- check
 420     std::map<label, LongList<labelledPoint> > check;
 421     forAll(procPoints, pI)
 422     {
 423         const label pointI = procPoints[pI];
 424
 425         forAllRow(pointAtProcs, pointI, i)
 426         {
 427             const label procI = pointAtProcs(pointI, i);
 428             if( procI == Pstream::myProcNo() )
 429                 continue;
 430             if( check.find(procI) == check.end() )
 431             {
 432                 check.insert(std::make_pair(procI, LongList<labelledPoint>()));
 433             }
 434
 435             LongList<labelledPoint>& data = check[procI];
 436             data.append(labelledPoint(globalPointLabel[pointI],points[pointI]));
 437         }
 438     }
 439
 440     LongList<labelledPoint> data;
 441     help::exchangeMap(check, data);
 442
 443     forAll(data, i)
 444     {
 445         const label pointI = globalToLocal[data[i].pointLabel()];
 446
 447         if( mag(points[pointI] - data[i].coordinates()) > SMALL )
 448             Pout << "Crap " << globalPointLabel[pointI] << " coordinates "
 449                 << points[pointI] << " point there " << data[i] << endl;
 450     }
 451     # endif
 452 }
 453
 454 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 455
 456 } // End namespace Foam
 457
 458 // ************************************************************************* //