src/mesh/cfMesh/utilities/smoothers/geometry/meshOptimizer/meshOptimizerOptimizePoint.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 #include "meshSurfaceEngine.H"
  32
  33 # ifdef USE_OMP
  34 #include <omp.h>
  35 # endif
  36
  37 //#define DEBUGSmooth
  38
  39 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  40
  41 namespace Foam
  42 {
  43
  44 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  45
  46 void meshOptimizer::laplaceSmoother::laplacian
  47 (
  48     const labelLongList& smoothPoints,
  49     const label nIterations
  50 )
  51 {
  52     const VRWGraph& pPoints = mesh_.addressingData().pointPoints();
  53     pointFieldPMG& points = mesh_.points();
  54
  55     for(label iterationI=0;iterationI<nIterations;++iterationI)
  56     {
  57         labelLongList procPoints;
  58
  59         forAll(smoothPoints, i)
  60         {
  61             const label pointI = smoothPoints[i];
  62
  63             if( vertexLocation_[pointI] & LOCKED )
  64                 continue;
  65
  66             if( vertexLocation_[pointI] & PARALLELBOUNDARY )
  67             {
  68                 procPoints.append(pointI);
  69
  70                 continue;
  71             }
  72
  73             vector newP(vector::zero);
  74
  75             const label nPointPoints = pPoints.sizeOfRow(pointI);
  76
  77             if( nPointPoints == 0 )
  78                 return;
  79
  80             for(label pI=0;pI<nPointPoints;++pI)
  81                 newP += points[pPoints(pointI, pI)];
  82
  83             newP /= pPoints.sizeOfRow(pointI);
  84             points[pointI] = newP;
  85         }
  86
  87         laplacianParallel(procPoints, false);
  88     }
  89
  90     updateMeshGeometry(smoothPoints);
  91 }
  92
  93 void meshOptimizer::laplaceSmoother::laplacianSurface
  94 (
  95     const labelLongList& smoothPoints,
  96     const label nIterations
  97 )
  98 {
  99     const VRWGraph& pPoints = mesh_.addressingData().pointPoints();
 100     pointFieldPMG& points = mesh_.points();
 101
 102     for(label iterationI=0;iterationI<nIterations;++iterationI)
 103     {
 104         labelLongList procPoints;
 105
 106         forAll(smoothPoints, i)
 107         {
 108             const label pointI = smoothPoints[i];
 109
 110             if( vertexLocation_[pointI] & LOCKED )
 111                 continue;
 112
 113             if( vertexLocation_[pointI] & PARALLELBOUNDARY )
 114             {
 115                 procPoints.append(pointI);
 116
 117                 continue;
 118             }
 119
 120             vector newP(vector::zero);
 121
 122             label counter(0);
 123             forAllRow(pPoints, pointI, pI)
 124             {
 125                 const label pLabel = pPoints(pointI, pI);
 126                 if( vertexLocation_[pLabel] & INSIDE )
 127                     continue;
 128
 129                 newP += points[pLabel];
 130                 ++counter;
 131             }
 132
 133             if( counter != 0 )
 134             {
 135                 newP /= counter;
 136                 points[pointI] = newP;
 137             }
 138         }
 139
 140         laplacianParallel(smoothPoints, true);
 141     }
 142
 143     updateMeshGeometry(smoothPoints);
 144 }
 145
 146 void meshOptimizer::laplaceSmoother::laplacianPC
 147 (
 148     const labelLongList& smoothPoints,
 149     const label nIterations
 150 )
 151 {
 152     const VRWGraph& pointCells = mesh_.addressingData().pointCells();
 153     const vectorField& centres = mesh_.addressingData().cellCentres();
 154     pointFieldPMG& points = mesh_.points();
 155
 156     for(label iterationI=0;iterationI<nIterations;++iterationI)
 157     {
 158         labelLongList procPoints;
 159
 160         # ifdef USE_OMP
 161         # pragma omp parallel for schedule(dynamic, 20)
 162         # endif
 163         forAll(smoothPoints, i)
 164         {
 165             const label pointI = smoothPoints[i];
 166
 167             if( vertexLocation_[pointI] & LOCKED )
 168                 continue;
 169
 170             if( pointCells.sizeOfRow(pointI) == 0 )
 171                 continue;
 172
 173             if( vertexLocation_[pointI] & PARALLELBOUNDARY )
 174             {
 175                 # ifdef USE_OMP
 176                 # pragma omp critical
 177                 # endif
 178                 procPoints.append(pointI);
 179
 180                 continue;
 181             }
 182
 183             point newP(vector::zero);
 184             forAllRow(pointCells, pointI, pcI)
 185                 newP += centres[pointCells(pointI, pcI)];
 186
 187             newP /= pointCells.sizeOfRow(pointI);
 188
 189             points[pointI] = newP;
 190         }
 191
 192         laplacianPCParallel(procPoints);
 193
 194         updateMeshGeometry(smoothPoints);
 195     }
 196 }
 197
 198 void meshOptimizer::laplaceSmoother::laplacianWPC
 199 (
 200     const labelLongList& smoothPoints,
 201     const label nIterations
 202 )
 203 {
 204     const VRWGraph& pointCells = mesh_.addressingData().pointCells();
 205     const vectorField& centres = mesh_.addressingData().cellCentres();
 206     const scalarField& volumes = mesh_.addressingData().cellVolumes();
 207
 208     pointFieldPMG& points = mesh_.points();
 209
 210     for(label iterationI=0;iterationI<nIterations;++iterationI)
 211     {
 212         labelLongList procPoints;
 213
 214         # ifdef USE_OMP
 215         # pragma omp parallel for schedule(dynamic, 20)
 216         # endif
 217         forAll(smoothPoints, i)
 218         {
 219             const label pointI = smoothPoints[i];
 220
 221             if( vertexLocation_[pointI] & LOCKED )
 222                 continue;
 223
 224             if( pointCells.sizeOfRow(pointI) == 0 )
 225                 continue;
 226
 227             if( vertexLocation_[pointI] & PARALLELBOUNDARY )
 228             {
 229                 # ifdef USE_OMP
 230                 # pragma omp critical
 231                 # endif
 232                 procPoints.append(pointI);
 233
 234                 continue;
 235             }
 236
 237             point newP(vector::zero);
 238             scalar sumWeights(0.0);
 239             forAllRow(pointCells, pointI, pcI)
 240             {
 241                 const label cellI = pointCells(pointI, pcI);
 242                 const scalar w = Foam::max(volumes[cellI], VSMALL);
 243                 newP += w * centres[cellI];
 244                 sumWeights += w;
 245             }
 246
 247             newP /= sumWeights;
 248             points[pointI] = newP;
 249         }
 250
 251         laplacianWPCParallel(procPoints);
 252
 253         updateMeshGeometry(smoothPoints);
 254     }
 255 }
 256
 257 void meshOptimizer::laplaceSmoother::updateMeshGeometry
 258 (
 259     const labelLongList& smoothPoints
 260 )
 261 {
 262     const cellListPMG& cells = mesh_.cells();
 263     const VRWGraph& pointCells = mesh_.addressingData().pointCells();
 264
 265     boolList chF(mesh_.faces().size(), false);
 266
 267     # ifdef USE_OMP
 268     # pragma omp parallel for if( smoothPoints.size() > 100 ) \
 269     schedule(dynamic, 20)
 270     # endif
 271     forAll(smoothPoints, i)
 272     {
 273         const label pointI = smoothPoints[i];
 274
 275         if( vertexLocation_[pointI] & LOCKED )
 276             continue;
 277
 278         forAllRow(pointCells, pointI, pcI)
 279         {
 280             const cell& c = cells[pointCells(pointI, pcI)];
 281
 282             forAll(c, fI)
 283                 chF[c[fI]] = true;
 284         }
 285     }
 286
 287     //- make sure that neighbouring processors get the same information
 288     const PtrList<processorBoundaryPatch>& pBnd = mesh_.procBoundaries();
 289     forAll(pBnd, patchI)
 290     {
 291         const label start = pBnd[patchI].patchStart();
 292         const label size = pBnd[patchI].patchSize();
 293
 294         labelLongList sendData;
 295         for(label faceI=0;faceI<size;++faceI)
 296         {
 297             if( chF[start+faceI] )
 298                 sendData.append(faceI);
 299         }
 300
 301         OPstream toOtherProc
 302         (
 303             Pstream::blocking,
 304             pBnd[patchI].neiProcNo(),
 305             sendData.byteSize()
 306         );
 307
 308         toOtherProc << sendData;
 309     }
 310
 311     forAll(pBnd, patchI)
 312     {
 313         labelList receivedData;
 314
 315         IPstream fromOtherProc
 316         (
 317             Pstream::blocking,
 318             pBnd[patchI].neiProcNo()
 319         );
 320
 321         fromOtherProc >> receivedData;
 322
 323         const label start = pBnd[patchI].patchStart();
 324         forAll(receivedData, i)
 325             chF[start+receivedData[i]] = true;
 326     }
 327
 328     //- update geometry information
 329     const_cast<polyMeshGenAddressing&>
 330     (
 331         mesh_.addressingData()
 332     ).updateGeometry(chF);
 333 }
 334
 335 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 336 // Constructor of laplaceSmoother
 337
 338 meshOptimizer::laplaceSmoother::laplaceSmoother
 339 (
 340     polyMeshGen& mesh,
 341     const List<direction>& vertexLocation
 342 )
 343 :
 344     mesh_(mesh),
 345     vertexLocation_(vertexLocation)
 346 {}
 347
 348 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 349
 350 meshOptimizer::laplaceSmoother::~laplaceSmoother()
 351 {}
 352
 353 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 354 // Member Functions
 355
 356 void meshOptimizer::laplaceSmoother::optimizeLaplacian(const label nIterations)
 357 {
 358     labelLongList smoothPoints;
 359
 360     forAll(vertexLocation_, pointI)
 361     {
 362         if( vertexLocation_[pointI] & INSIDE )
 363             smoothPoints.append(pointI);
 364     }
 365
 366     laplacian(smoothPoints, nIterations);
 367 }
 368
 369 void meshOptimizer::laplaceSmoother::optimizeLaplacian
 370 (
 371     const labelHashSet& badFaces,
 372     const label nIterations
 373 )
 374 {
 375     FatalError << "Not implemented " << exit(FatalError);
 376 }
 377
 378 void meshOptimizer::laplaceSmoother::optimizeSurfaceLaplacian
 379 (
 380     const labelHashSet& badFaces,
 381     const label nIterations
 382 )
 383 {
 384     FatalError << "Not implemented " << exit(FatalError);
 385 }
 386
 387 void meshOptimizer::laplaceSmoother::optimizeLaplacianPC
 388 (
 389     const label nIterations
 390 )
 391 {
 392     labelLongList smoothPoints;
 393
 394     forAll(vertexLocation_, pointI)
 395     {
 396         if( vertexLocation_[pointI] & INSIDE )
 397             smoothPoints.append(pointI);
 398     }
 399
 400     laplacianPC(smoothPoints, nIterations);
 401 }
 402
 403 void meshOptimizer::laplaceSmoother::optimizeLaplacianPC
 404 (
 405     const labelHashSet& badFaces,
 406     const label nIterations
 407 )
 408 {
 409     FatalError << "Not implemented " << exit(FatalError);
 410 }
 411
 412 void meshOptimizer::laplaceSmoother::optimizeLaplacianWPC
 413 (
 414     const label nIterations
 415 )
 416 {
 417     labelLongList smoothPoints;
 418
 419     forAll(vertexLocation_, pointI)
 420     {
 421         if( vertexLocation_[pointI] & INSIDE )
 422             smoothPoints.append(pointI);
 423     }
 424
 425     laplacianWPC(smoothPoints, nIterations);
 426 }
 427
 428 void meshOptimizer::laplaceSmoother::optimizeLaplacianWPC
 429 (
 430     const labelHashSet& badFaces,
 431     const label nIterations
 432 )
 433 {
 434     FatalError << "Not implemented " << exit(FatalError);
 435 }
 436
 437 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 438
 439 } // End namespace Foam
 440
 441 // ************************************************************************* //