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