src/mesh/cfMesh/utilities/smoothers/geometry/meshSurfaceOptimizer/meshSurfaceOptimizerOptimizeSurface.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 "meshSurfaceCheckInvertedVertices.H"
  32 #include "meshOctree.H"
  33 #include "triangle.H"
  34 #include "helperFunctionsPar.H"
  35 #include "meshSurfaceMapper.H"
  36 #include "meshSurfaceMapper2D.H"
  37 #include "polyMeshGen2DEngine.H"
  38 #include "polyMeshGenAddressing.H"
  39 #include "polyMeshGenChecks.H"
  40 #include "labelledPoint.H"
  41 #include "FIFOStack.H"
  42
  43 #include <map>
  44 #include <stdexcept>
  45
  46 # ifdef USE_OMP
  47 #include <omp.h>
  48 # endif
  49
  50 //#define DEBUGSmooth
  51
  52 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  53
  54 namespace Foam
  55 {
  56
  57 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  58
  59 label meshSurfaceOptimizer::findInvertedVertices
  60 (
  61     boolList& smoothVertex,
  62     const label nAdditionalLayers
  63 ) const
  64 {
  65     const labelList& bPoints = surfaceEngine_.boundaryPoints();
  66     const VRWGraph& pPoints = surfaceEngine_.pointPoints();
  67
  68     if( smoothVertex.size() != bPoints.size() )
  69     {
  70         smoothVertex.setSize(bPoints.size());
  71         smoothVertex = true;
  72     }
  73
  74     label nInvertedTria(0);
  75
  76     //- check the vertices at the surface
  77     //- mark the ones where the mesh is tangled
  78     meshSurfaceCheckInvertedVertices vrtCheck(*partitionerPtr_, smoothVertex);
  79     const labelHashSet& inverted = vrtCheck.invertedVertices();
  80
  81     smoothVertex = false;
  82     forAll(bPoints, bpI)
  83     {
  84         if( inverted.found(bPoints[bpI]) )
  85         {
  86             ++nInvertedTria;
  87             smoothVertex[bpI] = true;
  88         }
  89     }
  90
  91     if( Pstream::parRun() )
  92         reduce(nInvertedTria, sumOp<label>());
  93     Info << "Number of inverted boundary faces is " << nInvertedTria << endl;
  94
  95     if( nInvertedTria == 0 )
  96         return 0;
  97
  98     //- add additional layers around inverted points
  99     for(label i=0;i<nAdditionalLayers;++i)
 100     {
 101         boolList originallySelected = smoothVertex;
 102         forAll(smoothVertex, bpI)
 103             if( originallySelected[bpI] )
 104                 forAllRow(pPoints, bpI, ppI)
 105                     smoothVertex[pPoints(bpI, ppI)] = true;
 106
 107         if( Pstream::parRun() )
 108         {
 109             //- exchange global labels of inverted points
 110             const labelList& globalPointLabel =
 111                 surfaceEngine_.globalBoundaryPointLabel();
 112             const Map<label>& globalToLocal =
 113                 surfaceEngine_.globalToLocalBndPointAddressing();
 114             const VRWGraph& bpAtProcs = surfaceEngine_.bpAtProcs();
 115             const DynList<label>& neiProcs = surfaceEngine_.bpNeiProcs();
 116
 117             std::map<label, labelLongList> shareData;
 118             forAll(neiProcs, procI)
 119                 shareData.insert
 120                 (
 121                     std::make_pair(neiProcs[procI], labelLongList())
 122                 );
 123
 124             forAllConstIter(Map<label>, globalToLocal, iter)
 125             {
 126                 const label bpI = iter();
 127
 128                 if( !smoothVertex[bpI] )
 129                     continue;
 130
 131                 forAllRow(bpAtProcs, bpI, procI)
 132                 {
 133                     const label neiProc = bpAtProcs(bpI, procI);
 134
 135                     if( neiProc == Pstream::myProcNo() )
 136                         continue;
 137
 138                     shareData[neiProc].append(globalPointLabel[bpI]);
 139                 }
 140             }
 141
 142             //- exchange data with other processors
 143             labelLongList receivedData;
 144             help::exchangeMap(shareData, receivedData);
 145
 146             forAll(receivedData, j)
 147             {
 148                 const label bpI = globalToLocal[receivedData[j]];
 149
 150                 smoothVertex[bpI] = true;
 151             }
 152         }
 153     }
 154
 155     return nInvertedTria;
 156 }
 157
 158 void meshSurfaceOptimizer::smoothEdgePoints
 159 (
 160     const labelLongList& edgePoints,
 161     const labelLongList& procEdgePoints
 162 )
 163 {
 164     List<LongList<labelledPoint> > newPositions(1);
 165     # ifdef USE_OMP
 166     newPositions.setSize(omp_get_num_procs());
 167     # endif
 168
 169     //- smooth edge vertices
 170     # ifdef USE_OMP
 171     # pragma omp parallel num_threads(newPositions.size())
 172     # endif
 173     {
 174         # ifdef USE_OMP
 175         LongList<labelledPoint>& newPos =
 176             newPositions[omp_get_thread_num()];
 177         # else
 178         LongList<labelledPoint>& newPos = newPositions[0];
 179         # endif
 180
 181         # ifdef USE_OMP
 182         # pragma omp for schedule(dynamic, 40)
 183         # endif
 184         forAll(edgePoints, i)
 185         {
 186             const label bpI = edgePoints[i];
 187
 188             if( vertexType_[bpI] & (PROCBND | LOCKED) )
 189                 continue;
 190
 191             newPos.append(labelledPoint(bpI, newEdgePositionLaplacian(bpI)));
 192         }
 193     }
 194
 195     if( Pstream::parRun() )
 196         edgeNodeDisplacementParallel(procEdgePoints);
 197
 198     meshSurfaceEngineModifier surfaceModifier(surfaceEngine_);
 199     forAll(newPositions, threadI)
 200     {
 201         const LongList<labelledPoint>& newPos = newPositions[threadI];
 202
 203         forAll(newPos, i)
 204             surfaceModifier.moveBoundaryVertexNoUpdate
 205             (
 206                 newPos[i].pointLabel(),
 207                 newPos[i].coordinates()
 208             );
 209     }
 210
 211     surfaceModifier.updateGeometry(edgePoints);
 212 }
 213
 214 void meshSurfaceOptimizer::smoothLaplacianFC
 215 (
 216     const labelLongList& selectedPoints,
 217     const labelLongList& selectedProcPoints,
 218     const bool transform
 219 )
 220 {
 221     List<LongList<labelledPoint> > newPositions(1);
 222     # ifdef USE_OMP
 223     newPositions.setSize(omp_get_num_procs());
 224     # endif
 225
 226     # ifdef USE_OMP
 227     # pragma omp parallel num_threads(newPositions.size())
 228     # endif
 229     {
 230         # ifdef USE_OMP
 231         LongList<labelledPoint>& newPos =
 232             newPositions[omp_get_thread_num()];
 233         # else
 234         LongList<labelledPoint>& newPos = newPositions[0];
 235         # endif
 236
 237         # ifdef USE_OMP
 238         # pragma omp for schedule(dynamic, 40)
 239         # endif
 240         forAll(selectedPoints, i)
 241         {
 242             const label bpI = selectedPoints[i];
 243
 244             if( vertexType_[bpI] & (PROCBND | LOCKED) )
 245                 continue;
 246
 247             newPos.append
 248             (
 249                 labelledPoint(bpI, newPositionLaplacianFC(bpI, transform))
 250             );
 251         }
 252     }
 253
 254     if( Pstream::parRun() )
 255         nodeDisplacementLaplacianFCParallel(selectedProcPoints, transform);
 256
 257     meshSurfaceEngineModifier surfaceModifier(surfaceEngine_);
 258
 259     # ifdef USE_OMP
 260     # pragma omp parallel num_threads(newPositions.size())
 261     # endif
 262     {
 263         # ifdef USE_OMP
 264         const label threadI = omp_get_thread_num();
 265         # else
 266         const label threadI = 0;
 267         # endif
 268
 269         const LongList<labelledPoint>& newPos = newPositions[threadI];
 270
 271         forAll(newPos, i)
 272             surfaceModifier.moveBoundaryVertexNoUpdate
 273             (
 274                 newPos[i].pointLabel(),
 275                 newPos[i].coordinates()
 276             );
 277     }
 278
 279     surfaceModifier.updateGeometry(selectedPoints);
 280 }
 281
 282 void meshSurfaceOptimizer::smoothSurfaceOptimizer
 283 (
 284     const labelLongList& selectedPoints,
 285     const labelLongList& selectedProcPoints
 286 )
 287 {
 288     //- create partTriMesh is it is not yet present
 289     this->triMesh();
 290
 291     //- update coordinates of the triangulation
 292     updateTriMesh(selectedPoints);
 293
 294     pointField newPositions(selectedPoints.size());
 295
 296     # ifdef USE_OMP
 297     # pragma omp parallel for schedule(dynamic, 20)
 298     # endif
 299     forAll(selectedPoints, i)
 300     {
 301         const label bpI = selectedPoints[i];
 302
 303         newPositions[i] = newPositionSurfaceOptimizer(bpI);
 304     }
 305
 306     meshSurfaceEngineModifier surfaceModifier(surfaceEngine_);
 307
 308     # ifdef USE_OMP
 309     # pragma omp parallel for schedule(dynamic, 100)
 310     # endif
 311     forAll(newPositions, i)
 312     {
 313         const label bpI = selectedPoints[i];
 314
 315         surfaceModifier.moveBoundaryVertexNoUpdate(bpI, newPositions[i]);
 316     }
 317
 318     //- ensure that vertices at inter-processor boundaries are at the same
 319     //- location at all processors
 320     surfaceModifier.syncVerticesAtParallelBoundaries(selectedProcPoints);
 321
 322     //- update geometry addressing for moved points
 323     surfaceModifier.updateGeometry(selectedPoints);
 324 }
 325
 326 bool meshSurfaceOptimizer::untangleSurface
 327 (
 328     const labelLongList& selectedBoundaryPoints,
 329     const label nAdditionalLayers
 330 )
 331 {
 332     Info << "Starting untangling the surface of the volume mesh" << endl;
 333
 334     bool changed(false);
 335
 336     const labelList& bPoints = surfaceEngine_.boundaryPoints();
 337     const pointFieldPMG& points = surfaceEngine_.points();
 338     surfaceEngine_.pointFaces();
 339     surfaceEngine_.faceCentres();
 340     surfaceEngine_.pointPoints();
 341     surfaceEngine_.boundaryFacePatches();
 342     surfaceEngine_.pointNormals();
 343     surfaceEngine_.boundaryPointEdges();
 344
 345     if( Pstream::parRun() )
 346     {
 347         surfaceEngine_.bpAtProcs();
 348         surfaceEngine_.globalToLocalBndPointAddressing();
 349         surfaceEngine_.globalBoundaryPointLabel();
 350         surfaceEngine_.bpNeiProcs();
 351     }
 352
 353     boolList smoothVertex(bPoints.size(), false);
 354     # ifdef USE_OMP
 355     # pragma omp parallel for schedule(dynamic, 50)
 356     # endif
 357     forAll(selectedBoundaryPoints, i)
 358     {
 359         if( vertexType_[selectedBoundaryPoints[i]] & LOCKED )
 360             continue;
 361
 362         smoothVertex[selectedBoundaryPoints[i]] = true;
 363     }
 364
 365     meshSurfaceEngineModifier surfaceModifier(surfaceEngine_);
 366
 367     meshSurfaceMapper* mapperPtr = NULL;
 368     if( octreePtr_ )
 369         mapperPtr = new meshSurfaceMapper(*partitionerPtr_, *octreePtr_);
 370
 371     bool remapVertex(true);
 372     label nInvertedTria;
 373     label nGlobalIter(0);
 374
 375     labelLongList procBndPoints, movedPoints;
 376     labelLongList procEdgePoints, movedEdgePoints;
 377
 378     label minNumInverted(bPoints.size());
 379     FIFOStack<label> nInvertedHistory;
 380     pointField minInvertedPoints(bPoints.size());
 381
 382     do
 383     {
 384         label nIter(0), nAfterRefresh(0);
 385
 386         do
 387         {
 388             nInvertedTria =
 389                 findInvertedVertices(smoothVertex, nAdditionalLayers);
 390
 391             if( nInvertedTria == 0 )
 392             {
 393                 break;
 394             }
 395             else if( enforceConstraints_ && !remapVertex )
 396             {
 397                 polyMeshGen& mesh =
 398                     const_cast<polyMeshGen&>(surfaceEngine_.mesh());
 399
 400                 const label subsetId =
 401                     mesh.addPointSubset(badPointsSubsetName_);
 402
 403                 forAll(smoothVertex, bpI)
 404                     if( smoothVertex[bpI] )
 405                         mesh.addPointToSubset(subsetId, bPoints[bpI]);
 406
 407                 WarningIn
 408                 (
 409                     "bool meshSurfaceOptimizer::untangleSurface"
 410                     "(const labelLongList&, const label)"
 411                 ) << "Writing mesh with " << badPointsSubsetName_
 412                   << " subset. These points cannot be untangled"
 413                   << " without sacrificing geometry constraints. Exitting.."
 414                   << endl;
 415
 416                 returnReduce(1, sumOp<label>());
 417
 418                 throw std::logic_error
 419                 (
 420                     "bool meshSurfaceOptimizer::untangleSurface"
 421                     "(const labelLongList&, const label)"
 422                     "Cannot untangle mesh!!"
 423                 );
 424             }
 425
 426             //- find the min number of inverted points and
 427             //- add the last number to the stack
 428             if( nInvertedTria < minNumInverted )
 429             {
 430                 minNumInverted = nInvertedTria;
 431                 nAfterRefresh = 0;
 432
 433                 # ifdef USE_OMP
 434                 # pragma omp parallel for schedule(dynamic, 100)
 435                 # endif
 436                 forAll(bPoints, bpI)
 437                     minInvertedPoints[bpI] = points[bPoints[bpI]];
 438             }
 439
 440             //- count the number of iteration after the last minimum occurence
 441             ++nAfterRefresh;
 442
 443             //- update the stack
 444             nInvertedHistory.push(nInvertedTria);
 445             if( nInvertedHistory.size() > 2 )
 446                 nInvertedHistory.pop();
 447
 448             //- check if the number of inverted points reduces
 449             bool minimumInStack(false);
 450             forAllConstIter(FIFOStack<label>, nInvertedHistory, it)
 451                 if( it() == minNumInverted )
 452                     minimumInStack = true;
 453
 454             //- stop if the procedure does not minimise
 455             //- the number of inverted points
 456             if( !minimumInStack || (nAfterRefresh > 2) )
 457                 break;
 458
 459             //- find points which will be handled by the smoothers
 460             changed = true;
 461
 462             procBndPoints.clear();
 463             movedPoints.clear();
 464             procEdgePoints.clear();
 465             movedEdgePoints.clear();
 466
 467             forAll(bPoints, bpI)
 468             {
 469                 if( !smoothVertex[bpI] )
 470                     continue;
 471
 472                 if( vertexType_[bpI] & PARTITION )
 473                 {
 474                     movedPoints.append(bpI);
 475
 476                     if( vertexType_[bpI] & PROCBND )
 477                         procBndPoints.append(bpI);
 478                 }
 479                 else if( vertexType_[bpI] & EDGE )
 480                 {
 481                     movedEdgePoints.append(bpI);
 482
 483                     if( vertexType_[bpI] & PROCBND )
 484                         procEdgePoints.append(bpI);
 485                 }
 486             }
 487
 488             //- smooth edge vertices
 489             smoothEdgePoints(movedEdgePoints, procEdgePoints);
 490             if( remapVertex && mapperPtr )
 491                 mapperPtr->mapEdgeNodes(movedEdgePoints);
 492             surfaceModifier.updateGeometry(movedEdgePoints);
 493
 494             //- use laplacian smoothing
 495             smoothLaplacianFC(movedPoints, procBndPoints);
 496             surfaceModifier.updateGeometry(movedPoints);
 497
 498             //- use surface optimizer
 499             smoothSurfaceOptimizer(movedPoints, procBndPoints);
 500
 501             if( remapVertex && mapperPtr )
 502                 mapperPtr->mapVerticesOntoSurface(movedPoints);
 503
 504             //- update normals and other geometric data
 505             surfaceModifier.updateGeometry(movedPoints);
 506
 507         } while( nInvertedTria && (++nIter < 20) );
 508
 509         if( nInvertedTria > 0 )
 510         {
 511             //- use the combination with the minimum number of inverted points
 512             meshSurfaceEngineModifier sMod(surfaceEngine_);
 513             forAll(minInvertedPoints, bpI)
 514                 sMod.moveBoundaryVertexNoUpdate(bpI, minInvertedPoints[bpI]);
 515
 516             sMod.updateGeometry();
 517         }
 518
 519         if( nInvertedTria )
 520         {
 521             Info << "Smoothing remaining inverted vertices " << endl;
 522
 523             movedPoints.clear();
 524             procBndPoints.clear();
 525             forAll(smoothVertex, bpI)
 526                 if( smoothVertex[bpI] )
 527                 {
 528                     movedPoints.append(bpI);
 529
 530                     if( vertexType_[bpI] & PROCBND )
 531                         procBndPoints.append(bpI);
 532                 }
 533
 534             smoothLaplacianFC(movedPoints, procBndPoints, false);
 535
 536             if( remapVertex && mapperPtr )
 537                 mapperPtr->mapVerticesOntoSurface(movedPoints);
 538
 539             //- update normals and other geometric data
 540             surfaceModifier.updateGeometry(movedPoints);
 541
 542             if( nGlobalIter > 5 )
 543                 remapVertex = false;
 544         }
 545
 546     } while( nInvertedTria && (++nGlobalIter < 10) );
 547
 548     deleteDemandDrivenData(mapperPtr);
 549
 550     if( nInvertedTria != 0 )
 551     {
 552         //- the procedure has given up without success
 553         //- there exist some remaining inverted faces in the mesh
 554         polyMeshGen& mesh =
 555             const_cast<polyMeshGen&>(surfaceEngine_.mesh());
 556
 557         label subsetId = mesh.pointSubsetIndex(badPointsSubsetName_);
 558         if( subsetId >= 0 )
 559             mesh.removePointSubset(subsetId);
 560         subsetId = mesh.addPointSubset(badPointsSubsetName_);
 561
 562         forAll(smoothVertex, bpI)
 563             if( smoothVertex[bpI] )
 564                 mesh.addPointToSubset(subsetId, bPoints[bpI]);
 565     }
 566
 567     Info << "Finished untangling the surface of the volume mesh" << endl;
 568
 569     return changed;
 570 }
 571
 572 bool meshSurfaceOptimizer::untangleSurface(const label nAdditionalLayers)
 573 {
 574     labelLongList selectedPts(surfaceEngine_.boundaryPoints().size());
 575     forAll(selectedPts, i)
 576         selectedPts[i] = i;
 577
 578     return untangleSurface(selectedPts, nAdditionalLayers);
 579 }
 580
 581 void meshSurfaceOptimizer::optimizeSurface(const label nIterations)
 582 {
 583     const labelList& bPoints = surfaceEngine_.boundaryPoints();
 584
 585     //- needed for parallel execution
 586     surfaceEngine_.pointFaces();
 587     surfaceEngine_.faceCentres();
 588     surfaceEngine_.pointPoints();
 589     surfaceEngine_.boundaryPointEdges();
 590     surfaceEngine_.boundaryFacePatches();
 591     surfaceEngine_.pointNormals();
 592     surfaceEngine_.boundaryPointEdges();
 593
 594     meshSurfaceMapper* mapperPtr = NULL;
 595     if( octreePtr_ )
 596         mapperPtr = new meshSurfaceMapper(*partitionerPtr_, *octreePtr_);
 597
 598     labelLongList procBndPoints, edgePoints, partitionPoints, procPoints;
 599     forAll(bPoints, bpI)
 600     {
 601         if( vertexType_[bpI] & LOCKED )
 602             continue;
 603
 604         if( vertexType_[bpI] & EDGE )
 605         {
 606             edgePoints.append(bpI);
 607
 608             if( vertexType_[bpI] & PROCBND )
 609                 procBndPoints.append(bpI);
 610         }
 611         else if( vertexType_[bpI] & PARTITION )
 612         {
 613             partitionPoints.append(bpI);
 614
 615             if( vertexType_[bpI] & PROCBND )
 616                 procPoints.append(bpI);
 617         }
 618     }
 619
 620     //- optimize edge vertices
 621     Info << "Optimizing edges. Iteration:" << flush;
 622     for(label i=0;i<nIterations;++i)
 623     {
 624         Info << "." << flush;
 625
 626         meshSurfaceEngineModifier bMod(surfaceEngine_);
 627
 628         smoothEdgePoints(edgePoints, procBndPoints);
 629
 630         //- project vertices back onto the boundary
 631         if( mapperPtr )
 632             mapperPtr->mapEdgeNodes(edgePoints);
 633
 634         //- update the geometry information
 635         bMod.updateGeometry(edgePoints);
 636     }
 637     Info << endl;
 638
 639     //- delete the mapper
 640     deleteDemandDrivenData(mapperPtr);
 641
 642     //- optimize positions of surface vertices which are not on surface edges
 643     Info << "Optimizing surface vertices. Iteration:";
 644     for(label i=0;i<nIterations;++i)
 645     {
 646         smoothLaplacianFC(partitionPoints, procPoints, true);
 647
 648         smoothSurfaceOptimizer(partitionPoints, procPoints);
 649
 650         Info << "." << flush;
 651     }
 652
 653     Info << endl;
 654
 655     untangleSurface(0);
 656 }
 657
 658 void meshSurfaceOptimizer::optimizeSurface2D(const label nIterations)
 659 {
 660     const labelList& bPoints = surfaceEngine_.boundaryPoints();
 661     const edgeList& edges = surfaceEngine_.edges();
 662     const labelList& bp = surfaceEngine_.bp();
 663
 664     polyMeshGen2DEngine mesh2DEngine
 665     (
 666         const_cast<polyMeshGen&>(surfaceEngine_.mesh())
 667     );
 668     const boolList& zMinPoint = mesh2DEngine.zMinPoints();
 669
 670     //- needed for parallel execution
 671     surfaceEngine_.pointFaces();
 672     surfaceEngine_.faceCentres();
 673     surfaceEngine_.pointPoints();
 674     surfaceEngine_.boundaryPointEdges();
 675     surfaceEngine_.boundaryFacePatches();
 676     surfaceEngine_.pointNormals();
 677
 678     labelLongList procBndPoints, movedPoints, activeEdges, updatePoints;
 679     forAll(edges, beI)
 680     {
 681         const edge& e = edges[beI];
 682
 683         if( zMinPoint[e.start()] ^ zMinPoint[e.end()] )
 684         {
 685             label bpI = bp[e.start()];
 686             if( !zMinPoint[e.start()] )
 687                 bpI = bp[e.end()];
 688
 689             if( vertexType_[bpI] & EDGE )
 690             {
 691                 activeEdges.append(beI);
 692
 693                 updatePoints.append(bp[e.start()]);
 694                 updatePoints.append(bp[e.end()]);
 695
 696                 movedPoints.append(bpI);
 697
 698                 if( vertexType_[bpI] & PROCBND )
 699                     procBndPoints.append(bpI);
 700             }
 701         }
 702     }
 703
 704     meshSurfaceMapper2D* mapperPtr = NULL;
 705     if( octreePtr_ )
 706         mapperPtr = new meshSurfaceMapper2D(surfaceEngine_, *octreePtr_);
 707
 708     //- optimize edge vertices
 709     meshSurfaceEngineModifier bMod(surfaceEngine_);
 710
 711     Info << "Optimizing edges. Iteration:" << flush;
 712     for(label i=0;i<nIterations;++i)
 713     {
 714         Info << "." << flush;
 715
 716         smoothEdgePoints(movedPoints, procBndPoints);
 717
 718         //- move points with maximum z coordinate
 719         mesh2DEngine.correctPoints();
 720
 721         //- map boundary edges to the surface
 722         mapperPtr->mapVerticesOntoSurfacePatches(activeEdges);
 723
 724         //- update normal, centres, etc, after the surface has been modified
 725         bMod.updateGeometry(updatePoints);
 726     }
 727     Info << endl;
 728
 729     //- optimize Pts of surface vertices which are not on surface edges
 730     procBndPoints.clear();
 731     movedPoints.clear();
 732     forAll(bPoints, bpI)
 733         if( zMinPoint[bPoints[bpI]] && (vertexType_[bpI] & PARTITION) )
 734         {
 735             movedPoints.append(bpI);
 736
 737             if( vertexType_[bpI] & PROCBND )
 738                 procBndPoints.append(bpI);
 739         }
 740     Info << "Optimizing surface vertices. Iteration:";
 741     for(label i=0;i<nIterations;++i)
 742     {
 743         Info << "." << flush;
 744
 745         smoothSurfaceOptimizer(movedPoints, procBndPoints);
 746
 747         //- move the points which are not at minimum z coordinate
 748         mesh2DEngine.correctPoints();
 749
 750         //- update geometrical data due to movement of vertices
 751         bMod.updateGeometry();
 752     }
 753
 754     Info << endl;
 755
 756     deleteDemandDrivenData(mapperPtr);
 757 }
 758
 759 void meshSurfaceOptimizer::untangleSurface2D()
 760 {
 761     const polyMeshGen& mesh = surfaceEngine_.mesh();
 762     const faceListPMG& faces = mesh.faces();
 763     const VRWGraph& pointFaces = mesh.addressingData().pointFaces();
 764
 765     const labelList& bPoints = surfaceEngine_.boundaryPoints();
 766     const labelList& bp = surfaceEngine_.bp();
 767
 768     polyMeshGen2DEngine mesh2DEngine(const_cast<polyMeshGen&>(mesh));
 769     const boolList& zMinPoint = mesh2DEngine.zMinPoints();
 770     const boolList& activeFace = mesh2DEngine.activeFace();
 771
 772     //- needed for parallel execution
 773     surfaceEngine_.pointFaces();
 774     surfaceEngine_.faceCentres();
 775     surfaceEngine_.pointPoints();
 776     surfaceEngine_.boundaryPointEdges();
 777     surfaceEngine_.boundaryFacePatches();
 778     surfaceEngine_.pointNormals();
 779
 780     boolList activeBoundaryPoint(bPoints.size());
 781     boolList changedFace(activeFace.size(), true);
 782
 783     label iterationI(0);
 784     do
 785     {
 786         labelHashSet badFaces;
 787         const label nBadFaces =
 788             polyMeshGenChecks::findBadFaces
 789             (
 790                 mesh,
 791                 badFaces,
 792                 false,
 793                 &changedFace
 794             );
 795
 796         Info << "Iteration " << iterationI
 797              << ". Number of bad faces " << nBadFaces << endl;
 798
 799         if( nBadFaces == 0 )
 800             break;
 801
 802         //- update active points and faces affected by the movement
 803         //- of active points
 804         activeBoundaryPoint = false;
 805         changedFace = false;
 806         forAllConstIter(labelHashSet, badFaces, it)
 807         {
 808             const face& f = faces[it.key()];
 809
 810             forAll(f, pI)
 811             {
 812                 if( zMinPoint[f[pI]] )
 813                 {
 814                     activeBoundaryPoint[bp[f[pI]]] = true;
 815
 816                     forAllRow(pointFaces, f[pI], pfI)
 817                         changedFace[pointFaces(f[pI], pfI)] = true;
 818                 }
 819             }
 820         }
 821
 822         if( Pstream::parRun() )
 823         {
 824             const Map<label>& globalToLocal =
 825                 surfaceEngine_.globalToLocalBndPointAddressing();
 826             const DynList<label>& neiProcs = surfaceEngine_.bpNeiProcs();
 827             const VRWGraph& bpNeiProcs = surfaceEngine_.bpAtProcs();
 828
 829             std::map<label, labelLongList> exchangeData;
 830             forAll(neiProcs, i)
 831                 exchangeData[neiProcs[i]].clear();
 832
 833             //- collect active points at inter-processor boundaries
 834             forAllConstIter(Map<label>, globalToLocal, it)
 835             {
 836                 const label bpI = it();
 837
 838                 if( activeBoundaryPoint[bpI] )
 839                 {
 840                     forAllRow(bpNeiProcs, bpI, i)
 841                     {
 842                         const label neiProc = bpNeiProcs(bpI, i);
 843
 844                         if( neiProc == Pstream::myProcNo() )
 845                             continue;
 846
 847                         exchangeData[neiProc].append(it.key());
 848                     }
 849                 }
 850             }
 851
 852             //- exchange active points among the processors
 853             labelLongList receivedData;
 854             help::exchangeMap(exchangeData, receivedData);
 855
 856             //- ensure that all processors have the same Pts active
 857             forAll(receivedData, i)
 858             {
 859                 const label bpI = globalToLocal[receivedData[i]];
 860
 861                 //- activate this boundary point
 862                 activeBoundaryPoint[bpI] = true;
 863
 864                 //- set the changeFaces for the faces attached to this point
 865                 forAllRow(pointFaces, bPoints[bpI], pfI)
 866                     changedFace[pointFaces(bPoints[bpI], pfI)] = true;
 867             }
 868         }
 869
 870         //- apply smoothing to the activated points
 871         meshSurfaceEngineModifier bMod(surfaceEngine_);
 872
 873         labelLongList movedPts, procBndPts, edgePts, procEdgePts;
 874         forAll(bPoints, bpI)
 875         {
 876             if( !activeBoundaryPoint[bpI] )
 877                 continue;
 878
 879             if( vertexType_[bpI] & EDGE )
 880             {
 881                 edgePts.append(bpI);
 882
 883                 if( vertexType_[bpI] & PROCBND )
 884                     procEdgePts.append(bpI);
 885             }
 886             else if( vertexType_[bpI] & PARTITION )
 887             {
 888                 movedPts.append(bpI);
 889
 890                 if( vertexType_[bpI] & PROCBND )
 891                     procBndPts.append(bpI);
 892             }
 893         }
 894
 895         for(label i=0;i<5;++i)
 896         {
 897             smoothEdgePoints(edgePts, procEdgePts);
 898
 899             bMod.updateGeometry(edgePts);
 900
 901             smoothSurfaceOptimizer(movedPts, procBndPts);
 902
 903             bMod.updateGeometry(movedPts);
 904         }
 905
 906         //- move the points which are not at minimum z coordinate
 907         mesh2DEngine.correctPoints();
 908
 909         //- update geometrical data due to movement of vertices
 910         bMod.updateGeometry();
 911
 912         //- update cell centres and face centres
 913         const_cast<polyMeshGenAddressing&>
 914         (
 915             mesh.addressingData()
 916         ).updateGeometry(changedFace);
 917
 918     } while( ++iterationI < 20 );
 919
 920     //- delete invalid data
 921     mesh.clearAddressingData();
 922 }
 923
 924 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 925
 926 } // End namespace Foam
 927
 928 // ************************************************************************* //