src/meshLibrary/utilities/boundaryLayers/refineBoundaryLayersFaces.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 "demandDrivenData.H"
  31 #include "FixedList.H"
  32 #include "helperFunctions.H"
  33
  34 //#define DEBUGLayer
  35
  36 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  37
  38 namespace Foam
  39 {
  40
  41 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  42
  43 void refineBoundaryLayers::refineFace
  44 (
  45     const face& f,
  46     const FixedList<label, 2>& nLayersInDirection,
  47     DynList<DynList<label, 4>, 128>& newFaces
  48 )
  49 {
  50     //- this face must be a quad
  51     if( f.size() != 4 )
  52     {
  53         WarningIn
  54         (
  55             "void refineBoundaryLayers::refineFace(const face&,"
  56             " const FixedList<label, 2>&, DynList<DynList<label, 4> >&)"
  57         ) << "Face " << f << " is not a quad" << endl;
  58         return;
  59     }
  60
  61     //- direction 0 represents edges 0 and 2
  62     //- direction 1 represents edges 1 and 3
  63     if( (nLayersInDirection[0] <= 1) && (nLayersInDirection[1] <= 1) )
  64     {
  65         //- this face may comprise of some split edges
  66         DynList<label, 64> newF;
  67         forAll(f, eI)
  68         {
  69             const edge e = f.faceEdge(eI);
  70
  71             //- add the current point label
  72             newF.append(f[eI]);
  73
  74             //- check if a split edge matches this face edge
  75             forAllRow(splitEdgesAtPoint_, f[eI], peI)
  76             {
  77                 const label seI = splitEdgesAtPoint_(f[eI], peI);
  78                 const edge& se = splitEdges_[seI];
  79
  80                 if( e == se )
  81                 {
  82                     //- check the orientation and add new vertices created
  83                     //- on this edge
  84                     const label s = newVerticesForSplitEdge_.sizeOfRow(seI) - 1;
  85                     if( e.start() == se.start() )
  86                     {
  87                         for(label pI=1;pI<s;++pI)
  88                             newF.append(newVerticesForSplitEdge_(seI, pI));
  89                     }
  90                     else
  91                     {
  92                         for(label pI=s-1;pI>0;--pI)
  93                             newF.append(newVerticesForSplitEdge_(seI, pI));
  94                     }
  95                 }
  96             }
  97         }
  98
  99         newFaces.setSize(1);
 100         newFaces[0] = newF;
 101         return;
 102     }
 103
 104     //- check which face edge is a direction 0 and which one is a direction 1
 105     label dir0(-1), dir1(-1);
 106     labelPair dir0Edges(-1, -1), dir1Edges(-1, -1);
 107     forAll(f, eI)
 108     {
 109         const edge e = f.faceEdge(eI);
 110
 111         label ses(-1), see(-1);
 112         bool start(false), end(false);
 113         forAllRow(splitEdgesAtPoint_, e.start(), i)
 114         {
 115             const edge& se = splitEdges_[splitEdgesAtPoint_(e.start(), i)];
 116
 117             if( (se.start() == e.start()) && (se.end() == f.prevLabel(eI)) )
 118             {
 119                 ses = splitEdgesAtPoint_(e.start(), i);
 120                 start = true;
 121                 break;
 122             }
 123         }
 124
 125         forAllRow(splitEdgesAtPoint_, e.end(), i)
 126         {
 127             const edge& se = splitEdges_[splitEdgesAtPoint_(e.end(), i)];
 128
 129             if( (se.start() == e.end()) && (se.end() == f[(eI+2)%4]) )
 130             {
 131                 see = splitEdgesAtPoint_(e.end(), i);
 132                 end = true;
 133                 break;
 134             }
 135         }
 136
 137         if( start && end )
 138         {
 139             if( dir0 == -1 )
 140             {
 141                 dir0 = eI;
 142                 dir0Edges = labelPair(ses, see);
 143             }
 144             else if( dir1 == -1 )
 145             {
 146                 dir1 = eI;
 147                 dir1Edges = labelPair(ses, see);
 148             }
 149             else
 150             {
 151                 FatalErrorIn
 152                 (
 153                     "void refineBoundaryLayers::refineFace(const face&,"
 154                     " const FixedList<label, 2>&, DynList<DynList<label, 4> >&)"
 155                 ) << "More than two split directions for a face"
 156                   << abort(FatalError);
 157             }
 158         }
 159     }
 160
 161     # ifdef DEBUGLayer
 162     Pout << "Refining face " << f << endl;
 163     Pout << "Splits in direction " << nLayersInDirection << endl;
 164     Pout << "Here " << endl;
 165     Pout << "Dir0 " << dir0 << endl;
 166     Pout << "dir0Edges " << dir0Edges << endl;
 167     Pout << "Dir1 " << dir1 << endl;
 168     Pout << "dir1Edges " << dir1Edges << endl;
 169     # endif
 170
 171     if( (dir0 < 0) && (dir1 < 0) )
 172     {
 173         Pout << "Refining face " << f << endl;
 174         forAll(f, pI)
 175         {
 176             if( splitEdgesAtPoint_.size() >= f[pI] )
 177             Pout << "Split edges at point " << f[pI]
 178                  << " are " << splitEdgesAtPoint_[f[pI]] << endl;
 179         }
 180         Pout << "Splits in direction " << nLayersInDirection << endl;
 181         Pout << "Here " << endl;
 182         Pout << "Dir0 " << dir0 << endl;
 183         Pout << "dir0Edges " << dir0Edges << endl;
 184         Pout << "Dir1 " << dir1 << endl;
 185         Pout << "dir1Edges " << dir1Edges << endl;
 186
 187         FatalErrorIn
 188         (
 189             "void refineBoundaryLayers::refineFace(const face&,"
 190             " const FixedList<label, 2>&, DynList<DynList<label, 4> >&)"
 191         ) << "Cannot find split edges for a face" << abort(FatalError);
 192     }
 193
 194     //- in case of only one refinement direction, it must direction 0
 195     if( (dir1 != -1) && (dir0 == -1) )
 196     {
 197         dir0 = dir1;
 198         dir0Edges = dir1Edges;
 199         dir1 = -1;
 200     }
 201     else if( (dir0 != -1) && (dir1 != -1) && (dir1 != f.fcIndex(dir0)) )
 202     {
 203         //- alternate value to preserve correct face orientation
 204         const label add = dir0;
 205         dir0 = dir1;
 206         dir1 = add;
 207
 208         const labelPair lpAdd = dir0Edges;
 209         dir0Edges = dir1Edges;
 210         dir1Edges = lpAdd;
 211     }
 212
 213     //- permutate the number of refinements in each direction
 214     const label nLayersDir0 = dir0>=0?nLayersInDirection[dir0%2]:1;
 215     const label nLayersDir1 = dir1>=0?nLayersInDirection[dir1%2]:1;
 216
 217     # ifdef DEBUGLayer
 218     Pout << "Face has points " << f << endl;
 219     Pout << "dirEdges0 " << dir0Edges << endl;
 220     Pout << "dir1Edges " << dir1Edges << endl;
 221     if( dir0 >= 0 )
 222     {
 223         Pout << "Points on edge " << dir0Edges.first() << " with nodes "
 224              << splitEdges_[dir0Edges.first()]
 225              << " are " << newVerticesForSplitEdge_[dir0Edges.first()] << endl;
 226         Pout << "Points on edge " << dir0Edges.second() << " with nodes "
 227              << splitEdges_[dir0Edges.second()]
 228              << " are " << newVerticesForSplitEdge_[dir0Edges.second()] << endl;
 229     }
 230     if( dir1 >= 0 )
 231     {
 232         Pout << "Points on edge " << dir1Edges.first() << " with nodes "
 233              << splitEdges_[dir1Edges.first()]
 234              << " are " << newVerticesForSplitEdge_[dir1Edges.first()] << endl;
 235         Pout << "Points on edge " << dir1Edges.second() << " with nodes "
 236              << splitEdges_[dir1Edges.second()]
 237              << " are " << newVerticesForSplitEdge_[dir1Edges.second()] << endl;
 238     }
 239     Pout << "nLayersDir0 " << nLayersDir0 << endl;
 240     Pout << "nLayersDir1 " << nLayersDir1 << endl;
 241     # endif
 242
 243     //- map the face onto a matrix for easier orientation
 244     DynList<DynList<label> > facePoints;
 245     facePoints.setSize(nLayersDir0+1);
 246     forAll(facePoints, i)
 247     {
 248         facePoints[i].setSize(nLayersDir1+1);
 249         facePoints[i] = -1;
 250     }
 251
 252     //- add points in the matrix
 253     for(label i=0;i<nLayersDir0;++i)
 254     {
 255         facePoints[i][0] = newVerticesForSplitEdge_(dir0Edges.second(), i);
 256         facePoints[i][nLayersDir1] =
 257             newVerticesForSplitEdge_(dir0Edges.first(), i);
 258     }
 259     facePoints[nLayersDir0][0] = splitEdges_[dir0Edges.second()].end();
 260     facePoints[nLayersDir0][nLayersDir1] = splitEdges_[dir0Edges.first()].end();
 261
 262     for(label i=1;i<nLayersDir1;++i)
 263     {
 264         facePoints[0][i] = newVerticesForSplitEdge_(dir1Edges.first(), i);
 265         facePoints[nLayersDir0][i] =
 266             newVerticesForSplitEdge_(dir1Edges.second(), i);
 267     }
 268
 269     //- create missing vertices if there are any
 270     pointFieldPMG& points = mesh_.points();
 271     const point v00 = points[facePoints[0][0]];
 272     const point v10 = points[facePoints[nLayersDir0][0]];
 273     const point v01 = points[facePoints[0][nLayersDir1]];
 274     const point v11 = points[facePoints[nLayersDir0][nLayersDir1]];
 275
 276     # ifdef DEBUGLayer
 277     forAll(points, pointI)
 278         Pout << "Point " << pointI << " coordinates " << points[pointI] << endl;
 279     Pout << "v00 = " << v00 << endl;
 280     Pout << "v10 = " << v10 << endl;
 281     Pout << "v11 = " << v11 << endl;
 282     Pout << "v01 = " << v01 << endl;
 283     # endif
 284
 285     forAll(facePoints, i)
 286     {
 287         forAll(facePoints[i], j)
 288         {
 289             if( facePoints[i][j] < 0 )
 290             {
 291                 # ifdef DEBUGLayer
 292                 Pout << "Determining u " << facePoints[0][0]
 293                      << " coordinates " << points[facePoints[0][0]] << endl;
 294                 Pout << "Other point " << facePoints[i][0]
 295                      << " coordinates " << points[facePoints[i][0]] << endl;
 296                 Pout << "Points at aplit edge "
 297                      << newVerticesForSplitEdge_[dir0Edges.second()] << endl;}
 298                 # endif
 299
 300                 const scalar u
 301                 (
 302                     Foam::mag(points[facePoints[i][0]] - v00) /
 303                     splitEdges_[dir0Edges.second()].mag(points)
 304                 );
 305
 306                 # ifdef DEBUGLayer
 307                 Pout << "Determining v " << facePoints[0][0]
 308                      << " coordinates " << points[facePoints[0][0]] << endl;
 309                 Pout << "Other point " << facePoints[0][j]
 310                      << " coordinates " << points[facePoints[0][j]] << endl;
 311                 Pout << "Points at aplit edge "
 312                      << newVerticesForSplitEdge_[dir1Edges.first()] << endl;}
 313                 # endif
 314
 315                 const scalar v
 316                 (
 317                     Foam::mag(points[facePoints[0][j]] - v00) /
 318                     splitEdges_[dir1Edges.first()].mag(points)
 319                 );
 320
 321                 # ifdef DEBUGLayer
 322                 Pout << "Generating point of face " << endl;
 323                 Pout << "u = " << u << endl;
 324                 Pout << "v = " << v << endl;}
 325                 # endif
 326
 327                 //- calculate the coordinates of the missing point via
 328                 //- transfinite interpolation
 329                 const point newP
 330                 (
 331                     (1.0 - u) * (1.0 - v) * v00 +
 332                     u * (1.0 - v) * v10 +
 333                     u * v * v11 +
 334                     (1.0 - u) * v * v01
 335                 );
 336
 337                 # ifdef DEBUGLayer
 338                 Pout << "Point coordinate " << newP << endl;
 339                 # endif
 340
 341                 //- add the vertex to the mesh
 342                 facePoints[i][j] = points.size();
 343                 points.append(newP);
 344             }
 345         }
 346     }
 347
 348     # ifdef DEBUGLayer
 349     Pout << "Face points after creating vertices " << facePoints << endl;
 350     # endif
 351
 352     //- Finally, create the faces
 353     for(label j=0;j<nLayersDir1;++j)
 354     {
 355         for(label i=0;i<nLayersDir0;++i)
 356         {
 357             //- create quad face
 358             DynList<label, 4> f;
 359
 360             f.append(facePoints[i][j]);
 361
 362             if( (i == (nLayersDir0 - 1)) && (j == 0) )
 363             {
 364                 # ifdef DEBUGLayer
 365                 Pout << "1. Adding additional points on edge " << endl;
 366                 # endif
 367
 368                 //- add additional points on edge
 369                 const label eLabel = dir0Edges.second();
 370                 const label size =
 371                     newVerticesForSplitEdge_.sizeOfRow(eLabel) - 1;
 372
 373                 for(label index=i+1;index<size;++index)
 374                     f.append(newVerticesForSplitEdge_(eLabel, index));
 375             }
 376
 377             f.append(facePoints[i+1][j]);
 378
 379             if(
 380                 (dir1 != -1) &&
 381                 (i == (nLayersDir0 - 1)) &&
 382                 (j == (nLayersDir1 - 1))
 383             )
 384             {
 385                 # ifdef DEBUGLayer
 386                 Pout << "2. Adding additional points on edge " << endl;
 387                 # endif
 388
 389                 //- add additional points on edge
 390                 const label eLabel = dir1Edges.second();
 391                 const label size =
 392                     newVerticesForSplitEdge_.sizeOfRow(eLabel) - 1;
 393
 394                 for(label index=j+1;index<size;++index)
 395                     f.append(newVerticesForSplitEdge_(eLabel, index));
 396             }
 397
 398             f.append(facePoints[i+1][j+1]);
 399
 400             if( (i == (nLayersDir0 - 1)) && (j == (nLayersDir1 - 1)) )
 401             {
 402                 # ifdef DEBUGLayer
 403                 Pout << "3. Adding additional points on edge " << endl;
 404                 # endif
 405
 406                 const label eLabel = dir0Edges.first();
 407                 const label size =
 408                     newVerticesForSplitEdge_.sizeOfRow(eLabel) - 2;
 409                 for(label index=size;index>i;--index)
 410                     f.append(newVerticesForSplitEdge_(eLabel, index));
 411             }
 412
 413             f.append(facePoints[i][j+1]);
 414
 415             if( (dir1 != -1) && (i == 0) && (j == (nLayersDir1 - 1)) )
 416             {
 417                 # ifdef DEBUGLayer
 418                 Pout << "4. Adding additional points on edge " << endl;
 419                 # endif
 420
 421                 const label eLabel = dir1Edges.first();
 422                 const label size =
 423                     newVerticesForSplitEdge_.sizeOfRow(eLabel) - 2;
 424                 for(label index=size;index>j;--index)
 425                     f.append(newVerticesForSplitEdge_(eLabel, index));
 426             }
 427
 428             newFaces.append(f);
 429         }
 430     }
 431
 432     # ifdef DEBUGLayer
 433     Pout << "Input face " << f << endl;
 434     Pout << "Decomposed faces are " << newFaces << endl;
 435     //if( (nLayersInDirection[0] > 1) && (nLayersInDirection[1] > 1) )
 436     //::exit(1);
 437     # endif
 438 }
 439
 440 void refineBoundaryLayers::sortFacePoints
 441 (
 442     const label faceI,
 443     DynList<DynList<label> >& facePoints,
 444     const label transpose
 445 ) const
 446 {
 447     const faceListPMG& faces = mesh_.faces();
 448     const face& f = faces[faceI];
 449
 450     # ifdef DEBUGLayer
 451     Pout << "Creating matrix of points on a split face " << faceI << endl;
 452     Pout << "Face comprises of points " << f << endl;
 453     Pout << "New faces from face " << facesFromFace_.sizeOfRow(faceI) << endl;
 454     # endif
 455
 456     label procStart = mesh_.faces().size();
 457     const PtrList<processorBoundaryPatch>& procBoundaries = mesh_.procBoundaries();
 458     if( Pstream::parRun() )
 459         procStart = procBoundaries[0].patchStart();
 460
 461     if(
 462         (faceI < procStart) ||
 463         procBoundaries[mesh_.faceIsInProcPatch(faceI)].owner()
 464     )
 465     {
 466         //- orientation of new faces is the same as the face itself
 467         //- start the procedure by finding the number of splits in
 468         //- both i and j direction
 469         label numSplitsI(1);
 470
 471         const label pos = f.which(newFaces_(facesFromFace_(faceI, 0), 0));
 472
 473         forAllRow(facesFromFace_, faceI, i)
 474         {
 475             const label nfI = facesFromFace_(faceI, i);
 476
 477             if( (numSplitsI == 1) && newFaces_.contains(nfI, f.nextLabel(pos)) )
 478             {
 479                 numSplitsI = i + 1;
 480                 break;
 481             }
 482         }
 483
 484         const label numSplitsJ = (facesFromFace_.sizeOfRow(faceI) / numSplitsI);
 485
 486         # ifdef DEBUGLayer
 487         Pout << "Pos " << pos << endl;
 488         Pout << "Num splits in direction 0 " << numSplitsI << endl;
 489         Pout << "Num splits in direction 1 " << numSplitsJ << endl;
 490         # endif
 491
 492         facePoints.setSize(numSplitsI+1);
 493         forAll(facePoints, i)
 494             facePoints[i].setSize(numSplitsJ+1);
 495
 496         //- start filling in the matrix
 497         forAllRow(facesFromFace_, faceI, fI)
 498         {
 499             const label nfI = facesFromFace_(faceI, fI);
 500
 501             const label i = fI % numSplitsI;
 502             const label j = fI / numSplitsI;
 503
 504             # ifdef DEBUGLayer
 505             Pout << "New face " << fI << " is " << newFaces_[nfI] << endl;
 506             Pout << " i = " << i << endl;
 507             Pout << " j = " << j << endl;
 508             # endif
 509
 510             if( newFaces_.sizeOfRow(nfI) == 4 )
 511             {
 512                 facePoints[i][j] = newFaces_(nfI, 0);
 513                 facePoints[i+1][j] = newFaces_(nfI, 1);
 514                 facePoints[i+1][j+1] = newFaces_(nfI, 2);
 515                 facePoints[i][j+1] = newFaces_(nfI, 3);
 516             }
 517             else
 518             {
 519                 if( j == 0 )
 520                 {
 521                     forAllRow(newFaces_, nfI, pI)
 522                         if( f.which(newFaces_(nfI, pI)) >= 0 )
 523                         {
 524                             facePoints[i+1][0] = newFaces_(nfI, pI);
 525                             break;
 526                         }
 527                 }
 528                 else if( i == 0 )
 529                 {
 530                     forAllRow(newFaces_, nfI, pI)
 531                         if( f.which(newFaces_(nfI, pI)) >= 0 )
 532                         {
 533                             facePoints[0][j+1] = newFaces_(nfI, pI);
 534                             break;
 535                         }
 536                 }
 537                 else
 538                 {
 539                     forAllRow(newFaces_, nfI, pI)
 540                         if( f.which(newFaces_(nfI, pI)) >= 0 )
 541                         {
 542                             facePoints[i+1][j+1] = newFaces_(nfI, pI);
 543                             break;
 544                         }
 545                 }
 546             }
 547         }
 548
 549         # ifdef DEBUGLayer
 550         Pout << "Generated matrix of points on face " << faceI
 551              << " is " << facePoints << endl;
 552         # endif
 553     }
 554     else
 555     {
 556         //- this situation exists on inter-processor boundaries
 557         //- on neighbour processor. i and j coordinates are reversed
 558         label numSplitsJ(1);
 559
 560         const label pos = f.which(newFaces_(facesFromFace_(faceI, 0), 0));
 561
 562         forAllRow(facesFromFace_, faceI, j)
 563         {
 564             const label nfI = facesFromFace_(faceI, j);
 565
 566             if( (numSplitsJ == 1) && newFaces_.contains(nfI, f.prevLabel(pos)) )
 567             {
 568                 numSplitsJ = j + 1;
 569                 break;
 570             }
 571         }
 572
 573         const label numSplitsI = (facesFromFace_.sizeOfRow(faceI) / numSplitsJ);
 574
 575         # ifdef DEBUGLayer
 576         Pout << "2. Face comprises of points " << f << endl;
 577         Pout << "2. Num splits in direction 0 " << numSplitsI << endl;
 578         Pout << "2. Num splits in direction 1 " << numSplitsJ << endl;
 579         # endif
 580
 581         facePoints.setSize(numSplitsI+1);
 582         forAll(facePoints, i)
 583             facePoints[i].setSize(numSplitsJ+1);
 584
 585         //- start filling in the matrix
 586         forAllRow(facesFromFace_, faceI, fI)
 587         {
 588             const label nfI = facesFromFace_(faceI, fI);
 589
 590             const label i = fI / numSplitsJ;
 591             const label j = fI % numSplitsJ;
 592
 593             # ifdef DEBUGLayer
 594             Pout << "2. New face " << fI << " is " << newFaces_[nfI] << endl;
 595             Pout << "2. i = " << i << endl;
 596             Pout << "2. j = " << j << endl;
 597             # endif
 598
 599             if( newFaces_.sizeOfRow(nfI) == 4 )
 600             {
 601                 facePoints[i][j] = newFaces_(nfI, 0);
 602                 facePoints[i+1][j] = newFaces_(nfI, 1);
 603                 facePoints[i+1][j+1] = newFaces_(nfI, 2);
 604                 facePoints[i][j+1] = newFaces_(nfI, 3);
 605             }
 606             else
 607             {
 608                 if( i == 0 )
 609                 {
 610                     forAllRow(newFaces_, nfI, pI)
 611                         if( f.which(newFaces_(nfI, pI)) >= 0 )
 612                         {
 613                             facePoints[0][j+1] = newFaces_(nfI, pI);
 614                             break;
 615                         }
 616                 }
 617                 else if( j == 0 )
 618                 {
 619                     forAllRow(newFaces_, nfI, pI)
 620                         if( f.which(newFaces_(nfI, pI)) >= 0 )
 621                         {
 622                             facePoints[i+1][0] = newFaces_(nfI, pI);
 623                             break;
 624                         }
 625                 }
 626                 else
 627                 {
 628                     forAllRow(newFaces_, nfI, pI)
 629                         if( f.which(newFaces_(nfI, pI)) >= 0 )
 630                         {
 631                             facePoints[i+1][j+1] = newFaces_(nfI, pI);
 632                             break;
 633                         }
 634                 }
 635             }
 636         }
 637
 638         # ifdef DEBUGLayer
 639         Pout << "Generated matrix of points on processor face " << faceI
 640              << " is " << facePoints << endl;
 641         # endif
 642     }
 643
 644     if( transpose )
 645     {
 646         DynList<DynList<label> > transposedFacePoints;
 647         transposedFacePoints.setSize(facePoints[0].size());
 648         forAll(transposedFacePoints, j)
 649             transposedFacePoints[j].setSize(facePoints.size());
 650
 651         forAll(facePoints, i)
 652             forAll(facePoints[i], j)
 653                 transposedFacePoints[j][i] = facePoints[i][j];
 654
 655         facePoints = transposedFacePoints;
 656
 657         # ifdef DEBUGLayer
 658         Pout << "Transposed face points " << facePoints << endl;
 659         # endif
 660     }
 661 }
 662
 663 void refineBoundaryLayers::sortFaceFaces
 664 (
 665     const label faceI,
 666     DynList<DynList<label> >& faceFaces,
 667     const label transpose
 668 ) const
 669 {
 670     const faceListPMG& faces = mesh_.faces();
 671     const face& f = faces[faceI];
 672
 673     # ifdef DEBUGLayer
 674     Pout << "Creating matrix of faces on a split face " << faceI << endl;
 675     Pout << "Face comprises of points " << f << endl;
 676     # endif
 677
 678     label procStart = mesh_.faces().size();
 679     const PtrList<processorBoundaryPatch>& procBoundaries = mesh_.procBoundaries();
 680     if( Pstream::parRun() )
 681         procStart = procBoundaries[0].patchStart();
 682
 683     if(
 684         (faceI < procStart) ||
 685         procBoundaries[mesh_.faceIsInProcPatch(faceI)].owner()
 686     )
 687     {
 688         //- orientation of new faces is the same as the face itself
 689         //- start the procedure by finding the number of splits in
 690         //- both i and j direction
 691         label numSplitsI(1);
 692
 693         const label pos = f.which(newFaces_(facesFromFace_(faceI, 0), 0));
 694
 695         forAllRow(facesFromFace_, faceI, i)
 696         {
 697             const label nfI = facesFromFace_(faceI, i);
 698
 699             if( (numSplitsI == 1) && newFaces_.contains(nfI, f.nextLabel(pos)) )
 700             {
 701                 numSplitsI = i + 1;
 702                 break;
 703             }
 704         }
 705
 706         label numSplitsJ = (facesFromFace_.sizeOfRow(faceI) / numSplitsI);
 707
 708         # ifdef DEBUGLayer
 709         Pout << "3. Num splits in direction 0 " << numSplitsI << endl;
 710         Pout << "3. Num splits in direction 1 " << numSplitsJ << endl;
 711         # endif
 712
 713         faceFaces.setSize(numSplitsI);
 714         forAll(faceFaces, i)
 715             faceFaces[i].setSize(numSplitsJ);
 716
 717         //- start filling in the matrix
 718         forAllRow(facesFromFace_, faceI, fI)
 719         {
 720             const label nfI = facesFromFace_(faceI, fI);
 721
 722             const label i = fI % numSplitsI;
 723             const label j = fI / numSplitsI;
 724
 725             # ifdef DEBUGLayer
 726             Pout << "3. New face " << fI << " is " << newFaces_[nfI] << endl;
 727             Pout << "3. i = " << i << endl;
 728             Pout << "3. j = " << j << endl;
 729             # endif
 730
 731             faceFaces[i][j] = nfI;
 732         }
 733
 734         # ifdef DEBUGLayer
 735         Pout << "3. Generated matrix of points on face " << faceI
 736              << " is " << faceFaces << endl;
 737         # endif
 738     }
 739     else
 740     {
 741         //- this situation exists on inter-processor boundaries
 742         //- on neighbour processor. i and j coordinates are reversed
 743         label numSplitsJ(1);
 744
 745         const label pos = f.which(newFaces_(facesFromFace_(faceI, 0), 0));
 746
 747         forAllRow(facesFromFace_, faceI, j)
 748         {
 749             const label nfI = facesFromFace_(faceI, j);
 750
 751             if( (numSplitsJ == 1) && newFaces_.contains(nfI, f.prevLabel(pos)) )
 752             {
 753                 numSplitsJ = j + 1;
 754                 break;
 755             }
 756         }
 757
 758         const label numSplitsI = (facesFromFace_.sizeOfRow(faceI) / numSplitsJ);
 759
 760         # ifdef DEBUGLayer
 761         Pout << "4. Num splits in direction 0 " << numSplitsI << endl;
 762         Pout << "4. Num splits in direction 1 " << numSplitsJ << endl;
 763         # endif
 764
 765         faceFaces.setSize(numSplitsI);
 766         forAll(faceFaces, i)
 767             faceFaces[i].setSize(numSplitsJ);
 768
 769         //- start filling in the matrix
 770         forAllRow(facesFromFace_, faceI, fI)
 771         {
 772             const label nfI = facesFromFace_(faceI, fI);
 773
 774             const label i = fI / numSplitsJ;
 775             const label j = fI % numSplitsJ;
 776
 777             # ifdef DEBUGLayer
 778             Pout << "4. New face " << fI << " is " << newFaces_[nfI] << endl;
 779             Pout << "4. i = " << i << endl;
 780             Pout << "4. j = " << j << endl;
 781             # endif
 782
 783             faceFaces[i][j] = nfI;
 784         }
 785
 786         # ifdef DEBUGLayer
 787         Pout << "4. Generated matrix of faces on processor face " << faceI
 788              << " is " << faceFaces << endl;
 789         # endif
 790     }
 791
 792     if( transpose )
 793     {
 794         DynList<DynList<label> > transposedFaceFaces;
 795         transposedFaceFaces.setSize(faceFaces[0].size());
 796         forAll(transposedFaceFaces, j)
 797             transposedFaceFaces[j].setSize(faceFaces.size());
 798
 799         forAll(faceFaces, i)
 800             forAll(faceFaces[i], j)
 801                 transposedFaceFaces[j][i] = faceFaces[i][j];
 802
 803         faceFaces = transposedFaceFaces;
 804
 805         # ifdef DEBUGLayer
 806         Pout << "Transposed face faces " << faceFaces << endl;
 807         # endif
 808     }
 809 }
 810
 811 void refineBoundaryLayers::generateNewFaces()
 812 {
 813     //- generate new boundary and inter-processor faces
 814     const meshSurfaceEngine& mse = surfaceEngine();
 815     const faceList::subList& bFaces = mse.boundaryFaces();
 816     const labelList& facePatches = mse.boundaryFacePatches();
 817     const edgeList& edges = mse.edges();
 818     const labelList& bp = mse.bp();
 819     const VRWGraph& bfEdges = mse.faceEdges();
 820     const VRWGraph& bpEdges = mse.boundaryPointEdges();
 821     const VRWGraph& beFaces = mse.edgeFaces();
 822
 823     //- mesh data
 824     const label nInternalFaces = mesh_.nInternalFaces();
 825     const faceListPMG& faces = mesh_.faces();
 826
 827     //- container for faces
 828     facesFromFace_.setSize(faces.size());
 829     newFaces_.clear();
 830
 831     //- split internal faces
 832     for(label faceI=0;faceI<nInternalFaces;++faceI)
 833     {
 834         const face& f = faces[faceI];
 835
 836         //- only quad faces can be split
 837         if( f.size() != 4 )
 838         {
 839             facesFromFace_.append(faceI, newFaces_.size());
 840             newFaces_.appendList(f);
 841             continue;
 842         }
 843
 844         //- check if there exist an edge of the face at the boundary
 845         FixedList<label, 2> nRefinementInDirection(1);
 846
 847         forAll(f, eI)
 848         {
 849             const edge fe = f.faceEdge(eI);
 850
 851             const label bps = bp[fe.start()];
 852
 853             if( bps < 0 )
 854                 continue;
 855
 856             forAllRow(bpEdges, bps, bpsI)
 857             {
 858                 const label beI = bpEdges(bps, bpsI);
 859
 860                 if( edges[beI] == fe )
 861                 {
 862                     //- this edge is attached to the boundary
 863                     //- get the number of layers for neighbouring cells
 864                     const label nSplits0 = nLayersAtBndFace_[beFaces(beI, 0)];
 865                     const label nSplits1 = nLayersAtBndFace_[beFaces(beI, 1)];
 866
 867                     //- set the number of layers for the given direction
 868                     const label dir = eI % 2;
 869                     nRefinementInDirection[dir] =
 870                         Foam::max
 871                         (
 872                             nRefinementInDirection[dir],
 873                             Foam::max(nSplits0, nSplits1)
 874                         );
 875                 }
 876             }
 877         }
 878
 879         //- refine the face
 880         DynList<DynList<label, 4>, 128> newFacesForFace;
 881         refineFace(f, nRefinementInDirection, newFacesForFace);
 882
 883         //- store decomposed faces
 884         forAll(newFacesForFace, fI)
 885         {
 886             facesFromFace_.append(faceI, newFaces_.size());
 887             newFaces_.appendList(newFacesForFace[fI]);
 888         }
 889
 890         # ifdef DEBUGLayer
 891         Pout << "Internal face " << faceI << " with points " << f
 892              << " is refined " << endl;
 893         forAllRow(facesFromFace_, faceI, i)
 894             Pout << "New face " << i << " is "
 895                  << newFaces_[facesFromFace_(faceI, i)] << endl;
 896         DynList<DynList<label> > tralala;
 897         sortFacePoints(faceI, tralala);
 898         # endif
 899     }
 900
 901     //- refine boundary faces where needed
 902     //- it is required in locations where two or three layers intersect
 903     const label startingBoundaryFace = mesh_.boundaries()[0].patchStart();
 904     forAll(bFaces, bfI)
 905     {
 906         const face& bf = bFaces[bfI];
 907         const label faceI = startingBoundaryFace + bfI;
 908
 909         //- only quad faces can be split
 910         if( bf.size() != 4 )
 911         {
 912             facesFromFace_.append(faceI, newFaces_.size());
 913             newFaces_.appendList(bf);
 914             continue;
 915         }
 916
 917         //- check whether this face shall be refined and in which directions
 918         FixedList<label, 2> nRefinementInDirection(1);
 919
 920         forAll(bf, eI)
 921         {
 922             const label beI = bfEdges(bfI, eI);
 923
 924             //- get the neighbour face over the edge
 925             label neiFace = beFaces(beI, 0);
 926
 927             if( beFaces.sizeOfRow(beI) != 2 )
 928                 continue;
 929
 930             if( neiFace == bfI )
 931                 neiFace = beFaces(beI, 1);
 932
 933             //- faces cannot be in the same layer
 934             if(
 935                 layerAtPatch_[facePatches[neiFace]] ==
 936                 layerAtPatch_[facePatches[bfI]]
 937             )
 938                 continue;
 939
 940             //- set the refinement direction for this face
 941             nRefinementInDirection[eI%2] = nLayersAtBndFace_[neiFace];
 942         }
 943
 944         //- refine the face
 945         DynList<DynList<label, 4>, 128> newFacesForFace;
 946         refineFace(bf, nRefinementInDirection, newFacesForFace);
 947
 948         //- store the refined faces
 949         forAll(newFacesForFace, fI)
 950         {
 951             facesFromFace_.append(faceI, newFaces_.size());
 952             newFaces_.appendList(newFacesForFace[fI]);
 953         }
 954
 955         # ifdef DEBUGLayer
 956         Pout << "Boundary face " << faceI << " with points " << bf
 957              << " owner cell " << mesh_.owner()[faceI] << " is refined " << endl;
 958         forAllRow(facesFromFace_, faceI, i)
 959             Pout << "New face " << i << " is "
 960                  << newFaces_[facesFromFace_(faceI, i)] << endl;
 961         # endif
 962     }
 963
 964     if( Pstream::parRun() )
 965     {
 966         //- refine faces at interprocessor boundaries
 967         const PtrList<processorBoundaryPatch>& procBoundaries =
 968             mesh_.procBoundaries();
 969
 970         //- exchange information about the number of splits
 971         //- to other processors
 972         std::map<label, DynList<labelPair, 2> > localSplits;
 973         forAll(procBoundaries, patchI)
 974         {
 975             labelLongList sendData;
 976
 977             const label start = procBoundaries[patchI].patchStart();
 978             const label size = procBoundaries[patchI].patchSize();
 979
 980             for(label fI=0;fI<size;++fI)
 981             {
 982                 const label faceI = start + fI;
 983                 const face& f = faces[faceI];
 984
 985                 forAll(f, eI)
 986                 {
 987                     const edge fe = f.faceEdge(eI);
 988
 989                     const label bps = bp[fe.start()];
 990
 991                     if( bps < 0 )
 992                         continue;
 993
 994                     forAllRow(bpEdges, bps, bpeI)
 995                     {
 996                         const label beI = bpEdges(bps, bpeI);
 997
 998                         if( edges[beI] == fe )
 999                         {
1000                             //- this edge is attached to the boundary
1001                             //- get the number of layers for neighbouring cell
1002                             const label nSplits0 =
1003                                 nLayersAtBndFace_[beFaces(beI, 0)];
1004
1005                             //- add the data to the list for sending
1006                             const label dir = (eI % 2);
1007
1008                             # ifdef DEBUGLayer
1009                             Pout << "Face " << fI << " owner of proc patch "
1010                                  << procBoundaries[patchI].myProcNo()
1011                                  << " nei proc "
1012                                  << procBoundaries[patchI].neiProcNo()
1013                                  << " bnd face patch "
1014                                  << facePatches[beFaces(beI, 0)]
1015                                  << " direction " << dir
1016                                  << " nSplits " << nSplits0 << endl;
1017                             # endif
1018
1019                             //- add face label, direction
1020                             //- and the number of splits
1021                             sendData.append(fI);
1022                             sendData.append(dir);
1023                             sendData.append(nSplits0);
1024                             localSplits[faceI].append(labelPair(dir, nSplits0));
1025                         }
1026                     }
1027                 }
1028             }
1029
1030             OPstream toOtherProc
1031             (
1032                 Pstream::blocking,
1033                 procBoundaries[patchI].neiProcNo(),
1034                 sendData.byteSize()
1035             );
1036
1037             toOtherProc << sendData;
1038         }
1039
1040         //- receive data from other procesors
1041         forAll(procBoundaries, patchI)
1042         {
1043             //- get the data sent from the neighbour processor
1044             labelList receivedData;
1045
1046             IPstream fromOtherProc
1047             (
1048                 Pstream::blocking,
1049                 procBoundaries[patchI].neiProcNo()
1050             );
1051
1052             fromOtherProc >> receivedData;
1053
1054             const label start = procBoundaries[patchI].patchStart();
1055
1056             label counter(0);
1057             while( counter < receivedData.size() )
1058             {
1059                 const label fI = receivedData[counter++];
1060                 const label dir = ((receivedData[counter++] + 1) % 2);
1061                 const label nSplits = receivedData[counter++];
1062
1063                 DynList<labelPair, 2>& currentSplits = localSplits[start+fI];
1064                 forAll(currentSplits, i)
1065                 {
1066                     if( currentSplits[i].first() == dir )
1067                         currentSplits[i].second() =
1068                             Foam::max(currentSplits[i].second(), nSplits);
1069                 }
1070             }
1071         }
1072
1073         # ifdef DEBUGLayer
1074         returnReduce(1, sumOp<label>());
1075         Pout << "Starting splitting processor boundaries" << endl;
1076         # endif
1077
1078         //- perform splitting
1079         forAll(procBoundaries, patchI)
1080         {
1081             const label start = procBoundaries[patchI].patchStart();
1082             const label size = procBoundaries[patchI].patchSize();
1083
1084             for(label fI=0;fI<size;++fI)
1085             {
1086                 const label faceI = start + fI;
1087
1088                 std::map<label, DynList<labelPair, 2> >::const_iterator it =
1089                     localSplits.find(faceI);
1090
1091                 if( it == localSplits.end() )
1092                 {
1093                     //- this face is not split
1094                     facesFromFace_.append(faceI, newFaces_.size());
1095                     newFaces_.appendList(faces[faceI]);
1096                     continue;
1097                 }
1098
1099                 //- split the face and add the faces to the list
1100                 if( procBoundaries[patchI].owner() )
1101                 {
1102                     //- this processor owns this patch
1103                     FixedList<label, 2> nLayersInDirection(1);
1104                     const DynList<labelPair, 2>& dirSplits = it->second;
1105                     forAll(dirSplits, i)
1106                         nLayersInDirection[dirSplits[i].first()] =
1107                             dirSplits[i].second();
1108
1109                     # ifdef DEBUGLayer
1110                     Pout << "Face " << fI << " at owner processor "
1111                         << procBoundaries[patchI].myProcNo()
1112                         << " neighbour processor "
1113                         << procBoundaries[patchI].neiProcNo()
1114                         << " face " << faces[faceI] << " refinement direction "
1115                         << nLayersInDirection << endl;
1116                     # endif
1117
1118                     DynList<DynList<label, 4>, 128> facesFromFace;
1119                     refineFace(faces[faceI], nLayersInDirection, facesFromFace);
1120
1121                     //- add faces
1122                     forAll(facesFromFace, i)
1123                     {
1124                         facesFromFace_.append(faceI, newFaces_.size());
1125                         newFaces_.appendList(facesFromFace[i]);
1126                     }
1127                 }
1128                 else
1129                 {
1130                     //- reverse the face before splitting
1131                     FixedList<label, 2> nLayersInDirection(1);
1132                     const DynList<labelPair, 2>& dirSplits = it->second;
1133                     forAll(dirSplits, i)
1134                         nLayersInDirection[(dirSplits[i].first()+1)%2] =
1135                             dirSplits[i].second();
1136
1137                     const face rFace = faces[faceI].reverseFace();
1138
1139                     # ifdef DEBUGLayer
1140                     Pout << "Face " << fI << " at owner processor "
1141                         << procBoundaries[patchI].myProcNo()
1142                         << " neighbour processor "
1143                         << procBoundaries[patchI].neiProcNo()
1144                         << " face " << rFace << " refinement direction "
1145                         << nLayersInDirection << endl;
1146                     # endif
1147
1148                     DynList<DynList<label, 4>, 128> facesFromFace;
1149                     refineFace(rFace, nLayersInDirection, facesFromFace);
1150
1151                     forAll(facesFromFace, i)
1152                     {
1153                         const DynList<label, 4>& df = facesFromFace[i];
1154                         DynList<label, 4> rFace = help::reverseFace(df);
1155
1156                         facesFromFace_.append(faceI, newFaces_.size());
1157                         newFaces_.appendList(rFace);
1158                     }
1159                 }
1160             }
1161         }
1162
1163         # ifdef DEBUGLayer
1164         returnReduce(1, sumOp<label>());
1165         for(label procI=0;procI<Pstream::nProcs();++procI)
1166         {
1167             if( procI == Pstream::myProcNo() )
1168             {
1169                 forAll(procBoundaries, patchI)
1170                 {
1171                     const label start = procBoundaries[patchI].patchStart();
1172                     const label size = procBoundaries[patchI].patchSize();
1173
1174                     for(label fI=0;fI<size;++fI)
1175                     {
1176                         const label faceI = start + fI;
1177                         const face& f = faces[faceI];
1178                         Pout << "Face " << fI << " in patch "
1179                              << procBoundaries[patchI].patchName()
1180                              << " has nodes " << f
1181                              << " local splits " << localSplits[faceI]
1182                              << " new faces from face " << facesFromFace_[faceI]
1183                              << endl;
1184
1185                         Pout << " Face points ";
1186                         forAll(f, pI)
1187                             Pout << mesh_.points()[f[pI]] << " ";
1188                         Pout << endl;
1189
1190                         forAllRow(facesFromFace_, faceI, ffI)
1191                         {
1192                             const label nfI = facesFromFace_(faceI, ffI);
1193                             Pout << "New face " << ffI << " with label " << nfI
1194                                  << " consists of points ";
1195                             forAllRow(newFaces_, nfI, pI)
1196                                 Pout << mesh_.points()[newFaces_(nfI, pI)]
1197                                      << " ";
1198                             Pout << endl;
1199                         }
1200                     }
1201                 }
1202             }
1203
1204             returnReduce(1, sumOp<label>());
1205         }
1206
1207         returnReduce(1, sumOp<label>());
1208         //::exit(1);
1209         # endif
1210     }
1211
1212     # ifdef DEBUGLayer
1213     returnReduce(1, sumOp<label>());
1214
1215     for(label procI=0;procI<Pstream::nProcs();++procI)
1216     {
1217         if( procI == Pstream::myProcNo() )
1218         {
1219             Pout << "facesFromFace_ " << facesFromFace_ << endl;
1220             Pout << "newFaces_ " << newFaces_ << endl;
1221         }
1222
1223         returnReduce(1, sumOp<label>());
1224     }
1225     # endif
1226
1227     Info << "Finished refining boundary-layer faces " << endl;
1228 }
1229
1230 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
1231
1232 } // End namespace Foam
1233
1234 // ************************************************************************* //
1235