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