src/mesh/cfMesh/utilities/surfaceTools/edgeExtraction/edgeExtractor/edgeExtractor.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 "error.H"
  29 #include "objectRegistry.H"
  30 #include "foamTime.H"
  31 #include "polyMeshGenModifier.H"
  32 #include "edgeExtractor.H"
  33 #include "meshSurfaceEngine.H"
  34 #include "meshSurfaceEngineModifier.H"
  35 #include "meshSurfaceOptimizer.H"
  36 #include "meshOctree.H"
  37 #include "triSurf.H"
  38 #include "triSurfModifier.H"
  39 #include "helperFunctions.H"
  40 #include "DynList.H"
  41 #include "labelPair.H"
  42 #include "labelledScalar.H"
  43 #include "labelledPoint.H"
  44 #include "refLabelledPoint.H"
  45 #include "HashSet.H"
  46 #include "triSurfacePartitioner.H"
  47 #include "triSurfaceClassifyEdges.H"
  48 #include "meshSurfaceMapper.H"
  49 #include "meshSurfaceCheckInvertedVertices.H"
  50 #include "meshSurfaceCheckEdgeTypes.H"
  51
  52 # ifdef USE_OMP
  53 #include <omp.h>
  54 # endif
  55
  56 //#define DEBUGEdgeExtractor
  57
  58 namespace Foam
  59 {
  60
  61 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *//
  62 // Private member functions
  63
  64 void edgeExtractor::calculateValence()
  65 {
  66     const meshSurfaceEngine& mse = this->surfaceEngine();
  67     pointValence_.setSize(mse.boundaryPoints().size());
  68     pointValence_ = 0;
  69
  70     const faceList::subList& bFaces = mse.boundaryFaces();
  71     const labelList& bp = mse.bp();
  72
  73     forAll(bFaces, bfI)
  74     {
  75         const face& bf = bFaces[bfI];
  76
  77         forAll(bf, pI)
  78             ++pointValence_[bp[bf[pI]]];
  79     }
  80
  81     if( Pstream::parRun() )
  82     {
  83         const Map<label>& globalToLocal =
  84             mse.globalToLocalBndPointAddressing();
  85         const VRWGraph& bpAtProcs = mse.bpAtProcs();
  86         const DynList<label>& bpNeiProcs = mse.bpNeiProcs();
  87
  88         std::map<label, LongList<labelPair> > exchangeData;
  89         forAll(bpNeiProcs, i)
  90             exchangeData.insert
  91             (
  92                 std::make_pair(bpNeiProcs[i], LongList<labelPair>())
  93             );
  94
  95         forAllConstIter(Map<label>, globalToLocal, iter)
  96         {
  97             const label bpI = iter();
  98
  99             forAllRow(bpAtProcs, bpI, i)
 100             {
 101                 const label neiProc = bpAtProcs(bpI, i);
 102
 103                 if( neiProc == Pstream::myProcNo() )
 104                     continue;
 105
 106                 exchangeData[neiProc].append
 107                 (
 108                     labelPair(iter.key(), pointValence_[bpI])
 109                 );
 110             }
 111         }
 112
 113         LongList<labelPair> receivedData;
 114         help::exchangeMap(exchangeData, receivedData);
 115
 116         forAll(receivedData, i)
 117         {
 118             const labelPair& lp = receivedData[i];
 119
 120             pointValence_[globalToLocal[lp.first()]] += lp.second();
 121         }
 122     }
 123 }
 124
 125 void edgeExtractor::calculateSingleCellEdge()
 126 {
 127     const meshSurfaceEngine& mse = this->surfaceEngine();
 128     const edgeList& edges = mse.edges();
 129     const VRWGraph& bpEdges = mse.boundaryPointEdges();
 130     const VRWGraph& edgeFaces = mse.edgeFaces();
 131     const labelList& faceCells = mse.faceOwners();
 132
 133     //- find the number of boundary faces for each cell in the mesh
 134     edgeType_.setSize(edgeFaces.size());
 135     edgeType_ = NONE;
 136
 137     forAll(edgeFaces, eI)
 138     {
 139         if( edgeFaces.sizeOfRow(eI) == 2 )
 140         {
 141             const label c0 = faceCells[edgeFaces(eI, 0)];
 142             const label c1 = faceCells[edgeFaces(eI, 1)];
 143
 144             if( c0 == c1 )
 145                 edgeType_[eI] |= SINGLECELLEDGE;
 146         }
 147     }
 148
 149     //- calculate the number of cells attache to a boundary edge
 150     const labelList& bp = mse.bp();
 151     const cellListPMG& cells = mse.mesh().cells();
 152     const faceListPMG& faces = mse.faces();
 153
 154     nCellsAtEdge_.setSize(edgeFaces.size());
 155     nCellsAtEdge_ = 0;
 156
 157     # ifdef USE_OMP
 158     # pragma omp parallel for schedule(dynamic, 100)
 159     # endif
 160     forAll(cells, cellI)
 161     {
 162         const cell& c = cells[cellI];
 163
 164         DynList<edge> foundEdge;
 165
 166         forAll(c, fI)
 167         {
 168             const face& f = faces[c[fI]];
 169
 170             forAll(f, eI)
 171             {
 172                 const edge e = f.faceEdge(eI);
 173
 174                 const label bps = bp[e.start()];
 175
 176                 if( bps < 0 )
 177                     continue;
 178
 179                 forAllRow(bpEdges, bps, i)
 180                 {
 181                     const label beI = bpEdges(bps, i);
 182                     const edge& be = edges[beI];
 183
 184                     if( (e == be) && !foundEdge.contains(be) )
 185                     {
 186                         foundEdge.append(be);
 187
 188                         # ifdef USE_OMP
 189                         # pragma omp atomic
 190                         # endif
 191                         ++nCellsAtEdge_[beI];
 192                     }
 193                 }
 194             }
 195         }
 196     }
 197 }
 198
 199 void edgeExtractor::findPatchesNearSurfaceFace()
 200 {
 201     const meshSurfaceEngine& mse = this->surfaceEngine();
 202     const pointFieldPMG& points = mse.points();
 203     const faceList::subList& bFaces = mse.boundaryFaces();
 204     const triSurf& surface = meshOctree_.surface();
 205
 206     patchesNearFace_.setSize(bFaces.size());
 207     labelLongList nPatchesAtFace(bFaces.size());
 208
 209     # ifdef USE_OMP
 210     # pragma omp parallel
 211     # endif
 212     {
 213         labelLongList localData;
 214         DynList<label> nearFacets;
 215
 216         # ifdef USE_OMP
 217         # pragma omp for schedule(dynamic, 40)
 218         # endif
 219         forAll(bFaces, bfI)
 220         {
 221             const face& bf = bFaces[bfI];
 222
 223             const vector c = bf.centre(points);
 224
 225             // find a reasonable searching distance comparable to face size
 226             scalar d(0.0);
 227             forAll(bf, pI)
 228                 d = Foam::max(d, Foam::mag(c - points[bf[pI]]));
 229             d = 2.0 * d + VSMALL;
 230
 231             const boundBox bb(c - vector(d, d, d), c + vector(d, d, d));
 232
 233             //- get the patches near the current boundary face
 234             meshOctree_.findTrianglesInBox(bb, nearFacets);
 235             DynList<label> nearPatches;
 236             forAll(nearFacets, i)
 237                 nearPatches.appendIfNotIn(surface[nearFacets[i]].region());
 238
 239             localData.append(bfI);
 240             nPatchesAtFace[bfI] = nearPatches.size();
 241             forAll(nearPatches, i)
 242                 localData.append(nearPatches[i]);
 243         }
 244
 245         # ifdef USE_OMP
 246         # pragma omp barrier
 247
 248         # pragma omp master
 249         patchesNearFace_.setSizeAndRowSize(nPatchesAtFace);
 250
 251         # pragma omp barrier
 252         # else
 253         patchesNearFace_.setSizeAndRowSize(nPatchesAtFace);
 254         # endif
 255
 256         //- copy the data to the graph
 257         label counter(0);
 258         while( counter < localData.size() )
 259         {
 260             const label edgeI = localData[counter++];
 261
 262             const label size = nPatchesAtFace[edgeI];
 263
 264             for(label i=0;i<size;++i)
 265                 patchesNearFace_(edgeI, i) = localData[counter++];
 266         }
 267     }
 268 }
 269
 270 void edgeExtractor::findFeatureEdgesNearEdge()
 271 {
 272     const meshSurfaceEngine& mse = this->surfaceEngine();
 273     const pointFieldPMG& points = mse.points();
 274     const edgeList& edges = mse.edges();
 275
 276     featureEdgesNearEdge_.setSize(edges.size());
 277     labelLongList nFeatureEdgesAtEdge(edges.size());
 278
 279     # ifdef USE_OMP
 280     # pragma omp parallel
 281     # endif
 282     {
 283         labelLongList localData;
 284         DynList<label> nearEdges;
 285
 286         # ifdef USE_OMP
 287         # pragma omp for schedule(dynamic, 40)
 288         # endif
 289         forAll(edges, edgeI)
 290         {
 291             const edge& e = edges[edgeI];
 292             const vector c = e.centre(points);
 293             const scalar d = 1.5 * e.mag(points);
 294
 295             const boundBox bb(c - vector(d, d, d), c + vector(d, d, d));
 296
 297             //- get the edges near the current edge
 298             meshOctree_.findEdgesInBox(bb, nearEdges);
 299             forAllReverse(nearEdges, i)
 300             {
 301                 const label pos = nearEdges.containsAtPosition(nearEdges[i]);
 302
 303                 if( pos < i )
 304                     nearEdges.removeElement(i);
 305             }
 306
 307             localData.append(edgeI);
 308             nFeatureEdgesAtEdge[edgeI] = nearEdges.size();
 309             forAll(nearEdges, i)
 310                 localData.append(nearEdges[i]);
 311         }
 312
 313         # ifdef USE_OMP
 314         # pragma omp barrier
 315
 316         # pragma omp master
 317         featureEdgesNearEdge_.setSizeAndRowSize(nFeatureEdgesAtEdge);
 318
 319         # pragma omp barrier
 320         # else
 321         featureEdgesNearEdge_.setSizeAndRowSize(nFeatureEdgesAtEdge);
 322         # endif
 323
 324         //- copy the data to the graph
 325         label counter(0);
 326         while( counter < localData.size() )
 327         {
 328             const label edgeI = localData[counter++];
 329
 330             const label size = nFeatureEdgesAtEdge[edgeI];
 331
 332             for(label i=0;i<size;++i)
 333                 featureEdgesNearEdge_(edgeI, i) = localData[counter++];
 334         }
 335     }
 336 }
 337
 338 void edgeExtractor::markPatchPoints(boolList& patchPoint)
 339 {
 340     const meshSurfaceEngine& mse = this->surfaceEngine();
 341     const labelList& bPoints = mse.boundaryPoints();
 342     const edgeList& edges = mse.edges();
 343     const VRWGraph& edgeFaces = mse.edgeFaces();
 344     const labelList& bp = mse.bp();
 345
 346     patchPoint.setSize(bPoints.size());
 347     patchPoint = true;
 348
 349     std::map<label, label> otherProcPatch;
 350     if( Pstream::parRun() )
 351     {
 352         const Map<label>& otherProc = mse.otherEdgeFaceAtProc();
 353         const Map<label>& globalToLocal =
 354             mse.globalToLocalBndEdgeAddressing();
 355
 356         //- create communication matrix
 357         std::map<label, labelLongList> exchangeData;
 358         const DynList<label>& neiProcs = mse.beNeiProcs();
 359         forAll(neiProcs, procI)
 360             exchangeData.insert
 361             (
 362                 std::make_pair(neiProcs[procI], labelLongList())
 363             );
 364
 365         forAllConstIter(Map<label>, globalToLocal, it)
 366         {
 367             const label beI = it();
 368
 369             if( edgeFaces.sizeOfRow(beI) == 1 )
 370             {
 371                 labelLongList& dts = exchangeData[otherProc[beI]];
 372                 //- send data as follows:
 373                 //- 1. global edge label
 374                 //- 2. patch of the attached boundary face
 375                 dts.append(it.key());
 376                 dts.append(facePatch_[edgeFaces(beI, 0)]);
 377             }
 378         }
 379
 380         labelLongList receivedData;
 381         help::exchangeMap(exchangeData, receivedData);
 382
 383         label counter(0);
 384         while( counter < receivedData.size() )
 385         {
 386             const label beI = globalToLocal[receivedData[counter++]];
 387             const label fPatch = receivedData[counter++];
 388
 389             otherProcPatch[beI] = fPatch;
 390         }
 391     }
 392
 393     //- set the patchPoint to false for all vertices at feature edges
 394     # ifdef USE_OMP
 395     # pragma omp parallel for schedule(dynamic, 40)
 396     # endif
 397     forAll(edgeFaces, beI)
 398     {
 399         if( edgeFaces.sizeOfRow(beI) == 2 )
 400         {
 401             //- an ordinary edge
 402             if( facePatch_[edgeFaces(beI, 0)] != facePatch_[edgeFaces(beI, 1)] )
 403             {
 404                 const edge& e = edges[beI];
 405                 patchPoint[bp[e.start()]] = false;
 406                 patchPoint[bp[e.end()]] = false;
 407             }
 408         }
 409         else if( edgeFaces.sizeOfRow(beI) == 1 )
 410         {
 411             //- an edge at a parallel interface
 412             const label otherPatch = otherProcPatch[beI];
 413
 414             if( facePatch_[edgeFaces(beI, 0)] != otherPatch )
 415             {
 416                 const edge& e = edges[beI];
 417                 patchPoint[bp[e.start()]] = false;
 418                 patchPoint[bp[e.end()]] = false;
 419             }
 420         }
 421         else
 422         {
 423             //- this is a non-manifold edge
 424             const edge& e = edges[beI];
 425             patchPoint[bp[e.start()]] = false;
 426             patchPoint[bp[e.end()]] = false;
 427         }
 428     }
 429
 430     if( Pstream::parRun() )
 431     {
 432         //- make sure that the information is spread to all processors
 433         const VRWGraph& bpAtProcs = mse.bpAtProcs();
 434         const DynList<label>& neiProcs = mse.bpNeiProcs();
 435         const labelList& globalPointLabel =
 436             mse.globalBoundaryPointLabel();
 437         const Map<label>& globalToLocal =
 438             mse.globalToLocalBndPointAddressing();
 439
 440
 441         std::map<label, labelLongList> sendData;
 442         forAll(neiProcs, i)
 443             sendData.insert(std::make_pair(neiProcs[i], labelLongList()));
 444
 445         forAll(bpAtProcs, bpI)
 446         {
 447             forAllRow(bpAtProcs, bpI, i)
 448             {
 449                 const label neiProc = bpAtProcs(bpI, i);
 450
 451                 if( neiProc != Pstream::myProcNo() )
 452                     sendData[neiProc].append(globalPointLabel[bpI]);
 453             }
 454         }
 455
 456         labelLongList receivedData;
 457         help::exchangeMap(sendData, receivedData);
 458
 459         forAll(receivedData, i)
 460                 patchPoint[globalToLocal[receivedData[i]]] = false;
 461     }
 462 }
 463
 464 const meshSurfaceEngine& edgeExtractor::surfaceEngine() const
 465 {
 466     if( !surfaceEnginePtr_ )
 467     {
 468         # ifdef USE_OMP
 469         # pragma omp critical
 470         # endif
 471         {
 472             if( !surfaceEnginePtr_ )
 473             {
 474                 surfaceEnginePtr_ = new meshSurfaceEngine(mesh_);
 475             }
 476         }
 477     }
 478
 479     return *surfaceEnginePtr_;
 480 }
 481
 482 const triSurfacePartitioner& edgeExtractor::partitioner() const
 483 {
 484     if( !surfPartitionerPtr_ )
 485     {
 486         # ifdef USE_OMP
 487         # pragma omp critical
 488         # endif
 489         {
 490             if( !surfPartitionerPtr_ )
 491             {
 492                 surfPartitionerPtr_ =
 493                     new triSurfacePartitioner(meshOctree_.surface());
 494             }
 495         }
 496     }
 497
 498     return *surfPartitionerPtr_;
 499 }
 500
 501 const triSurfaceClassifyEdges& edgeExtractor::edgeClassifier() const
 502 {
 503     if( !surfEdgeClassificationPtr_ )
 504     {
 505         surfEdgeClassificationPtr_ =
 506             new triSurfaceClassifyEdges(meshOctree_);
 507     }
 508
 509     return *surfEdgeClassificationPtr_;
 510 }
 511
 512 void edgeExtractor::findFaceCandidates
 513 (
 514     labelLongList& faceCandidates,
 515     const labelList* facePatchPtr,
 516     const Map<label>* otherFacePatchPtr
 517 ) const
 518 {
 519     faceCandidates.clear();
 520     if( !facePatchPtr )
 521         facePatchPtr = &facePatch_;
 522
 523     const labelList& fPatches = *facePatchPtr;
 524
 525     if( !otherFacePatchPtr )
 526     {
 527         Map<label> otherFacePatch;
 528         findOtherFacePatchesParallel(otherFacePatch, &fPatches);
 529
 530         otherFacePatchPtr = &otherFacePatch;
 531     }
 532
 533     const Map<label>& otherFacePatch = *otherFacePatchPtr;
 534
 535     const meshSurfaceEngine& mse = surfaceEngine();
 536     const VRWGraph& faceEdges = mse.faceEdges();
 537     const VRWGraph& edgeFaces = mse.edgeFaces();
 538
 539     # ifdef USE_OMP
 540     # pragma omp parallel if( faceEdges.size() > 1000 )
 541     # endif
 542     {
 543         # ifdef USE_OMP
 544         labelLongList procCandidates;
 545         # pragma omp for schedule(dynamic, 40)
 546         # endif
 547         forAll(faceEdges, bfI)
 548         {
 549             DynList<label> allNeiPatches;
 550             forAllRow(faceEdges, bfI, eI)
 551             {
 552                 const label beI = faceEdges(bfI, eI);
 553
 554                 if( edgeFaces.sizeOfRow(beI) == 2 )
 555                 {
 556                     label fNei = edgeFaces(beI, 0);
 557                     if( fNei == bfI )
 558                         fNei = edgeFaces(faceEdges(bfI, eI), 1);
 559
 560                     allNeiPatches.appendIfNotIn(fPatches[fNei]);
 561                 }
 562                 else if( edgeFaces.sizeOfRow(beI) == 1 )
 563                 {
 564                     allNeiPatches.appendIfNotIn(otherFacePatch[beI]);
 565                 }
 566             }
 567
 568             if( allNeiPatches.size() > 1 )
 569             {
 570                 //- this face is probably near some feature edge
 571                 # ifdef USE_OMP
 572                 procCandidates.append(bfI);
 573                 # else
 574                 faceCandidates.append(bfI);
 575                 # endif
 576             }
 577         }
 578
 579         # ifdef USE_OMP
 580         # pragma omp critical
 581         {
 582             forAll(procCandidates, i)
 583                 faceCandidates.append(procCandidates[i]);
 584         }
 585         # endif
 586     }
 587 }
 588
 589 void edgeExtractor::findOtherFacePatchesParallel
 590 (
 591     Map<label>& otherFacePatch,
 592     const labelList* facePatchPtr
 593 ) const
 594 {
 595     otherFacePatch.clear();
 596
 597     if( !facePatchPtr )
 598         facePatchPtr = &facePatch_;
 599
 600     const labelList& fPatches = *facePatchPtr;
 601
 602     if( Pstream::parRun() )
 603     {
 604         const meshSurfaceEngine& mse = this->surfaceEngine();
 605         const VRWGraph& edgeFaces = mse.edgeFaces();
 606         const Map<label>& otherProc = mse.otherEdgeFaceAtProc();
 607         const Map<label>& globalToLocal =
 608             mse.globalToLocalBndEdgeAddressing();
 609
 610         //- create communication matrix
 611         std::map<label, labelLongList> exchangeData;
 612         const DynList<label>& neiProcs = mse.beNeiProcs();
 613         forAll(neiProcs, procI)
 614             exchangeData.insert
 615             (
 616                 std::make_pair(neiProcs[procI], labelLongList())
 617             );
 618
 619         forAllConstIter(Map<label>, globalToLocal, it)
 620         {
 621             const label beI = it();
 622
 623             if( edgeFaces.sizeOfRow(beI) == 1 )
 624             {
 625                 labelLongList& dts = exchangeData[otherProc[beI]];
 626                 //- send data as follows:
 627                 //- 1. global edge label
 628                 //- 2. patch of the attached boundary face
 629                 dts.append(it.key());
 630                 dts.append(fPatches[edgeFaces(beI, 0)]);
 631             }
 632         }
 633
 634         labelLongList receivedData;
 635         help::exchangeMap(exchangeData, receivedData);
 636
 637         label counter(0);
 638         while( counter < receivedData.size() )
 639         {
 640             const label beI = globalToLocal[receivedData[counter++]];
 641             const label fPatch = receivedData[counter++];
 642
 643             otherFacePatch.insert(beI, fPatch);
 644         }
 645     }
 646 }
 647
 648 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *//
 649
 650 edgeExtractor::edgeExtractor
 651 (
 652     polyMeshGen& mesh,
 653     const meshOctree& octree
 654 )
 655 :
 656     mesh_(mesh),
 657     surfaceEnginePtr_(NULL),
 658     meshOctree_(octree),
 659     surfPartitionerPtr_(NULL),
 660     surfEdgeClassificationPtr_(NULL),
 661     pointValence_(),
 662     pointPatch_(),
 663     facePatch_(),
 664     nCellsAtEdge_(),
 665     edgeType_(),
 666     patchesNearFace_(),
 667     featureEdgesNearEdge_()
 668 {
 669     calculateValence();
 670
 671     calculateSingleCellEdge();
 672
 673     findPatchesNearSurfaceFace();
 674
 675     findFeatureEdgesNearEdge();
 676 }
 677
 678 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *//
 679 // Destructor
 680
 681 edgeExtractor::~edgeExtractor()
 682 {
 683     deleteDemandDrivenData(surfaceEnginePtr_);
 684     deleteDemandDrivenData(surfPartitionerPtr_);
 685     deleteDemandDrivenData(surfEdgeClassificationPtr_);
 686 }
 687
 688 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *//
 689
 690 void edgeExtractor::moveVerticesTowardsDiscontinuities(const label nIterations)
 691 {
 692     Info << "Reducing Hausdorff distance:" << flush;
 693
 694     const meshSurfaceEngine& mse = this->surfaceEngine();
 695     const labelList& bPoints = mse.boundaryPoints();
 696     const VRWGraph& pointFaces = mse.pointFaces();
 697     const pointFieldPMG& points = mse.points();
 698     const faceList::subList& bFaces = mse.boundaryFaces();
 699
 700     meshSurfaceEngineModifier modifier(mse);
 701
 702     vectorField faceCentreDisplacement(bFaces.size());
 703     List<labelledPoint> pointDisplacements(bPoints.size());
 704
 705     for(label iterI=0;iterI<nIterations;++iterI)
 706     {
 707         # ifdef USE_OMP
 708         # pragma omp parallel
 709         # endif
 710         {
 711             //- find displacements of face centres
 712             # ifdef USE_OMP
 713             # pragma omp for schedule(dynamic, 40)
 714             # endif
 715             forAll(bFaces, bfI)
 716             {
 717                 const vector centre = bFaces[bfI].centre(points);
 718
 719                 point newP;
 720                 scalar distSq;
 721                 label patchI, nearestTri;
 722                 meshOctree_.findNearestSurfacePoint
 723                 (
 724                     newP,
 725                     distSq,
 726                     nearestTri,
 727                     patchI,
 728                     centre
 729                 );
 730
 731                 faceCentreDisplacement[bfI] = newP - centre;
 732             }
 733
 734             //- initialise displacements
 735             # ifdef USE_OMP
 736             # pragma omp for schedule(dynamic, 40)
 737             # endif
 738             forAll(pointDisplacements, bpI)
 739                 pointDisplacements[bpI] = labelledPoint(0, vector::zero);
 740
 741             # ifdef USE_OMP
 742             # pragma omp barrier
 743             # endif
 744
 745             //- calculate displacements of boundary points as the average
 746             //- of face centre displacements
 747             # ifdef USE_OMP
 748             # pragma omp for schedule(dynamic, 40)
 749             # endif
 750             forAll(pointFaces, bpI)
 751             {
 752                 forAllRow(pointFaces, bpI, pfI)
 753                 {
 754                     pointDisplacements[bpI].coordinates() +=
 755                         faceCentreDisplacement[pointFaces(bpI, pfI)];
 756                     ++pointDisplacements[bpI].pointLabel();
 757                 }
 758             }
 759         }
 760
 761         if( Pstream::parRun() )
 762         {
 763             const Map<label>& globalToLocal =
 764                 mse.globalToLocalBndPointAddressing();
 765             const DynList<label>& neiProcs = mse.bpNeiProcs();
 766             const VRWGraph& bpAtProcs = mse.bpAtProcs();
 767
 768             std::map<label, LongList<refLabelledPoint> > exchangeData;
 769             forAll(neiProcs, i)
 770                 exchangeData[i] = LongList<refLabelledPoint>();
 771
 772             forAllConstIter(Map<label>, globalToLocal, iter)
 773             {
 774                 const label bpI = iter();
 775
 776                 forAllRow(bpAtProcs, bpI, i)
 777                 {
 778                     const label neiProc = bpAtProcs(bpI, i);
 779
 780                     if( neiProc == Pstream::myProcNo() )
 781                         continue;
 782
 783                     exchangeData[neiProc].append
 784                     (
 785                         refLabelledPoint(iter.key(), pointDisplacements[bpI])
 786                     );
 787                 }
 788             }
 789
 790             LongList<refLabelledPoint> receivedData;
 791             help::exchangeMap(exchangeData, receivedData);
 792
 793             forAll(receivedData, i)
 794             {
 795                 const label globalLabel = receivedData[i].objectLabel();
 796                 const labelledPoint& lp = receivedData[i].lPoint();
 797
 798                 const label bpI = globalToLocal[globalLabel];
 799
 800                 pointDisplacements[bpI].coordinates() += lp.coordinates();
 801                 pointDisplacements[bpI].pointLabel() += lp.pointLabel();
 802             }
 803         }
 804
 805         # ifdef USE_OMP
 806         # pragma omp parallel for schedule(dynamic, 40)
 807         # endif
 808         forAll(pointDisplacements, bpI)
 809         {
 810             const labelledPoint& lp = pointDisplacements[bpI];
 811             const point mp =
 812                 points[bPoints[bpI]] + lp.coordinates() / lp.pointLabel();
 813
 814             //Info << "Original point " << bPoints[bpI] << " has coordinates "
 815             //        << points[bPoints[bpI]] << endl;
 816             //Info << "Displacement vector " << lp.coordinates() / lp.pointLabel() << endl;
 817             //Info << "Moved point " << mp << endl;
 818
 819             point newPoint;
 820             label patchI, nt;
 821             scalar distSq;
 822
 823             meshOctree_.findNearestSurfacePoint(newPoint, distSq, nt, patchI, mp);
 824
 825             //Info << "Mapped point " << newPoint << nl << endl;
 826
 827             modifier.moveBoundaryVertexNoUpdate(bpI, newPoint);
 828         }
 829
 830         //- update geometry
 831         modifier.updateGeometry();
 832
 833         Info << '.' << flush;
 834     }
 835
 836     Info << endl;
 837 }
 838
 839 void edgeExtractor::distributeBoundaryFaces()
 840 {
 841     const meshSurfaceEngine& mse = this->surfaceEngine();
 842     const labelList& bPoints = mse.boundaryPoints();
 843     const faceList::subList& bFaces = mse.boundaryFaces();
 844     const pointFieldPMG& points = mse.points();
 845
 846     //- set the size of the facePatch list
 847     facePatch_.setSize(bFaces.size());
 848
 849     //- check if the mesh already has patches
 850     if( mesh_.boundaries().size() > 1 )
 851         Warning << "Mesh patches are already assigned!" << endl;
 852
 853     //- set size of patchNames, newBoundaryFaces_ and newBoundaryOwners_
 854     const triSurf& surface = meshOctree_.surface();
 855     const label nPatches = surface.patches().size();
 856
 857     //- find patches to which the surface points are mapped to
 858     pointPatch_.setSize(bPoints.size());
 859
 860     # ifdef USE_OMP
 861     # pragma omp parallel for schedule(dynamic, 40)
 862     # endif
 863     forAll(bPoints, bpI)
 864     {
 865         const point& bp = points[bPoints[bpI]];
 866
 867         label fPatch, nTri;
 868         point p;
 869         scalar distSq;
 870
 871         meshOctree_.findNearestSurfacePoint(p, distSq, nTri, fPatch, bp);
 872
 873         if( (fPatch > -1) && (fPatch < nPatches) )
 874         {
 875             pointPatch_[bpI] = fPatch;
 876         }
 877         else
 878         {
 879             pointPatch_[bpI] = nPatches;
 880             FatalErrorIn
 881             (
 882                 "void meshSurfaceEdgeExtractorNonTopo::"
 883                 "distributeBoundaryFaces()"
 884             ) << "Cannot distribute a boundary points " << bPoints[bpI]
 885                  << " into any surface patch!. Exiting.." << exit(FatalError);
 886         }
 887     }
 888
 889     //- find the patch for face by finding the patch nearest
 890     //- to the face centre
 891     # ifdef USE_OMP
 892     # pragma omp parallel for schedule(dynamic, 40)
 893     # endif
 894     forAll(bFaces, bfI)
 895     {
 896         const face& bf = bFaces[bfI];
 897
 898         const point c = bf.centre(points);
 899
 900         //- find the nearest surface patch to face centre
 901         label fPatch, nTri;
 902         point p;
 903         scalar distSq;
 904
 905         meshOctree_.findNearestSurfacePoint(p, distSq, nTri, fPatch, c);
 906
 907         if( (fPatch > -1) && (fPatch < nPatches) )
 908         {
 909             facePatch_[bfI] = fPatch;
 910         }
 911         else
 912         {
 913             facePatch_[bfI] = nPatches;
 914
 915             FatalErrorIn
 916             (
 917                 "void meshSurfaceEdgeExtractorNonTopo::"
 918                 "distributeBoundaryFaces()"
 919             ) << "Cannot distribute a face " << bFaces[bfI] << " into any "
 920                 << "surface patch!. Exiting.." << exit(FatalError);
 921         }
 922     }
 923 }
 924
 925 bool edgeExtractor::distributeBoundaryFacesNormalAlignment()
 926 {
 927     bool changed(false);
 928
 929     const pointFieldPMG& points = mesh_.points();
 930     const meshSurfaceEngine& mse = this->surfaceEngine();
 931     const faceList::subList& bFaces = mse.boundaryFaces();
 932     const VRWGraph& faceEdges = mse.faceEdges();
 933     const VRWGraph& edgeFaces = mse.edgeFaces();
 934
 935     const triSurf& surf = meshOctree_.surface();
 936     const pointField& sPoints = surf.points();
 937
 938     label nCorrected, nIterations(0);
 939     Map<label> otherProcNewPatch;
 940
 941     do
 942     {
 943         nCorrected = 0;
 944
 945         //- allocate a copy of boundary patches
 946         labelList newBoundaryPatches(facePatch_);
 947
 948         //- check whether there exist situations where a boundary face
 949         //- is surrounded by more faces in different patches than the
 950         //- faces in the current patch
 951         if( Pstream::parRun() )
 952         {
 953             findOtherFacePatchesParallel
 954             (
 955                 otherProcNewPatch,
 956                 &facePatch_
 957             );
 958         }
 959
 960         //- find the faces which have neighbouring faces in other patches
 961         labelLongList candidates;
 962         findFaceCandidates(candidates, &facePatch_, &otherProcNewPatch);
 963
 964         //- go through the list of faces and check if they shall remain
 965         //- in the current patch
 966         # ifdef USE_OMP
 967         # pragma omp parallel for schedule(dynamic, 40) \
 968         reduction(+ : nCorrected)
 969         # endif
 970         forAll(candidates, i)
 971         {
 972             const label bfI = candidates[i];
 973             const face& bf = bFaces[bfI];
 974
 975             DynList<label> allNeiPatches;
 976             DynList<label> neiPatches;
 977             neiPatches.setSize(faceEdges.sizeOfRow(bfI));
 978
 979             forAllRow(faceEdges, bfI, eI)
 980             {
 981                 const label beI = faceEdges(bfI, eI);
 982
 983                 if( edgeFaces.sizeOfRow(beI) == 2 )
 984                 {
 985                     label fNei = edgeFaces(beI, 0);
 986                     if( fNei == bfI )
 987                         fNei = edgeFaces(faceEdges(bfI, eI), 1);
 988
 989                     allNeiPatches.appendIfNotIn(facePatch_[fNei]);
 990                     neiPatches[eI] = facePatch_[fNei];
 991                 }
 992                 else if( edgeFaces.sizeOfRow(beI) == 1 )
 993                 {
 994                     allNeiPatches.appendIfNotIn(otherProcNewPatch[beI]);
 995                     neiPatches[eI] = otherProcNewPatch[beI];
 996                 }
 997             }
 998
 999             //- do not modify faces with all neighbours in the same patch
1000             if
1001             (
1002                 (allNeiPatches.size() == 1) &&
1003                 (allNeiPatches[0] == facePatch_[bfI])
1004             )
1005                 continue;
1006
1007             //- check whether there exist edges which are more suitable for
1008             //- projection onto feature edges than the currently selected ones
1009             label newPatch(-1);
1010             DynList<scalar> normalAlignment(allNeiPatches.size());
1011             DynList<scalar> distanceSq(allNeiPatches.size());
1012             scalar maxDSq(0.0);
1013             forAll(allNeiPatches, i)
1014             {
1015                 point pMap;
1016                 scalar dSq(VGREAT);
1017                 label nearestTriangle;
1018
1019                 point p = bf.centre(points);
1020                 meshOctree_.findNearestSurfacePointInRegion
1021                 (
1022                     pMap,
1023                     dSq,
1024                     nearestTriangle,
1025                     allNeiPatches[i],
1026                     p
1027                 );
1028
1029                 maxDSq = Foam::max(dSq, maxDSq);
1030
1031                 //- calculate normal vectors
1032                 vector tn = surf[nearestTriangle].normal(sPoints);
1033                 tn /= (mag(tn) + VSMALL);
1034                 vector fn = bf.normal(points);
1035                 fn /= (mag(fn) + SMALL);
1036
1037                 //- calculate alignment
1038                 normalAlignment[i] = mag(tn & fn);
1039                 distanceSq[i] = dSq;
1040             }
1041
1042             scalar maxAlignment(0.0);
1043             forAll(normalAlignment, i)
1044             {
1045                 const scalar metric
1046                 (
1047                     sqrt(maxDSq / (distanceSq[i] + VSMALL)) * normalAlignment[i]
1048                 );
1049
1050                 if( metric > maxAlignment )
1051                 {
1052                     maxAlignment = metric;
1053                     newPatch = allNeiPatches[i];
1054                 }
1055             }
1056
1057             if( (newPatch >= 0) && (newPatch != facePatch_[bfI]) )
1058             {
1059                 newBoundaryPatches[bfI] = newPatch;
1060                 ++nCorrected;
1061             }
1062         }
1063
1064         reduce(nCorrected, sumOp<label>());
1065
1066         if( nCorrected )
1067         {
1068             changed = true;
1069
1070             //- transfer the new patches back
1071             facePatch_.transfer(newBoundaryPatches);
1072         }
1073     } while( (nCorrected != 0) && (++nIterations < 5) );
1074
1075     return changed;
1076 }
1077
1078 void edgeExtractor::findEdgeCandidates()
1079 {
1080     const triSurf& surface = meshOctree_.surface();
1081     const vectorField& sp = surface.points();
1082     const VRWGraph& facetEdges = surface.facetEdges();
1083     const VRWGraph& edgeFacets = surface.edgeFacets();
1084
1085     const triSurfacePartitioner& partitioner = this->partitioner();
1086
1087     const meshSurfaceEngine& mse = this->surfaceEngine();
1088     const pointFieldPMG& points = mse.points();
1089     const labelList& bPoints = mse.boundaryPoints();
1090     const labelList& bp = mse.bp();
1091     const VRWGraph& faceEdges = mse.faceEdges();
1092
1093     Map<label> otherFacePatch;
1094     findOtherFacePatchesParallel(otherFacePatch, &facePatch_);
1095     labelLongList faceCandidates;
1096     findFaceCandidates(faceCandidates, &facePatch_, &otherFacePatch);
1097
1098     # ifdef USE_OMP
1099     # pragma omp parallel for schedule(dynamic, 40) \
1100     if( faceCandidates.size() > 100 )
1101     # endif
1102     forAll(faceCandidates, fcI)
1103     {
1104         const label bfI = faceCandidates[fcI];
1105
1106         forAllRow(faceEdges, bfI, i)
1107         {
1108             const label eI = faceEdges(bfI, i);
1109             edgeType_[eI] |= CANDIDATE;
1110         }
1111     }
1112
1113     //- find distances of all vertices supporting CANDIDATE edges
1114     //- from feature edges separating various patches
1115     const VRWGraph& pEdges = mse.boundaryPointEdges();
1116     const edgeList& edges = mse.edges();
1117
1118     List<List<labelledScalar> > featureEdgesNearPoint(bPoints.size());
1119
1120     DynList<label> containedTriangles;
1121     # ifdef USE_OMP
1122     # pragma omp parallel for schedule(dynamic, 40) private(containedTriangles)
1123     # endif
1124     forAll(pEdges, bpI)
1125     {
1126         // TODO rewrite for execution on distributed machines
1127         bool check(false);
1128         forAllRow(pEdges, bpI, peI)
1129         {
1130             const label eI = pEdges(bpI, peI);
1131
1132             if( edgeType_[eI] & CANDIDATE )
1133             {
1134                 check = true;
1135                 break;
1136             }
1137         }
1138
1139         if( check )
1140         {
1141             //- check the squared distance from the nearest feature edge
1142             scalar rSq(0.0);
1143             forAllRow(pEdges, bpI, peI)
1144             {
1145                 const label eI = pEdges(bpI, peI);
1146                 const edge& e = edges[eI];
1147                 const scalar dSq = magSqr(points[e.start()] - points[e.end()]);
1148
1149                 rSq = Foam::max(rSq, dSq);
1150             }
1151
1152             rSq *= 2.0;
1153             const scalar r = Foam::sqrt(rSq);
1154
1155             //- create a boundaing box used for searching neighbour edges
1156             const point& p = points[bPoints[bpI]];
1157             boundBox bb(p - point(r, r, r), p + point(r, r, r));
1158
1159             //- find the surface triangles in the vicinity of the point
1160             //- check for potential feature edges
1161             containedTriangles.clear();
1162             meshOctree_.findTrianglesInBox(bb, containedTriangles);
1163
1164             DynList<label> featureEdgeCandidates;
1165
1166             forAll(containedTriangles, ctI)
1167             {
1168                 const label tI = containedTriangles[ctI];
1169
1170                 forAllRow(facetEdges, tI, feI)
1171                 {
1172                     const label seI = facetEdges(tI, feI);
1173
1174                     if( edgeFacets.sizeOfRow(seI) == 2 )
1175                     {
1176                         const label p0 = surface[edgeFacets(seI, 0)].region();
1177                         const label p1 = surface[edgeFacets(seI, 1)].region();
1178
1179                         if( p0 != p1 )
1180                         {
1181                             featureEdgeCandidates.appendIfNotIn(seI);
1182                         }
1183                     }
1184                     else
1185                     {
1186                         featureEdgeCandidates.appendIfNotIn(seI);
1187                     }
1188                 }
1189             }
1190
1191             //- check the distance of the vertex from the candidates
1192             List<labelledScalar>& featureEdgeDistances =
1193                 featureEdgesNearPoint[bpI];
1194             featureEdgeDistances.setSize(featureEdgeCandidates.size());
1195             forAll(featureEdgeCandidates, i)
1196             {
1197                 const label seI = featureEdgeCandidates[i];
1198
1199                 const point s = sp[edges[seI].start()];
1200                 const point e = sp[edges[seI].end()];
1201                 const point np = help::nearestPointOnTheEdgeExact(s, e, p);
1202
1203                 featureEdgeDistances[i] = labelledScalar(seI, magSqr(np - p));
1204             }
1205
1206             //- find nearest edges
1207             sort(featureEdgeDistances);
1208         }
1209     }
1210
1211     //- start post-processing gathered data
1212     const labelList& edgeGroup = partitioner.edgeGroups();
1213
1214     List<List<labelledScalar> > edgeGroupAndWeights(edges.size());
1215
1216     # ifdef USE_OMP
1217     # pragma omp parallel for schedule(dynamic, 40) \
1218     if( edges.size() > 1000 )
1219     # endif
1220     forAll(edgeType_, edgeI)
1221     {
1222         if( edgeType_[edgeI] & CANDIDATE )
1223         {
1224             const edge& e = edges[edgeI];
1225
1226             const List<labelledScalar>& sc =
1227                 featureEdgesNearPoint[bp[e.start()]];
1228             const List<labelledScalar>& ec =
1229                 featureEdgesNearPoint[bp[e.end()]];
1230
1231             //- find the feature-edge partition for which the sum of
1232             //- node weights is minimal.
1233             DynList<labelledScalar> weights;
1234             forAll(sc, i)
1235             {
1236                 const label sPart = edgeGroup[sc[i].scalarLabel()];
1237
1238                 forAll(ec, j)
1239                 {
1240                     const label ePart = edgeGroup[ec[j].scalarLabel()];
1241
1242                     if( (sPart >= 0) && (sPart == ePart) )
1243                     {
1244                         weights.append
1245                         (
1246                             labelledScalar
1247                             (
1248                                 sPart,
1249                                 sc[i].value() + ec[j].value()
1250                             )
1251                         );
1252                     }
1253                 }
1254             }
1255
1256             //- store the data
1257             edgeGroupAndWeights[edgeI].setSize(weights.size());
1258             forAll(edgeGroupAndWeights[edgeI], epI)
1259                 edgeGroupAndWeights[edgeI][epI] = weights[epI];
1260
1261             //- sort the data according to the weights
1262             stableSort(edgeGroupAndWeights[edgeI]);
1263         }
1264     }
1265
1266     Info << "Edge partitions and weights " << edgeGroupAndWeights << endl;
1267 }
1268
1269 void edgeExtractor::findNeiPatches
1270 (
1271     const label bfI,
1272     const Map<label>& otherProcPatch,
1273     DynList<label>& neiPatches
1274 ) const
1275 {
1276     const meshSurfaceEngine& mse = surfaceEngine();
1277
1278     const VRWGraph& faceEdges = mse.faceEdges();
1279     const VRWGraph& edgeFaces = mse.edgeFaces();
1280
1281     neiPatches.setSize(faceEdges.sizeOfRow(bfI));
1282
1283     forAllRow(faceEdges, bfI, feI)
1284     {
1285         const label beI = faceEdges(bfI, feI);
1286
1287         if( edgeFaces.sizeOfRow(beI) == 2 )
1288         {
1289             label nei = edgeFaces(beI, 0);
1290             if( nei == bfI )
1291                 nei = edgeFaces(beI, 1);
1292
1293             neiPatches[feI] = facePatch_[nei];
1294         }
1295         else if( edgeFaces.sizeOfRow(beI) == 1 )
1296         {
1297             neiPatches[feI] = otherProcPatch[beI];
1298         }
1299     }
1300 }
1301
1302 scalar edgeExtractor::calculateAlignmentForEdge
1303 (
1304     const label beI,
1305     const label patch0,
1306     const label patch1
1307 ) const
1308 {
1309     scalar val(0.0);
1310
1311     DynList<label> patches(2);
1312     patches[0] = patch0;
1313     patches[1] = patch1;
1314
1315     const pointFieldPMG& points = surfaceEnginePtr_->mesh().points();
1316
1317     const edge& e = surfaceEnginePtr_->edges()[beI];
1318     const point& ps = points[e.start()];
1319     const point& pe = points[e.end()];
1320
1321     vector ev = e.vec(points);
1322     const scalar magE = mag(ev) + VSMALL;
1323     ev /= magE;
1324
1325     point mps, mpe;
1326     scalar dSqS, dSqE;
1327
1328     meshOctree_.findNearestPointToPatches(mps, dSqS, ps, patches);
1329     meshOctree_.findNearestPointToPatches(mpe, dSqE, pe, patches);
1330
1331     vector fv = mpe - mps;
1332     fv /= (mag(fv) + VSMALL);
1333
1334     val = 0.5 * (1.0 + (ev & fv));
1335
1336     val = min(val, 1.0);
1337     val = max(val, 0.0);
1338
1339     return val;
1340 }
1341
1342 scalar edgeExtractor::calculateDeformationMetricForEdge
1343 (
1344     const label beI,
1345     const label patch0,
1346     const label patch1
1347 ) const
1348 {
1349     scalar val(0.0);
1350
1351     DynList<label> patches(2);
1352     patches[0] = patch0;
1353     patches[1] = patch1;
1354
1355     const pointFieldPMG& points = surfaceEnginePtr_->mesh().points();
1356
1357     const edge& e = surfaceEnginePtr_->edges()[beI];
1358     const point& ps = points[e.start()];
1359     const point& pe = points[e.end()];
1360
1361     vector ev = e.vec(points);
1362     const scalar magE = mag(ev) + VSMALL;
1363     ev /= magE;
1364
1365     point mps, mpe;
1366     scalar dSqS, dSqE;
1367
1368     meshOctree_.findNearestPointToPatches(mps, dSqS, ps, patches);
1369     meshOctree_.findNearestPointToPatches(mpe, dSqE, pe, patches);
1370
1371     vector fv = mpe - mps;
1372     fv /= (mag(fv) + VSMALL);
1373
1374     scalar c = min(fv & ev, 1.0);
1375     c = max(-1.0, c);
1376     const scalar angle = acos(c);
1377
1378     val = 0.5 * (sqrt(dSqS) + sqrt(dSqE)) + magE * angle;
1379
1380     return val;
1381 }
1382
1383 scalar edgeExtractor::calculateDeformationMetricForFace
1384 (
1385     const label bfI,
1386     const DynList<label>& neiPatches,
1387     const label facePatch
1388 ) const
1389 {
1390     scalar Enorm(0.0);
1391
1392     const VRWGraph& faceEdges = surfaceEnginePtr_->faceEdges();
1393
1394     if( neiPatches.size() != faceEdges.sizeOfRow(bfI) )
1395     {
1396         FatalErrorIn
1397         (
1398             "scalar edgeExtractor::calculateDeformationMetricForFace"
1399             "(const label, const DynList<label>&, const label) const"
1400         ) << "Number of neiPatches and faceEdge does not match for face "
1401           << bfI << abort(FatalError);
1402     }
1403
1404     const label patch0 = facePatch == -1 ? facePatch_[bfI] : facePatch;
1405
1406     forAllRow(faceEdges, bfI, i)
1407     {
1408         const label beI = faceEdges(bfI, i);
1409
1410         if( neiPatches[i] == patch0 )
1411             continue;
1412
1413         Enorm += calculateDeformationMetricForEdge(beI, patch0, neiPatches[i]);
1414     }
1415
1416     return Enorm;
1417 }
1418
1419 bool edgeExtractor::checkConcaveEdgeCells()
1420 {
1421     bool changed(false);
1422
1423     const triSurf& surf = meshOctree_.surface();
1424     const VRWGraph& edgeFacets = surf.edgeFacets();
1425
1426     const pointFieldPMG& points = mesh_.points();
1427     const faceListPMG& faces = mesh_.faces();
1428     const cellListPMG& cells = mesh_.cells();
1429     const label bndStartFace = mesh_.boundaries()[0].patchStart();
1430
1431     const meshSurfaceEngine& mse = this->surfaceEngine();
1432     const faceList::subList& bFaces = mse.boundaryFaces();
1433     const labelList& bp = mse.bp();
1434     const labelList& faceCells = mse.faceOwners();
1435     const VRWGraph& edgeFaces = mse.edgeFaces();
1436
1437     //- analyse the surface mesh and find out which edges are concave or convex
1438     const triSurfaceClassifyEdges& edgeClassifier = this->edgeClassifier();
1439     const List<direction>& edgeType = edgeClassifier.edgeTypes();
1440
1441     //- create a copy of facePatch array for local modifications
1442     labelList newBoundaryPatches(facePatch_);
1443
1444     //- start checking the surface of the mesh
1445     label nChanged;
1446
1447     boolList patchPoint(mse.boundaryPoints().size(), false);
1448
1449     do
1450     {
1451         nChanged = 0;
1452
1453         //- check which surface points are surrounded by boundary faces
1454         //- in the same surface patch
1455         markPatchPoints(patchPoint);
1456
1457         //- check whether exist edges of a single cell which shall be projected
1458         //- onto a concave edge
1459         # ifdef USE_OMP
1460         # pragma omp parallel for schedule(dynamic, 40) reduction(+ : nChanged)
1461         # endif
1462         forAll(edgeType_, beI)
1463         {
1464             if( !(edgeType_[beI] & SINGLECELLEDGE) )
1465                 continue;
1466
1467             //- check if all faces are assigned to the same patch
1468             const label firstPatch = newBoundaryPatches[edgeFaces(beI, 0)];
1469             const label secondPatch = newBoundaryPatches[edgeFaces(beI, 1)];
1470
1471             if( firstPatch == secondPatch )
1472                 continue;
1473
1474             const cell& c = cells[faceCells[edgeFaces(beI, 0)]];
1475
1476             //- find edges within the bounding box determined by the cell
1477             point pMin(VGREAT, VGREAT, VGREAT);
1478             point pMax(-VGREAT, -VGREAT, -VGREAT);
1479             forAll(c, fI)
1480             {
1481                 const face& f = faces[c[fI]];
1482
1483                 forAll(f, pI)
1484                 {
1485                     pMin = Foam::min(pMin, points[f[pI]]);
1486                     pMax = Foam::max(pMax, points[f[pI]]);
1487                 }
1488             }
1489
1490             const point cc = 0.5 * (pMin + pMax);
1491             const point diff = pMax - pMin;
1492             boundBox bb(cc-diff, cc+diff);
1493             DynList<label> containedEdges;
1494             meshOctree_.findEdgesInBox(bb, containedEdges);
1495
1496             //- check if there exists concave edges boundaing patches
1497             //- assigned to boundary faces of the current cell
1498             forAll(containedEdges, ceI)
1499             {
1500                 const label eI = containedEdges[ceI];
1501
1502                 if( edgeFacets.sizeOfRow(eI) != 2 )
1503                     continue;
1504                 if( !(edgeType[eI] & triSurfaceClassifyEdges::FEATUREEDGE) )
1505                     continue;
1506
1507                 if( edgeType[eI] & triSurfaceClassifyEdges::CONCAVEEDGE )
1508                 {
1509                     const label patch0 = surf[edgeFacets(eI, 0)].region();
1510                     const label patch1 = surf[edgeFacets(eI, 1)].region();
1511
1512                     if
1513                     (
1514                         ((firstPatch == patch0) && (secondPatch == patch1)) ||
1515                         ((firstPatch == patch1) && (secondPatch == patch0))
1516                     )
1517                     {
1518                         DynList<DynList<label>, 2> facesInPatch;
1519                         facesInPatch.setSize(2);
1520
1521                         DynList<label, 2> nFacesInPatch;
1522                         nFacesInPatch.setSize(2);
1523                         nFacesInPatch = 0;
1524
1525                         DynList<bool, 2> hasPatchPoints;
1526                         hasPatchPoints.setSize(2);
1527                         hasPatchPoints = false;
1528
1529                         forAll(c, fI)
1530                         {
1531                             if( c[fI] < bndStartFace )
1532                                 continue;
1533
1534                             const label bfI = c[fI] - bndStartFace;
1535                             const face& bf = bFaces[bfI];
1536
1537                             if( newBoundaryPatches[bfI] == patch1 )
1538                             {
1539                                 facesInPatch[1].append(bfI);
1540                                 ++nFacesInPatch[1];
1541
1542                                 forAll(bf, pI)
1543                                 {
1544                                     if( patchPoint[bp[bf[pI]]] )
1545                                         hasPatchPoints[1] = true;
1546                                 }
1547                             }
1548                             else if( newBoundaryPatches[bfI] == patch0 )
1549                             {
1550                                 facesInPatch[0].append(bfI);
1551                                 ++nFacesInPatch[0];
1552
1553                                 forAll(bf, pI)
1554                                 {
1555                                     if( patchPoint[bp[bf[pI]]] )
1556                                         hasPatchPoints[0] = true;
1557                                 }
1558                             }
1559                         }
1560
1561                         if( nFacesInPatch[1] > nFacesInPatch[0] )
1562                         {
1563                             //- there exist more faces in patch 1
1564                             //- assign all boundary faces to the same patch
1565                             forAll(facesInPatch[0], i)
1566                                 newBoundaryPatches[facesInPatch[0][i]] = patch1;
1567                             ++nChanged;
1568                             break;
1569                         }
1570                         else if( nFacesInPatch[0] > nFacesInPatch[1] )
1571                         {
1572                             //- there exist more faces in patch 0
1573                             //- assign all boundary faces to the same patch
1574                             forAll(facesInPatch[1], i)
1575                                 newBoundaryPatches[facesInPatch[1][i]] = patch0;
1576                             ++nChanged;
1577                             break;
1578                         }
1579                         else
1580                         {
1581                             if( hasPatchPoints[0] && !hasPatchPoints[1] )
1582                             {
1583                                 //- transfer all faces to patch 1
1584                                 forAll(facesInPatch[0], i)
1585                                     newBoundaryPatches[facesInPatch[0][i]] =
1586                                         patch1;
1587                                 ++nChanged;
1588                                 break;
1589                             }
1590                             else if( !hasPatchPoints[0] && hasPatchPoints[1] )
1591                             {
1592                                 //- transfer all faces to patch 0
1593                                 forAll(facesInPatch[1], i)
1594                                     newBoundaryPatches[facesInPatch[1][i]] =
1595                                         patch0;
1596                                 ++nChanged;
1597                                 break;
1598                             }
1599                             else
1600                             {
1601                                 //- just transfer all faces to the same patch
1602                                 forAll(facesInPatch[1], i)
1603                                     newBoundaryPatches[facesInPatch[1][i]] =
1604                                         patch0;
1605                                 ++nChanged;
1606                                 break;
1607                             }
1608                         }
1609                     }
1610                 }
1611             }
1612         }
1613
1614         if( Pstream::parRun() )
1615             reduce(nChanged, sumOp<label>());
1616
1617         if( nChanged )
1618             changed = true;
1619
1620     } while( nChanged != 0 );
1621
1622     //- transfer the information back to facePatch
1623     facePatch_.transfer(newBoundaryPatches);
1624
1625     return changed;
1626 }
1627
1628 bool edgeExtractor::checkFacePatchesTopology()
1629 {
1630     bool changed(false);
1631
1632     const meshSurfaceEngine& mse = this->surfaceEngine();
1633     const faceList::subList& bFaces = mse.boundaryFaces();
1634     const labelList& bp = mse.bp();
1635     const VRWGraph& faceEdges = mse.faceEdges();
1636     const VRWGraph& edgeFaces = mse.edgeFaces();
1637
1638     label nCorrected;
1639     Map<label> otherProcNewPatch;
1640
1641     label nIter(0);
1642     do
1643     {
1644         # ifdef DEBUGEdgeExtractor
1645         {
1646             const triSurf* surfPtr = surfaceWithPatches();
1647             surfPtr->writeSurface
1648             (
1649                 "surfaceTopologyIter_"+help::scalarToText(nIter)+".stl"
1650             );
1651             delete surfPtr;
1652         }
1653         # endif
1654
1655         nCorrected = 0;
1656
1657         //- allocate a copy of boundary patches
1658         labelList newBoundaryPatches(facePatch_);
1659
1660         //- check whether there exist situations where a boundary face
1661         //- is surrounded by more faces in different patches than the
1662         //- faces in the current patch
1663         if( Pstream::parRun() )
1664         {
1665             findOtherFacePatchesParallel
1666             (
1667                 otherProcNewPatch,
1668                 &facePatch_
1669             );
1670         }
1671
1672         //- find the faces which have neighbouring faces in other patches
1673         labelLongList candidates;
1674         findFaceCandidates(candidates, &facePatch_, &otherProcNewPatch);
1675
1676         //- go through the list of faces and check if they shall remain
1677         //- in the current patch
1678         # ifdef USE_OMP
1679         # pragma omp parallel for schedule(dynamic, 40) \
1680         reduction(+ : nCorrected)
1681         # endif
1682         forAll(candidates, i)
1683         {
1684             const label bfI = candidates[i];
1685             const face& bf = bFaces[bfI];
1686
1687             //- do not change patches of faces where all points are mapped
1688             //- onto the same patch
1689             bool allInSamePatch(true);
1690             forAll(bf, pI)
1691                 if( pointPatch_[bp[bf[pI]]] != facePatch_[bfI] )
1692                 {
1693                     allInSamePatch = false;
1694                     break;
1695                 }
1696
1697             if( allInSamePatch )
1698                 continue;
1699
1700             DynList<label> allNeiPatches;
1701             DynList<label> neiPatches;
1702             neiPatches.setSize(faceEdges.sizeOfRow(bfI));
1703
1704             forAllRow(faceEdges, bfI, eI)
1705             {
1706                 const label beI = faceEdges(bfI, eI);
1707
1708                 if( edgeFaces.sizeOfRow(beI) == 2 )
1709                 {
1710                     label fNei = edgeFaces(beI, 0);
1711                     if( fNei == bfI )
1712                         fNei = edgeFaces(faceEdges(bfI, eI), 1);
1713
1714                     allNeiPatches.appendIfNotIn(facePatch_[fNei]);
1715                     neiPatches[eI] = facePatch_[fNei];
1716                 }
1717                 else if( edgeFaces.sizeOfRow(beI) == 1 )
1718                 {
1719                     allNeiPatches.appendIfNotIn(otherProcNewPatch[beI]);
1720                     neiPatches[eI] = otherProcNewPatch[beI];
1721                 }
1722             }
1723
1724             //- do not modify faces with all neighbours in the same patch
1725             if
1726             (
1727                 (allNeiPatches.size() == 1) &&
1728                 (allNeiPatches[0] == facePatch_[bfI])
1729             )
1730                 continue;
1731
1732             //- check whether there exist edges which are more suitable for
1733             //- projection onto feature edges than the currently selected ones
1734             label newPatch(-1);
1735
1736             //- check if some faces have to be distributed to another patch
1737             //- in order to reduce the number of feature edges
1738             Map<label> nNeiInPatch(allNeiPatches.size());
1739             forAll(allNeiPatches, i)
1740                 nNeiInPatch.insert(allNeiPatches[i], 0);
1741             forAll(neiPatches, eI)
1742                 ++nNeiInPatch[neiPatches[eI]];
1743
1744             newPatch = -1;
1745             label nNeiEdges(0);
1746             forAllConstIter(Map<label>, nNeiInPatch, it)
1747             {
1748                 if( it() > nNeiEdges )
1749                 {
1750                     newPatch = it.key();
1751                     nNeiEdges = it();
1752                 }
1753                 else if
1754                 (
1755                     (it() == nNeiEdges) && (it.key() == facePatch_[bfI])
1756                 )
1757                 {
1758                     newPatch = it.key();
1759                 }
1760             }
1761
1762             //- do not swap in case the
1763             if( (newPatch < 0) || (newPatch == facePatch_[bfI]) )
1764                 continue;
1765
1766             //- check whether the edges shared ith the neighbour patch form
1767             //- a singly linked chain
1768             DynList<bool> sharedEdge;
1769             sharedEdge.setSize(bFaces[bfI].size());
1770             sharedEdge = false;
1771
1772             forAll(neiPatches, eI)
1773                 if( neiPatches[eI] == newPatch )
1774                     sharedEdge[eI] = true;
1775
1776             if( help::areElementsInChain(sharedEdge) )
1777             {
1778                 //- change the patch to the newPatch
1779                 ++nCorrected;
1780                 newBoundaryPatches[bfI] = newPatch;
1781             }
1782         }
1783
1784         //- evaluate the new situation and ensure that no oscillation occur
1785         reduce(nCorrected, sumOp<label>());
1786         if( nCorrected )
1787         {
1788             faceEvaluator faceEvaluator(*this);
1789
1790             faceEvaluator.setNewBoundaryPatches(newBoundaryPatches);
1791
1792             # ifdef USE_OMP
1793             # pragma omp parallel for schedule(dynamic, 50)
1794             # endif
1795             forAll(candidates, i)
1796             {
1797                 const label bfI = candidates[i];
1798
1799                 const label bestPatch =
1800                     faceEvaluator.bestPatchAfterModification(bfI);
1801
1802                 newBoundaryPatches[bfI] = bestPatch;
1803             }
1804         }
1805
1806         if( nCorrected )
1807         {
1808             changed = true;
1809
1810             //- transfer the new patches back
1811             facePatch_.transfer(newBoundaryPatches);
1812         }
1813
1814     } while( nCorrected != 0 && (nIter++ < 3) );
1815
1816     return changed;
1817 }
1818
1819 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *//
1820
1821 namespace featureEdgeHelpers
1822 {
1823
1824 class featureEdgesNeiOp
1825 {
1826     // Private data
1827         //- reference to meshSurfaceEngine
1828         const meshSurfaceEngine& mse_;
1829
1830         //- refence to a list holding information which edges are feature edges
1831         const boolList& isFeatureEdge_;
1832
1833         //- number of feature edges at a surface point
1834         labelList nFeatureEdgesAtPoint_;
1835
1836     // Private member functions
1837         //- calculate the number of feature edges connected to a surface vertex
1838         void calculateNumberOfEdgesAtPoint()
1839         {
1840             const labelList& bp = mse_.bp();
1841             const edgeList& edges = mse_.edges();
1842
1843             nFeatureEdgesAtPoint_.setSize(mse_.boundaryPoints().size());
1844             nFeatureEdgesAtPoint_ = 0;
1845
1846             forAll(isFeatureEdge_, edgeI)
1847             {
1848                 if( !isFeatureEdge_[edgeI] )
1849                     continue;
1850
1851                 const edge& e = edges[edgeI];
1852                 ++nFeatureEdgesAtPoint_[bp[e.start()]];
1853                 ++nFeatureEdgesAtPoint_[bp[e.end()]];
1854             }
1855
1856             if( Pstream::parRun() )
1857             {
1858                 const Map<label>& globalToLocal =
1859                     mse_.globalToLocalBndPointAddressing();
1860                 const DynList<label>& neiProcs = mse_.bpNeiProcs();
1861                 const VRWGraph& bpAtProcs = mse_.bpAtProcs();
1862
1863                 std::map<label, DynList<labelPair> > exchangeData;
1864                 forAll(neiProcs, i)
1865                     exchangeData[neiProcs[i]].clear();
1866
1867                 //- fill the data from sending
1868                 forAllConstIter(Map<label>, globalToLocal, it)
1869                 {
1870                     const label bpI = it();
1871
1872                     forAllRow(bpAtProcs, bpI, i)
1873                     {
1874                         const label neiProc = bpAtProcs(bpI, i);
1875
1876                         if( neiProc == Pstream::myProcNo() )
1877                             continue;
1878
1879                         exchangeData[neiProc].append
1880                         (
1881                             labelPair(it.key(), nFeatureEdgesAtPoint_[bpI])
1882                         );
1883                     }
1884                 }
1885
1886                 //- exchange the data between the procesors
1887                 LongList<labelPair> receivedData;
1888                 help::exchangeMap(exchangeData, receivedData);
1889
1890                 forAll(receivedData, i)
1891                 {
1892                     const labelPair& lp = receivedData[i];
1893
1894                     nFeatureEdgesAtPoint_[globalToLocal[lp.first()]] +=
1895                         lp.second();
1896                 }
1897             }
1898         }
1899
1900 public:
1901
1902     featureEdgesNeiOp
1903     (
1904         const meshSurfaceEngine& mse,
1905         const boolList& isFeatureEdge
1906     )
1907     :
1908         mse_(mse),
1909         isFeatureEdge_(isFeatureEdge),
1910         nFeatureEdgesAtPoint_()
1911     {
1912         calculateNumberOfEdgesAtPoint();
1913     }
1914
1915     label size() const
1916     {
1917         return isFeatureEdge_.size();
1918     }
1919
1920     void operator()(const label beI, DynList<label>& neighbourEdges) const
1921     {
1922         neighbourEdges.clear();
1923
1924         const VRWGraph& bpEdges = mse_.boundaryPointEdges();
1925         const labelList& bp = mse_.bp();
1926         const edgeList& edges = mse_.edges();
1927
1928         const edge& e = edges[beI];
1929
1930         const label bps = bp[e.start()];
1931         const label bpe = bp[e.end()];
1932
1933         if( nFeatureEdgesAtPoint_[bps] == 2 )
1934         {
1935             forAllRow(bpEdges, bps, peI)
1936             {
1937                 const label beJ = bpEdges(bps, peI);
1938
1939                 if( (beJ == beI) || !isFeatureEdge_[beJ] )
1940                     continue;
1941
1942                 neighbourEdges.append(beJ);
1943             }
1944         }
1945
1946         if( nFeatureEdgesAtPoint_[bpe] == 2 )
1947         {
1948             forAllRow(bpEdges, bpe, peI)
1949             {
1950                 const label beJ = bpEdges(bpe, peI);
1951
1952                 if( (beJ == beI) || !isFeatureEdge_[beJ] )
1953                     continue;
1954
1955                 neighbourEdges.append(beJ);
1956             }
1957         }
1958     }
1959
1960     template<class labelListType>
1961     void collectGroups
1962     (
1963         std::map<label, DynList<label> >& neiGroups,
1964         const labelListType& elementInGroup,
1965         const DynList<label>& localGroupLabel
1966     ) const
1967     {
1968         const Map<label>& globalToLocal = mse_.globalToLocalBndPointAddressing();
1969         const VRWGraph& bpAtProcs = mse_.bpAtProcs();
1970         const VRWGraph& bpEdges = mse_.boundaryPointEdges();
1971
1972         const DynList<label>& neiProcs = mse_.beNeiProcs();
1973
1974         std::map<label, DynList<labelPair> > exchangeData;
1975         forAll(neiProcs, i)
1976             exchangeData[neiProcs[i]].clear();
1977
1978         forAllConstIter(Map<label>, globalToLocal, it)
1979         {
1980             const label bpI = it();
1981
1982             if( nFeatureEdgesAtPoint_[bpI] != 2 )
1983                 continue;
1984
1985             forAllRow(bpEdges, bpI, i)
1986             {
1987                 const label beI = bpEdges(bpI, i);
1988
1989                 if( !isFeatureEdge_[beI] )
1990                     continue;
1991
1992                 const label groupI = elementInGroup[beI];
1993
1994                 forAllRow(bpAtProcs, bpI, ppI)
1995                 {
1996                     const label neiProc = bpAtProcs(bpI, ppI);
1997
1998                     if( neiProc == Pstream::myProcNo() )
1999                         continue;
2000
2001                     exchangeData[neiProc].append
2002                     (
2003                         labelPair(it.key(), localGroupLabel[groupI])
2004                     );
2005                 }
2006             }
2007         }
2008
2009         LongList<labelPair> receivedData;
2010         help::exchangeMap(exchangeData, receivedData);
2011
2012         forAll(receivedData, i)
2013         {
2014             const labelPair& lp = receivedData[i];
2015             const label groupI = elementInGroup[globalToLocal[lp.first()]];
2016
2017             DynList<label>& ng = neiGroups[localGroupLabel[groupI]];
2018
2019             //- store the connection over the inter-processor boundary
2020             ng.appendIfNotIn(lp.second());
2021         }
2022     }
2023 };
2024
2025 class featureEdgesSelOp
2026 {
2027     // Private data
2028         //- reference to a list holding information which edge is afeture edge
2029         const boolList& isFeatureEdge_;
2030
2031 public:
2032
2033     featureEdgesSelOp(const boolList& isFeatureEdge)
2034     :
2035         isFeatureEdge_(isFeatureEdge)
2036     {}
2037
2038     bool operator()(const label beI) const
2039     {
2040         return isFeatureEdge_[beI];
2041     }
2042 };
2043
2044 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *//
2045
2046 } // End namespace featureEdgeHelpers
2047
2048 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
2049
2050 bool edgeExtractor::checkFacePatchesGeometry()
2051 {
2052     bool changed(false);
2053
2054     const meshSurfaceEngine& mse = this->surfaceEngine();
2055     const labelList& bPoints = mse.boundaryPoints();
2056     const faceList::subList& bFaces = mse.boundaryFaces();
2057     const labelList& bp = mse.bp();
2058
2059     //- allocate a copy of boundary patches
2060     labelList newBoundaryPatches(facePatch_.size());
2061
2062     label nCorrected;
2063     Map<label> otherProcNewPatch;
2064
2065     boolList activePoints(bPoints.size(), true);
2066     labelLongList activePointLabel(bPoints.size());
2067     forAll(bPoints, bpI)
2068         activePointLabel[bpI] = bpI;
2069
2070     label iter(0);
2071
2072     do
2073     {
2074         # ifdef DEBUGEdgeExtractor
2075         {
2076             const triSurf* surfPtr = surfaceWithPatches();
2077             surfPtr->writeSurface
2078             (
2079                 "surfaceIter_"+help::scalarToText(iter)+".stl"
2080             );
2081             delete surfPtr;
2082         }
2083         # endif
2084
2085         //- create feature edges and corners information
2086         meshSurfacePartitioner mPart(mse, facePatch_);
2087
2088         //- project vertices onto the surface mesh
2089         meshSurfaceMapper(mPart, meshOctree_).mapVerticesOntoSurfacePatches
2090         (
2091             activePointLabel
2092         );
2093
2094         //- stop after a certain number of iterations
2095         if( ++iter > 3 )
2096             break;
2097
2098         //- update surface geometry data
2099         meshSurfaceEngineModifier(mse).updateGeometry();
2100
2101         //- check if there exist any inverted faces
2102         meshSurfaceCheckInvertedVertices surfCheck(mPart, activePoints);
2103         const labelHashSet& invertedPoints = surfCheck.invertedVertices();
2104
2105         if( returnReduce(invertedPoints.size(), sumOp<label>()) == 0 )
2106             return false;
2107
2108         WarningIn
2109         (
2110             "void edgeExtractor::extractEdges()"
2111         ) << "Found " << invertedPoints.size()
2112           << " points with inverted surface normals. Getting rid of them..."
2113           << endl;
2114
2115         //- untangle the surface
2116         activePointLabel.clear();
2117         activePoints = false;
2118         forAllConstIter(labelHashSet, invertedPoints, it)
2119         {
2120             activePointLabel.append(bp[it.key()]);
2121             activePoints[bp[it.key()]] = true;
2122         }
2123
2124         //- untangle the surface
2125         meshSurfaceOptimizer mso(mPart, meshOctree_);
2126         mso.untangleSurface(activePointLabel, 1);
2127
2128         nCorrected = 0;
2129         newBoundaryPatches = facePatch_;
2130
2131         //- check whether there exist situations where a boundary face
2132         //- is surrounded by more faces in different patches than the
2133         //- faces in the current patch
2134         if( Pstream::parRun() )
2135         {
2136             findOtherFacePatchesParallel
2137             (
2138                 otherProcNewPatch,
2139                 &facePatch_
2140             );
2141         }
2142
2143         //- find the faces which have neighbouring faces in other patches
2144         labelLongList candidates;
2145         forAll(bFaces, bfI)
2146         {
2147             const face& bf = bFaces[bfI];
2148
2149             forAll(bf, pI)
2150                 if( invertedPoints.found(bf[pI]) )
2151                 {
2152                     candidates.append(bfI);
2153                     break;
2154                 }
2155         }
2156
2157         # ifdef USE_OMP
2158         # pragma omp parallel for schedule(dynamic, 5) reduction(+ : nCorrected)
2159         # endif
2160         forAll(candidates, i)
2161         {
2162             const label bfI = candidates[i];
2163
2164             DynList<label> neiPatches;
2165             findNeiPatches(bfI, otherProcNewPatch, neiPatches);
2166
2167             DynList<label> allNeiPatches;
2168             forAll(neiPatches, i)
2169                 allNeiPatches.appendIfNotIn(neiPatches[i]);
2170
2171             //- check the deformation energy and find the minimum energy which
2172             //- can be achieved by switching face patch
2173             scalar minE(VGREAT);
2174             label minEPatch(-1);
2175             DynList<scalar> Enorm(allNeiPatches.size());
2176             forAll(allNeiPatches, i)
2177             {
2178                 Enorm[i] =
2179                     calculateDeformationMetricForFace
2180                     (
2181                         bfI,
2182                         neiPatches,
2183                         allNeiPatches[i]
2184                     );
2185
2186                 if( Enorm[i] < minE )
2187                 {
2188                     minE = Enorm[i];
2189                     minEPatch = allNeiPatches[i];
2190                 }
2191             }
2192
2193             if( minEPatch != facePatch_[bfI] )
2194             {
2195                 newBoundaryPatches[bfI] = minEPatch;
2196                 ++nCorrected;
2197             }
2198         }
2199
2200         //- check if any faces are re-assigned to some other patch
2201         reduce(nCorrected, sumOp<label>());
2202         if( nCorrected == 0 )
2203             break;
2204
2205         //- update faceEvaluator with information after patches have been
2206         //- altered. It blocks chaning of patches if it causes oscillations
2207         faceEvaluator facePatchEvaluator(*this);
2208         facePatchEvaluator.setNewBoundaryPatches(newBoundaryPatches);
2209
2210         //- compare face patches before and after
2211         //- disallow modification which may trigger oscillating behaviour
2212         labelLongList changedFaces;
2213         forAll(newBoundaryPatches, bfI)
2214         {
2215             if( newBoundaryPatches[bfI] != facePatch_[bfI] )
2216             {
2217                 const label patchI =
2218                     facePatchEvaluator.bestPatchAfterModification(bfI);
2219                 newBoundaryPatches[bfI] = patchI;
2220
2221                 if( patchI != facePatch_[bfI] )
2222                     changedFaces.append(bfI);
2223             }
2224         }
2225
2226         nCorrected = changedFaces.size();
2227
2228         reduce(nCorrected, sumOp<label>());
2229
2230         if( nCorrected )
2231         {
2232             changed = true;
2233             facePatch_ = newBoundaryPatches;
2234         }
2235
2236     } while( nCorrected != 0 );
2237
2238     return changed;
2239 }
2240
2241 void edgeExtractor::projectDeterminedFeatureVertices()
2242 {
2243     List<DynList<label, 5> > pointPatches;
2244     pointPatches.setSize(pointValence_.size());
2245
2246     const meshSurfaceEngine& mse = surfaceEngine();
2247     const pointFieldPMG& points = mse.mesh().points();
2248     const labelList& bPoints = mse.boundaryPoints();
2249     const labelList& bp = mse.bp();
2250     const faceList::subList& bFaces = mse.boundaryFaces();
2251     meshOctree_.surface().pointEdges();
2252
2253     //- calculate patches for each point
2254     forAll(bFaces, bfI)
2255     {
2256         const face& bf = bFaces[bfI];
2257
2258         forAll(bf, pI)
2259             pointPatches[bp[bf[pI]]].appendIfNotIn(facePatch_[bfI]);
2260     }
2261
2262     if( Pstream::parRun() )
2263     {
2264         const Map<label>& globalToLocal =
2265             mse.globalToLocalBndPointAddressing();
2266         const VRWGraph& bpAtProcs = mse.bpAtProcs();
2267         const DynList<label>& bpNeiProcs = mse.bpNeiProcs();
2268
2269         std::map<label, LongList<labelPair> > exchangeData;
2270         forAll(bpNeiProcs, i)
2271             exchangeData.insert
2272             (
2273                 std::make_pair(bpNeiProcs[i], LongList<labelPair>())
2274             );
2275
2276         //- collect the data distributed to others
2277         forAllConstIter(Map<label>, globalToLocal, it)
2278         {
2279             const label bpI = it();
2280
2281             const DynList<label, 5>& pPatches = pointPatches[bpI];
2282
2283             forAllRow(bpAtProcs, bpI, i)
2284             {
2285                 const label neiProc = bpAtProcs(bpI, i);
2286
2287                 if( neiProc == Pstream::myProcNo() )
2288                     continue;
2289
2290                 LongList<labelPair>& data = exchangeData[neiProc];
2291
2292                 forAll(pPatches, ppI)
2293                     data.append(labelPair(it.key(), pPatches[ppI]));
2294             }
2295         }
2296
2297         //- exchange information
2298         LongList<labelPair> receivedData;
2299         help::exchangeMap(exchangeData, receivedData);
2300
2301         //- unify the data
2302         forAll(receivedData, i)
2303         {
2304             const labelPair& lp = receivedData[i];
2305
2306             pointPatches[globalToLocal[lp.first()]].appendIfNotIn(lp.second());
2307         }
2308     }
2309
2310     meshSurfaceEngineModifier surfMod(mse);
2311
2312     # ifdef USE_OMP
2313     # pragma omp parallel for schedule(dynamic, 10)
2314     # endif
2315     forAll(pointPatches, bpI)
2316     {
2317         if( pointPatches[bpI].size() < 2 )
2318             continue;
2319
2320         const DynList<label> pPatches = pointPatches[bpI];
2321
2322         const point& p = points[bPoints[bpI]];
2323
2324         //- find the nearest object on the surface mesh
2325         point newP;
2326         scalar dSqExact;
2327         if( pPatches.size() == 2 )
2328         {
2329             label nse;
2330             meshOctree_.findNearestEdgePoint(newP, dSqExact, nse, p, pPatches);
2331         }
2332         else
2333         {
2334             label nsp;
2335             meshOctree_.findNearestCorner(newP, dSqExact, nsp, p, pPatches);
2336         }
2337
2338         //- find the nearest object in an iterative procedure
2339         point pp(p);
2340         for(label iterI=0;iterI<20;++iterI)
2341         {
2342             point inp(vector::zero);
2343
2344             forAll(pPatches, i)
2345             {
2346                 point np;
2347                 scalar dSq;
2348                 label nt;
2349
2350                 meshOctree_.findNearestSurfacePointInRegion
2351                 (
2352                     np,
2353                     dSq,
2354                     nt,
2355                     pPatches[i],
2356                     pp
2357                 );
2358
2359                 inp += np;
2360             }
2361
2362             inp /= pPatches.size();
2363             const scalar currDSq = magSqr(inp - pp);
2364             pp = inp;
2365
2366             if( currDSq < 1e-2 * dSqExact )
2367                 break;
2368         }
2369
2370         //- check if the exact position of the corner is further away
2371         //- than the iteratively found object
2372         if( dSqExact > 1.1 * magSqr(pp - p) )
2373             newP = pp;
2374
2375         surfMod.moveBoundaryVertexNoUpdate(bpI, newP);
2376     }
2377
2378     surfMod.syncVerticesAtParallelBoundaries();
2379     surfMod.updateGeometry();
2380 }
2381
2382 bool edgeExtractor::untangleSurface()
2383 {
2384     bool changed(false);
2385
2386     meshSurfaceEngine& mse =
2387         const_cast<meshSurfaceEngine&>(this->surfaceEngine());
2388     meshSurfaceOptimizer optimizer(mse, meshOctree_);
2389     changed = optimizer.untangleSurface();
2390
2391     return changed;
2392 }
2393
2394 void edgeExtractor::extractEdges()
2395 {
2396     distributeBoundaryFaces();
2397
2398     distributeBoundaryFacesNormalAlignment();
2399
2400     # ifdef DEBUGEdgeExtractor
2401     const triSurf* sPtr = surfaceWithPatches();
2402     sPtr->writeSurface("initialDistributionOfPatches.stl");
2403     deleteDemandDrivenData(sPtr);
2404     # endif
2405
2406     Info << "Starting topological adjustment of patches" << endl;
2407     if( checkFacePatchesTopology() )
2408     {
2409         Info << "Finished topological adjustment of patches" << endl;
2410
2411         # ifdef DEBUGEdgeExtractor
2412         Info << "Changes due to face patches" << endl;
2413         fileName sName("checkFacePatches"+help::scalarToText(nIter)+".stl");
2414         sPtr = surfaceWithPatches();
2415         sPtr->writeSurface(sName);
2416         deleteDemandDrivenData(sPtr);
2417         # endif
2418     }
2419     else
2420     {
2421         Info << "No topological adjustment was needed" << endl;
2422     }
2423
2424     Info << "Starting geometrical adjustment of patches" << endl;
2425     if( checkFacePatchesGeometry() )
2426     {
2427         Info << "Finished geometrical adjustment of patches" << endl;
2428     }
2429     else
2430     {
2431         Info << "No geometrical adjustment was needed" << endl;
2432     }
2433
2434     # ifdef DEBUGEdgeExtractor
2435     const triSurf* sPtr = surfaceWithPatches();
2436     sPtr->writeSurface("finalDistributionOfPatches.stl");
2437     deleteDemandDrivenData(sPtr);
2438     # endif
2439 }
2440
2441 const triSurf* edgeExtractor::surfaceWithPatches() const
2442 {
2443     //- allocate the memory for the surface mesh
2444     triSurf* surfPtr = new triSurf();
2445
2446     //- surface of the volume mesh
2447     const meshSurfaceEngine& mse = surfaceEngine();
2448     const faceList::subList& bFaces = mse.boundaryFaces();
2449     const labelList& bp = mse.bp();
2450     const pointFieldPMG& points = mesh_.points();
2451
2452     //- modifier of the new surface mesh
2453     triSurfModifier surfModifier(*surfPtr);
2454     surfModifier.patchesAccess() = meshOctree_.surface().patches();
2455     pointField& sPts = surfModifier.pointsAccess();
2456     sPts.setSize(mse.boundaryPoints().size());
2457
2458     //- copy points
2459     forAll(bp, pointI)
2460     {
2461         if( bp[pointI] < 0 )
2462             continue;
2463
2464         sPts[bp[pointI]] = points[pointI];
2465     }
2466
2467     //- create the triangulation of the volume mesh surface
2468     forAll(bFaces, bfI)
2469     {
2470         const face& bf = bFaces[bfI];
2471
2472         labelledTri tri;
2473         tri.region() = facePatch_[bfI];
2474         tri[0] = bp[bf[0]];
2475
2476         for(label i=bf.size()-2;i>0;--i)
2477         {
2478             tri[1] = bp[bf[i]];
2479             tri[2] = bp[bf[i+1]];
2480
2481             surfPtr->appendTriangle(tri);
2482         }
2483     }
2484
2485     return surfPtr;
2486 }
2487
2488 const triSurf* edgeExtractor::surfaceWithPatches(const label bpI) const
2489 {
2490     //- allocate the memory for the surface mesh
2491     triSurf* surfPtr = new triSurf();
2492
2493     //- surface of the volume mesh
2494     const meshSurfaceEngine& mse = surfaceEngine();
2495     const faceList::subList& bFaces = mse.boundaryFaces();
2496     const VRWGraph& pFaces = mse.pointFaces();
2497     const pointFieldPMG& points = mesh_.points();
2498
2499     //- modifier of the new surface mesh
2500     triSurfModifier surfModifier(*surfPtr);
2501     surfModifier.patchesAccess() = meshOctree_.surface().patches();
2502     pointField& sPts = surfModifier.pointsAccess();
2503
2504     //- create the triangulation of the volume mesh surface
2505     labelLongList newPointLabel(points.size(), -1);
2506     label nPoints(0);
2507     forAllRow(pFaces, bpI, pfI)
2508     {
2509         const label bfI = pFaces(bpI, pfI);
2510         const face& bf = bFaces[bfI];
2511
2512         forAll(bf, pI)
2513             if( newPointLabel[bf[pI]] == -1 )
2514                 newPointLabel[bf[pI]] = nPoints++;
2515
2516         labelledTri tri;
2517         tri.region() = facePatch_[bfI];
2518         tri[0] = newPointLabel[bf[0]];
2519
2520         for(label i=bf.size()-2;i>0;--i)
2521         {
2522             tri[1] = newPointLabel[bf[i]];
2523             tri[2] = newPointLabel[bf[i+1]];
2524
2525             surfPtr->appendTriangle(tri);
2526         }
2527     }
2528
2529     //- copy points
2530     sPts.setSize(nPoints);
2531     forAll(newPointLabel, pointI)
2532         if( newPointLabel[pointI] != -1 )
2533         {
2534             sPts[newPointLabel[pointI]] = points[pointI];
2535         }
2536
2537     return surfPtr;
2538 }
2539
2540 void edgeExtractor::updateMeshPatches()
2541 {
2542     const triSurf& surface = meshOctree_.surface();
2543     const geometricSurfacePatchList& surfPatches = surface.patches();
2544     const label nPatches = surfPatches.size();
2545
2546     const meshSurfaceEngine& mse = this->surfaceEngine();
2547     const faceList::subList& bFaces = mse.boundaryFaces();
2548     const labelList& faceOwner = mse.faceOwners();
2549
2550     wordList patchNames(nPatches);
2551     VRWGraph newBoundaryFaces;
2552     labelLongList newBoundaryOwners(bFaces.size());
2553     labelLongList newBoundaryPatches(bFaces.size());
2554
2555     //- set patchNames
2556     forAll(surface.patches(), patchI)
2557         patchNames[patchI] = surface.patches()[patchI].name();
2558
2559     //- append boundary faces
2560     forAll(bFaces, bfI)
2561     {
2562         newBoundaryFaces.appendList(bFaces[bfI]);
2563         newBoundaryOwners[bfI] = faceOwner[bfI];
2564         newBoundaryPatches[bfI] = facePatch_[bfI];
2565     }
2566
2567     //- replace the boundary with the new patches
2568     polyMeshGenModifier meshModifier(mesh_);
2569     meshModifier.replaceBoundary
2570     (
2571         patchNames,
2572         newBoundaryFaces,
2573         newBoundaryOwners,
2574         newBoundaryPatches
2575     );
2576
2577     //- set the new patch types
2578     PtrList<boundaryPatch>& boundaries = meshModifier.boundariesAccess();
2579
2580     forAll(surfPatches, patchI)
2581         boundaries[patchI].patchType() = surfPatches[patchI].geometricType();
2582 }
2583
2584 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *//
2585
2586 } // End namespace Foam
2587
2588 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //