src/mesh/cfMesh/utilities/meshes/partTetMesh/partTetMesh.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | cfMesh: A library for mesh generation
   4    \\    /   O peration     |
   5     \\  /    A nd           | Author: Franjo Juretic (franjo.juretic@c-fields.com)
   6      \\/     M anipulation  | Copyright (C) Creative Fields, Ltd.
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of cfMesh.
  10
  11     cfMesh is free software; you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by the
  13     Free Software Foundation; either version 3 of the License, or (at your
  14     option) any later version.
  15
  16     cfMesh is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with cfMesh.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 Description
  25
  26 \*---------------------------------------------------------------------------*/
  27
  28 #include "demandDrivenData.H"
  29 #include "partTetMesh.H"
  30 #include "polyMeshGenModifier.H"
  31 #include "VRWGraphList.H"
  32 #include "polyMeshGenAddressing.H"
  33 #include "helperFunctions.H"
  34
  35 #include <map>
  36
  37 # ifdef USE_OMP
  38 #include <omp.h>
  39 # endif
  40
  41 //#define DEBUGSmooth
  42
  43 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  44
  45 namespace Foam
  46 {
  47
  48 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  49
  50 partTetMesh::partTetMesh(polyMeshGen& mesh, const labelLongList& lockedPoints)
  51 :
  52     origMesh_(mesh),
  53     points_(),
  54     tets_(),
  55     nodeLabelInOrigMesh_(),
  56     smoothVertex_(),
  57     pointTets_(),
  58     internalPointsOrderPtr_(NULL),
  59     boundaryPointsOrderPtr_(NULL),
  60     globalPointLabelPtr_(NULL),
  61     pAtProcsPtr_(NULL),
  62     globalToLocalPointAddressingPtr_(NULL),
  63     neiProcsPtr_(NULL),
  64     pAtParallelBoundariesPtr_(NULL),
  65     pAtBufferLayersPtr_(NULL)
  66 {
  67     List<direction> useCell(mesh.cells().size(), direction(1));
  68
  69     boolList lockedPoint(mesh.points().size(), false);
  70     forAll(lockedPoints, i)
  71         lockedPoint[lockedPoints[i]] = true;
  72
  73     createPointsAndTets(useCell, lockedPoint);
  74 }
  75
  76 partTetMesh::partTetMesh
  77 (
  78     polyMeshGen& mesh,
  79     const labelLongList& lockedPoints,
  80     const direction nLayers
  81 )
  82 :
  83     origMesh_(mesh),
  84     points_(),
  85     tets_(),
  86     nodeLabelInOrigMesh_(),
  87     smoothVertex_(),
  88     pointTets_(),
  89     internalPointsOrderPtr_(NULL),
  90     boundaryPointsOrderPtr_(NULL),
  91     globalPointLabelPtr_(NULL),
  92     pAtProcsPtr_(NULL),
  93     globalToLocalPointAddressingPtr_(NULL),
  94     neiProcsPtr_(NULL),
  95     pAtParallelBoundariesPtr_(NULL),
  96     pAtBufferLayersPtr_(NULL)
  97 {
  98     const faceListPMG& faces = mesh.faces();
  99     const cellListPMG& cells = mesh.cells();
 100     const labelList& owner = mesh.owner();
 101     const PtrList<boundaryPatch>& boundaries = mesh.boundaries();
 102     const VRWGraph& pointCells = mesh.addressingData().pointCells();
 103
 104     List<direction> useCell(cells.size(), direction(0));
 105
 106     //- select cells containing at least one vertex of the bad faces
 107     forAll(boundaries, patchI)
 108     {
 109         const label start = boundaries[patchI].patchStart();
 110         const label size = boundaries[patchI].patchSize();
 111
 112         for(label fI=0;fI<size;++fI)
 113         {
 114             useCell[owner[start+fI]] = 1;
 115         }
 116     }
 117
 118     //- add additional layer of cells
 119     for(direction layerI=1;layerI<(nLayers+1);++layerI)
 120     {
 121         forAll(useCell, cI)
 122             if( useCell[cI] == layerI )
 123             {
 124                 const cell& c = cells[cI];
 125
 126                 forAll(c, fI)
 127                 {
 128                     const face& f = faces[c[fI]];
 129
 130                     forAll(f, pI)
 131                     {
 132                         forAllRow(pointCells, f[pI], pcI)
 133                         {
 134                             const label cLabel = pointCells(f[pI], pcI);
 135                             if( !useCell[cLabel] )
 136                                 useCell[cLabel] = layerI + 1;
 137                         }
 138                     }
 139                 }
 140             }
 141
 142         if( Pstream::parRun() )
 143         {
 144             const labelLongList& globalPointLabel =
 145                 mesh.addressingData().globalPointLabel();
 146             const VRWGraph& pProcs = mesh.addressingData().pointAtProcs();
 147             const Map<label>& globalToLocal =
 148                 mesh.addressingData().globalToLocalPointAddressing();
 149
 150             std::map<label, LongList<label> > eData;
 151             forAllConstIter(Map<label>, globalToLocal, iter)
 152             {
 153                 const label pointI = iter();
 154
 155                 forAllRow(pProcs, pointI, procI)
 156                 {
 157                     const label neiProc = pProcs(pointI, procI);
 158                     if( neiProc == Pstream::myProcNo() )
 159                         continue;
 160
 161                     if( eData.find(neiProc) == eData.end() )
 162                     {
 163                         eData.insert
 164                         (
 165                             std::make_pair(neiProc, LongList<label>())
 166                         );
 167                     }
 168
 169                     forAllRow(pointCells, pointI, pcI)
 170                         if( useCell[pointCells(pointI, pcI)] == layerI )
 171                         {
 172                             eData[neiProc].append(globalPointLabel[pointI]);
 173                             break;
 174                         }
 175                 }
 176             }
 177
 178             //- exchange data with other processors
 179             labelLongList receivedData;
 180             help::exchangeMap(eData, receivedData);
 181
 182             forAll(receivedData, i)
 183             {
 184                 const label pointI = globalToLocal[receivedData[i]];
 185
 186                 forAllRow(pointCells, pointI, pcI)
 187                 {
 188                     const label cLabel = pointCells(pointI, pcI);
 189                     if( !useCell[cLabel] )
 190                         useCell[cLabel] = layerI + 1;
 191                 }
 192             }
 193         }
 194     }
 195
 196     boolList lockedPoint(mesh.points().size(), false);
 197     forAll(lockedPoints, i)
 198         lockedPoint[lockedPoints[i]] = true;
 199
 200     createPointsAndTets(useCell, lockedPoint);
 201 }
 202
 203 partTetMesh::partTetMesh
 204 (
 205     polyMeshGen& mesh,
 206     const labelLongList& lockedPoints,
 207     labelHashSet& badFaces,
 208     const direction additionalLayers
 209 )
 210 :
 211     origMesh_(mesh),
 212     points_(),
 213     tets_(),
 214     nodeLabelInOrigMesh_(),
 215     smoothVertex_(),
 216     pointTets_(),
 217     internalPointsOrderPtr_(NULL),
 218     boundaryPointsOrderPtr_(NULL),
 219     globalPointLabelPtr_(NULL),
 220     pAtProcsPtr_(NULL),
 221     globalToLocalPointAddressingPtr_(NULL),
 222     neiProcsPtr_(NULL),
 223     pAtParallelBoundariesPtr_(NULL),
 224     pAtBufferLayersPtr_(NULL)
 225 {
 226     const faceListPMG& faces = mesh.faces();
 227     const cellListPMG& cells = mesh.cells();
 228     const VRWGraph& pointCells = mesh.addressingData().pointCells();
 229
 230     List<direction> useCell(cells.size(), direction(0));
 231
 232     //- select cells containing at least one vertex of the bad faces
 233     forAll(faces, faceI)
 234         if( badFaces.found(faceI) )
 235         {
 236             const face& f = faces[faceI];
 237
 238             forAll(f, pI)
 239             {
 240                 forAllRow(pointCells, f[pI], pcI)
 241                     useCell[pointCells(f[pI], pcI)] = 1;
 242             }
 243         }
 244
 245     //- add additional layer of cells
 246     for(direction layerI=1;layerI<(additionalLayers+1);++layerI)
 247     {
 248         forAll(useCell, cI)
 249             if( useCell[cI] == layerI )
 250             {
 251                 const cell& c = cells[cI];
 252
 253                 forAll(c, fI)
 254                 {
 255                     const face& f = faces[c[fI]];
 256
 257                     forAll(f, pI)
 258                     {
 259                         forAllRow(pointCells, f[pI], pcI)
 260                         {
 261                             const label cLabel = pointCells(f[pI], pcI);
 262                             if( !useCell[cLabel] )
 263                                 useCell[cLabel] = layerI + 1;
 264                         }
 265                     }
 266                 }
 267             }
 268
 269         if( Pstream::parRun() )
 270         {
 271             const labelLongList& globalPointLabel =
 272                 mesh.addressingData().globalPointLabel();
 273             const VRWGraph& pProcs = mesh.addressingData().pointAtProcs();
 274             const Map<label>& globalToLocal =
 275                 mesh.addressingData().globalToLocalPointAddressing();
 276
 277             std::map<label, LongList<label> > eData;
 278             forAllConstIter(Map<label>, globalToLocal, iter)
 279             {
 280                 const label pointI = iter();
 281
 282                 forAllRow(pProcs, pointI, procI)
 283                 {
 284                     const label neiProc = pProcs(pointI, procI);
 285                     if( neiProc == Pstream::myProcNo() )
 286                         continue;
 287
 288                     if( eData.find(neiProc) == eData.end() )
 289                     {
 290                         eData.insert
 291                         (
 292                             std::make_pair(neiProc, LongList<label>())
 293                         );
 294                     }
 295
 296                     forAllRow(pointCells, pointI, pcI)
 297                         if( useCell[pointCells(pointI, pcI)] == layerI )
 298                         {
 299                             eData[neiProc].append(globalPointLabel[pointI]);
 300                             break;
 301                         }
 302                 }
 303             }
 304
 305             //- exchange data with other processors
 306             labelLongList receivedData;
 307             help::exchangeMap(eData, receivedData);
 308
 309             forAll(receivedData, i)
 310             {
 311                 const label pointI = globalToLocal[receivedData[i]];
 312
 313                 forAllRow(pointCells, pointI, pcI)
 314                 {
 315                     const label cLabel = pointCells(pointI, pcI);
 316                     if( !useCell[cLabel] )
 317                         useCell[cLabel] = layerI + 1;
 318                 }
 319             }
 320         }
 321     }
 322
 323     boolList lockedPoint(mesh.points().size(), false);
 324     forAll(lockedPoints, i)
 325         lockedPoint[lockedPoints[i]] = true;
 326
 327     createPointsAndTets(useCell, lockedPoint);
 328 }
 329
 330 partTetMesh::~partTetMesh()
 331 {
 332     deleteDemandDrivenData(internalPointsOrderPtr_);
 333     deleteDemandDrivenData(boundaryPointsOrderPtr_);
 334     deleteDemandDrivenData(globalPointLabelPtr_);
 335     deleteDemandDrivenData(pAtProcsPtr_);
 336     deleteDemandDrivenData(globalToLocalPointAddressingPtr_);
 337     deleteDemandDrivenData(neiProcsPtr_);
 338     deleteDemandDrivenData(pAtParallelBoundariesPtr_);
 339     deleteDemandDrivenData(pAtBufferLayersPtr_);
 340 }
 341
 342 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 343
 344 const VRWGraph& partTetMesh::internalPointOrdering() const
 345 {
 346     # ifdef USE_OMP
 347     if( omp_in_parallel() )
 348     {
 349         FatalErrorIn
 350         (
 351             "const VRWGraph& partTetMesh::internalPointOrdering() const"
 352         ) << "Calculating addressing inside a parallel region."
 353           << " This is not thread safe" << exit(FatalError);
 354     }
 355     # endif
 356
 357     if( !internalPointsOrderPtr_ )
 358         createSMOOTHPointsOrdering();
 359
 360     return *internalPointsOrderPtr_;
 361 }
 362
 363 const VRWGraph& partTetMesh::boundaryPointOrdering() const
 364 {
 365     # ifdef USE_OMP
 366     if( omp_in_parallel() )
 367     {
 368         FatalErrorIn
 369         (
 370             "const VRWGraph& partTetMesh::boundaryPointOrdering() const"
 371         ) << "Calculating addressing inside a parallel region."
 372           << " This is not thread safe" << exit(FatalError);
 373     }
 374     # endif
 375
 376     if( !boundaryPointsOrderPtr_ )
 377         createBOUNDARYPointsOrdering();
 378
 379     return *boundaryPointsOrderPtr_;
 380 }
 381 void partTetMesh::updateVertex(const label pointI, const point& newP)
 382 {
 383     points_[pointI] = newP;
 384
 385     if( smoothVertex_[pointI] & (FACECENTRE+CELLCENTRE) )
 386     {
 387         Warning << "Smoothing auxiliary vertex."
 388             << " This has no effect on the original mesh" << endl;
 389         return;
 390     }
 391
 392     //- find face centres attached
 393     DynList<label, 64> helper;
 394     forAllRow(pointTets_, pointI, ptI)
 395     {
 396         const label centreI = tets_[pointTets_(pointI, ptI)][2];
 397         if( smoothVertex_[centreI] & FACECENTRE )
 398             helper.appendIfNotIn(centreI);
 399     }
 400
 401     //- update coordinates of FACECENTRE vertices
 402     forAll(helper, i)
 403     {
 404         const label centreI = helper[i];
 405
 406         point centre(vector::zero);
 407         scalar faceArea(0.0);
 408         forAllRow(pointTets_, centreI, ptI)
 409         {
 410             const partTet& tet = tets_[pointTets_(centreI, ptI)];
 411             point c(vector::zero);
 412             for(label i=0;i<3;++i)
 413                 c += points_[tet[i]];
 414             c /= 3;
 415             const scalar area = Foam::mag(tet.Sd(points_)) + VSMALL;
 416
 417             centre += c * area;
 418             faceArea += area;
 419         }
 420
 421         points_[centreI] = centre / faceArea;
 422     }
 423
 424     //- find cell centres attached
 425     helper.clear();
 426     forAllRow(pointTets_, pointI, ptI)
 427     {
 428         const label centreI = tets_[pointTets_(pointI, ptI)][3];
 429         if( smoothVertex_[centreI] & CELLCENTRE )
 430             helper.appendIfNotIn(centreI);
 431     }
 432
 433     //- update coordinates of CELLCENTRE vertices
 434     forAll(helper, i)
 435     {
 436         const label centreI = helper[i];
 437
 438         point centre(vector::zero);
 439         scalar cellVol(0.0);
 440         forAllRow(pointTets_, centreI, ptI)
 441         {
 442             const partTet& tet = tets_[pointTets_(centreI, ptI)];
 443             const point c = tet.centroid(points_);
 444             const scalar vol = Foam::mag(tet.mag(points_)) + VSMALL;
 445
 446             centre += c * vol;
 447             cellVol += vol;
 448         }
 449
 450         points_[centreI] = centre / cellVol;
 451     }
 452 }
 453
 454 void partTetMesh::updateVerticesSMP(const List<LongList<labelledPoint> >& np)
 455 {
 456     List<direction> updateType(points_.size(), direction(0));
 457
 458     # ifdef USE_OMP
 459     # pragma omp parallel num_threads(np.size())
 460     # endif
 461     {
 462         # ifdef USE_OMP
 463         const LongList<labelledPoint>& newPoints = np[omp_get_thread_num()];
 464         # else
 465         const LongList<labelledPoint>& newPoints = np[0];
 466         # endif
 467
 468         forAll(newPoints, i)
 469         {
 470             const labelledPoint& lp = newPoints[i];
 471             const label pointI = lp.pointLabel();
 472
 473             points_[pointI] = lp.coordinates();
 474
 475             forAllRow(pointTets_, pointI, ptI)
 476             {
 477                 const partTet& pt = tets_[pointTets_(pointI, ptI)];
 478
 479                 if( smoothVertex_[pt[3]] & CELLCENTRE )
 480                     updateType[pt[3]] |= CELLCENTRE;
 481                 if( smoothVertex_[pt[2]] & FACECENTRE )
 482                     updateType[pt[2]] |= FACECENTRE;
 483             }
 484         }
 485
 486         # ifdef USE_OMP
 487         # pragma omp barrier
 488
 489         # pragma omp flush(updateType)
 490
 491         //- update coordinates of CELLCENTRE and FACECENTRE vertices
 492         # pragma omp for schedule(dynamic, 20)
 493         # endif
 494         forAll(updateType, pI)
 495         {
 496             if( updateType[pI] & CELLCENTRE )
 497             {
 498                 point centre(vector::zero);
 499                 scalar cellVol(0.0);
 500                 forAllRow(pointTets_, pI, ptI)
 501                 {
 502                     const partTet& tet = tets_[pointTets_(pI, ptI)];
 503                     const point c = tet.centroid(points_);
 504                     const scalar vol = Foam::mag(tet.mag(points_)) + VSMALL;
 505
 506                     centre += c * vol;
 507                     cellVol += vol;
 508                 }
 509
 510                 points_[pI] = centre / cellVol;
 511             }
 512             else if( updateType[pI] & FACECENTRE )
 513             {
 514                 point centre(vector::zero);
 515                 scalar faceArea(0.0);
 516                 forAllRow(pointTets_, pI, ptI)
 517                 {
 518                     const partTet& tet = tets_[pointTets_(pI, ptI)];
 519                     point c(vector::zero);
 520                     for(label i=0;i<3;++i)
 521                         c += points_[tet[i]];
 522                     c /= 3;
 523                     const scalar area = Foam::mag(tet.Sd(points_)) + VSMALL;
 524
 525                     centre += c * area;
 526                     faceArea += area;
 527                 }
 528
 529                 points_[pI] = centre / faceArea;
 530             }
 531         }
 532     }
 533 }
 534
 535 void partTetMesh::updateOrigMesh(boolList* changedFacePtr)
 536 {
 537     pointFieldPMG& pts = origMesh_.points();
 538
 539     boolList changedNode(pts.size(), false);
 540
 541     # ifdef USE_OMP
 542     # pragma omp parallel for if( pts.size() > 1000 ) \
 543     schedule(guided, 10)
 544     # endif
 545     forAll(nodeLabelInOrigMesh_, pI)
 546         if( nodeLabelInOrigMesh_[pI] != -1 )
 547         {
 548             changedNode[nodeLabelInOrigMesh_[pI]] = true;
 549             pts[nodeLabelInOrigMesh_[pI]] = points_[pI];
 550         }
 551
 552     if( changedFacePtr )
 553     {
 554         boolList& chF = *changedFacePtr;
 555         chF = false;
 556
 557         const cellListPMG& cells = origMesh_.cells();
 558         const VRWGraph& pointCells = origMesh_.addressingData().pointCells();
 559
 560         # ifdef USE_OMP
 561         # pragma omp parallel for if( pointCells.size() > 100 ) \
 562         schedule(dynamic, 20)
 563         # endif
 564         forAll(pointCells, pointI)
 565         {
 566             if( changedNode[pointI] )
 567             {
 568                 forAllRow(pointCells, pointI, pcI)
 569                 {
 570                     const cell& c = cells[pointCells(pointI, pcI)];
 571
 572                     forAll(c, fI)
 573                         chF[c[fI]] = true;
 574                 }
 575             }
 576         }
 577
 578         //- make sure that neighbouring processors get the same information
 579         const PtrList<processorBoundaryPatch>& pBnd = origMesh_.procBoundaries();
 580         forAll(pBnd, patchI)
 581         {
 582             const label start = pBnd[patchI].patchStart();
 583             const label size = pBnd[patchI].patchSize();
 584
 585             labelLongList sendData;
 586             for(label faceI=0;faceI<size;++faceI)
 587             {
 588                 if( chF[start+faceI] )
 589                     sendData.append(faceI);
 590             }
 591
 592             OPstream toOtherProc
 593             (
 594                 Pstream::blocking,
 595                 pBnd[patchI].neiProcNo(),
 596                 sendData.byteSize()
 597             );
 598
 599             toOtherProc << sendData;
 600         }
 601
 602         forAll(pBnd, patchI)
 603         {
 604             labelList receivedData;
 605
 606             IPstream fromOtherProc
 607             (
 608                 Pstream::blocking,
 609                 pBnd[patchI].neiProcNo()
 610             );
 611
 612             fromOtherProc >> receivedData;
 613
 614             const label start = pBnd[patchI].patchStart();
 615             forAll(receivedData, i)
 616                 chF[start+receivedData[i]] = true;
 617         }
 618
 619         //- update geometry information
 620         const_cast<polyMeshGenAddressing&>
 621         (
 622             origMesh_.addressingData()
 623         ).updateGeometry(chF);
 624     }
 625     else
 626     {
 627         const_cast<polyMeshGenAddressing&>
 628         (
 629             origMesh_.addressingData()
 630         ).clearGeom();
 631     }
 632 }
 633
 634 void partTetMesh::createPolyMesh(polyMeshGen& pmg) const
 635 {
 636     polyMeshGenModifier meshModifier(pmg);
 637
 638     pointFieldPMG& pAccess = meshModifier.pointsAccess();
 639     pAccess.setSize(points_.size());
 640     forAll(points_, pI)
 641         pAccess[pI] = points_[pI];
 642
 643     VRWGraphList cellFaces;
 644
 645     forAll(tets_, tetI)
 646     {
 647         const partTet& tet = tets_[tetI];
 648
 649         FixedList<FixedList<label, 3>, 4> tetFaces;
 650
 651         //- face 0
 652         tetFaces[0][0] = tet[0];
 653         tetFaces[0][1] = tet[2];
 654         tetFaces[0][2] = tet[1];
 655
 656         //- face 1
 657         tetFaces[1][0] = tet[0];
 658         tetFaces[1][1] = tet[1];
 659         tetFaces[1][2] = tet[3];
 660
 661         //- face 2
 662         tetFaces[2][0] = tet[0];
 663         tetFaces[2][1] = tet[3];
 664         tetFaces[2][2] = tet[2];
 665
 666         //- face 3
 667         tetFaces[3][0] = tet[1];
 668         tetFaces[3][1] = tet[2];
 669         tetFaces[3][2] = tet[3];
 670
 671         cellFaces.appendGraph(tetFaces);
 672     }
 673
 674     meshModifier.addCells(cellFaces);
 675     meshModifier.reorderBoundaryFaces();
 676
 677     //- store points into subsets
 678     const label bndPointID = pmg.addPointSubset("boundaryPoints");
 679     const label smoothPointID = pmg.addPointSubset("smoothPoints");
 680     const label faceCentreID = pmg.addPointSubset("faceCentres");
 681     const label cellCentreID = pmg.addPointSubset("cellCentres");
 682
 683     forAll(smoothVertex_, pointI)
 684     {
 685         if( smoothVertex_[pointI] & SMOOTH )
 686             pmg.addPointToSubset(smoothPointID, pointI);
 687         if( smoothVertex_[pointI] & BOUNDARY )
 688             pmg.addPointToSubset(bndPointID, pointI);
 689         if( smoothVertex_[pointI] & FACECENTRE )
 690             pmg.addPointToSubset(faceCentreID, pointI);
 691         if( smoothVertex_[pointI] & CELLCENTRE )
 692             pmg.addPointToSubset(cellCentreID, pointI);
 693     }
 694
 695     const VRWGraph& internalOrdering = internalPointOrdering();
 696     const VRWGraph& boundaryOrdering = boundaryPointOrdering();
 697
 698     forAll(internalOrdering, i)
 699     {
 700         const word name = "smoothPoints_"+help::scalarToText(i);
 701         const label orderID = pmg.addPointSubset(name);
 702
 703         forAllRow(internalOrdering, i, nI)
 704             pmg.addPointToSubset(orderID, internalOrdering(i, nI));
 705     }
 706
 707     forAll(boundaryOrdering, i)
 708     {
 709         const word name = "boundaryPoints_"+help::scalarToText(i);
 710         const label orderID = pmg.addPointSubset(name);
 711
 712         forAllRow(boundaryOrdering, i, nI)
 713             pmg.addPointToSubset(orderID, boundaryOrdering(i, nI));
 714     }
 715 }
 716
 717 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 718
 719 } // End namespace Foam
 720
 721 // ************************************************************************* //