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