src/mesh/cfMesh/utilities/boundaryLayers/extrudeLayer/extrudeLayer.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 "extrudeLayer.H"
  29 #include "helperFunctions.H"
  30 #include "polyMeshGenAddressing.H"
  31 #include "meshSurfaceEngine.H"
  32 #include "meshSurfacePartitioner.H"
  33 #include "labelledPointScalar.H"
  34
  35 # ifdef USE_OMP
  36 #include <omp.h>
  37 # endif
  38
  39 #ifdef DEBUGExtrudeLayer
  40 #include "polyMeshGenChecks.H"
  41 #endif
  42
  43 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  44
  45 namespace Foam
  46 {
  47
  48 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  49 // Nested class definition
  50
  51 extrudeLayer::addressingCalculator::addressingCalculator
  52 (
  53     const faceListPMG& faces,
  54     const LongList<labelPair>& extrudedFaces,
  55     const LongList<bool>& pairOrientation,
  56     const VRWGraph& pointExtruded
  57 )
  58 :
  59     faces_(faces),
  60     extrudedFaces_(extrudedFaces),
  61     pairOrientation_(pairOrientation),
  62     pointExtruded_(pointExtruded)
  63 {}
  64
  65 extrudeLayer::addressingCalculator::~addressingCalculator()
  66 {}
  67
  68 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  69
  70 void extrudeLayer::createDuplicateFrontFaces(const LongList<labelPair>& front)
  71 {
  72     polyMeshGenModifier meshModifier(mesh_);
  73
  74     const labelList& owner = mesh_.owner();
  75     const labelList& neighbour = mesh_.neighbour();
  76     faceListPMG& faces = meshModifier.facesAccess();
  77     cellListPMG& cells = meshModifier.cellsAccess();
  78
  79     labelList faceInFront(faces.size(), -1);
  80     LongList<labelPair> newFaceLabels;
  81     label counter(0);
  82     forAll(front, ffI)
  83     {
  84         const labelPair& lp = front[ffI];
  85
  86         if( faceInFront[lp.first()] == -1 )
  87         {
  88             faceInFront[lp.first()] = newFaceLabels.size();
  89             newFaceLabels.append(labelPair(-1, -1));
  90         }
  91
  92         labelPair& cf = newFaceLabels[faceInFront[lp.first()]];
  93
  94         if( (lp.second() == owner[lp.first()]) && (cf.first() == -1) )
  95         {
  96             cf.first() = counter;
  97             ++counter;
  98         }
  99         else if( (lp.second() == neighbour[lp.first()]) && (cf.second() == -1) )
 100         {
 101             cf.second() = counter;
 102             ++counter;
 103         }
 104     }
 105
 106     //- create a copy of the faces in the extrusion front
 107     //- to create space for the layer
 108     faces.setSize(nOrigFaces_+counter);
 109     extrudedFaces_.setSize(counter);
 110     pairOrientation_.setSize(counter);
 111
 112     # ifdef USE_OMP
 113     # pragma omp parallel for if( faceInFront.size() > 100 ) schedule(guided)
 114     # endif
 115     forAll(faceInFront, faceI)
 116     {
 117         if( faceInFront[faceI] < 0 )
 118             continue;
 119
 120         const label fOwn = newFaceLabels[faceInFront[faceI]].first();
 121         const label fNei = newFaceLabels[faceInFront[faceI]].second();
 122
 123         if( neighbour[faceI] < 0 )
 124         {
 125             //- boundary faces
 126             if( fOwn != -1 )
 127             {
 128                 faces[nOrigFaces_+fOwn] = faces[faceI];
 129                 extrudedFaces_[fOwn] = labelPair(nOrigFaces_+fOwn, faceI);
 130                 pairOrientation_[fOwn] = true;
 131             }
 132         }
 133         else
 134         {
 135             //- internal face
 136             if( fOwn != -1 && fNei != -1 )
 137             {
 138                 if( fOwn < fNei )
 139                 {
 140                     faces[nOrigFaces_+fOwn] = faces[faceI];
 141                     extrudedFaces_[fOwn] = labelPair(nOrigFaces_+fOwn, faceI);
 142                     pairOrientation_[fOwn] = true;
 143
 144                     faces[nOrigFaces_+fNei] = faces[faceI].reverseFace();
 145                     extrudedFaces_[fNei] = labelPair(nOrigFaces_+fNei, faceI);
 146                     pairOrientation_[fNei] = false;
 147                 }
 148                 else
 149                 {
 150                     faces[nOrigFaces_+fOwn].transfer(faces[faceI]);
 151                     extrudedFaces_[fOwn] = labelPair(nOrigFaces_+fOwn, faceI);
 152                     pairOrientation_[fOwn] = false;
 153
 154                     faces[faceI] = faces[nOrigFaces_+fOwn].reverseFace();
 155
 156                     faces[nOrigFaces_+fNei] = faces[faceI];
 157                     extrudedFaces_[fNei] = labelPair(nOrigFaces_+fNei, faceI);
 158                     pairOrientation_[fNei] = true;
 159                 }
 160             }
 161             else if( fOwn != -1 )
 162             {
 163                 faces[nOrigFaces_+fOwn].transfer(faces[faceI]);
 164                 faces[faceI] = faces[nOrigFaces_+fOwn].reverseFace();
 165                 extrudedFaces_[fOwn] = labelPair(nOrigFaces_+fOwn, faceI);
 166                 pairOrientation_[fOwn] = false;
 167             }
 168             else if( fNei != -1 )
 169             {
 170                 faces[nOrigFaces_+fNei] = faces[faceI].reverseFace();
 171                 extrudedFaces_[fNei] = labelPair(nOrigFaces_+fNei, faceI);
 172                 pairOrientation_[fNei] = false;
 173             }
 174         }
 175     }
 176
 177     //- renumber the cells
 178     # ifdef USE_OMP
 179     # pragma omp parallel for if( faceInFront.size() > 100 ) schedule(guided)
 180     # endif
 181     forAll(faceInFront, faceI)
 182     {
 183         if( faceInFront[faceI] < 0 )
 184             continue;
 185
 186         const labelPair& lp = newFaceLabels[faceInFront[faceI]];
 187
 188         if( lp.first() != -1 )
 189         {
 190             cell& c = cells[owner[faceI]];
 191
 192             forAll(c, fI)
 193                 if( c[fI] == faceI )
 194                     c[fI] = nOrigFaces_ + lp.first();
 195         }
 196         if( lp.second() != -1 )
 197         {
 198             cell& c = cells[neighbour[faceI]];
 199
 200             forAll(c, fI)
 201                 if( c[fI] == faceI )
 202                     c[fI] = nOrigFaces_ + lp.second();
 203         }
 204     }
 205
 206     meshModifier.clearAll();
 207 }
 208
 209 void extrudeLayer::createNewVertices()
 210 {
 211     polyMeshGenModifier meshModifier(mesh_);
 212     meshModifier.clearAll();
 213     pointFieldPMG& points = meshModifier.pointsAccess();
 214     faceListPMG& faces = meshModifier.facesAccess();
 215
 216     const labelList& owner = mesh_.owner();
 217     const labelList& neighbour = mesh_.neighbour();
 218
 219     //- find the points in the marked front
 220     List<direction> frontPoints(points.size(), NONE);
 221
 222     # ifdef USE_OMP
 223     # pragma omp parallel for if( points.size() > 1000 ) schedule(guided)
 224     # endif
 225     forAll(extrudedFaces_, efI)
 226     {
 227         const face& f = faces[extrudedFaces_[efI].first()];
 228
 229         forAll(f, pI)
 230             frontPoints[f[pI]] |= FRONTVERTEX;
 231     }
 232
 233     //- propagate this information to other processors in case of a parallel run
 234     if( Pstream::parRun() )
 235     {
 236         const polyMeshGenAddressing& addr = mesh_.addressingData();
 237         const VRWGraph& pAtProcs = addr.pointAtProcs();
 238         const Map<label>& globalToLocal = addr.globalToLocalPointAddressing();
 239         const DynList<label>& pProcs = addr.pointNeiProcs();
 240
 241         //- allocate the map
 242         std::map<label, labelLongList> exchangeData;
 243         forAll(pProcs, i)
 244             exchangeData.insert(std::make_pair(pProcs[i], labelLongList()));
 245
 246         //- collect the information about markes points at processor boundaries
 247         forAllConstIter(Map<label>, globalToLocal, it)
 248             if( frontPoints[it()] & FRONTVERTEX )
 249             {
 250                 //- mark points at processor boundaries
 251                 frontPoints[it()] |= FRONTVERTEXPROCBND;
 252
 253                 forAllRow(pAtProcs, it(), i)
 254                 {
 255                     const label neiProc = pAtProcs(it(), i);
 256
 257                     if( neiProc == Pstream::myProcNo() )
 258                         continue;
 259
 260                     exchangeData[neiProc].append(it.key());
 261                 }
 262             }
 263
 264         LongList<label> receivedData;
 265         help::exchangeMap(exchangeData, receivedData);
 266
 267         # ifdef USE_OMP
 268         # pragma omp parallel for if( receivedData.size() > 1000 ) \
 269         schedule(guided)
 270         # endif
 271         forAll(receivedData, i)
 272         {
 273             frontPoints[globalToLocal[receivedData[i]]] =
 274                 FRONTVERTEX+FRONTVERTEXPROCBND;
 275         }
 276     }
 277
 278     //- create a new vertex for each group of faces
 279     //- faces are grouped such that the faces belonging to the same group
 280     //- can be visited over cells without crossing any front faces
 281     VRWGraph pointFaces;
 282     pointFaces.reverseAddressing(points.size(), faces);
 283
 284     if( Pstream::parRun() )
 285     {
 286         Pout << "Creating new points at processor boundaries" << endl;
 287         for(label procI=0;procI<Pstream::nProcs();++procI)
 288         {
 289             if( Pstream::myProcNo() == procI )
 290             {
 291                Pout << "Front points are " << frontPoints << endl;
 292             }
 293
 294             returnReduce(1, sumOp<label>());
 295         }
 296
 297         //- create new vertices at processor boundaries
 298         const PtrList<processorBoundaryPatch>& procBoundaries =
 299             mesh_.procBoundaries();
 300         const polyMeshGenAddressing& addr = mesh_.addressingData();
 301         const VRWGraph& pAtProcs = addr.pointAtProcs();
 302         const labelLongList& globalPointLabel = addr.globalPointLabel();
 303         const Map<label>& globalToLocal = addr.globalToLocalPointAddressing();
 304         const DynList<label>& pProcs = addr.pointNeiProcs();
 305         const labelLongList& globalCellLabel = addr.globalCellLabel();
 306
 307         //- create the information which faces are attached to points
 308         //- at parallel boundaries in dual form where each edge represents
 309         //- global labels of cells sharing a face
 310         typedef std::map<label, DynList<edge> > dualEdgesMap;
 311         dualEdgesMap procPointsDual;
 312
 313         //- fill in local data
 314         forAllConstIter(Map<label>, globalToLocal, it)
 315         {
 316             if( frontPoints[it()] & FRONTVERTEXPROCBND )
 317             {
 318                 const label pointI = it();
 319                 DynList<edge>& pEdges = procPointsDual[pointI];
 320
 321                 forAllRow(pointFaces, pointI, pfI)
 322                 {
 323                     const label faceI = pointFaces(pointI, pfI);
 324
 325                     //- do not store faces at processor boundaries
 326                     //- and newly created faces
 327                     if( (neighbour[faceI] < 0) || (owner[faceI] < 0) )
 328                         continue;
 329
 330                     const edge e
 331                     (
 332                         globalCellLabel[owner[faceI]],
 333                         globalCellLabel[neighbour[faceI]]
 334                     );
 335
 336                     pEdges.append(e);
 337                 }
 338             }
 339         }
 340
 341         for(label procI=0;procI<Pstream::nProcs();++procI)
 342         {
 343             if( Pstream::myProcNo() == procI )
 344             {
 345                forAllConstIter(dualEdgesMap, procPointsDual, it)
 346                     Pout << "Point " << it->first << " local dual edges " << it->second << endl;
 347             }
 348
 349             returnReduce(1, sumOp<label>());
 350         }
 351
 352         //- fill-in with data at processor boundaries. Store edges
 353         //- on the processor with the lower label not to duplicate the data
 354         returnReduce(1, sumOp<label>());
 355         Pout << "Adding data from processor boundaries" << endl;
 356         forAll(procBoundaries, patchI)
 357         {
 358             if( procBoundaries[patchI].owner() )
 359                 continue;
 360
 361             const label start = procBoundaries[patchI].patchStart();
 362             labelList globalLabels(procBoundaries[patchI].patchSize());
 363             forAll(globalLabels, fI)
 364             {
 365                 if( owner[start+fI] < 0 )
 366                 {
 367                     globalLabels[fI] = -1;
 368                 }
 369                 else
 370                 {
 371                     globalLabels[fI] = globalCellLabel[owner[start+fI]];
 372                 }
 373             }
 374
 375             OPstream toOtherProc
 376             (
 377                 Pstream::blocking,
 378                 procBoundaries[patchI].neiProcNo(),
 379                 globalLabels.byteSize()
 380             );
 381
 382             toOtherProc << globalLabels;
 383         }
 384
 385         forAll(procBoundaries, patchI)
 386         {
 387             if( !procBoundaries[patchI].owner() )
 388                 continue;
 389
 390             labelList receivedData;
 391             IPstream fromOtherProc
 392             (
 393                 Pstream::blocking,
 394                 procBoundaries[patchI].neiProcNo()
 395             );
 396
 397             fromOtherProc >> receivedData;
 398
 399             const label start = procBoundaries[patchI].patchStart();
 400
 401             forAll(receivedData, i)
 402             {
 403                 if( owner[start+i] < 0 )
 404                     continue;
 405
 406                 const face& f = faces[start+i];
 407
 408                 forAll(f, pI)
 409                 {
 410                     if( frontPoints[f[pI]] & FRONTVERTEXPROCBND )
 411                     {
 412                         dualEdgesMap::iterator dIter =
 413                             procPointsDual.find(f[pI]);
 414                         const label cOwn = globalCellLabel[owner[start+i]];
 415                         const label cNei = receivedData[i];
 416                         dIter->second.append(edge(cOwn, cNei));
 417                     }
 418                 }
 419             }
 420         }
 421
 422         //- exchange this information with neighbouring processors
 423         returnReduce(1, sumOp<label>());
 424         Pout << "Exchanging data with other processors" << endl;
 425
 426         std::map<label, labelLongList> exchangeData;
 427         forAll(pProcs, i)
 428             exchangeData.insert(std::make_pair(pProcs[i], labelLongList()));
 429
 430         //- fill in the exchangeData map
 431         forAllConstIter(dualEdgesMap, procPointsDual, dIter)
 432         {
 433             const label pointI = dIter->first;
 434
 435             forAllRow(pAtProcs, pointI, i)
 436             {
 437                 const label neiProc = pAtProcs(pointI, i);
 438
 439                 if( neiProc == Pstream::myProcNo() )
 440                     continue;
 441
 442                 labelLongList& dts = exchangeData[neiProc];
 443
 444                 dts.append(globalPointLabel[pointI]);
 445
 446                 const DynList<edge>& dualEdges = dIter->second;
 447                 dts.append(dualEdges.size());
 448                 forAll(dualEdges, edgeI)
 449                 {
 450                     const edge& e = dualEdges[edgeI];
 451                     dts.append(e.start());
 452                     dts.append(e.end());
 453                 }
 454             }
 455         }
 456
 457         //- exchange data with other processors
 458         labelLongList receivedData;
 459         help::exchangeMap(exchangeData, receivedData);
 460
 461         //- update local data
 462         label counter(0);
 463         while( counter < receivedData.size() )
 464         {
 465             const label pointI = globalToLocal[receivedData[counter++]];
 466
 467             DynList<edge>& dualEdges = procPointsDual[pointI];
 468
 469             const label nDualEdges = receivedData[counter++];
 470             for(label eI=0;eI<nDualEdges;++eI)
 471             {
 472                 edge e;
 473                 e.start() = receivedData[counter++];
 474                 e.end() = receivedData[counter++];
 475
 476                 dualEdges.append(e);
 477             }
 478         }
 479
 480         for(label procI=0;procI<Pstream::nProcs();++procI)
 481         {
 482             if( Pstream::myProcNo() == procI )
 483             {
 484                forAllConstIter(dualEdgesMap, procPointsDual, it)
 485                     Pout << "Point " << it->first << " dual edges " << it->second << endl;
 486             }
 487
 488             returnReduce(1, sumOp<label>());
 489         }
 490
 491         //- Finally, find groups of faces and create new vertices
 492         returnReduce(1, sumOp<label>());
 493         Pout << "Finding groups of edges at vertex" << endl;
 494         forAllConstIter(dualEdgesMap, procPointsDual, dIter)
 495         {
 496             const label pointI = dIter->first;
 497             const DynList<edge>& dEdges = dIter->second;
 498
 499             DynList<label> edgeGroup;
 500             edgeGroup.setSize(dEdges.size());
 501             edgeGroup = -1;
 502
 503             //- check edge connections and store all edges which can be reached
 504             //- over other edges into the same group
 505             label group(0);
 506
 507             forAll(dEdges, eI)
 508             {
 509                 if( edgeGroup[eI] != -1 )
 510                     continue;
 511
 512                 edgeGroup[eI] = group;
 513                 DynList<label> front;
 514                 front.append(eI);
 515
 516                 while( front.size() != 0 )
 517                 {
 518                     const label eLabel = front.removeLastElement();
 519
 520                     forAll(dEdges, eJ)
 521                     {
 522                         if( edgeGroup[eJ] != -1 )
 523                             continue;
 524                         if( eJ == eLabel )
 525                             continue;
 526
 527                         if( dEdges[eLabel].commonVertex(dEdges[eJ]) != -1 )
 528                         {
 529                             front.append(eJ);
 530                             edgeGroup[eJ] = group;
 531                         }
 532                     }
 533                 }
 534
 535                 ++group;
 536             }
 537
 538             //- find face groups from the groups assigned to dual edges
 539             DynList<label> faceGroups;
 540             faceGroups.setSize(pointFaces.sizeOfRow(pointI));
 541             faceGroups = -1;
 542
 543             forAllRow(pointFaces, pointI, pfI)
 544             {
 545                 const label faceI = pointFaces(pointI, pfI);
 546
 547                 if( owner[faceI] < 0 )
 548                     continue;
 549
 550                 const label own = globalCellLabel[owner[faceI]];
 551
 552                 forAll(dEdges, eI)
 553                 {
 554                     const edge& dualEdge = dEdges[eI];
 555
 556                     if( (own == dualEdge.start()) || (own == dualEdge.end()) )
 557                     {
 558                         faceGroups[pfI] = edgeGroup[eI];
 559                         break;
 560                     }
 561                 }
 562             }
 563
 564             Info << "Edge groups for point " << pointI << " are " << edgeGroup << endl;
 565             Info << "Face groups at point " << pointI << " are " << faceGroups
 566                 << " point faces " << pointFaces[pointI] << endl;
 567
 568             //- stop in case there is only one group
 569             //- of faces attached to this point
 570             if( group == 1 )
 571             {
 572                 bool problem(true);
 573                 forAll(faceGroups, i)
 574                     if( faceGroups[i] != 0 )
 575                     {
 576                         problem = false;
 577                         break;
 578                     }
 579
 580                 if( problem )
 581                     FatalErrorIn("void extrudeLayer::createNewVertices()")
 582                         << "Front is not defined at point " << pointI
 583                         << ". Cannot continue.." << exit(FatalError);
 584             }
 585
 586             //- create new vertices
 587             const label currentNumPoints = points.size();
 588             for(label i=0;i<group;++i)
 589             {
 590                 const point p = points[pointI];
 591                 origPointLabel_.append(pointI);
 592                 points.append(p);
 593             }
 594
 595             //- renumber faces
 596             forAllRow(pointFaces, pointI, pfI)
 597             {
 598                 if( faceGroups[pfI] == -1 )
 599                     continue;
 600
 601                 face& f = faces[pointFaces(pointI, pfI)];
 602
 603                 const label pos = f.which(pointI);
 604                 f[pos] = currentNumPoints + faceGroups[pfI];
 605             }
 606         }
 607
 608         Pout << "Finishing creating new vertices at paralle boundaries" << endl;
 609         returnReduce(1, sumOp<label>());
 610     }
 611
 612     //- treat local points
 613     forAll(pointFaces, pointI)
 614     {
 615         if( frontPoints[pointI] != FRONTVERTEX )
 616             continue;
 617
 618         //- assign groups to faces and cells
 619         DynList<label> faceGroup;
 620         faceGroup.setSize(pointFaces.sizeOfRow(pointI));
 621         faceGroup = -1;
 622
 623         label group(0);
 624
 625         forAllRow(pointFaces, pointI, pfI)
 626         {
 627             if( pointFaces(pointI, pfI) < nOrigFaces_ )
 628                 continue;
 629             if( faceGroup[pfI] != -1 )
 630                 continue;
 631
 632             DynList<label> front;
 633             front.append(pfI);
 634             faceGroup[pfI] = group;
 635
 636             while( front.size() )
 637             {
 638                 const label fLabel = front.removeLastElement();
 639
 640                 const label own = owner[pointFaces(pointI, fLabel)];
 641                 const label nei = neighbour[pointFaces(pointI, fLabel)];
 642
 643                 forAllRow(pointFaces, pointI, pfJ)
 644                 {
 645                     if( faceGroup[pfJ] != -1 )
 646                         continue;
 647
 648                     const label faceJ = pointFaces(pointI, pfJ);
 649
 650                     if( owner[faceJ] < 0 )
 651                         continue;
 652
 653                     if( owner[faceJ] == own || owner[faceJ] == nei )
 654                     {
 655                         faceGroup[pfJ] = group;
 656                         front.append(pfJ);
 657                     }
 658
 659                     if( neighbour[faceJ] < 0 )
 660                         continue;
 661
 662                     if( neighbour[faceJ] == own || neighbour[faceJ] == nei )
 663                     {
 664                         faceGroup[pfJ] = group;
 665                         front.append(pfJ);
 666                     }
 667                 }
 668             }
 669
 670             ++group;
 671         }
 672
 673         //- stop in case there is only one group of faces attached to this point
 674         if( group == 1 )
 675         {
 676             bool problem(true);
 677             forAll(faceGroup, i)
 678                 if( faceGroup[i] != 0 )
 679                 {
 680                     problem = false;
 681                     break;
 682                 }
 683
 684             if( problem )
 685                 FatalErrorIn("void extrudeLayer::createNewVertices()")
 686                     << "Front is not defined at point " << pointI
 687                     << ". Cannot continue.." << exit(FatalError);
 688         }
 689
 690         //- create new vertices
 691         const label currentNumPoints = points.size();
 692         for(label i=0;i<group;++i)
 693         {
 694             const point p = points[pointI];
 695             origPointLabel_.append(pointI);
 696             points.append(p);
 697         }
 698
 699         //- renumber faces
 700         forAllRow(pointFaces, pointI, pfI)
 701         {
 702             if( faceGroup[pfI] == -1 )
 703                 continue;
 704
 705             face& f = faces[pointFaces(pointI, pfI)];
 706
 707             const label pos = f.which(pointI);
 708             f[pos] = currentNumPoints + faceGroup[pfI];
 709         }
 710     }
 711
 712     polyMeshGenModifier(mesh_).clearOut();
 713
 714     # ifdef DEBUGExtrudeLayer
 715     const label procID = mesh_.addPointSubset("parPoints");
 716     const label frontID = mesh_.addPointSubset("frontPoints");
 717     forAll(frontPoints, pI)
 718     {
 719         if( frontPoints[pI] & FRONTVERTEXPROCBND )
 720             mesh_.addPointToSubset(procID, pI);
 721         if( frontPoints[pI] & FRONTVERTEX )
 722             mesh_.addPointToSubset(frontID, pI);
 723     }
 724
 725     returnReduce(1, sumOp<label>());
 726     //::exit(1);
 727     # endif
 728 }
 729
 730 void extrudeLayer::movePoints()
 731 {
 732     pointFieldPMG& points = mesh_.points();
 733     const faceListPMG& faces = mesh_.faces();
 734
 735     vectorField displacements(points.size()-nOrigPoints_);
 736     boolList pointAtProcBnd(displacements.size(), false);
 737
 738     VRWGraph pointFaces;
 739     pointFaces.reverseAddressing(points.size(), mesh_.faces());
 740
 741     if( Pstream::parRun() )
 742     {
 743         const polyMeshGenAddressing& addr = mesh_.addressingData();
 744         const Map<label>& globalToLocal = addr.globalToLocalPointAddressing();
 745         const VRWGraph& pAtProcs = addr.pointAtProcs();
 746         const DynList<label>& pProcs = addr.pointNeiProcs();
 747
 748         std::map<label, vector> normals;
 749         std::map<label, scalar> distances;
 750
 751         //- create a map for exchanging data
 752         std::map<label, LongList<labelledPointScalar> > exchangeData;
 753         forAll(pProcs, i)
 754             exchangeData.insert
 755             (
 756                 std::make_pair(pProcs[i], LongList<labelledPointScalar>())
 757             );
 758
 759         //- create displacements from local data
 760         forAllConstIter(Map<label>, globalToLocal, it)
 761         {
 762             if( it() >= nOrigPoints_ )
 763             {
 764                 const label pointI = it();
 765
 766                 pointAtProcBnd[pointI-nOrigPoints_] = true;
 767
 768                 //- create local part of the displacement vector
 769                 labelledPointScalar lps;
 770                 lps.pointLabel() = it.key();
 771                 lps.coordinates() = vector::zero;
 772                 lps.scalarValue() = VGREAT;
 773
 774                 forAllRow(pointFaces, pointI, pfI)
 775                 {
 776                     const label faceI = pointFaces(pointI, pfI);
 777
 778                     if( faceI >= nOrigFaces_ )
 779                     {
 780                         const face& f = faces[faceI];
 781
 782                         lps.coordinates() -= f.normal(points);
 783
 784                         if( thickness_ < 0.0 )
 785                         {
 786                             const vector c = f.centre(points);
 787                             scalar d(VGREAT);
 788
 789                             forAll(f, pI)
 790                                 d = Foam::min(d, Foam::mag(points[f[pI]] - c));
 791
 792                             d *= 0.4;
 793                             lps.scalarValue() = Foam::min(lps.scalarValue(), d);
 794                         }
 795                         else
 796                         {
 797                             lps.scalarValue() = thickness_;
 798                         }
 799                     }
 800                 }
 801
 802                 normals[pointI] = lps.coordinates();
 803                 distances[pointI] = lps.scalarValue();
 804
 805                 //- store data in the exchangeData map
 806                 forAllRow(pAtProcs, pointI, i)
 807                 {
 808                     const label neiProc = pAtProcs(pointI, i);
 809
 810                     if( neiProc == Pstream::myProcNo() )
 811                         continue;
 812
 813                     exchangeData[neiProc].append(lps);
 814                 }
 815             }
 816         }
 817
 818         LongList<labelledPointScalar> receivedData;
 819         help::exchangeMap(exchangeData, receivedData);
 820
 821         forAll(receivedData, i)
 822         {
 823             const labelledPointScalar& lps = receivedData[i];
 824             const label pointI = globalToLocal[lps.pointLabel()];
 825
 826             normals[pointI] += lps.coordinates();
 827             scalar& d = distances[pointI];
 828             d = Foam::min(d, lps.scalarValue());
 829         }
 830
 831         //- calculate displacements of vertices at processor boundaries
 832         for
 833         (
 834             std::map<label, vector>::const_iterator it=normals.begin();
 835             it!=normals.end();
 836             ++it
 837         )
 838         {
 839             vector n = it->second;
 840             if( mag(n) > VSMALL )
 841                 n /= mag(n);
 842             displacements[it->first-nOrigPoints_] = n * distances[it->first];
 843         }
 844     }
 845
 846     # ifdef USE_OMP
 847     # pragma omp parallel if( displacements.size() > 100 )
 848     # endif
 849     {
 850         //- find displacement vectors
 851         # ifdef USE_OMP
 852         # pragma omp for schedule(guided)
 853         # endif
 854         forAll(displacements, pI)
 855         {
 856             if( pointAtProcBnd[pI] )
 857                 continue;
 858
 859             const label pointI = nOrigPoints_ + pI;
 860
 861             vector normal(vector::zero);
 862             scalar thickness(VGREAT);
 863
 864             forAllRow(pointFaces, pointI, pfI)
 865             {
 866                 const label faceI = pointFaces(pointI, pfI);
 867
 868                 if( faceI >= nOrigFaces_ )
 869                 {
 870                     const face& f = faces[faceI];
 871
 872                     normal -= f.normal(points);
 873
 874                     if( thickness_ < 0.0 )
 875                     {
 876                         const vector c = f.centre(points);
 877                         scalar d(VGREAT);
 878
 879                         forAll(f, pI)
 880                             d = Foam::min(d, Foam::mag(points[f[pI]] - c));
 881
 882                         thickness = Foam::min(thickness, d);
 883                     }
 884                 }
 885             }
 886
 887             if( thickness_ >= 0.0 )
 888             {
 889                 thickness = thickness_;
 890             }
 891             else
 892             {
 893                 thickness *= 0.4;
 894             }
 895
 896             const scalar d = mag(normal);
 897             if( d > VSMALL )
 898                 normal /= d;
 899
 900             displacements[pI] = normal * thickness;
 901         }
 902
 903         # ifdef USE_OMP
 904         # pragma omp barrier
 905         # endif
 906
 907         # ifdef USE_OMP
 908         # pragma omp for schedule(guided)
 909         # endif
 910         forAll(displacements, pI)
 911             points[nOrigPoints_+pI] += displacements[pI];
 912     }
 913
 914     # ifdef DEBUGExtrudeLayer
 915     mesh_.write();
 916     returnReduce(1, sumOp<label>());
 917     //::exit(1);
 918     # endif
 919 }
 920
 921 void extrudeLayer::createNewFacesParallel()
 922 {
 923     if( !Pstream::parRun() )
 924         return;
 925
 926     VRWGraph newProcFaces;
 927     labelLongList faceProcPatch;
 928
 929     //- add faces into the mesh
 930     polyMeshGenModifier(mesh_).addProcessorFaces(newProcFaces, faceProcPatch);
 931 }
 932
 933 void extrudeLayer::createLayerCells()
 934 {
 935     const faceListPMG& faces = mesh_.faces();
 936
 937     VRWGraphList newCells;
 938
 939     //- create cells from corners and edges
 940     faceList::subList newFaces(faces, faces.size()-nOrigFaces_, nOrigFaces_);
 941     VRWGraph pointFaces;
 942     pointFaces.reverseAddressing(mesh_.points().size(), newFaces);
 943
 944     //- create new cells extruded from the faces in the selected front
 945     forAll(extrudedFaces_, extrudedI)
 946     {
 947         const face& f = faces[extrudedFaces_[extrudedI].first()];
 948         const face& of = faces[extrudedFaces_[extrudedI].second()];
 949
 950         //- create new cell from the front pair
 951         DynList<DynList<label, 4> > newCell;
 952         newCell.setSize(f.size()+2);
 953
 954         newCell[0] = f;
 955         newCell[1] = of;
 956
 957         if( pairOrientation_[extrudedI] )
 958         {
 959             //- create a new cell. Faces are of the same orientation
 960             forAll(f, pI)
 961             {
 962                 DynList<label, 4>& ef = newCell[pI+2];
 963
 964                 ef.setSize(4);
 965                 ef[0] = f[pI];
 966                 ef[1] = f.nextLabel(pI);
 967                 ef[2] = of[f.fcIndex(pI)];
 968                 ef[3] = of[pI];
 969             }
 970         }
 971         else
 972         {
 973             //- create a new cell. Faces are of the opposite orientation
 974             forAll(f, pI)
 975             {
 976                 DynList<label, 4>& ef = newCell[pI+2];
 977
 978                 ef.setSize(4);
 979                 ef[0] = f[pI];
 980                 ef[1] = f.nextLabel(pI);
 981                 ef[2] = of[(f.size()-f.fcIndex(pI))%f.size()];
 982                 ef[3] = of[(f.size()-pI)%f.size()];
 983             }
 984         }
 985
 986         newCells.appendGraph(newCell);
 987     }
 988
 989     addressingCalculator adc
 990     (
 991         faces,
 992         extrudedFaces_,
 993         pairOrientation_,
 994         pointFaces
 995     );
 996
 997     //- create cells created as a consequence of self-intersection
 998     //- over edges. And edge is transformed into a hex cell
 999     forAll(extrudedFaces_, extrudedI)
1000     {
1001         const face& f = faces[extrudedFaces_[extrudedI].first()];
1002
1003         forAll(f, eI)
1004         {
1005             const label pointI = f[eI];
1006             if( pointFaces.sizeOfRow(pointI) == 0 )
1007                 continue;
1008             const label nextI = f.nextLabel(eI);
1009             if( pointFaces.sizeOfRow(nextI) == 0 )
1010                 continue;
1011             if( pointI < nextI )
1012                 continue;
1013
1014             //- find point labels of the edge which is the origin
1015             //- of the edge (pointI, nextI) with respect to face extrudedI
1016             const label origFacePointI = adc.origPoint(extrudedI, pointI);
1017             const label origFaceNextI = adc.origPoint(extrudedI, nextI);
1018
1019             //- points of the current edge and of the original edge must have
1020             //- the same point as their origin
1021             if( origPointLabel_[pointI] != origPointLabel_[origFacePointI] )
1022                 continue;
1023             if( origPointLabel_[nextI] != origPointLabel_[origFaceNextI] )
1024                 continue;
1025
1026             //- Finally, start creating a cell from this edge
1027             //- find the other extruded face which shares this edge
1028             const label otherExtI = adc.faceSharingEdge(extrudedI, eI);
1029
1030             //- find point labels of the edge which is the origin
1031             //- of the edge (pointI, nextI)
1032             //- with respect to face otherExtrudedFace
1033             const label otherOrigFacePointI = adc.origPoint(otherExtI, pointI);
1034             const label otherOrigFaceNextI = adc.origPoint(otherExtI, nextI);
1035
1036             //- find points of the edge opposite to the edge (pointI, nextI)
1037             //- these points make the original edge which is blown into
1038             //- a hex cell
1039             label origPointI(-1), origNextI(-1);
1040
1041             if(
1042                 (origFacePointI >= nOrigPoints_) &&
1043                 (origFaceNextI >= nOrigPoints_)
1044             )
1045             {
1046                 //- find an extruded face attached to edge
1047                 //- (origFacePointI, origFaceNextI). The must exist only one
1048                 //- such face
1049                 DynList<label> fc;
1050                 adc.facesSharingEdge(origFacePointI, origFaceNextI, fc);
1051
1052                 if( fc.size() != 1 )
1053                     FatalErrorIn("void extrudeLayer::createLayerCells()")
1054                         << "Cannot find searched face" << abort(FatalError);
1055
1056                 origPointI = adc.origPoint(fc[0], origFacePointI);
1057                 origNextI = adc.origPoint(fc[0], origFaceNextI);
1058             }
1059             else
1060             {
1061                 FatalErrorIn("void extrudeLayer::createLayerCells()")
1062                     << "Cannot find points" << abort(FatalError);
1063             }
1064
1065             //- create new cell and add it to the list
1066             FixedList<FixedList<label, 4>, 6> newCell;
1067
1068             //- face 0
1069             newCell[0][0] = pointI;
1070             newCell[0][1] = origFacePointI;
1071             newCell[0][2] = origFaceNextI;
1072             newCell[0][3] = nextI;
1073
1074             //- face 1
1075             newCell[1][0] = pointI;
1076             newCell[1][1] = nextI;
1077             newCell[1][2] = otherOrigFaceNextI;
1078             newCell[1][3] = otherOrigFacePointI;
1079
1080             //- face 2
1081             newCell[2][0] = otherOrigFacePointI;
1082             newCell[2][1] = otherOrigFaceNextI;
1083             newCell[2][2] = origNextI;
1084             newCell[2][3] = origPointI;
1085
1086             //- face 3
1087             newCell[3][0] = origFacePointI;
1088             newCell[3][1] = origPointI;
1089             newCell[3][2] = origNextI;
1090             newCell[3][3] = origFaceNextI;
1091
1092             //- face 4
1093             newCell[4][0] = pointI;
1094             newCell[4][1] = otherOrigFacePointI;
1095             newCell[4][2] = origPointI;
1096             newCell[4][3] = origFacePointI;
1097
1098             //- face 5
1099             newCell[5][0] = nextI;
1100             newCell[5][1] = origFaceNextI;
1101             newCell[5][2] = origNextI;
1102             newCell[5][3] = otherOrigFaceNextI;
1103
1104             newCells.appendGraph(newCell);
1105         }
1106     }
1107
1108     //- create cells at points where three or more self-intersecting layers meet
1109     forAll(pointFaces, pointI)
1110     {
1111         if( pointFaces.sizeOfRow(pointI) < 3 )
1112             continue;
1113
1114         //- find labels of points
1115         DynList<label> origFacePoints;
1116         origFacePoints.setSize(pointFaces.sizeOfRow(pointI));
1117         origFacePoints = -1;
1118
1119         forAllRow(pointFaces, pointI, pfI)
1120         {
1121             const label extI = pointFaces(pointI, pfI);
1122
1123             origFacePoints[pfI] = adc.origPoint(extI, pointI);
1124         }
1125
1126         //- cell shall be created only if all points are different
1127         bool createCell(true);
1128         for(label i=origFacePoints.size()-1;i>0;--i)
1129         {
1130             for(label j=i-1;j>=0;--j)
1131                 if( origFacePoints[i] == origFacePoints[j] )
1132                 {
1133                     createCell = false;
1134                     break;
1135                 }
1136
1137             if( !createCell )
1138                 break;
1139         }
1140
1141         if( !createCell )
1142             continue;
1143
1144         DynList<FixedList<label, 4>, 6> newCell;
1145         newCell.setSize(2*origFacePoints.size());
1146
1147         //- start creating faces attached to the pointI
1148         //- this is performed by finding original points with respect to
1149         //- the face extI and the face which shares edge eI of the face extI
1150         //- this is needed to ensure correct normal orientation
1151         forAllRow(pointFaces, pointI, pfI)
1152         {
1153             const label extI = pointFaces(pointI, pfI);
1154
1155             //- find original labels of points making a forward circular edge
1156             //- with respect to pointI
1157             const face& f = faces[extrudedFaces_[extI].first()];
1158             const label eI = f.which(pointI);
1159             const label origFacePointI = origFacePoints[pfI];
1160             const label origFaceNextI = adc.origPointLabel(extI, f.fcIndex(eI));
1161
1162             //- find another face attached to the edge eI
1163             const label fFace = adc.faceSharingEdge(extI, eI);
1164             const label pos = pointFaces.containsAtPosition(pointI, fFace);
1165
1166             //- find an extrusion face attached to the edge consisting of points
1167             //- origFacePointI and origFaceNextI
1168             DynList<label> fe;
1169             adc.facesSharingEdge(origFacePointI, origFaceNextI, fe);
1170             const label origPointI = adc.origPoint(fe[0], origFacePointI);
1171
1172             //- create a face attached to pointI
1173             FixedList<label, 4>& cf = newCell[pfI];
1174             cf[0] = pointI;
1175             cf[1] = origFacePoints[pfI];
1176             cf[2] = origPointI;
1177             cf[3] = origFacePoints[pos];
1178         }
1179
1180         //- close the cell by creating new faces from the existing
1181         //- faces which obey pre-determined order. If a face contains
1182         //- a point in the origFacePoints list at the second position, then
1183         //- the point after shall be the previous point of the face. If a face
1184         //- contains such point at the last position then the point before it
1185         //- shall be the next point in the face. The last point of the face
1186         //- is the original points from which all these points were generated.
1187         forAll(origFacePoints, opI)
1188         {
1189             FixedList<label, 4>& cf = newCell[opI+origFacePoints.size()];
1190
1191             //- find previous and next point
1192             label prev(-1), next(-1);
1193             forAll(origFacePoints, fJ)
1194             {
1195                 if( newCell[fJ][1] == origFacePoints[opI] )
1196                     prev = newCell[fJ][2];
1197                 if( newCell[fJ][3] == origFacePoints[opI] )
1198                     next = newCell[fJ][2];
1199             }
1200
1201             if( (prev < 0) || (next < 0) )
1202                 FatalErrorIn("void extrudeLayer::createLayerCells()")
1203                     << "Corner cell is invalid" << abort(FatalError);
1204
1205             cf[0] = prev;
1206             cf[1] = origFacePoints[opI];
1207             cf[2] = next;
1208             cf[3] = origPointLabel_[pointI];
1209         }
1210
1211         newCells.appendGraph(newCell);
1212     }
1213
1214     //- create new faces at parallel boundaries
1215     createNewFacesParallel();
1216
1217     //- add cells into the mesh
1218     polyMeshGenModifier(mesh_).addCells(newCells);
1219 }
1220
1221 void extrudeLayer::updateBoundary()
1222 {
1223     wordList patchNames(mesh_.boundaries().size());
1224     wordList patchTypes(mesh_.boundaries().size());
1225     forAll(patchNames, patchI)
1226     {
1227         patchNames[patchI] = mesh_.boundaries()[patchI].patchName();
1228         patchTypes[patchI] = mesh_.boundaries()[patchI].patchType();
1229     }
1230
1231     VRWGraph newBoundaryFaces;
1232     labelLongList newBoundaryOwners;
1233     labelLongList newBoundaryPatches;
1234
1235     meshSurfaceEngine mse(mesh_);
1236     const faceList::subList& bFaces = mse.boundaryFaces();
1237     const labelList& bfOwner = mse.faceOwners();
1238     const labelList& facePatch = mse.boundaryFacePatches();
1239     const labelList& bp = mse.bp();
1240
1241     meshSurfacePartitioner mPart(mse);
1242     const VRWGraph& pointPatches = mPart.pointPatches();
1243
1244     //- store existing boundary faces. They remain in the mesh
1245     forAll(bFaces, bfI)
1246     {
1247         newBoundaryFaces.appendList(bFaces[bfI]);
1248         newBoundaryOwners.append(bfOwner[bfI]);
1249         newBoundaryPatches.append(facePatch[bfI]);
1250     }
1251
1252     //- find and store boundary faces which have been generated as a consequence
1253     //- of layer insertion
1254     const faceListPMG& faces = mesh_.faces();
1255     const labelList& owner = mesh_.owner();
1256     const labelList& neighbour = mesh_.neighbour();
1257
1258     for(label faceI=nOrigFaces_;faceI<faces.size();++faceI)
1259     {
1260         if( neighbour[faceI] >= 0 )
1261             continue;
1262
1263         const face& f = faces[faceI];
1264
1265         if( f.size() != 4 )
1266             continue;
1267
1268         DynList<label> origPts, newPts;
1269         forAll(f, pI)
1270             if( origPointLabel_[f[pI]] < 0 )
1271             {
1272                 origPts.appendIfNotIn(f[pI]);
1273             }
1274             else
1275             {
1276                 origPts.appendIfNotIn(origPointLabel_[f[pI]]);
1277                 newPts.append(f[pI]);
1278             }
1279
1280         if( origPts.size() > 2 )
1281             continue;
1282         if( newPts.size() == 0 )
1283             continue;
1284
1285         //- this face belongs to the extruded layer
1286         //- find patches of boundary faces attached to the newly created points
1287         DynList<label> patches;
1288         DynList<label> nAppearances;
1289         forAll(newPts, npI)
1290         {
1291             const label bpI = bp[newPts[npI]];
1292
1293             if( bpI < 0 )
1294                 continue;
1295
1296             forAllRow(pointPatches, bpI, ppI)
1297             {
1298                 const label patchI = pointPatches(bpI, ppI);
1299                 if( patches.contains(patchI) )
1300                 {
1301                     ++nAppearances[patches.containsAtPosition(patchI)];
1302                 }
1303                 else
1304                 {
1305                     patches.append(patchI);
1306                     nAppearances.append(1);
1307                 }
1308             }
1309         }
1310
1311         label patch(-1);
1312
1313         forAll(patches, i)
1314         {
1315             if( nAppearances[i] == newPts.size() )
1316             {
1317                 if( patch != -1 )
1318                     FatalErrorIn("void extrudeLayer::updateBoundary()")
1319                         << "Found more than one patch!" << abort(FatalError);
1320
1321                 patch = patches[i];
1322             }
1323         }
1324
1325         if( patch != -1 )
1326         {
1327             //- append boundary face
1328             newBoundaryFaces.appendList(f);
1329             newBoundaryOwners.append(owner[faceI]);
1330             newBoundaryPatches.append(patch);
1331         }
1332     }
1333
1334     polyMeshGenModifier meshModifier(mesh_);
1335     meshModifier.reorderBoundaryFaces();
1336     meshModifier.replaceBoundary
1337     (
1338         patchNames,
1339         newBoundaryFaces,
1340         newBoundaryOwners,
1341         newBoundaryPatches
1342     );
1343
1344     PtrList<boundaryPatch>& boundaries = meshModifier.boundariesAccess();
1345     forAll(boundaries, patchI)
1346         boundaries[patchI].patchType() = patchTypes[patchI];
1347
1348     meshModifier.clearAll();
1349 }
1350
1351 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
1352
1353 // Construct from mesh reference
1354 extrudeLayer::extrudeLayer
1355 (
1356     polyMeshGen& mesh,
1357     const LongList<labelPair>& extrusionFront,
1358     const scalar thickness
1359 )
1360 :
1361     mesh_(mesh),
1362     thickness_(thickness),
1363     nOrigPoints_(mesh.points().size()),
1364     nOrigFaces_(mesh.faces().size()),
1365     nOrigCells_(mesh.cells().size()),
1366     extrudedFaces_(),
1367     pairOrientation_(),
1368     origPointLabel_(nOrigPoints_, -1)
1369 {
1370     createDuplicateFrontFaces(extrusionFront);
1371
1372     createNewVertices();
1373
1374     movePoints();
1375
1376     createLayerCells();
1377
1378     updateBoundary();
1379
1380     mesh_.clearAddressingData();
1381
1382     # ifdef DEBUGExtrudeLayer
1383     polyMeshGenChecks::checkMesh(mesh_, true);
1384     # endif
1385 }
1386
1387 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
1388
1389 extrudeLayer::~extrudeLayer()
1390 {
1391     mesh_.clearAddressingData();
1392 }
1393
1394 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
1395
1396 } // End namespace Foam
1397
1398 // ************************************************************************* //