src/mesh/cfMesh/utilities/boundaryLayers/refineBoundaryLayers/refineBoundaryLayersFunctions.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | cfMesh: A library for mesh generation
   4    \\    /   O peration     |
   5     \\  /    A nd           | Author: Franjo Juretic (franjo.juretic@c-fields.com)
   6      \\/     M anipulation  | Copyright (C) Creative Fields, Ltd.
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of cfMesh.
  10
  11     cfMesh is free software; you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by the
  13     Free Software Foundation; either version 3 of the License, or (at your
  14     option) any later version.
  15
  16     cfMesh is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with cfMesh.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 Description
  25
  26 \*---------------------------------------------------------------------------*/
  27
  28 #include "refineBoundaryLayers.H"
  29 #include "meshSurfaceEngine.H"
  30 #include "helperFunctions.H"
  31 #include "polyMeshGenAddressing.H"
  32 #include "polyMeshGen2DEngine.H"
  33 #include "VRWGraphList.H"
  34 #include "meshSurfacePartitioner.H"
  35 #include "detectBoundaryLayers.H"
  36
  37 #include "labelledPair.H"
  38 #include "labelledScalar.H"
  39
  40 # ifdef USE_OMP
  41 #include <omp.h>
  42 # endif
  43
  44 //#define DEBUGLayer
  45
  46 # ifdef DEBUGLayer
  47 #include "OFstream.H"
  48 # endif
  49
  50 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  51
  52 namespace Foam
  53 {
  54
  55 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  56
  57 bool refineBoundaryLayers::analyseLayers()
  58 {
  59     const meshSurfaceEngine& mse = surfaceEngine();
  60     const faceList::subList& bFaces = mse.boundaryFaces();
  61     const labelList& facePatch = mse.boundaryFacePatches();
  62
  63     meshSurfacePartitioner mPart(mse);
  64     detectBoundaryLayers dbl(mPart, is2DMesh_);
  65
  66     const label nGroups = dbl.nDistinctLayers();
  67     const labelList& faceInLayer = dbl.faceInLayer();
  68
  69     //- get the hair edges
  70     splitEdges_ = dbl.hairEdges();
  71
  72     # ifdef DEBUGLayer
  73     OFstream file("hairEdges.vtk");
  74
  75     //- write the header
  76     file << "# vtk DataFile Version 3.0\n";
  77     file << "vtk output\n";
  78     file << "ASCII\n";
  79     file << "DATASET POLYDATA\n";
  80
  81     //- write points
  82     file << "POINTS " << 2*splitEdges_.size() << " float\n";
  83     forAll(splitEdges_, seI)
  84     {
  85         const point& p = mse.mesh().points()[splitEdges_[seI].start()];
  86
  87         file << p.x() << ' ' << p.y() << ' ' << p.z() << nl;
  88
  89         const point op = mse.mesh().points()[splitEdges_[seI].end()];
  90
  91         file << op.x() << ' ' << op.y() << ' ' << op.z() << nl;
  92     }
  93
  94     //- write lines
  95     file << "\nLINES " << splitEdges_.size()
  96          << " " << 3*splitEdges_.size() << nl;
  97     forAll(splitEdges_, eI)
  98     {
  99         file << 2 << " " << 2*eI << " " << (2*eI+1) << nl;
 100     }
 101
 102     file << "\n";
 103     # endif
 104
 105     //- create point to split edges addressing
 106     splitEdgesAtPoint_.reverseAddressing(splitEdges_);
 107
 108     //- check if the layer is valid
 109     bool validLayer(true);
 110     # ifdef USE_OMP
 111     # pragma omp parallel for schedule(dynamic, 40)
 112     # endif
 113     forAll(faceInLayer, bfI)
 114     {
 115         if( faceInLayer[bfI] < 0 )
 116             continue;
 117
 118         const face& bf = bFaces[bfI];
 119
 120         forAll(bf, pI)
 121             if( splitEdgesAtPoint_.sizeOfRow(bf[pI]) == 0 )
 122                 validLayer = false;
 123     }
 124
 125     # ifdef DEBUGLayer
 126     Info << "Number of independent layers in the mesh is " << nGroups << endl;
 127     Info << "Is valid layer " << validLayer << endl;
 128     # endif
 129
 130     const PtrList<boundaryPatch>& boundaries = mesh_.boundaries();
 131
 132     //- create patch name to index addressing
 133     std::map<word, label> patchNameToIndex;
 134     forAll(boundaries, patchI)
 135         patchNameToIndex[boundaries[patchI].patchName()] = patchI;
 136
 137     //- check layer labels over a patch
 138     layerAtPatch_.setSize(boundaries.size());
 139     forAll(layerAtPatch_, i)
 140         layerAtPatch_[i].clear();
 141     List<DynList<label> > groupsAtPatch(boundaries.size());
 142     forAll(faceInLayer, bfI)
 143         groupsAtPatch[facePatch[bfI]].appendIfNotIn(faceInLayer[bfI]);
 144
 145     //- set the information which patches have an extruded layer
 146     forAll(groupsAtPatch, patchI)
 147     {
 148         const DynList<label>& layers = groupsAtPatch[patchI];
 149
 150         forAll(layers, i)
 151         {
 152             if( layers[i] < 0 )
 153             {
 154                 layerAtPatch_[patchI].clear();
 155                 break;
 156             }
 157             else
 158             {
 159                 layerAtPatch_[patchI].append(layers[i]);
 160             }
 161         }
 162     }
 163
 164     # ifdef DEBUGLayer
 165     Info << "Layer at patch " << layerAtPatch_ << endl;
 166     # endif
 167
 168     //- set the information which patches are a single boundary layer face
 169     patchesInLayer_.setSize(nGroups);
 170     forAll(layerAtPatch_, patchI)
 171     {
 172         const DynList<label>& layers = layerAtPatch_[patchI];
 173
 174         forAll(layers, i)
 175             patchesInLayer_[layers[i]].append
 176             (
 177                 boundaries[patchI].patchName()
 178             );
 179     }
 180
 181     # ifdef DEBUGLayer
 182     Info << "Patches in layer " << patchesInLayer_ << endl;
 183     # endif
 184
 185     //- set the number of boundary layers for each patch
 186     labelList nLayersAtPatch(layerAtPatch_.size(), -1);
 187     boolList protectedValue(layerAtPatch_.size(), false);
 188
 189     forAll(patchesInLayer_, layerI)
 190     {
 191         const DynList<word>& layerPatches = patchesInLayer_[layerI];
 192
 193         label maxNumLayers(1);
 194         bool hasLocalValue(false);
 195
 196         //- find the maximum requested number of layers over the layer
 197         forAll(layerPatches, lpI)
 198         {
 199             const word pName = layerPatches[lpI];
 200
 201             std::map<word, label>::const_iterator it =
 202                 numLayersForPatch_.find(pName);
 203
 204             if( it != numLayersForPatch_.end() )
 205             {
 206                 //- check if the layer is interrupted at this patch
 207                 if(
 208                     discontinuousLayersForPatch_.find(pName) !=
 209                     discontinuousLayersForPatch_.end()
 210                 )
 211                 {
 212                     //- set the number of layers and lock this location
 213                     nLayersAtPatch[patchNameToIndex[pName]] = it->second;
 214                     protectedValue[patchNameToIndex[pName]] = true;
 215                     hasLocalValue = true;
 216                 }
 217                 else
 218                 {
 219                     //- take the maximum number of layers
 220                     maxNumLayers = Foam::max(maxNumLayers, it->second);
 221                     hasLocalValue = true;
 222                 }
 223             }
 224         }
 225
 226         //- apply the global value if no local values exist
 227         if( !hasLocalValue )
 228             maxNumLayers = globalNumLayers_;
 229
 230         //- apply the maximum number of ayer of all unprotected patches
 231         forAll(layerPatches, lpI)
 232         {
 233             const label ptchI = patchNameToIndex[layerPatches[lpI]];
 234
 235             if( !protectedValue[ptchI] )
 236                 nLayersAtPatch[ptchI] = maxNumLayers;
 237         }
 238     }
 239
 240     if( is2DMesh_ )
 241     {
 242         polyMeshGen2DEngine mesh2DEngine(mesh_);
 243         const boolList& zMinPoint = mesh2DEngine.zMinPoints();
 244         const boolList& zMaxPoint = mesh2DEngine.zMaxPoints();
 245
 246         const faceList::subList& bFaces = mse.boundaryFaces();
 247
 248         boolList allZMax(mesh_.boundaries().size(), true);
 249         boolList allZMin(mesh_.boundaries().size(), true);
 250
 251         # ifdef USE_OMP
 252         # pragma omp parallel for schedule(dynamic, 50)
 253         # endif
 254         forAll(bFaces, bfI)
 255         {
 256             const face& bf = bFaces[bfI];
 257
 258             forAll(bf, pI)
 259             {
 260                 if( !zMinPoint[bf[pI]] )
 261                     allZMin[facePatch[bfI]] = false;
 262                 if( !zMaxPoint[bf[pI]] )
 263                     allZMax[facePatch[bfI]] = false;
 264             }
 265         }
 266
 267         //- mark empty patches as already used
 268         forAll(allZMin, patchI)
 269         {
 270             if( allZMin[patchI] ^ allZMax[patchI] )
 271             {
 272                 nLayersAtPatch[patchI] = -1;
 273                 layerAtPatch_[patchI].clear();
 274             }
 275         }
 276     }
 277
 278     //- perform reduction over all processors
 279     reduce(nLayersAtPatch, maxOp<labelList>());
 280
 281     # ifdef DEBUGLayer
 282     Pout << "nLayersAtPatch " << nLayersAtPatch << endl;
 283     # endif
 284
 285     //- set the number of boundary layers which shall be generated above
 286     //- each boundary face
 287     nLayersAtBndFace_.setSize(facePatch.size());
 288     nLayersAtBndFace_ = globalNumLayers_;
 289
 290     # ifdef USE_OMP
 291     # pragma omp parallel for schedule(dynamic, 50)
 292     # endif
 293     forAll(nLayersAtBndFace_, bfI)
 294     {
 295         const label patchI = facePatch[bfI];
 296
 297         if( nLayersAtPatch[patchI] < 0 )
 298         {
 299             nLayersAtBndFace_[bfI] = 1;
 300         }
 301         else
 302         {
 303             nLayersAtBndFace_[bfI] = nLayersAtPatch[patchI];
 304
 305             if( specialMode_ )
 306             {
 307                 ++nLayersAtBndFace_[bfI];
 308             }
 309         }
 310     }
 311
 312     # ifdef DEBUGLayer
 313     forAll(nLayersAtBndFace_, bfI)
 314     Pout << "Boundary face " << bfI << " in patch "
 315         << facePatch[bfI] << " num layers " << nLayersAtBndFace_[bfI] << endl;
 316     //::exit(1);
 317     # endif
 318
 319     return validLayer;
 320 }
 321
 322 void refineBoundaryLayers::generateNewVertices()
 323 {
 324     const PtrList<boundaryPatch>& boundaries = mesh_.boundaries();
 325     pointFieldPMG& points = mesh_.points();
 326
 327     const meshSurfaceEngine& mse = surfaceEngine();
 328     const faceList::subList& bFaces = mse.boundaryFaces();
 329     const VRWGraph& pointFaces = mse.pointFaces();
 330     const labelList& facePatch = mse.boundaryFacePatches();
 331     const labelList& bp = mse.bp();
 332
 333     //- allocate the data from storing parameters applying to a split edge
 334     LongList<scalar> firstLayerThickness(splitEdges_.size());
 335     LongList<scalar> thicknessRatio(splitEdges_.size());
 336     labelLongList nNodesAtEdge(splitEdges_.size());
 337     labelLongList nLayersAtEdge(splitEdges_.size());
 338
 339     //- count the number of vertices for each split edge
 340     # ifdef USE_OMP
 341     const label nThreads = 3 * omp_get_num_procs();
 342     # else
 343     const label nThreads = 1;
 344     # endif
 345
 346     # ifdef USE_OMP
 347     # pragma omp parallel num_threads(nThreads)
 348     # endif
 349     {
 350         //- start counting vertices at each thread
 351         # ifdef USE_OMP
 352         # pragma omp for schedule(static, 1)
 353         # endif
 354         forAll(splitEdges_, seI)
 355         {
 356             const edge& e = splitEdges_[seI];
 357
 358             //- get the requested number of boundary layers
 359             label nLayers(1);
 360             scalar ratio(globalThicknessRatio_);
 361             scalar thickness(globalMaxThicknessFirstLayer_);
 362             bool overridenThickness(false);
 363
 364             const label bpI = bp[e.start()];
 365
 366             forAllRow(pointFaces, bpI, pfI)
 367             {
 368                 const label bfI = pointFaces(bpI, pfI);
 369                 const label pos = help::positionOfEdgeInFace(e, bFaces[bfI]);
 370                 if( pos >= 0 )
 371                     continue;
 372
 373                 const word& patchName =
 374                     boundaries[facePatch[bfI]].patchName();
 375
 376                 //- overrride the global value with the maximum number of layers
 377                 //- at this edge
 378                 nLayers = Foam::max(nLayers, nLayersAtBndFace_[bfI]);
 379
 380                 //- override with the maximum ratio
 381                 const std::map<word, scalar>::const_iterator rIt =
 382                     thicknessRatioForPatch_.find(patchName);
 383                 if( rIt != thicknessRatioForPatch_.end() )
 384                 {
 385                     ratio = rIt->second;
 386                 }
 387
 388                 //- override with the minimum thickness set for this edge
 389                 const std::map<word, scalar>::const_iterator tIt =
 390                     maxThicknessForPatch_.find(patchName);
 391                 if( tIt != maxThicknessForPatch_.end() )
 392                 {
 393                     if( overridenThickness )
 394                     {
 395                         thickness = Foam::min(thickness, tIt->second);
 396                     }
 397                     else
 398                     {
 399                         thickness = tIt->second;
 400                         overridenThickness = true;
 401                     }
 402                 }
 403             }
 404
 405             //- store the information
 406             firstLayerThickness[seI] = thickness;
 407             thicknessRatio[seI] = ratio;
 408             nLayersAtEdge[seI] = nLayers;
 409
 410             if( !specialMode_ )
 411             {
 412                 nNodesAtEdge[seI] = nLayers + 1;
 413             }
 414             else
 415             {
 416                 nNodesAtEdge[seI] = 3;
 417             }
 418         }
 419     }
 420
 421     if( Pstream::parRun() )
 422     {
 423         //- transfer the information over all processor for edges
 424         //- at inter-processor boundaries
 425         const labelLongList& globalEdgeLabel =
 426             mesh_.addressingData().globalEdgeLabel();
 427         const VRWGraph& edgeAtProcs = mesh_.addressingData().edgeAtProcs();
 428         const Map<label>& globalToLocal =
 429             mesh_.addressingData().globalToLocalEdgeAddressing();
 430         const DynList<label>& neiProcs = mesh_.addressingData().edgeNeiProcs();
 431         const edgeList& edges = mesh_.addressingData().edges();
 432         const VRWGraph& pointEdges = mesh_.addressingData().pointEdges();
 433
 434         //- exchange point number of layers
 435         std::map<label, LongList<labelPair> > exchangeNumLayers;
 436         std::map<label, LongList<labelPair> > exchangeNumNodesAtEdge;
 437         std::map<label, LongList<labelledScalar> > exchangeThickness;
 438         std::map<label, LongList<labelledScalar> > exchangeRatio;
 439         forAll(neiProcs, i)
 440         {
 441             exchangeNumNodesAtEdge.insert
 442             (
 443                 std::make_pair(neiProcs[i], LongList<labelPair>())
 444             );
 445             exchangeNumLayers.insert
 446             (
 447                 std::make_pair(neiProcs[i], LongList<labelPair>())
 448             );
 449             exchangeThickness.insert
 450             (
 451                 std::make_pair(neiProcs[i], LongList<labelledScalar>())
 452             );
 453             exchangeRatio.insert
 454             (
 455                 std::make_pair(neiProcs[i], LongList<labelledScalar>())
 456             );
 457         }
 458
 459         //- exchange the number of layers
 460         forAll(splitEdges_, seI)
 461         {
 462             const edge& se = splitEdges_[seI];
 463
 464             const label s = se.start();
 465             label edgeI(-1);
 466             forAllRow(pointEdges, s, peI)
 467             {
 468                 const label eI = pointEdges(s, peI);
 469
 470                 if( edges[eI] == se )
 471                 {
 472                     edgeI = eI;
 473                     break;
 474                 }
 475             }
 476
 477             const label geI = globalEdgeLabel[edgeI];
 478
 479             if( globalToLocal.found(geI) )
 480             {
 481                 forAllRow(edgeAtProcs, edgeI, i)
 482                 {
 483                     const label neiProc = edgeAtProcs(edgeI, i);
 484
 485                     if( neiProc == Pstream::myProcNo() )
 486                         continue;
 487
 488                     exchangeNumNodesAtEdge[neiProc].append
 489                     (
 490                         labelPair(geI, nNodesAtEdge[seI])
 491                     );
 492                     exchangeNumLayers[neiProc].append
 493                     (
 494                         labelPair(geI, nLayersAtEdge[seI])
 495                     );
 496                     exchangeThickness[neiProc].append
 497                     (
 498                         labelledScalar(geI, firstLayerThickness[seI])
 499                     );
 500                     exchangeRatio[neiProc].append
 501                     (
 502                         labelledScalar(geI, thicknessRatio[seI])
 503                     );
 504                 }
 505             }
 506         }
 507
 508         //- exchange number of nodes at split edge
 509         LongList<labelPair> receivedNumLayers;
 510         help::exchangeMap(exchangeNumNodesAtEdge, receivedNumLayers);
 511
 512         forAll(receivedNumLayers, i)
 513         {
 514             const labelPair& lp = receivedNumLayers[i];
 515             const label eI = globalToLocal[lp.first()];
 516             const edge& e = edges[eI];
 517             label seI(-1);
 518             forAllRow(splitEdgesAtPoint_, e.start(), i)
 519             {
 520                 const label seJ = splitEdgesAtPoint_(e.start(), i);
 521                 if( splitEdges_[seJ] == e )
 522                 {
 523                     seI = seJ;
 524                     break;
 525                 }
 526             }
 527             nNodesAtEdge[seI] = std::max(nNodesAtEdge[seI], lp.second());
 528         }
 529
 530         //- exchange number of layers
 531         receivedNumLayers.clear();
 532         help::exchangeMap(exchangeNumLayers, receivedNumLayers);
 533
 534         forAll(receivedNumLayers, i)
 535         {
 536             const labelPair& lp = receivedNumLayers[i];
 537             const label eI = globalToLocal[lp.first()];
 538             const edge& e = edges[eI];
 539             label seI(-1);
 540             forAllRow(splitEdgesAtPoint_, e.start(), i)
 541             {
 542                 const label seJ = splitEdgesAtPoint_(e.start(), i);
 543                 if( splitEdges_[seJ] == e )
 544                 {
 545                     seI = seJ;
 546                     break;
 547                 }
 548             }
 549             nLayersAtEdge[seI] = std::max(nLayersAtEdge[seI], lp.second());
 550         }
 551
 552         //- exchange thickness ratio
 553         LongList<labelledScalar> receivedScalar;
 554         help::exchangeMap(exchangeRatio, receivedScalar);
 555
 556         forAll(receivedScalar, i)
 557         {
 558             const labelledScalar& ls = receivedScalar[i];
 559             const label eI = globalToLocal[ls.scalarLabel()];
 560             const edge& e = edges[eI];
 561             label seI(-1);
 562             forAllRow(splitEdgesAtPoint_, e.start(), i)
 563             {
 564                 const label seJ = splitEdgesAtPoint_(e.start(), i);
 565                 if( splitEdges_[seJ] == e )
 566                 {
 567                     seI = seJ;
 568                     break;
 569                 }
 570             }
 571             thicknessRatio[seI] = std::max(thicknessRatio[seI], ls.value());
 572         }
 573
 574         //- exchange maximum thickness of the first layer
 575         receivedScalar.clear();
 576         help::exchangeMap(exchangeThickness, receivedScalar);
 577
 578         forAll(receivedScalar, i)
 579         {
 580             const labelledScalar& ls = receivedScalar[i];
 581             const label eI = globalToLocal[ls.scalarLabel()];
 582             const edge& e = edges[eI];
 583             label seI(-1);
 584             forAllRow(splitEdgesAtPoint_, e.start(), i)
 585             {
 586                 const label seJ = splitEdgesAtPoint_(e.start(), i);
 587                 if( splitEdges_[seJ] == e )
 588                 {
 589                     seI = seJ;
 590                     break;
 591                 }
 592             }
 593             firstLayerThickness[seI] =
 594                 std::min(firstLayerThickness[seI], ls.value());
 595         }
 596     }
 597
 598     //- calculate the number of additional vertices which will be generated
 599     //- on edges of the mesh
 600     DynList<label> numPointsAtThread;
 601     numPointsAtThread.setSize(nThreads);
 602     numPointsAtThread = 0;
 603
 604     # ifdef USE_OMP
 605     # pragma omp parallel for num_threads(nThreads) schedule(static, 1)
 606     # endif
 607     forAll(nNodesAtEdge, seI)
 608     {
 609         # ifdef USE_OMP
 610         const label threadI = omp_get_thread_num();
 611         # else
 612         const label threadI(0);
 613         # endif
 614
 615         numPointsAtThread[threadI] += nNodesAtEdge[seI] - 2;
 616     }
 617
 618     //- allocate the space in a graph storing ids of points on a split edge
 619     newVerticesForSplitEdge_.setSizeAndRowSize(nNodesAtEdge);
 620
 621     //- calculate the number of points which will be generated
 622     //- on split edges
 623     label numPoints = points.size();
 624     forAll(numPointsAtThread, threadI)
 625     {
 626         const label nPts = numPointsAtThread[threadI];
 627         numPointsAtThread[threadI] = numPoints;
 628         numPoints += nPts;
 629     }
 630
 631     points.setSize(numPoints);
 632
 633     # ifdef DEBUGLayer
 634     Info << "Generating split vertices" << endl;
 635     # endif
 636
 637     //- generate vertices on split edges
 638     # ifdef USE_OMP
 639     # pragma omp parallel num_threads(nThreads)
 640     # endif
 641     {
 642         # ifdef USE_OMP
 643         const label threadI = omp_get_thread_num();
 644         # else
 645         const label threadI(0);
 646         # endif
 647
 648         label& nPoints = numPointsAtThread[threadI];
 649
 650         # ifdef USE_OMP
 651         # pragma omp for schedule(static, 1)
 652         # endif
 653         forAll(splitEdges_, seI)
 654         {
 655             const edge& e = splitEdges_[seI];
 656
 657             const vector v = e.vec(points);
 658             const scalar magv = mag(v);
 659
 660             const label nLayers = newVerticesForSplitEdge_.sizeOfRow(seI) - 1;
 661
 662             scalar firstThickness = magv / nLayersAtEdge[seI];
 663             if( thicknessRatio[seI] > (1. + SMALL) )
 664             {
 665                 firstThickness =
 666                     magv /
 667                     (
 668                         (1 - Foam::pow(thicknessRatio[seI], nLayersAtEdge[seI]))
 669                         / (1.0 - thicknessRatio[seI])
 670                     );
 671
 672                 # ifdef DEBUGLayer
 673                 Pout << "Thread " << threadI << endl;
 674                 Pout << "Generating vertices at split edge "
 675                      << " start point " << points[e.start()]
 676                      << " end point " << points[e.end()] << endl;
 677                 Pout << "Edge length " << magv << endl;
 678                 Pout << "Thickness of the first layer "
 679                      << firstThickness << endl;
 680                 # endif
 681             }
 682
 683             firstThickness =
 684                 Foam::min
 685                 (
 686                     Foam::max(firstLayerThickness[seI], SMALL),
 687                     firstThickness
 688                 );
 689
 690             if( specialMode_ )
 691             {
 692                 scalar t = firstThickness;
 693
 694                 for(label i=1;i<nLayersAtEdge[seI]-1;++i)
 695                     t += firstThickness * Foam::pow(thicknessRatio[seI], i);
 696
 697                 firstThickness = t;
 698             }
 699
 700             //- generate vertices for this edge
 701             newVerticesForSplitEdge_(seI, 0) = e.start();
 702
 703             scalar param = firstThickness;
 704             const vector vec = v / (magv + VSMALL);
 705
 706             for(label pI=1;pI<nLayers;++pI)
 707             {
 708                 //- generate the new vertex
 709                 const point newP = points[e.start()] + param * vec;
 710
 711                 # ifdef DEBUGLayer
 712                 Pout << "Split edge " << seI << " edge points " << e
 713                     << " start point " << points[e.start()]
 714                     << " end point " << points[e.end()]
 715                     << " param " << param
 716                     << " new point " << nPoints
 717                     << " has coordinates " << newP << endl;
 718                 # endif
 719
 720                 param += firstThickness * Foam::pow(thicknessRatio[seI], pI);
 721
 722                 newVerticesForSplitEdge_(seI, pI) = nPoints;
 723                 points[nPoints++] = newP;
 724             }
 725
 726             newVerticesForSplitEdge_(seI, nLayers) = e.end();
 727         }
 728     }
 729
 730     if( specialMode_ )
 731     {
 732         //- set the number of layers to 2
 733         forAll(nLayersAtBndFace_, bfI)
 734             if( nLayersAtBndFace_[bfI] > 1 )
 735                 nLayersAtBndFace_[bfI] = 2;
 736     }
 737
 738     # ifdef DEBUGLayer
 739     for(label procI=0;procI<Pstream::nProcs();++procI)
 740     {
 741         if( procI == Pstream::myProcNo() )
 742         {
 743             forAll(splitEdges_, seI)
 744             {
 745                 Pout << "\nSplit edge " << seI << " nodes " << splitEdges_[seI]
 746                     << " coordinates " << points[splitEdges_[seI][0]]
 747                     << " " << points[splitEdges_[seI][1]]
 748                     << " has new points "
 749                     << newVerticesForSplitEdge_[seI] << endl;
 750
 751                 forAllRow(newVerticesForSplitEdge_, seI, i)
 752                     Pout << "Point " << i << " on edge ha coordinates "
 753                          << points[newVerticesForSplitEdge_(seI, i)] << endl;
 754             }
 755         }
 756
 757         returnReduce(1, sumOp<label>());
 758     }
 759
 760     Info << "Finished generating vertices at split edges" << endl;
 761     //::exit(1);
 762     # endif
 763 }
 764
 765 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 766
 767 } // End namespace Foam
 768
 769 // ************************************************************************* //