src/mesh/cfMesh/utilities/boundaryLayers/refineBoundaryLayers/refineBoundaryLayersCells.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 "refineBoundaryLayers.H"
  29 #include "meshSurfaceEngine.H"
  30 #include "helperFunctions.H"
  31 #include "demandDrivenData.H"
  32
  33 //#define DEBUGLayer
  34
  35 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  36
  37 namespace Foam
  38 {
  39
  40 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  41
  42 void refineBoundaryLayers::generateNewCellsPrism
  43 (
  44     const label cellI,
  45     DynList<DynList<DynList<label, 8>, 10>, 64>& cellsFromCell
  46 ) const
  47 {
  48     cellsFromCell.clear();
  49
  50     const cell& c = mesh_.cells()[cellI];
  51     const labelList& owner = mesh_.owner();
  52
  53     # ifdef DEBUGLayer
  54     Pout << "New cells from cell " << cellI << endl;
  55     # endif
  56
  57     const label startBoundary = mesh_.boundaries()[0].patchStart();
  58
  59     //- find the number of lyers for this cell
  60     label nLayers(1), baseFace(-1);
  61     forAll(c, fI)
  62     {
  63         const label bfI = c[fI] - startBoundary;
  64
  65         if( (bfI < 0) || (bfI >= nLayersAtBndFace_.size()) )
  66             continue;
  67
  68         if( nLayersAtBndFace_[bfI] < 2 )
  69             continue;
  70
  71         # ifdef DEBUGLayer
  72         Pout << "Boundary face " << bfI << endl;
  73         # endif
  74
  75         nLayers = nLayersAtBndFace_[bfI];
  76         baseFace = fI;
  77     }
  78
  79     # ifdef DEBUGLayer
  80     Pout << "Number of layers " << nLayers << endl;
  81     Pout << "Base face " << baseFace << " has points "
  82          << mesh_.faces()[c[baseFace]] << endl;
  83     forAll(c, fI)
  84     {
  85         Pout << "Faces from face " << fI << " are "
  86              << facesFromFace_[c[fI]] << endl;
  87
  88         forAllRow(facesFromFace_, c[fI], i)
  89             Pout << "Face " << facesFromFace_(c[fI], i)
  90                  << " is " << newFaces_[facesFromFace_(c[fI], i)] << endl;
  91     }
  92     # endif
  93
  94     //- set the number of layers
  95     cellsFromCell.setSize(nLayers);
  96
  97     //- distribute existing faces into new cells
  98     label otherBaseFace(-1);
  99     forAll(c, fI)
 100     {
 101         if( fI == baseFace )
 102         {
 103             const label faceI = facesFromFace_(c[fI], 0);
 104             DynList<label, 8> f;
 105             f = newFaces_[faceI];
 106             cellsFromCell[nLayers-1].append(f);
 107         }
 108         else if( facesFromFace_.sizeOfRow(c[fI]) == 1 )
 109         {
 110             const label faceI = facesFromFace_(c[fI], 0);
 111             otherBaseFace = fI;
 112             DynList<label, 8> f;
 113             f = newFaces_[faceI];
 114             cellsFromCell[0].append(f);
 115         }
 116         else
 117         {
 118             forAllRow(facesFromFace_, c[fI], cfI)
 119             {
 120                 const label nfI = facesFromFace_(c[fI], cfI);
 121
 122                 DynList<label, 8> cf;
 123                 cf = newFaces_[nfI];
 124
 125                 if( owner[c[fI]] != cellI )
 126                     cf = help::reverseFace(cf);
 127
 128                 cellsFromCell[Foam::max(nLayers-1-cfI, 0)].append(cf);
 129             }
 130         }
 131     }
 132
 133     //- generate missing faces
 134     const faceListPMG& faces = mesh_.faces();
 135     const face& bf = faces[c[baseFace]];
 136     const face& obf = faces[c[otherBaseFace]];
 137     for(label layerI=1;layerI<nLayers;++layerI)
 138     {
 139         //- create new face from points at the same height
 140         DynList<label, 8> cf;
 141         forAll(bf, pI)
 142         {
 143             const label pointI = bf[pI];
 144
 145             # ifdef DEBUGLayer
 146             Pout << "Split edges at point " << pointI << " are "
 147                  << splitEdgesAtPoint_[pointI] << endl;
 148             # endif
 149
 150             label seI(-1);
 151             if( splitEdgesAtPoint_.sizeOfRow(pointI) == 1 )
 152             {
 153                 seI = splitEdgesAtPoint_(pointI, 0);
 154             }
 155             else
 156             {
 157                 forAllRow(splitEdgesAtPoint_, pointI, sepI)
 158                 {
 159                     const label seJ = splitEdgesAtPoint_(pointI, sepI);
 160                     const edge& se = splitEdges_[seJ];
 161
 162                     if( obf.which(se.end()) >= 0 || obf.which(se.start()) >= 0 )
 163                     {
 164                         seI = seJ;
 165                         break;
 166                     }
 167                 }
 168             }
 169
 170             cf.append(newVerticesForSplitEdge_(seI, layerI));
 171         }
 172
 173         //- add faces to cells
 174         cellsFromCell[nLayers-layerI].append(cf);
 175         cellsFromCell[nLayers-1-layerI].append(cf);
 176     }
 177
 178     # ifdef DEBUGLayer
 179     Pout << "New cells from cell " << cellI << " are " << cellsFromCell << endl;
 180     //::exit(1);
 181
 182     Pout << "1. Newly generated cells " << cellsFromCell << endl;
 183
 184     //- check if all generated cells are topologically closed
 185     forAll(cellsFromCell, cI)
 186     {
 187         const DynList<DynList<label, 8>, 10>& cellFaces = cellsFromCell[cI];
 188
 189         DynList<edge, 12> edges;
 190         DynList<label, 12> nAppearances;
 191
 192         forAll(cellFaces, fI)
 193         {
 194             const DynList<label, 8>& f = cellFaces[fI];
 195
 196             forAll(f, eI)
 197             {
 198                 const edge e(f[eI], f.fcElement(eI));
 199
 200                 const label pos = edges.containsAtPosition(e);
 201
 202                 if( pos < 0 )
 203                 {
 204                     edges.append(e);
 205                     nAppearances.append(1);
 206                 }
 207                 else
 208                 {
 209                     ++nAppearances[pos];
 210                 }
 211             }
 212         }
 213
 214         forAll(nAppearances, eI)
 215             if( nAppearances[eI] != 2 )
 216             {
 217                 Pout << "Prism cell " << cI << " edge " << edges[eI]
 218                     << " is present " << nAppearances[eI] << " times!" << endl;
 219                 abort(FatalError);
 220             }
 221     }
 222     # endif
 223 }
 224
 225 void refineBoundaryLayers::storeFacesIntoCells
 226 (
 227     const label faceI,
 228     const bool reverseOrientation,
 229     const label normalDirection,
 230     const bool maxCoordinate,
 231     const label nLayersI,
 232     const label nLayersJ,
 233     const label nLayersK,
 234     DynList<DynList<DynList<label, 4>, 6>, 256>& cellsFromCell
 235 ) const
 236 {
 237     DynList<DynList<label> > faceFaces;
 238     sortFaceFaces(faceI, faceFaces, reverseOrientation);
 239
 240     const label maxI = nLayersI - 1;
 241     const label maxJ = nLayersJ - 1;
 242     const label maxK = nLayersK - 1;
 243
 244     # ifdef DEBUGLayer
 245     Pout << "Storing new faces from face " << faceI
 246          << " reverseOrientation = " << reverseOrientation
 247          << " normal direction " << normalDirection
 248          << " maxCoordinate " << maxCoordinate << endl;
 249     Pout << "faceFaces " << faceFaces << endl;
 250     # endif
 251
 252     label i(-1), j(-1), k(-1);
 253
 254     forAll(faceFaces, nI)
 255     {
 256         forAll(faceFaces[nI], nJ)
 257         {
 258             const label nfI = faceFaces[nI][nJ];
 259
 260             # ifdef DEBUGLayer
 261             Pout << "nI = " << nI << " nJ = " << nJ << endl;
 262             # endif
 263
 264             if( normalDirection == 0 )
 265             {
 266                 //- k is const
 267                 i = Foam::min(nI, maxI);
 268                 j = Foam::min(nJ, maxJ);
 269                 k = maxCoordinate?maxK:0;
 270             }
 271             else if( normalDirection == 1 )
 272             {
 273                 //- j is const
 274                 i = Foam::min(nJ, maxI);
 275                 j = maxCoordinate?maxJ:0;
 276                 k = Foam::min(nI, maxK);
 277             }
 278             else if( normalDirection == 2 )
 279             {
 280                 //- i is const
 281                 i = maxCoordinate?maxI:0;
 282                 j = Foam::min(nI, maxJ);
 283                 k = Foam::min(nJ, maxK);
 284             }
 285
 286             //- store the face into a new cell
 287             const label cI
 288             (
 289                 j * nLayersI +
 290                 k * nLayersI * nLayersJ +
 291                 i
 292             );
 293
 294             # ifdef DEBUGLayer
 295             Pout << "Storing face " << newFaces_[nfI]
 296                  << " i = " << i << " j = " << j << " k = " << k
 297                  << "\n cell label " << cI << endl;
 298             # endif
 299
 300             cellsFromCell[cI].append(newFaces_[nfI]);
 301         }
 302     }
 303 }
 304
 305 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 306
 307 void refineBoundaryLayers::refineEdgeHexCell::determineFacesInDirections()
 308 {
 309     const labelList& nLayersAtBndFace = bndLayers_.nLayersAtBndFace_;
 310     const polyMeshGen& mesh = bndLayers_.mesh_;
 311     const faceListPMG& faces = mesh.faces();
 312     const cell& c = mesh.cells()[cellI_];
 313
 314     # ifdef DEBUGLayer
 315     Pout << "Generating new cells from edge cell " << cellI_ << endl;
 316     # endif
 317
 318     const PtrList<boundaryPatch>& bnd = mesh.boundaries();
 319     const label startBoundary = bnd[0].patchStart();
 320
 321     //- find the number of layers for this cell
 322     FixedList<label, 2> layersInDirection(-1), dirFace;
 323     label currDir(0);
 324
 325     FixedList<bool, 6> determinedFace(false);
 326
 327     forAll(c, fI)
 328     {
 329         const label bfI = c[fI] - startBoundary;
 330
 331         if( (bfI < 0) || (bfI >= nLayersAtBndFace.size()) )
 332             continue;
 333
 334         # ifdef DEBUGLayer
 335         Pout << "Boundary face " << bfI << endl;
 336         # endif
 337
 338         if( nLayersAtBndFace[bfI] < 2 )
 339             continue;
 340
 341         layersInDirection[currDir] = nLayersAtBndFace[bfI];
 342         dirFace[currDir] = fI;
 343         ++currDir;
 344     }
 345
 346     //- set the number of newly create cells
 347     nLayersI_ = layersInDirection[0];
 348     nLayersJ_ = layersInDirection[1];
 349     cellsFromCell_.setSize(nLayersI_ * nLayersJ_);
 350
 351     //- find the shared edge between the boundary faces
 352     const edge commonEdge =
 353         help::sharedEdge(faces[c[dirFace[0]]], faces[c[dirFace[1]]]);
 354
 355     //- faces at i = const in the local coordinate system
 356     faceInDirection_[4] = dirFace[0];
 357     determinedFace[dirFace[0]] = true;
 358     forAll(c, fI)
 359     {
 360         if( determinedFace[fI] )
 361             continue;
 362
 363         if( !help::shareAnEdge(faces[c[dirFace[0]]], faces[c[fI]]) )
 364         {
 365             faceInDirection_[5] = fI;
 366             determinedFace[fI] = true;
 367             break;
 368         }
 369     }
 370
 371     //- faces k = const in the local coordinate system
 372     faceInDirection_[2] = dirFace[1];
 373     determinedFace[dirFace[1]] = true;
 374     forAll(c, fI)
 375     {
 376         if( determinedFace[fI] )
 377             continue;
 378
 379         if( !help::shareAnEdge(faces[c[dirFace[1]]], faces[c[fI]]) )
 380         {
 381             faceInDirection_[3] = fI;
 382             determinedFace[fI] = true;
 383             break;
 384         }
 385     }
 386
 387     # ifdef DEBUGLayer
 388     Pout << "Common edge " << commonEdge << endl;
 389     Pout << "Donor face " << dirFace[0] << endl;
 390     Pout << "Donor face points " << faces[c[dirFace[0]]] << endl;
 391     # endif
 392
 393     //- find the face attached to the starting point of the edge and
 394     //- the face attached to the end point of the edge
 395     forAll(c, fI)
 396     {
 397         if( determinedFace[fI] )
 398             continue;
 399
 400         if(
 401             (faces[c[fI]].which(commonEdge.start()) >= 0) &&
 402             (help::positionOfEdgeInFace(commonEdge, faces[c[fI]]) < 0)
 403         )
 404             faceInDirection_[0] = fI;
 405
 406         if(
 407             (faces[c[fI]].which(commonEdge.end()) >= 0) &&
 408             (help::positionOfEdgeInFace(commonEdge, faces[c[fI]]) < 0)
 409         )
 410             faceInDirection_[1] = fI;
 411     }
 412
 413     //- check the orientation of faces
 414     const labelList& owner = mesh.owner();
 415
 416     //- checking face at direction k = 0
 417     faceOrientation_[0] = owner[c[faceInDirection_[0]]] == cellI_?true:false;
 418
 419     //- checking face in direction k = 1
 420     faceOrientation_[1] = owner[c[faceInDirection_[1]]] == cellI_?false:true;
 421
 422     //- set orientation flag for face in direction j = 0
 423     faceOrientation_[2] = true;
 424
 425     //- checking face in direction j = nLayersJ_
 426     faceOrientation_[3] = owner[c[faceInDirection_[3]]] == cellI_?false:true;
 427
 428     //- set orientation flag for face in direction i = 0
 429     faceOrientation_[4] = true;
 430
 431     //- checking face in direction i = nLayersI_
 432     faceOrientation_[5] = owner[c[faceInDirection_[5]]] == cellI_?false:true;
 433
 434     # ifdef DEBUGLayer
 435     Pout << "Face at start " << faces[c[faceInDirection_[0]]] << endl;
 436     Pout << "Face at end " << faces[c[faceInDirection_[1]]] << endl;
 437     forAll(faceInDirection_, i)
 438         Pout << "Face in direction " << i << " is "
 439              << faces[c[faceInDirection_[i]]]
 440              << " orientation " << faceOrientation_[i] << endl;
 441     # endif
 442 }
 443
 444 void refineBoundaryLayers::refineEdgeHexCell::populateExistingFaces()
 445 {
 446     const cell& c = bndLayers_.mesh_.cells()[cellI_];
 447     const VRWGraph& facesFromFace = bndLayers_.facesFromFace_;
 448     const VRWGraph& newFaces = bndLayers_.newFaces_;
 449
 450     cellsFromCell_.setSize(nLayersI_ * nLayersJ_);
 451     forAll(cellsFromCell_, cI)
 452         cellsFromCell_[cI].clear();
 453
 454     //- store new faces at k = 0
 455     bndLayers_.storeFacesIntoCells
 456     (
 457         c[faceInDirection_[0]], faceOrientation_[0],
 458         0, 0,
 459         nLayersI_, nLayersJ_, 1,
 460         cellsFromCell_
 461     );
 462
 463     //- store new faces at k = 1
 464     bndLayers_.storeFacesIntoCells
 465     (
 466         c[faceInDirection_[1]], faceOrientation_[1],
 467         0, 1,
 468         nLayersI_, nLayersJ_, 1,
 469         cellsFromCell_
 470     );
 471
 472     //- store new faces at j = 0
 473     forAllRow(facesFromFace, c[faceInDirection_[2]], i)
 474     {
 475         const label faceI = facesFromFace(c[faceInDirection_[2]], i);
 476         cellsFromCell_[i].append(newFaces[faceI]);
 477     }
 478
 479     //- store faces at j = nLayersJ
 480     const label maxJ = nLayersJ_ - 1;
 481     forAllRow(facesFromFace, c[faceInDirection_[3]], i)
 482     {
 483         const label faceI = facesFromFace(c[faceInDirection_[3]], i);
 484         cellsFromCell_[i + maxJ * nLayersI_].append(newFaces[faceI]);
 485     }
 486
 487     //- store new faces at i = 0
 488     forAllRow(facesFromFace, c[faceInDirection_[4]], j)
 489     {
 490         const label faceI = facesFromFace(c[faceInDirection_[4]], j);
 491         cellsFromCell_[j * nLayersI_].append(newFaces[faceI]);
 492     }
 493
 494     //- store new faces at i = nLayersI
 495     const label maxI = nLayersI_ - 1;
 496     forAllRow(facesFromFace, c[faceInDirection_[5]], j)
 497     {
 498         const label faceI = facesFromFace(c[faceInDirection_[5]], j);
 499         cellsFromCell_[j * nLayersI_ + maxI].append(newFaces[faceI]);
 500     }
 501
 502     # ifdef DEBUGLayer
 503     Pout << "New cells after populating existing faces "
 504          << cellsFromCell_ << endl;
 505     # endif
 506 }
 507
 508 void refineBoundaryLayers::refineEdgeHexCell::generateMissingFaces()
 509 {
 510     const cell& c = bndLayers_.mesh_.cells()[cellI_];
 511
 512     //- fill up the matrix of points for this cell
 513     //- the matrix is used for generation of new cells
 514     FixedList<DynList<DynList<label> >, 2> cellPoints;
 515
 516     //- fill in the data for a cross-split faces
 517     bndLayers_.sortFacePoints
 518     (
 519         c[faceInDirection_[0]],
 520         cellPoints[0],
 521         faceOrientation_[0]
 522     );
 523     bndLayers_.sortFacePoints
 524     (
 525         c[faceInDirection_[1]],
 526         cellPoints[1],
 527         faceOrientation_[1]
 528     );
 529
 530     //- generate new internal faces for this cell
 531     //- generate faces with normal in the i direction
 532     const label maxI = nLayersI_ - 1;
 533     const label maxJ = nLayersJ_ - 1;
 534
 535     for(label i=1;i<nLayersI_;++i)
 536     {
 537         for(label j=0;j<nLayersJ_;++j)
 538         {
 539             const label own = j * nLayersI_ + i - 1;
 540             const label nei = own + 1;
 541
 542             if( j < maxJ )
 543             {
 544                 //- generate a quad face
 545                 FixedList<label, 4> mf;
 546
 547                 //- populate the points form cellPoints
 548                 mf[0] = cellPoints[0][i][j];
 549                 mf[1] = cellPoints[0][i][j+1];
 550                 mf[2] = cellPoints[1][i][j+1];
 551                 mf[3] = cellPoints[1][i][j];
 552
 553                 # ifdef DEBUGLayer
 554                 Pout << "1. Adding missing face " << mf
 555                      << " to cells " << own << " and " << nei << endl;
 556                 # endif
 557
 558                 cellsFromCell_[own].append(mf);
 559                 cellsFromCell_[nei].append(help::reverseFace(mf));
 560             }
 561             else
 562             {
 563                 DynList<label> mf;
 564                 for(label index=j;index<cellPoints[0][i].size();++index)
 565                     mf.append(cellPoints[0][i][index]);
 566                 for(label index=cellPoints[1][i].size()-1;index>=j;--index)
 567                     mf.append(cellPoints[1][i][index]);
 568
 569                 # ifdef DEBUGLayer
 570                 Pout << "2. Adding missing face " << mf
 571                      << " to cells " << own << " and " << nei << endl;
 572                 # endif
 573
 574                 cellsFromCell_[own].append(mf);
 575                 cellsFromCell_[nei].append(help::reverseFace(mf));
 576             };
 577         }
 578     }
 579
 580     //- generate faces with the normal in j direction
 581     for(label i=0;i<nLayersI_;++i)
 582     {
 583         for(label j=1;j<nLayersJ_;++j)
 584         {
 585             const label nei = j * nLayersI_ + i;
 586             const label own = (j - 1) * nLayersI_ + i;
 587
 588             if( i < maxI )
 589             {
 590                 //- generate a quad face
 591                 FixedList<label, 4> mf;
 592
 593                 //- populate the points form cellPoints
 594                 mf[0] = cellPoints[0][i][j];
 595                 mf[1] = cellPoints[1][i][j];
 596                 mf[2] = cellPoints[1][i+1][j];
 597                 mf[3] = cellPoints[0][i+1][j];
 598
 599                 # ifdef DEBUGLayer
 600                 Pout << "3. Adding missing face " << mf
 601                      << " to cells " << own << " and " << nei << endl;
 602                 # endif
 603
 604                 cellsFromCell_[own].append(mf);
 605                 cellsFromCell_[nei].append(help::reverseFace(mf));
 606             }
 607             else
 608             {
 609                 DynList<label> mf;
 610                 for(label index=i;index<cellPoints[1].size();++index)
 611                     mf.append(cellPoints[1][index][j]);
 612                 for(label index=cellPoints[0].size()-1;index>=i;--index)
 613                     mf.append(cellPoints[0][index][j]);
 614
 615                 # ifdef DEBUGLayer
 616                 Pout << "4. Adding missing face " << mf
 617                      << " to cells " << own << " and " << nei << endl;
 618                 # endif
 619
 620                 cellsFromCell_[own].append(mf);
 621                 cellsFromCell_[nei].append(help::reverseFace(mf));
 622             };
 623         }
 624     }
 625
 626     # ifdef DEBUGLayer
 627     Pout << "Cell " << cellI_ << " new cells are " << cellsFromCell_ << endl;
 628     //::exit(1);
 629     # endif
 630 }
 631
 632 refineBoundaryLayers::refineEdgeHexCell::refineEdgeHexCell
 633 (
 634     const label cellI,
 635     const refineBoundaryLayers& ref
 636 )
 637 :
 638     cellI_(cellI),
 639     nLayersI_(),
 640     nLayersJ_(),
 641     cellsFromCell_(),
 642     bndLayers_(ref),
 643     faceInDirection_(),
 644     faceOrientation_(),
 645     cellPoints_()
 646 {
 647     determineFacesInDirections();
 648
 649     populateExistingFaces();
 650
 651     generateMissingFaces();
 652 }
 653
 654 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 655
 656 void refineBoundaryLayers::refineCornerHexCell::determineFacesInDirections()
 657 {
 658     const polyMeshGen& mesh = bndLayers_.mesh_;
 659     const cell& c = mesh.cells()[cellI_];
 660     const faceListPMG& faces = mesh.faces();
 661     const labelList& nLayersAtBndFace = bndLayers_.nLayersAtBndFace_;
 662
 663     # ifdef DEBUGLayer
 664     Pout << "Generating new cells from corner hex cell " << cellI_ << endl;
 665     Pout << "Cell faces " << c << endl;
 666     # endif
 667
 668     const label startBoundary = mesh.boundaries()[0].patchStart();
 669
 670     //- find the number of layers for this cell
 671     FixedList<label, 3> layersInDirection(-1), dirFace;
 672     FixedList<bool, 6> usedDirection(false);
 673     label currDir(0);
 674
 675     forAll(c, fI)
 676     {
 677         const label bfI = c[fI] - startBoundary;
 678
 679         if( (bfI < 0) || (bfI >= nLayersAtBndFace.size()) )
 680             continue;
 681
 682         # ifdef DEBUGLayer
 683         Pout << "Boundary face " << bfI << endl;
 684         # endif
 685
 686         if( nLayersAtBndFace[bfI] < 2 )
 687             continue;
 688
 689         usedDirection[fI] = true;
 690         layersInDirection[currDir] = nLayersAtBndFace[bfI];
 691         dirFace[currDir] = fI;
 692         ++currDir;
 693     }
 694
 695     //- find a common point for all three boundary faces
 696     FixedList<DynList<label, 4>, 3> bndFaces;
 697     forAll(dirFace, i)
 698     {
 699         bndFaces[i] = faces[c[dirFace[i]]];
 700     }
 701
 702     const label commonPoint = help::sharedVertex(bndFaces);
 703
 704     # ifdef DEBUGLayer
 705     Pout << "Used directions " << usedDirection << endl;
 706     Pout << "Layers in direction " << layersInDirection << endl;
 707     Pout << "dirFace " << dirFace << endl;
 708     Pout << "Common point " << commonPoint << endl;
 709
 710     forAll(dirFace, i)
 711         Pout << "bnd face " << i << " is " << faces[c[dirFace[i]]] << endl;
 712     # endif
 713
 714     //- find the position of the common point in each boundary face
 715     const edgeLongList& splitEdges = bndLayers_.splitEdges_;
 716     const VRWGraph& splitEdgesAtPoint = bndLayers_.splitEdgesAtPoint_;
 717
 718     const face& baseFace = faces[c[dirFace[0]]];
 719     const label posInBndFace = baseFace.which(commonPoint);
 720
 721     //- find split edges starting at the commonPoints
 722     forAllRow(splitEdgesAtPoint, commonPoint, i)
 723     {
 724         const edge& se = splitEdges[splitEdgesAtPoint(commonPoint, i)];
 725
 726         if( se == baseFace.faceEdge(posInBndFace) )
 727         {
 728             //- this edge is in j direction
 729             splitEdgeInDirection_[1] = splitEdgesAtPoint(commonPoint, i);
 730         }
 731         else if( se == baseFace.faceEdge(baseFace.rcIndex(posInBndFace)) )
 732         {
 733             //- this edge is in i diretion
 734             splitEdgeInDirection_[0] = splitEdgesAtPoint(commonPoint, i);
 735         }
 736         else if( splitEdgesAtPoint.sizeOfRow(commonPoint) == 3 )
 737         {
 738             //- this point is in k direction
 739             splitEdgeInDirection_[2] = splitEdgesAtPoint(commonPoint, i);
 740         }
 741         else
 742         {
 743             //- this situation is not allowed
 744             FatalErrorIn
 745             (
 746                 "void refineBoundaryLayers::refineCornerHexCell::"
 747                 "determineFacesInDirections()"
 748             ) << "Cannot refine layer for cell " << cellI_ << abort(FatalError);
 749         }
 750     }
 751
 752     # ifdef DEBUGLayer
 753     const VRWGraph& newVerticesForSplitEdge =
 754         bndLayers_.newVerticesForSplitEdge_;
 755     forAll(splitEdgeInDirection_, i)
 756         Pout << "Split edge in direction " << i << " has nodes "
 757              << splitEdges[splitEdgeInDirection_[i]]
 758              << " number of points on split edge "
 759              << newVerticesForSplitEdge.sizeOfRow(splitEdgeInDirection_[i])
 760              << endl;
 761     # endif
 762
 763     //- find the direction od other boundary faces
 764     //- in the local coordinate system
 765     FixedList<label, 3> permutation;
 766     permutation[0] = 0;
 767
 768     label helper = help::positionOfEdgeInFace
 769     (
 770         baseFace.faceEdge(baseFace.rcIndex(posInBndFace)),
 771         faces[c[dirFace[1]]]
 772     );
 773
 774     if( helper >= 0 )
 775     {
 776         permutation[1] = 1;
 777         permutation[2] = 2;
 778     }
 779     else
 780     {
 781         permutation[1] = 2;
 782         permutation[2] = 1;
 783     }
 784
 785     //- find the number of layers and a split in each direction
 786     nLayersI_ = layersInDirection[permutation[2]];
 787     nLayersJ_ = layersInDirection[permutation[1]];
 788     nLayersK_ = layersInDirection[permutation[0]];
 789
 790     //- determine the directions of cell faces
 791     //- store boundary faces first. Their normals point in the wrong direction
 792     //- face at k = 0
 793     faceInDirection_[0] = dirFace[permutation[0]];
 794     faceOrientation_[0] = true;
 795     //- face at j = 0
 796     faceInDirection_[2] = dirFace[permutation[1]];
 797     faceOrientation_[2] = true;
 798     //- face at i = 0
 799     faceInDirection_[4] = dirFace[permutation[2]];
 800     faceOrientation_[4] = true;
 801
 802     //- find directions of other faces and thrie orientation
 803     const labelList& owner = mesh.owner();
 804     forAll(c, fI)
 805     {
 806         if( usedDirection[fI] )
 807             continue;
 808
 809         const bool orientation = owner[c[fI]]==cellI_?false:true;
 810
 811         if( !help::shareAnEdge(faces[c[fI]], faces[c[faceInDirection_[0]]]) )
 812         {
 813             //- face at k = nLayersK_
 814             faceInDirection_[1] = fI;
 815             faceOrientation_[1] = orientation;
 816         }
 817         else if
 818         (
 819             !help::shareAnEdge(faces[c[fI]], faces[c[faceInDirection_[2]]])
 820         )
 821         {
 822             //- face at j = nLayersJ_
 823             faceInDirection_[3] = fI;
 824             faceOrientation_[3] = orientation;
 825         }
 826         else if
 827         (
 828             !help::shareAnEdge(faces[c[fI]], faces[c[faceInDirection_[4]]])
 829         )
 830         {
 831             //- face at i = nLayersI_
 832             faceInDirection_[5] = fI;
 833             faceOrientation_[5] = orientation;
 834         }
 835     }
 836
 837     # ifdef DEBUGLayer
 838     forAll(faceInDirection_, i)
 839         Pout << "Face in direction " << i
 840              << " is " << faces[c[faceInDirection_[i]]]
 841              << " orientation " << faceOrientation_[i] << endl;
 842     Pout << "nLayersI = " << nLayersI_
 843          << " nLayersJ = " << nLayersJ_
 844          << " nLayersK = " << nLayersK_ << endl;
 845     # endif
 846 }
 847
 848 void refineBoundaryLayers::refineCornerHexCell::populateExistingFaces()
 849 {
 850     const cell& c = bndLayers_.mesh_.cells()[cellI_];
 851
 852     //- set the number of cells
 853     cellsFromCell_.setSize(nLayersI_ * nLayersJ_ * nLayersK_);
 854     forAll(cellsFromCell_, i)
 855         cellsFromCell_[i].clear();
 856
 857     //- add new faces from existing faces into new cells
 858     forAll(faceInDirection_, dirI)
 859     {
 860         bndLayers_.storeFacesIntoCells
 861         (
 862             c[faceInDirection_[dirI]], faceOrientation_[dirI],
 863             dirI / 2, dirI % 2,
 864             nLayersI_, nLayersJ_, nLayersK_,
 865             cellsFromCell_
 866         );
 867     }
 868
 869     # ifdef DEBUGLayer
 870     Pout << "cellsFromCell_ before new faces " << cellsFromCell_ << endl;
 871     //::exit(1);
 872     # endif
 873 }
 874
 875 void refineBoundaryLayers::refineCornerHexCell::generateNewPoints()
 876 {
 877     const cell& c = bndLayers_.mesh_.cells()[cellI_];
 878
 879     //- allocate space for points generated inside the cell
 880     cellPoints_.setSize(nLayersI_+1);
 881     forAll(cellPoints_, i)
 882     {
 883         cellPoints_[i].setSize(nLayersJ_+1);
 884
 885         forAll(cellPoints_[i], j)
 886         {
 887             cellPoints_[i][j].setSize(nLayersK_+1);
 888             cellPoints_[i][j] = -1;
 889         }
 890     }
 891
 892     //- collect information about points generated on faces of the cell
 893     forAll(faceInDirection_, dirI)
 894     {
 895         bndLayers_.sortFacePoints
 896         (
 897             c[faceInDirection_[dirI]],
 898             facePoints_[dirI],
 899             faceOrientation_[dirI]
 900         );
 901     }
 902
 903     # ifdef DEBUGLayer
 904     Pout << "Face points " << facePoints_ << endl;
 905     # endif
 906
 907     //- fill in cellPoints at the boundary
 908     forAll(cellPoints_, i)
 909     {
 910         forAll(cellPoints_[i], j)
 911         {
 912             cellPoints_[i][j][0] = facePoints_[0][i][j];
 913             cellPoints_[i][j][nLayersK_] = facePoints_[1][i][j];
 914         }
 915     }
 916
 917     forAll(cellPoints_, i)
 918     {
 919         forAll(cellPoints_[i][0], k)
 920         {
 921             cellPoints_[i][0][k] = facePoints_[2][k][i];
 922             cellPoints_[i][nLayersJ_][k] = facePoints_[3][k][i];
 923         }
 924     }
 925
 926     forAll(cellPoints_[0], j)
 927     {
 928         forAll(cellPoints_[0][j], k)
 929         {
 930             cellPoints_[0][j][k] = facePoints_[4][j][k];
 931             cellPoints_[nLayersI_][j][k] = facePoints_[5][j][k];
 932         }
 933     }
 934
 935     //- useful data for generating missing points
 936     const edgeLongList& splitEdges = bndLayers_.splitEdges_;
 937     const edge& seDirI = splitEdges[splitEdgeInDirection_[0]];
 938     const edge& seDirJ = splitEdges[splitEdgeInDirection_[1]];
 939     const edge& seDirK = splitEdges[splitEdgeInDirection_[2]];
 940     const VRWGraph& ptsAtEdge = bndLayers_.newVerticesForSplitEdge_;
 941
 942     //- const references to vertices of the cell ordered in a local
 943     //- i, j, k coordinate system
 944     pointFieldPMG& points = bndLayers_.mesh_.points();
 945     const point v000 = points[seDirI.start()];
 946     const point v100 = points[seDirI.end()];
 947     const point v110 = points[facePoints_[0].lastElement().lastElement()];
 948     const point v010 = points[seDirJ.end()];
 949     const point v001 = points[seDirK.end()];
 950     const point v101 = points[facePoints_[1].lastElement()[0]];
 951     const point v111 = points[facePoints_[1].lastElement().lastElement()];
 952     const point v011 = points[facePoints_[1][0].lastElement()];
 953
 954     for(label i=1;i<nLayersI_;++i)
 955     {
 956         const scalar u
 957         (
 958             Foam::mag
 959             (
 960                 points[ptsAtEdge(splitEdgeInDirection_[0], i)] -
 961                 points[seDirI.start()]
 962             ) /
 963             seDirI.mag(points)
 964         );
 965
 966         for(label j=1;j<nLayersJ_;++j)
 967         {
 968             const scalar v
 969             (
 970                 Foam::mag
 971                 (
 972                     points[ptsAtEdge(splitEdgeInDirection_[1], j)] -
 973                     points[seDirJ.start()]
 974                 ) /
 975                 seDirJ.mag(points)
 976             );
 977
 978             for(label k=1;k<nLayersK_;++k)
 979             {
 980                 const scalar w
 981                 (
 982                     Foam::mag
 983                     (
 984                         points[ptsAtEdge(splitEdgeInDirection_[2], k)] -
 985                         points[seDirK.start()]
 986                     ) /
 987                     seDirK.mag(points)
 988                 );
 989
 990                 # ifdef DEBUGLayer
 991                 Pout << "Generating point in corner cell local coordinates "
 992                      << "u = " << u << " v = " << v << " w = " << w << endl;
 993                 # endif
 994
 995                 //- calculate coordinates of the new vertex
 996                 const point newP =
 997                     (1.0 - u) * (1.0 - v) * (1.0 - w) * v000 +
 998                     u * (1.0 - v) * (1.0 - w) * v100 +
 999                     u * v * (1.0 - w) * v110 +
1000                     (1.0 - u) * v * (1.0 - w) * v010 +
1001                     (1.0 - u) * (1.0 - v) * w * v001 +
1002                     u * (1.0 - v) * w * v101 +
1003                     u * v * w * v111 +
1004                     (1.0 - u) * v * w * v011;
1005
1006                 # ifdef DEBUGLayer
1007                 Pout << "New point " << points.size() << " in corner hex "
1008                     << "has coordinates " << newP << endl;
1009                 # endif
1010
1011                 //- add the point to the mesh
1012                 cellPoints_[i][j][k] = points.size();
1013                 points.append(newP);
1014             }
1015         }
1016     }
1017
1018     # ifdef DEBUGLayer
1019     Pout << "New cell points " << cellPoints_ << endl;
1020     //::exit(1);
1021     # endif
1022 }
1023
1024 void refineBoundaryLayers::refineCornerHexCell::generateMissingFaces()
1025 {
1026     //- generate face in direction i
1027     for(label i=1;i<nLayersI_;++i)
1028     {
1029         //- generate quad faces
1030         for(label j=0;j<nLayersJ_;++j)
1031         {
1032             for(label k=0;k<nLayersK_;++k)
1033             {
1034                 //- skip generating last face because it might not be a quad
1035                 if( (j == (nLayersJ_-1)) && (k == (nLayersK_-1)) )
1036                     continue;
1037
1038                 const label own
1039                 (
1040                     k * nLayersI_ * nLayersJ_ +
1041                     j * nLayersI_ +
1042                     i - 1
1043                 );
1044                 const label nei = own + 1;
1045
1046                 FixedList<label, 4> mf;
1047
1048                 mf[0] = cellPoints_[i][j][k];
1049                 mf[1] = cellPoints_[i][j+1][k];
1050                 mf[2] = cellPoints_[i][j+1][k+1];
1051                 mf[3] = cellPoints_[i][j][k+1];
1052
1053                 cellsFromCell_[own].append(mf);
1054                 cellsFromCell_[nei].append(help::reverseFace(mf));
1055             }
1056         }
1057
1058         //- generate faces which might not be a quads
1059         DynList<label> mf;
1060
1061         mf.append(cellPoints_[i][nLayersJ_-1][nLayersK_-1]);
1062
1063         //- this face might not be a quad
1064         //- add points fom the last face in direction j
1065         const DynList<DynList<label> >& f3 = facePoints_[3];
1066         for(label index=nLayersK_-1;index<f3.size()-1;++index)
1067             mf.append(f3[index][i]);
1068
1069         //- add points from the last face in direction k
1070         const DynList<DynList<label> >& f1 = facePoints_[1];
1071         for(label index=f1[i].size()-1;index>=nLayersJ_-1;--index)
1072             mf.append(f1[i][index]);
1073
1074         const label own
1075         (
1076             (nLayersK_-1) * nLayersI_ * nLayersJ_ +
1077             (nLayersJ_-1) * nLayersI_ +
1078             i - 1
1079         );
1080
1081         const label nei = own + 1;
1082
1083         # ifdef DEBUGLayer
1084         Pout << "Additional face in direction i = " << i
1085              << " j = " << (nLayersJ_-1)
1086              << " has owner " << own
1087              << " neighbour " << nei << " with nodes " << mf << endl;
1088         # endif
1089
1090         cellsFromCell_[own].append(mf);
1091         cellsFromCell_[nei].append(help::reverseFace(mf));
1092     }
1093
1094     //- generate faces in direction j
1095     for(label j=1;j<nLayersJ_;++j)
1096     {
1097         //- generate quad faces
1098         for(label i=0;i<nLayersI_;++i)
1099         {
1100             for(label k=0;k<nLayersK_;++k)
1101             {
1102                 //- skip generating late face because it might not be a quad
1103                 if( (i == (nLayersI_-1)) && (k == (nLayersK_-1)) )
1104                     continue;
1105
1106                 const label own
1107                 (
1108                     k * nLayersI_ * nLayersJ_ +
1109                     (j-1) * nLayersI_ +
1110                     i
1111                 );
1112
1113                 const label nei
1114                 (
1115                     k * nLayersI_ * nLayersJ_ +
1116                     j * nLayersI_ +
1117                     i
1118                 );
1119
1120                 FixedList<label, 4> mf;
1121
1122                 mf[0] = cellPoints_[i][j][k];
1123                 mf[1] = cellPoints_[i][j][k+1];
1124                 mf[2] = cellPoints_[i+1][j][k+1];
1125                 mf[3] = cellPoints_[i+1][j][k];
1126
1127                 cellsFromCell_[own].append(mf);
1128                 cellsFromCell_[nei].append(help::reverseFace(mf));
1129             }
1130         }
1131
1132         //- generate a face which might not be a quad
1133         DynList<label> mf;
1134
1135         mf.append(cellPoints_[nLayersI_-1][j][nLayersK_-1]);
1136
1137         //- add points from the last face in direction k
1138         const DynList<DynList<label> >& fp1 = facePoints_[1];
1139         for(label index=nLayersI_-1;index<fp1.size()-1;++index)
1140             mf.append(fp1[index][j]);
1141
1142         //- add points from the last face in direction i
1143         const DynList<DynList<label> >& fp5 = facePoints_[5];
1144         for(label index=fp5[j].size()-1;index>=nLayersK_-1;--index)
1145             mf.append(fp5[j][index]);
1146
1147         const label own
1148         (
1149             (nLayersK_-1) * nLayersI_ * nLayersJ_ +
1150             (j-1) * nLayersI_ +
1151             (nLayersI_ - 1)
1152         );
1153
1154         const label nei
1155         (
1156             (nLayersK_-1) * nLayersI_ * nLayersJ_ +
1157             j * nLayersI_ +
1158             (nLayersI_ - 1)
1159         );
1160
1161         # ifdef DEBUGLayer
1162         Pout << "Additional face at i = " << (nLayersI_-1)
1163              << " j = " << j << " k = " << (nLayersK_-1)
1164              << " has owner " << own
1165              << " neighbour " << nei << " with nodes " << mf << endl;
1166         # endif
1167
1168         cellsFromCell_[own].append(mf);
1169         cellsFromCell_[nei].append(help::reverseFace(mf));
1170     }
1171
1172     //- generate faces in direction k
1173     for(label k=1;k<nLayersK_;++k)
1174     {
1175         //- generate quad faces
1176         for(label i=0;i<nLayersI_;++i)
1177         {
1178             for(label j=0;j<nLayersJ_;++j)
1179             {
1180                 //- skip the last face because it might not be a quad
1181                 if( (i == (nLayersI_-1)) && (j == (nLayersJ_-1)) )
1182                     continue;
1183
1184                 const label own
1185                 (
1186                     (k-1) * nLayersI_ * nLayersJ_ +
1187                     j * nLayersI_ +
1188                     i
1189                 );
1190
1191                 const label nei
1192                 (
1193                     k * nLayersI_ * nLayersJ_ +
1194                     j * nLayersI_ +
1195                     i
1196                 );
1197
1198                 FixedList<label, 4> mf;
1199
1200                 mf[0] = cellPoints_[i][j][k];
1201                 mf[1] = cellPoints_[i+1][j][k];
1202                 mf[2] = cellPoints_[i+1][j+1][k];
1203                 mf[3] = cellPoints_[i][j+1][k];
1204
1205                 cellsFromCell_[own].append(mf);
1206                 cellsFromCell_[nei].append(help::reverseFace(mf));
1207             }
1208         }
1209
1210         //- generate a face which might not be a quad
1211         DynList<label> mf;
1212
1213         mf.append(cellPoints_[nLayersI_-1][nLayersJ_-1][k]);
1214
1215         //- this face might not be a quad
1216         //- add points from the last face in direction i
1217         const DynList<DynList<label> >& fp5 = facePoints_[5];
1218         for(label index=nLayersJ_-1;index<fp5.size()-1;++index)
1219             mf.append(fp5[index][k]);
1220
1221         //- add points from the last face in direction j
1222         const DynList<DynList<label> >& fp3 = facePoints_[3];
1223         for(label index=fp3[k].size()-1;index>=nLayersI_-1;--index)
1224             mf.append(fp3[k][index]);
1225
1226         const label own
1227         (
1228             (k-1) * nLayersI_ * nLayersJ_ +
1229             (nLayersJ_-1) * nLayersI_ +
1230             (nLayersI_ - 1)
1231         );
1232
1233         const label nei
1234         (
1235             k * nLayersI_ * nLayersJ_ +
1236             (nLayersJ_-1) * nLayersI_ +
1237             (nLayersI_ - 1)
1238         );
1239
1240         # ifdef DEBUGLayer
1241         Pout << "Additional face at position i = " << (nLayersI_-1)
1242              << " j = " << (nLayersJ_-1) << " k = " << k
1243              << " has owner " << own
1244              << " neighbour " << nei << " with nodes " << mf << endl;
1245         # endif
1246
1247         cellsFromCell_[own].append(mf);
1248         cellsFromCell_[nei].append(help::reverseFace(mf));
1249     }
1250
1251     # ifdef DEBUGLayer
1252     Pout << "Generated cells " << cellsFromCell_ << endl;
1253
1254     forAll(cellsFromCell_, cI)
1255     {
1256         const DynList<DynList<label, 4>, 6>& cellFaces = cellsFromCell_[cI];
1257
1258         DynList<edge, 12> edges;
1259         DynList<label, 12> nAppearances;
1260
1261         forAll(cellFaces, fI)
1262         {
1263             const DynList<label, 4>& f = cellFaces[fI];
1264
1265             forAll(f, eI)
1266             {
1267                 const edge e(f[eI], f.fcElement(eI));
1268
1269                 const label pos = edges.containsAtPosition(e);
1270
1271                 if( pos < 0 )
1272                 {
1273                     edges.append(e);
1274                     nAppearances.append(1);
1275                 }
1276                 else
1277                 {
1278                     ++nAppearances[pos];
1279                 }
1280             }
1281         }
1282
1283         forAll(nAppearances, eI)
1284             if( nAppearances[eI] != 2 )
1285             {
1286                 Pout << "Edge hex cell " << cI << " edge " << edges[eI]
1287                     << " is present " << nAppearances[eI] << " times!" << endl;
1288                 abort(FatalError);
1289             }
1290     }
1291
1292     //::exit(1);
1293     # endif
1294 }
1295
1296 refineBoundaryLayers::refineCornerHexCell::refineCornerHexCell
1297 (
1298     const label cellI,
1299     const refineBoundaryLayers& ref
1300 )
1301 :
1302     cellI_(cellI),
1303     nLayersI_(),
1304     nLayersJ_(),
1305     nLayersK_(),
1306     splitEdgeInDirection_(),
1307     cellsFromCell_(),
1308     bndLayers_(ref),
1309     faceInDirection_(),
1310     faceOrientation_(),
1311     facePoints_(),
1312     cellPoints_()
1313 {
1314     determineFacesInDirections();
1315
1316     populateExistingFaces();
1317
1318     generateNewPoints();
1319
1320     generateMissingFaces();
1321 }
1322
1323 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
1324
1325 void refineBoundaryLayers::generateNewCells()
1326 {
1327     labelList nCellsFromCell(mesh_.cells().size(), 1);
1328     labelList refType(mesh_.cells().size(), 0);
1329
1330     const meshSurfaceEngine& mse = surfaceEngine();
1331     const labelList& faceOwners = mse.faceOwners();
1332
1333     //- calculate the number new cells generated from a cell
1334     forAll(faceOwners, bfI)
1335     {
1336         const label cellI = faceOwners[bfI];
1337
1338         nCellsFromCell[cellI] *= nLayersAtBndFace_[bfI];
1339
1340         if( nLayersAtBndFace_[bfI] > 1 )
1341             ++refType[cellI];
1342     }
1343
1344     //- add cells which shall be refined in a subset
1345     if( cellSubsetName_ != "" )
1346     {
1347         label subsetI = mesh_.cellSubsetIndex(cellSubsetName_);
1348         if( subsetI >= 0 )
1349             Warning << "The subset with name " << cellSubsetName_
1350                     << " already exists. Skipping generation of a new subset"
1351                     << endl;
1352
1353         subsetI = mesh_.addCellSubset(cellSubsetName_);
1354
1355         forAll(nCellsFromCell, cI)
1356             if( nCellsFromCell[cI] > 1 )
1357                 mesh_.addCellToSubset(subsetI, cI);
1358     }
1359
1360     //- check the number of cells which will be generated
1361     label nNewCells(0);
1362     forAll(nCellsFromCell, cellI)
1363         nNewCells += (nCellsFromCell[cellI] - 1);
1364
1365     # ifdef DEBUGLayer
1366     forAll(nCellsFromCell, cellI)
1367     {
1368         Pout << "\nCell " << cellI << endl;
1369         Pout << "nCellsFromCell " << nCellsFromCell[cellI] << endl;
1370         Pout << "Ref type " << refType[cellI] << endl;
1371     }
1372     #  endif
1373
1374     const label totalNumNewCells = returnReduce(nNewCells, sumOp<label>());
1375     Info << "Number of newly generated cells " << totalNumNewCells << endl;
1376
1377     //- create mesh modifier
1378     polyMeshGenModifier meshModifier(mesh_);
1379     faceListPMG& faces = meshModifier.facesAccess();
1380
1381     const label numFacesBefore = newFaces_.size();
1382
1383     //- set the number of cells to the new value
1384     cellListPMG& cells = meshModifier.cellsAccess();
1385     label nCells = cells.size();
1386     cells.setSize(nCells+nNewCells);
1387
1388     //- start creating new cells
1389     //- store the information which new cells were generated from
1390     //- an existing cell
1391     VRWGraph newCellsFromCell(refType.size());
1392
1393     VRWGraph pointNewFaces;
1394     pointNewFaces.reverseAddressing(newFaces_);
1395
1396     forAll(nCellsFromCell, cellI)
1397     {
1398         if( refType[cellI] == 0 )
1399         {
1400             //- this cell is not refined
1401             //- update face labels
1402             newCellsFromCell.append(cellI, cellI);
1403
1404             cell& c = cells[cellI];
1405
1406             //- copy the new faces of this cell
1407             DynList<label, 64> newC;
1408             forAll(c, fI)
1409             {
1410                 forAllRow(facesFromFace_, c[fI], cfI)
1411                     newC.append(facesFromFace_(c[fI], cfI));
1412             }
1413
1414             //- update the cell
1415             c.setSize(newC.size());
1416             forAll(c, fI)
1417                 c[fI] = newC[fI];
1418         }
1419         else if( refType[cellI] == 1 )
1420         {
1421             //- generate new cells from this prism refined in one direction
1422             DynList<DynList<DynList<label, 8>, 10>, 64> cellsFromCell;
1423             generateNewCellsPrism(cellI, cellsFromCell);
1424
1425             forAll(cellsFromCell, cI)
1426             {
1427                 const DynList<DynList<label, 8>, 10>& nc = cellsFromCell[cI];
1428
1429                 const label newCellI = cI==0?cellI:nCells++;
1430
1431                 newCellsFromCell.append(cellI, newCellI);
1432
1433                 cell& c = cells[newCellI];
1434                 c.setSize(nc.size());
1435
1436                 //- find face labels for this cell
1437                 forAll(nc, fI)
1438                 {
1439                     const DynList<label, 8>& nf = nc[fI];
1440
1441                     label faceLabel(-1);
1442                     forAllRow(pointNewFaces, nf[0], pfI)
1443                     {
1444                         const label nfI = pointNewFaces(nf[0], pfI);
1445
1446                         if( help::areFacesEqual(nf, newFaces_[nfI]) )
1447                         {
1448                             c[fI] = nfI;
1449                             faceLabel = nfI;
1450                             break;
1451                         }
1452                     }
1453
1454                     if( faceLabel < 0 )
1455                     {
1456                         forAll(nf, pI)
1457                             pointNewFaces.append(nf[pI], newFaces_.size());
1458                         c[fI] = newFaces_.size();
1459                         newFaces_.appendList(nf);
1460                     }
1461                 }
1462             }
1463         }
1464         else if( refType[cellI] == 2 )
1465         {
1466             //- generate new cell from a hex cell where two layers intersect
1467             //- generate mostly hex cells;
1468             refineEdgeHexCell refEdgeHex(cellI, *this);
1469             const DynList<DynList<DynList<label, 4>, 6>, 256>& cellsFromCell =
1470                 refEdgeHex.newCells();
1471
1472             forAll(cellsFromCell, cI)
1473             {
1474                 const DynList<DynList<label, 4>, 6>& nc = cellsFromCell[cI];
1475
1476                 # ifdef DEBUGLayer
1477                 Pout << "Adding cell " << (cI==0?cellI:nCells)
1478                      << " originating from cell " << cellI << endl;
1479                 # endif
1480
1481                 const label newCellI = cI==0?cellI:nCells++;
1482
1483                 newCellsFromCell.append(cellI, newCellI);
1484
1485                 cell& c = cells[newCellI];
1486                 c.setSize(nc.size());
1487
1488                 //- find face labels for this cell
1489                 forAll(nc, fI)
1490                 {
1491                     const DynList<label, 4>& nf = nc[fI];
1492
1493                     label faceLabel(-1);
1494                     forAllRow(pointNewFaces, nf[0], pfI)
1495                     {
1496                         const label nfI = pointNewFaces(nf[0], pfI);
1497
1498                         if( help::areFacesEqual(nf, newFaces_[nfI]) )
1499                         {
1500                             c[fI] = nfI;
1501                             faceLabel = nfI;
1502                             break;
1503                         }
1504                     }
1505
1506                     if( faceLabel < 0 )
1507                     {
1508                         forAll(nf, pI)
1509                             pointNewFaces.append(nf[pI], newFaces_.size());
1510                         c[fI] = newFaces_.size();
1511                         newFaces_.appendList(nf);
1512                     }
1513                 }
1514             }
1515         }
1516         else if( refType[cellI] == 3 )
1517         {
1518             //- generate new cells from a hex at a corner where three
1519             //- layers intersect
1520             //- generate mostly hex cells
1521             refineCornerHexCell refCell(cellI, *this);
1522             const DynList<DynList<DynList<label, 4>, 6>, 256>& cellsFromCell =
1523                 refCell.newCells();
1524
1525             //- new points have been generated
1526             pointNewFaces.setSize(mesh_.points().size());
1527
1528             //- recognise face cells in the graph of new faces
1529             forAll(cellsFromCell, cI)
1530             {
1531                 const DynList<DynList<label, 4>, 6>& nc = cellsFromCell[cI];
1532
1533                 const label newCellI = cI==0?cellI:nCells++;
1534
1535                 newCellsFromCell.append(cellI, newCellI);
1536
1537                 cell& c = cells[newCellI];
1538                 c.setSize(nc.size());
1539
1540                 //- find face labels for this cell
1541                 forAll(nc, fI)
1542                 {
1543                     const DynList<label, 4>& nf = nc[fI];
1544
1545                     label faceLabel(-1);
1546                     forAllRow(pointNewFaces, nf[0], pfI)
1547                     {
1548                         const label nfI = pointNewFaces(nf[0], pfI);
1549
1550                         if( help::areFacesEqual(nf, newFaces_[nfI]) )
1551                         {
1552                             c[fI] = nfI;
1553                             faceLabel = nfI;
1554                             break;
1555                         }
1556                     }
1557
1558                     if( faceLabel < 0 )
1559                     {
1560                         forAll(nf, pI)
1561                             pointNewFaces.append(nf[pI], newFaces_.size());
1562                         c[fI] = newFaces_.size();
1563                         newFaces_.appendList(nf);
1564                     }
1565                 }
1566             }
1567         }
1568         else
1569         {
1570             FatalErrorIn
1571             (
1572                 "void refineBoundaryLayers::generateNewCells()"
1573             ) << "Cannot refine boundary layer for cell "
1574               << cellI << abort(FatalError);
1575         }
1576     }
1577
1578     //- check the orientation of new faces, because owner and neighbour cells
1579     //- may require a face to be flipped
1580     const label nOrigInternalFaces = mesh_.nInternalFaces();
1581
1582     # ifdef USE_OMP
1583     # pragma omp parallel
1584     # endif
1585     {
1586         const labelList& owner = mesh_.owner();
1587         const labelList& neighbour = mesh_.neighbour();
1588
1589         # ifdef USE_OMP
1590         # pragma omp for schedule(dynamic, 40)
1591         # endif
1592         for(label fI=0;fI<nOrigInternalFaces;++fI)
1593         {
1594             const label own = owner[fI];
1595             const label nei = neighbour[fI];
1596
1597             if( facesFromFace_.sizeOfRow(fI) == 1 )
1598                 continue;
1599
1600             forAllRow(facesFromFace_, fI, cfI)
1601             {
1602                 const label nfI = facesFromFace_(fI, cfI);
1603
1604                 //- find the new owner and neighbour cells of the new face
1605                 label newOwner(-1), newNeighbour(-1);
1606                 forAllRow(newCellsFromCell, own, cI)
1607                 {
1608                     const cell& cOwn = cells[newCellsFromCell(own, cI)];
1609
1610                     const label pos = help::positionInList(nfI, cOwn);
1611
1612                     if( pos >= 0 )
1613                     {
1614                         newOwner = newCellsFromCell(own, cI);
1615                         break;
1616                     }
1617                 }
1618
1619                 forAllRow(newCellsFromCell, nei, cI)
1620                 {
1621                     const cell& cNei = cells[newCellsFromCell(nei, cI)];
1622
1623                     const label pos = help::positionInList(nfI, cNei);
1624
1625                     if( pos >= 0 )
1626                     {
1627                         newNeighbour = newCellsFromCell(nei, cI);
1628                         break;
1629                     }
1630                 }
1631
1632                 if( newOwner > newNeighbour )
1633                 {
1634                     DynList<label> rf;
1635                     rf.setSize(newFaces_.sizeOfRow(nfI));
1636
1637                     forAll(rf, i)
1638                         rf[i] = newFaces_(nfI, i);
1639
1640                     rf = help::reverseFace(rf);
1641
1642                     newFaces_.setRow(nfI, rf);
1643                 }
1644             }
1645         }
1646     }
1647
1648     //- update cell sets
1649     mesh_.updateCellSubsets(newCellsFromCell);
1650     newCellsFromCell.setSize(0);
1651
1652     //- point-faces addressing is not needed any more
1653     pointNewFaces.setSize(0);
1654
1655     //- copy newFaces to the mesh
1656     # ifdef DEBUGLayer
1657     Pout << "Copying internal faces " << endl;
1658     Pout << "Original number of internal faces " << nOrigInternalFaces << endl;
1659     # endif
1660
1661     //- store internal faces originating from existing faces
1662     labelLongList newFaceLabel(newFaces_.size());
1663     faces.setSize(newFaces_.size());
1664
1665     label currFace = 0;
1666     for(label faceI=0;faceI<nOrigInternalFaces;++faceI)
1667     {
1668         forAllRow(facesFromFace_, faceI, ffI)
1669         {
1670             face& f = faces[currFace];
1671             newFaceLabel[currFace] = currFace;
1672             ++currFace;
1673
1674             const label newFaceI = facesFromFace_(faceI, ffI);
1675
1676             f.setSize(newFaces_.sizeOfRow(newFaceI));
1677
1678             forAll(f, pI)
1679                 f[pI] = newFaces_(newFaceI, pI);
1680         }
1681     }
1682
1683     //- store newly-generated internal faces
1684     # ifdef DEBUGLayer
1685     Pout << "Copying newly generated internal faces" << endl;
1686     Pout << "nNewInternalFaces " << currFace << endl;
1687     Pout << "numFacesBefore " << numFacesBefore << endl;
1688     Pout << "Total number of faces " << newFaces_.size() << endl;
1689     # endif
1690
1691     for(label faceI=numFacesBefore;faceI<newFaces_.size();++faceI)
1692     {
1693         newFaceLabel[faceI] = currFace;
1694         face& f = faces[currFace];
1695         ++currFace;
1696
1697         f.setSize(newFaces_.sizeOfRow(faceI));
1698
1699         forAll(f, pI)
1700             f[pI] = newFaces_(faceI, pI);
1701     }
1702
1703     //- store new boundary faces
1704     # ifdef DEBUGLayer
1705     Pout << "Copying boundary faces " << endl;
1706     Pout << "currFace " << currFace << endl;
1707     Pout << "Faces size " << faces.size() << endl;
1708     Pout << "Initial number of faces " << facesFromFace_.size() << endl;
1709     # endif
1710
1711     PtrList<boundaryPatch>& boundaries = meshModifier.boundariesAccess();
1712     forAll(boundaries, patchI)
1713     {
1714         const label start = boundaries[patchI].patchStart();
1715         const label size = boundaries[patchI].patchSize();
1716
1717         const label newStart = currFace;
1718         label nNewFacesInPatch(0);
1719         for(label fI=0;fI<size;++fI)
1720         {
1721             const label faceI = start + fI;
1722
1723             forAllRow(facesFromFace_, faceI, nfI)
1724             {
1725                 face& f = faces[currFace];
1726
1727                 //- update the new label
1728                 const label origFaceI = facesFromFace_(faceI, nfI);
1729                 newFaceLabel[origFaceI] = currFace;
1730                 facesFromFace_(faceI, nfI) = currFace;
1731                 ++currFace;
1732
1733                 //- copy the face into the mesh
1734                 f.setSize(newFaces_.sizeOfRow(origFaceI));
1735                 forAll(f, pI)
1736                     f[pI] = newFaces_(origFaceI, pI);
1737
1738                 ++nNewFacesInPatch;
1739             }
1740         }
1741
1742         //- update patch
1743         boundaries[patchI].patchStart() = newStart;
1744         boundaries[patchI].patchSize() = nNewFacesInPatch;
1745     }
1746
1747     if( Pstream::parRun() )
1748     {
1749         # ifdef DEBUGLayer
1750         Pout << "Copying processor faces" << endl;
1751         # endif
1752
1753         //- copy faces at inter-processor boundaries
1754         PtrList<processorBoundaryPatch>& procBoundaries =
1755             meshModifier.procBoundariesAccess();
1756
1757         forAll(procBoundaries, patchI)
1758         {
1759             const label start = procBoundaries[patchI].patchStart();
1760             const label size = procBoundaries[patchI].patchSize();
1761
1762             const label newStart = currFace;
1763             label nNewFacesInPatch(0);
1764             for(label fI=0;fI<size;++fI)
1765             {
1766                 const label faceI = start + fI;
1767                 forAllRow(facesFromFace_, faceI, nfI)
1768                 {
1769                     face& f = faces[currFace];
1770
1771                     //- update the new label
1772                     const label origFaceI = facesFromFace_(faceI, nfI);
1773                     newFaceLabel[origFaceI] = currFace;
1774                     facesFromFace_(faceI, nfI) = currFace;
1775                     ++currFace;
1776
1777                     //- copy the face into the mesh
1778                     f.setSize(newFaces_.sizeOfRow(origFaceI));
1779                     forAll(f, pI)
1780                         f[pI] = newFaces_(origFaceI, pI);
1781
1782                     ++nNewFacesInPatch;
1783                 }
1784             }
1785
1786             //- update patch
1787             procBoundaries[patchI].patchStart() = newStart;
1788             procBoundaries[patchI].patchSize() = nNewFacesInPatch;
1789         }
1790     }
1791
1792     # ifdef DEBUGLayer
1793     Pout << "Faces after refinement " << faces << endl;
1794     Pout << "newFaceLabel " << newFaceLabel << endl;
1795     # endif
1796
1797     //- update face subsets
1798     mesh_.updateFaceSubsets(facesFromFace_);
1799     facesFromFace_.setSize(0);
1800     newFaces_.setSize(0);
1801
1802     //- update cells to match the faces
1803     # ifdef DEBUGLayer
1804     Pout << "Updating cells to match new faces" << endl;
1805     # endif
1806
1807     forAll(cells, cellI)
1808     {
1809         cell& c = cells[cellI];
1810
1811         forAll(c, fI)
1812             c[fI] = newFaceLabel[c[fI]];
1813     }
1814
1815     # ifdef DEBUGLayer
1816     Pout << "Cleaning mesh " << endl;
1817     # endif
1818
1819     //- delete all adressing which is no longer up-to-date
1820     meshModifier.clearAll();
1821     deleteDemandDrivenData(msePtr_);
1822
1823     # ifdef DEBUGLayer
1824     for(label procI=0;procI<Pstream::nProcs();++procI)
1825     {
1826         if( procI == Pstream::myProcNo() )
1827         {
1828             forAll(cells, cellI)
1829             {
1830                 const cell& c = cells[cellI];
1831
1832                 DynList<edge> edges;
1833                 DynList<label> nAppearances;
1834                 forAll(c, fI)
1835                 {
1836                     const face& f = faces[c[fI]];
1837
1838                     forAll(f, eI)
1839                     {
1840                         const edge e = f.faceEdge(eI);
1841
1842                         const label pos = edges.containsAtPosition(e);
1843
1844                         if( pos < 0 )
1845                         {
1846                             edges.append(e);
1847                             nAppearances.append(1);
1848                         }
1849                         else
1850                         {
1851                             ++nAppearances[pos];
1852                         }
1853                     }
1854                 }
1855
1856                 bool badCell(false);
1857                 forAll(nAppearances, i)
1858                     if( nAppearances[i] != 2 )
1859                     {
1860                         badCell = true;
1861                         break;
1862
1863                     }
1864
1865                 if( badCell )
1866                 {
1867                     Pout << "Cell " << cellI
1868                          << " is not topologically closed" << endl;
1869
1870                     forAll(c, fI)
1871                         Pout << "Face " << c[fI] << " with points "
1872                              << faces[c[fI]] << endl;
1873                     Pout << "Cell edges " << edges << endl;
1874                     Pout << "nAppearances " << nAppearances << endl;
1875                     ::exit(1);
1876                 }
1877             }
1878         }
1879
1880         returnReduce(1, sumOp<label>());
1881     }
1882
1883     const labelList& owner = mesh_.owner();
1884     const labelList& neighbour = mesh_.neighbour();
1885     const label nInternalFaces = mesh_.nInternalFaces();
1886
1887     for(label procI=0;procI<Pstream::nProcs();++procI)
1888     {
1889         if( procI == Pstream::myProcNo() )
1890         {
1891             forAll(faces, faceI)
1892             {
1893                 if( faceI < nInternalFaces && neighbour[faceI] < 0 )
1894                 {
1895                     Pout << "Num interface faces " << nInternalFaces
1896                          << " current face " << faceI
1897                          << " face points " << faces[faceI] << endl;
1898                     ::exit(1);
1899                 }
1900                 Pout << "Face " << faceI << " owner " << owner[faceI]
1901                      << " neighbour " << neighbour[faceI]
1902                      << " face points " << faces[faceI] << endl;
1903             }
1904
1905             forAll(cells, cellI)
1906                 Pout << "Cell " << cellI << " has faces "
1907                      << cells[cellI] << endl;
1908         }
1909
1910         returnReduce(procI, maxOp<label>());
1911     }
1912     # endif
1913
1914     Info << "Finished generating new cells " << endl;
1915 }
1916
1917 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
1918
1919 } // End namespace Foam
1920
1921 // ************************************************************************* //