src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of OpenFOAM.
  10
  11     OpenFOAM 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 2 of the License, or (at your
  14     option) any later version.
  15
  16     OpenFOAM 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 OpenFOAM; if not, write to the Free Software Foundation,
  23     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
  24
  25 \*---------------------------------------------------------------------------*/
  26
  27 #include "addPatchCellLayer.H"
  28 #include "polyMesh.H"
  29 #include "polyTopoChange.H"
  30 #include "meshTools.H"
  31 #include "mapPolyMesh.H"
  32 #include "syncTools.H"
  33 #include "polyAddPoint.H"
  34 #include "polyAddFace.H"
  35 #include "polyModifyFace.H"
  36 #include "polyAddCell.H"
  37 #include "wallPoint.H"
  38 #include "globalIndex.H"
  39
  40 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  41
  42 defineTypeNameAndDebug(Foam::addPatchCellLayer, 0);
  43
  44
  45 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  46
  47 // Calculate global faces per pp edge.
  48 Foam::labelListList Foam::addPatchCellLayer::calcGlobalEdgeFaces
  49 (
  50     const polyMesh& mesh,
  51     const globalIndex& globalFaces,
  52     const indirectPrimitivePatch& pp,
  53     const labelList& meshEdges
  54 )
  55 {
  56     //// Determine coupled edges just so we don't have to have storage
  57     //// for all non-coupled edges.
  58     //
  59     //PackedBoolList isCoupledEdge(mesh.nEdges());
  60     //
  61     //const polyBoundaryMesh& patches = mesh.boundaryMesh();
  62     //
  63     //forAll(patches, patchI)
  64     //{
  65     //    const polyPatch& pp = patches[patchI];
  66     //
  67     //    if (pp.coupled())
  68     //    {
  69     //        const labelList& meshEdges = pp.meshEdges();
  70     //
  71     //        forAll(meshEdges, i)
  72     //        {
  73     //            isCoupledEdge.set(meshEdges[i], 1);
  74     //        }
  75     //    }
  76     //}
  77
  78     // From mesh edge to global face labels. Sized only for pp edges.
  79     labelListList globalEdgeFaces(mesh.nEdges());
  80
  81     const labelListList& edgeFaces = pp.edgeFaces();
  82
  83     forAll(edgeFaces, edgeI)
  84     {
  85         label meshEdgeI = meshEdges[edgeI];
  86
  87         //if (isCoupledEdge.get(meshEdgeI) == 1)
  88         {
  89             const labelList& eFaces = edgeFaces[edgeI];
  90
  91             // Store face and processor as unique tag.
  92             labelList& globalEFaces = globalEdgeFaces[meshEdgeI];
  93             globalEFaces.setSize(eFaces.size());
  94             forAll(eFaces, i)
  95             {
  96                 globalEFaces[i] =
  97                     globalFaces.toGlobal(pp.addressing()[eFaces[i]]);
  98             }
  99         }
 100     }
 101
 102     // Synchronise across coupled edges.
 103     syncTools::syncEdgeList
 104     (
 105         mesh,
 106         globalEdgeFaces,
 107         uniqueEqOp(),
 108         labelList(),    // null value
 109         false           // no separation
 110     );
 111
 112     // Extract pp part
 113     return UIndirectList<labelList>(globalEdgeFaces, meshEdges);
 114 }
 115
 116
 117 Foam::label Foam::addPatchCellLayer::nbrFace
 118 (
 119     const labelListList& edgeFaces,
 120     const label edgeI,
 121     const label faceI
 122 )
 123 {
 124     const labelList& eFaces = edgeFaces[edgeI];
 125
 126     if (eFaces.size() == 2)
 127     {
 128         return (eFaces[0] != faceI ? eFaces[0] : eFaces[1]);
 129     }
 130     else
 131     {
 132         return -1;
 133     }
 134 }
 135
 136
 137 void Foam::addPatchCellLayer::addVertex
 138 (
 139     const label pointI,
 140     face& f,
 141     label& fp
 142 )
 143 {
 144     if (fp == 0)
 145     {
 146         f[fp++] = pointI;
 147     }
 148     else
 149     {
 150         if (f[fp-1] != pointI && f[0] != pointI)
 151         {
 152             f[fp++] = pointI;
 153         }
 154     }
 155
 156     // Check for duplicates.
 157     if (debug)
 158     {
 159         label n = 0;
 160         for (label i = 0; i < fp; i++)
 161         {
 162             if (f[i] == pointI)
 163             {
 164                 n++;
 165
 166                 if (n == 2)
 167                 {
 168                     f.setSize(fp);
 169                     FatalErrorIn
 170                     (
 171                         "addPatchCellLayer::addVertex(const label, face&"
 172                         ", label&)"
 173                     )   << "Point " << pointI << " already present in face "
 174                         << f << abort(FatalError);
 175                 }
 176             }
 177         }
 178     }
 179 }
 180
 181
 182 // Is edge to the same neighbour? (and needs extrusion and has not been
 183 // dealt with already)
 184 bool Foam::addPatchCellLayer::sameEdgeNeighbour
 185 (
 186     const indirectPrimitivePatch& pp,
 187     const labelListList& globalEdgeFaces,
 188     const boolList& doneEdge,
 189     const label thisGlobalFaceI,
 190     const label nbrGlobalFaceI,
 191     const label edgeI
 192 ) const
 193 {
 194     const edge& e = pp.edges()[edgeI];
 195
 196     return
 197         !doneEdge[edgeI]                            // not yet handled
 198      && (
 199             addedPoints_[e[0]].size()               // is extruded
 200          || addedPoints_[e[1]].size()
 201         )
 202      && (
 203             nbrFace(globalEdgeFaces, edgeI, thisGlobalFaceI)
 204          == nbrGlobalFaceI  // is to same neighbour
 205         );
 206 }
 207
 208
 209 // Collect consecutive string of edges that connects the same two
 210 // (possibly coupled) faces. Returns -1 if no unvisited edge can be found.
 211 // Otherwise returns start and end index in face.
 212 Foam::labelPair Foam::addPatchCellLayer::getEdgeString
 213 (
 214     const indirectPrimitivePatch& pp,
 215     const labelListList& globalEdgeFaces,
 216     const boolList& doneEdge,
 217     const label patchFaceI,
 218     const label globalFaceI
 219 ) const
 220 {
 221     const labelList& fEdges = pp.faceEdges()[patchFaceI];
 222
 223     label startFp = -1;
 224     label endFp = -1;
 225
 226     // Get edge that hasn't been done yet but needs extrusion
 227     forAll(fEdges, fp)
 228     {
 229         label edgeI = fEdges[fp];
 230         const edge& e = pp.edges()[edgeI];
 231
 232         if
 233         (
 234             !doneEdge[edgeI]
 235          && ( addedPoints_[e[0]].size() || addedPoints_[e[1]].size() )
 236         )
 237         {
 238             startFp = fp;
 239             break;
 240         }
 241     }
 242
 243     if (startFp != -1)
 244     {
 245         // We found an edge that needs extruding but hasn't been done yet.
 246         // Now find the face on the other side
 247         label nbrGlobalFaceI = nbrFace
 248         (
 249             globalEdgeFaces,
 250             fEdges[startFp],
 251             globalFaceI
 252         );
 253
 254         if (nbrGlobalFaceI == -1)
 255         {
 256             // Proper boundary edge. Only extrude single edge.
 257             endFp = startFp;
 258         }
 259         else
 260         {
 261             // Search back for edge
 262             // - which hasn't been handled yet
 263             // - with same neighbour
 264             // - that needs extrusion
 265             while(true)
 266             {
 267                 label prevFp = fEdges.rcIndex(startFp);
 268
 269                 if
 270                 (
 271                     !sameEdgeNeighbour
 272                     (
 273                         pp,
 274                         globalEdgeFaces,
 275                         doneEdge,
 276                         globalFaceI,
 277                         nbrGlobalFaceI,
 278                         fEdges[prevFp]
 279                     )
 280                 )
 281                 {
 282                     break;
 283                 }
 284                 startFp = prevFp;
 285             }
 286
 287             // Search forward for end of string
 288             endFp = startFp;
 289             while(true)
 290             {
 291                 label nextFp = fEdges.fcIndex(endFp);
 292
 293                 if
 294                 (
 295                     !sameEdgeNeighbour
 296                     (
 297                         pp,
 298                         globalEdgeFaces,
 299                         doneEdge,
 300                         globalFaceI,
 301                         nbrGlobalFaceI,
 302                         fEdges[nextFp]
 303                     )
 304                 )
 305                 {
 306                     break;
 307                 }
 308                 endFp = nextFp;
 309             }
 310         }
 311     }
 312
 313     return labelPair(startFp, endFp);
 314 }
 315
 316
 317 // Adds a side face i.e. extrudes a patch edge.
 318 Foam::label Foam::addPatchCellLayer::addSideFace
 319 (
 320     const indirectPrimitivePatch& pp,
 321     const labelList& patchID,           // prestored patch per pp face
 322     const labelListList& addedCells,    // per pp face the new extruded cell
 323     const face& newFace,
 324     const label ownFaceI,               // pp face that provides owner
 325     const label nbrFaceI,
 326     const label patchEdgeI,             // edge to add to
 327     const label meshEdgeI,              // corresponding mesh edge
 328     const label layerI,                 // layer
 329     const label numEdgeFaces,           // number of layers for edge
 330     const labelList& meshFaces,         // precalculated edgeFaces
 331     polyTopoChange& meshMod
 332 ) const
 333 {
 334     // Edge to 'inflate' from
 335     label inflateEdgeI = -1;
 336
 337     // Check mesh faces using edge
 338     forAll(meshFaces, i)
 339     {
 340         if (mesh_.isInternalFace(meshFaces[i]))
 341         {
 342             // meshEdge uses internal faces so ok to inflate from it
 343             inflateEdgeI = meshEdgeI;
 344             break;
 345         }
 346     }
 347
 348
 349     // Get my mesh face and its zone.
 350     label meshFaceI = pp.addressing()[ownFaceI];
 351     // Zone info comes from any side patch face. Otherwise -1 since we
 352     // don't know what to put it in - inherit from the extruded faces?
 353     label zoneI = -1;   //mesh_.faceZones().whichZone(meshFaceI);
 354     bool flip = false;
 355
 356     label addedFaceI = -1;
 357
 358     // Is patch edge external edge of indirectPrimitivePatch?
 359     if (nbrFaceI == -1)
 360     {
 361         // External edge so external face. Patch id is obtained from
 362         // any other patch connected to edge.
 363
 364         const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 365
 366         // Loop over all faces connected to edge to inflate and
 367         // see if any boundary face (but not meshFaceI)
 368         label otherPatchID = patchID[ownFaceI];
 369
 370         forAll(meshFaces, k)
 371         {
 372             label faceI = meshFaces[k];
 373
 374             if
 375             (
 376                 faceI != meshFaceI
 377              && !mesh_.isInternalFace(faceI)
 378             )
 379             {
 380                 otherPatchID = patches.whichPatch(faceI);
 381                 zoneI = mesh_.faceZones().whichZone(faceI);
 382                 if (zoneI != -1)
 383                 {
 384                     label index = mesh_.faceZones()[zoneI].whichFace(faceI);
 385                     flip = mesh_.faceZones()[zoneI].flipMap()[index];
 386                 }
 387                 break;
 388             }
 389         }
 390
 391         // Determine if different number of layer on owner and neighbour side
 392         // (relevant only for coupled faces). See section for internal edge
 393         // below.
 394
 395         label layerOwn;
 396
 397         if (addedCells[ownFaceI].size() < numEdgeFaces)
 398         {
 399             label offset = numEdgeFaces - addedCells[ownFaceI].size();
 400             if (layerI <= offset)
 401             {
 402                 layerOwn = 0;
 403             }
 404             else
 405             {
 406                 layerOwn = layerI - offset;
 407             }
 408         }
 409         else
 410         {
 411             layerOwn = layerI;
 412         }
 413
 414
 415         //Pout<< "Added boundary face:" << newFace
 416         //    << " own:" << addedCells[ownFaceI][layerOwn]
 417         //    << " patch:" << otherPatchID
 418         //    << endl;
 419
 420         addedFaceI = meshMod.setAction
 421         (
 422             polyAddFace
 423             (
 424                 newFace,                    // face
 425                 addedCells[ownFaceI][layerOwn],   // owner
 426                 -1,                         // neighbour
 427                 -1,                         // master point
 428                 inflateEdgeI,               // master edge
 429                 -1,                         // master face
 430                 false,                      // flux flip
 431                 otherPatchID,               // patch for face
 432                 zoneI,                      // zone for face
 433                 flip                        // face zone flip
 434             )
 435         );
 436     }
 437     else
 438     {
 439         // When adding side faces we need to modify neighbour and owners
 440         // in region where layer mesh is stopped. Determine which side
 441         // has max number of faces and make sure layers match closest to
 442         // original pp if there are different number of layers.
 443
 444         label layerNbr;
 445         label layerOwn;
 446
 447         if (addedCells[ownFaceI].size() > addedCells[nbrFaceI].size())
 448         {
 449             label offset =
 450                 addedCells[ownFaceI].size() - addedCells[nbrFaceI].size();
 451
 452             layerOwn = layerI;
 453
 454             if (layerI <= offset)
 455             {
 456                 layerNbr = 0;
 457             }
 458             else
 459             {
 460                 layerNbr = layerI - offset;
 461             }
 462         }
 463         else if (addedCells[nbrFaceI].size() > addedCells[ownFaceI].size())
 464         {
 465             label offset =
 466                 addedCells[nbrFaceI].size() - addedCells[ownFaceI].size();
 467
 468             layerNbr = layerI;
 469
 470             if (layerI <= offset)
 471             {
 472                 layerOwn = 0;
 473             }
 474             else
 475             {
 476                 layerOwn = layerI - offset;
 477             }
 478         }
 479         else
 480         {
 481             // Same number of layers on both sides.
 482             layerNbr = layerI;
 483             layerOwn = layerI;
 484         }
 485
 486         addedFaceI = meshMod.setAction
 487         (
 488             polyAddFace
 489             (
 490                 newFace,                    // face
 491                 addedCells[ownFaceI][layerOwn],   // owner
 492                 addedCells[nbrFaceI][layerNbr],   // neighbour
 493                 -1,                         // master point
 494                 inflateEdgeI,               // master edge
 495                 -1,                         // master face
 496                 false,                      // flux flip
 497                 -1,                         // patch for face
 498                 zoneI,                      // zone for face
 499                 flip                        // face zone flip
 500             )
 501         );
 502
 503        //Pout<< "Added internal face:" << newFace
 504         //    << " own:" << addedCells[ownFaceI][layerOwn]
 505         //    << " nei:" << addedCells[nbrFaceI][layerNbr]
 506         //    << endl;
 507     }
 508
 509     return addedFaceI;
 510 }
 511
 512
 513 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 514
 515 // Construct from mesh
 516 Foam::addPatchCellLayer::addPatchCellLayer(const polyMesh& mesh)
 517 :
 518     mesh_(mesh),
 519     addedPoints_(0),
 520     layerFaces_(0)
 521 {}
 522
 523
 524 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 525
 526 Foam::labelListList Foam::addPatchCellLayer::addedCells
 527 (
 528     const polyMesh& mesh,
 529     const labelListList& layerFaces
 530 )
 531 {
 532     labelListList layerCells(layerFaces.size());
 533
 534     forAll(layerFaces, patchFaceI)
 535     {
 536         const labelList& faceLabels = layerFaces[patchFaceI];
 537
 538         if (faceLabels.size())
 539         {
 540             labelList& added = layerCells[patchFaceI];
 541             added.setSize(faceLabels.size()-1);
 542
 543             for (label i = 0; i < faceLabels.size()-1; i++)
 544             {
 545                 added[i] = mesh.faceNeighbour()[faceLabels[i]];
 546             }
 547         }
 548     }
 549     return layerCells;
 550 }
 551
 552
 553 Foam::labelListList Foam::addPatchCellLayer::addedCells() const
 554 {
 555     return addedCells(mesh_, layerFaces_);
 556 }
 557
 558
 559 void Foam::addPatchCellLayer::setRefinement
 560 (
 561     const scalarField& expansionRatio,
 562     const indirectPrimitivePatch& pp,
 563     const labelList& nFaceLayers,
 564     const labelList& nPointLayers,
 565     const vectorField& firstLayerDisp,
 566     polyTopoChange& meshMod
 567 )
 568 {
 569     if (debug)
 570     {
 571         Pout<< "addPatchCellLayer::setRefinement : Adding up to "
 572             << max(nPointLayers)
 573             << " layers of cells to indirectPrimitivePatch with "
 574             << pp.nPoints() << " points" << endl;
 575     }
 576
 577     if
 578     (
 579         pp.nPoints() != firstLayerDisp.size()
 580      || pp.nPoints() != nPointLayers.size()
 581      || pp.size() != nFaceLayers.size()
 582     )
 583     {
 584         FatalErrorIn
 585         (
 586             "addPatchCellLayer::setRefinement"
 587             "(const scalar, const indirectPrimitivePatch&"
 588             ", const labelList&, const vectorField&, polyTopoChange&)"
 589         )   << "Size of new points is not same as number of points used by"
 590             << " the face subset" << endl
 591             << "  patch.nPoints:" << pp.nPoints()
 592             << "  displacement:" << firstLayerDisp.size()
 593             << "  nPointLayers:" << nPointLayers.size() << nl
 594             << " patch.nFaces:" << pp.size()
 595             << "  nFaceLayers:" << nFaceLayers.size()
 596             << abort(FatalError);
 597     }
 598
 599     forAll(nPointLayers, i)
 600     {
 601         if (nPointLayers[i] < 0)
 602         {
 603             FatalErrorIn
 604             (
 605                 "addPatchCellLayer::setRefinement"
 606                 "(const scalar, const indirectPrimitivePatch&"
 607                 ", const labelList&, const vectorField&, polyTopoChange&)"
 608             )   << "Illegal number of layers " << nPointLayers[i]
 609                 << " at patch point " << i << abort(FatalError);
 610         }
 611     }
 612     forAll(nFaceLayers, i)
 613     {
 614         if (nFaceLayers[i] < 0)
 615         {
 616             FatalErrorIn
 617             (
 618                 "addPatchCellLayer::setRefinement"
 619                 "(const scalar, const indirectPrimitivePatch&"
 620                 ", const labelList&, const vectorField&, polyTopoChange&)"
 621             )   << "Illegal number of layers " << nFaceLayers[i]
 622                 << " at patch face " << i << abort(FatalError);
 623         }
 624     }
 625
 626     const labelList& meshPoints = pp.meshPoints();
 627
 628     // Some storage for edge-face-addressing.
 629     DynamicList<label> ef;
 630
 631     // Precalculate mesh edges for pp.edges.
 632     labelList meshEdges(calcMeshEdges(mesh_, pp));
 633
 634     if (debug)
 635     {
 636         // Check synchronisation
 637         // ~~~~~~~~~~~~~~~~~~~~~
 638
 639         {
 640             labelList n(mesh_.nPoints(), 0);
 641             UIndirectList<label>(n, meshPoints) = nPointLayers;
 642             syncTools::syncPointList(mesh_, n, maxEqOp<label>(), 0, false);
 643
 644             // Non-synced
 645             forAll(meshPoints, i)
 646             {
 647                 label meshPointI = meshPoints[i];
 648
 649                 if (n[meshPointI] != nPointLayers[i])
 650                 {
 651                     FatalErrorIn
 652                     (
 653                         "addPatchCellLayer::setRefinement"
 654                         "(const scalar, const indirectPrimitivePatch&"
 655                         ", const labelList&, const vectorField&"
 656                         ", polyTopoChange&)"
 657                     )   << "At mesh point:" << meshPointI
 658                         << " coordinate:" << mesh_.points()[meshPointI]
 659                         << " specified nLayers:" << nPointLayers[i] << endl
 660                         << "On coupled point a different nLayers:"
 661                         << n[meshPointI] << " was specified."
 662                         << abort(FatalError);
 663                 }
 664             }
 665
 666
 667             // Check that nPointLayers equals the max layers of connected faces
 668             // (or 0). Anything else makes no sense.
 669             labelList nFromFace(mesh_.nPoints(), 0);
 670             forAll(nFaceLayers, i)
 671             {
 672                 const face& f = pp[i];
 673
 674                 forAll(f, fp)
 675                 {
 676                     label pointI = f[fp];
 677
 678                     nFromFace[pointI] = max(nFromFace[pointI], nFaceLayers[i]);
 679                 }
 680             }
 681             syncTools::syncPointList
 682             (
 683                 mesh_,
 684                 nFromFace,
 685                 maxEqOp<label>(),
 686                 0,
 687                 false
 688             );
 689
 690             forAll(nPointLayers, i)
 691             {
 692                 label meshPointI = meshPoints[i];
 693
 694                 if
 695                 (
 696                     nPointLayers[i] > 0
 697                  && nPointLayers[i] != nFromFace[meshPointI]
 698                 )
 699                 {
 700                     FatalErrorIn
 701                     (
 702                         "addPatchCellLayer::setRefinement"
 703                         "(const scalar, const indirectPrimitivePatch&"
 704                         ", const labelList&, const vectorField&"
 705                         ", polyTopoChange&)"
 706                     )   << "At mesh point:" << meshPointI
 707                         << " coordinate:" << mesh_.points()[meshPointI]
 708                         << " specified nLayers:" << nPointLayers[i] << endl
 709                         << "but the max nLayers of surrounding faces is:"
 710                         << nFromFace[meshPointI]
 711                         << abort(FatalError);
 712                 }
 713             }
 714         }
 715
 716         {
 717             pointField d(mesh_.nPoints(), wallPoint::greatPoint);
 718             UIndirectList<point>(d, meshPoints) = firstLayerDisp;
 719             syncTools::syncPointList
 720             (
 721                 mesh_,
 722                 d,
 723                 minEqOp<vector>(),
 724                 wallPoint::greatPoint,
 725                 false
 726             );
 727
 728             forAll(meshPoints, i)
 729             {
 730                 label meshPointI = meshPoints[i];
 731
 732                 if (mag(d[meshPointI] - firstLayerDisp[i]) > SMALL)
 733                 {
 734                     FatalErrorIn
 735                     (
 736                         "addPatchCellLayer::setRefinement"
 737                         "(const scalar, const indirectPrimitivePatch&"
 738                         ", const labelList&, const vectorField&"
 739                         ", polyTopoChange&)"
 740                     )   << "At mesh point:" << meshPointI
 741                         << " coordinate:" << mesh_.points()[meshPointI]
 742                         << " specified displacement:" << firstLayerDisp[i]
 743                         << endl
 744                         << "On coupled point a different displacement:"
 745                         << d[meshPointI] << " was specified."
 746                         << abort(FatalError);
 747                 }
 748             }
 749         }
 750
 751         // Check that edges of pp (so ones that become boundary faces)
 752         // connect to only one boundary face. Guarantees uniqueness of
 753         // patch that they go into so if this is a coupled patch both
 754         // sides decide the same.
 755         // ~~~~~~~~~~~~~~~~~~~~~~
 756
 757         for (label edgeI = pp.nInternalEdges(); edgeI < pp.nEdges(); edgeI++)
 758         {
 759             const edge& e = pp.edges()[edgeI];
 760
 761             if (nPointLayers[e[0]] > 0 || nPointLayers[e[1]] > 0)
 762             {
 763                 // Edge is to become a face
 764
 765                 const labelList& eFaces = pp.edgeFaces()[edgeI];
 766
 767                 // First check: pp should be single connected.
 768                 if (eFaces.size() != 1)
 769                 {
 770                     FatalErrorIn
 771                     (
 772                         "addPatchCellLayer::setRefinement"
 773                         "(const scalar, const indirectPrimitivePatch&"
 774                         ", const labelList&, const vectorField&"
 775                         ", polyTopoChange&)"
 776                     )   << "boundary-edge-to-be-extruded:"
 777                         << pp.points()[meshPoints[e[0]]]
 778                         << pp.points()[meshPoints[e[1]]]
 779                         << " has more than two faces using it:" << eFaces
 780                         << abort(FatalError);
 781                 }
 782
 783                 label myFaceI = pp.addressing()[eFaces[0]];
 784
 785                 label meshEdgeI = meshEdges[edgeI];
 786
 787                 // Mesh faces using edge
 788
 789                 // Mesh faces using edge
 790                 const labelList& meshFaces = mesh_.edgeFaces(meshEdgeI, ef);
 791
 792                 // Check that there is only one patchface using edge.
 793                 const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 794
 795                 label bFaceI = -1;
 796
 797                 forAll(meshFaces, i)
 798                 {
 799                     label faceI = meshFaces[i];
 800
 801                     if (faceI != myFaceI)
 802                     {
 803                         if (!mesh_.isInternalFace(faceI))
 804                         {
 805                             if (bFaceI == -1)
 806                             {
 807                                 bFaceI = faceI;
 808                             }
 809                             else
 810                             {
 811                                 FatalErrorIn
 812                                 (
 813                                     "addPatchCellLayer::setRefinement"
 814                                     "(const scalar"
 815                                     ", const indirectPrimitivePatch&"
 816                                     ", const labelList&, const vectorField&"
 817                                     ", polyTopoChange&)"
 818                                 )   << "boundary-edge-to-be-extruded:"
 819                                     << pp.points()[meshPoints[e[0]]]
 820                                     << pp.points()[meshPoints[e[1]]]
 821                                     << " has more than two boundary faces"
 822                                     << " using it:"
 823                                     << bFaceI << " fc:"
 824                                     << mesh_.faceCentres()[bFaceI]
 825                                     << " patch:" << patches.whichPatch(bFaceI)
 826                                     << " and " << faceI << " fc:"
 827                                     << mesh_.faceCentres()[faceI]
 828                                     << " patch:" << patches.whichPatch(faceI)
 829                                     << abort(FatalError);
 830                             }
 831                         }
 832                     }
 833                 }
 834             }
 835         }
 836     }
 837
 838
 839     // From master point (in patch point label) to added points (in mesh point
 840     // label)
 841     addedPoints_.setSize(pp.nPoints());
 842
 843     // Mark points that do not get extruded by setting size of addedPoints_ to 0
 844     label nTruncated = 0;
 845
 846     forAll(nPointLayers, patchPointI)
 847     {
 848         if (nPointLayers[patchPointI] > 0)
 849         {
 850             addedPoints_[patchPointI].setSize(nPointLayers[patchPointI]);
 851         }
 852         else
 853         {
 854             nTruncated++;
 855         }
 856     }
 857
 858     if (debug)
 859     {
 860         Pout<< "Not adding points at " << nTruncated << " out of "
 861             << pp.nPoints() << " points" << endl;
 862     }
 863
 864
 865     //
 866     // Create new points
 867     //
 868
 869     forAll(firstLayerDisp, patchPointI)
 870     {
 871         if (addedPoints_[patchPointI].size())
 872         {
 873             label meshPointI = meshPoints[patchPointI];
 874
 875             label zoneI = mesh_.pointZones().whichZone(meshPointI);
 876
 877             point pt = mesh_.points()[meshPointI];
 878
 879             vector disp = firstLayerDisp[patchPointI];
 880
 881             for (label i = 0; i < addedPoints_[patchPointI].size(); i++)
 882             {
 883                 pt += disp;
 884
 885                 label addedVertI = meshMod.setAction
 886                 (
 887                     polyAddPoint
 888                     (
 889                         pt,         // point
 890                         meshPointI, // master point
 891                         zoneI,      // zone for point
 892                         true        // supports a cell
 893                     )
 894                 );
 895
 896                 addedPoints_[patchPointI][i] = addedVertI;
 897
 898                 disp *= expansionRatio[patchPointI];
 899             }
 900         }
 901     }
 902
 903
 904     //
 905     // Add cells to all boundaryFaces
 906     //
 907
 908     labelListList addedCells(pp.size());
 909
 910     forAll(pp, patchFaceI)
 911     {
 912         if (nFaceLayers[patchFaceI] > 0)
 913         {
 914             addedCells[patchFaceI].setSize(nFaceLayers[patchFaceI]);
 915
 916             label meshFaceI = pp.addressing()[patchFaceI];
 917
 918             label ownZoneI = mesh_.cellZones().whichZone
 919             (
 920                 mesh_.faceOwner()[meshFaceI]
 921             );
 922
 923             for (label i = 0; i < nFaceLayers[patchFaceI]; i++)
 924             {
 925                 // Note: add from cell (owner of patch face) or from face?
 926                 // for now add from cell so we can map easily.
 927                 addedCells[patchFaceI][i] = meshMod.setAction
 928                 (
 929                     polyAddCell
 930                     (
 931                         -1,             // master point
 932                         -1,             // master edge
 933                         -1,             // master face
 934                         mesh_.faceOwner()[meshFaceI],   // master cell id
 935                         ownZoneI        // zone for cell
 936                     )
 937                 );
 938
 939                 //Pout<< "For patchFace:" << patchFaceI
 940                 //    << " meshFace:" << pp.addressing()[patchFaceI]
 941                 //    << " layer:" << i << " added cell:"
 942                 //    << addedCells[patchFaceI][i]
 943                 //    << endl;
 944             }
 945         }
 946     }
 947
 948
 949     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 950
 951     // Precalculated patchID for each patch face
 952     labelList patchID(pp.size());
 953
 954     forAll(pp, patchFaceI)
 955     {
 956         label meshFaceI = pp.addressing()[patchFaceI];
 957
 958         patchID[patchFaceI] = patches.whichPatch(meshFaceI);
 959     }
 960
 961
 962
 963     // Create faces on top of the original patch faces.
 964     // These faces are created from original patch faces outwards so in order
 965     // of increasing cell number. So orientation should be same as original
 966     // patch face for them to have owner<neighbour.
 967
 968     layerFaces_.setSize(pp.size());
 969
 970     forAll(pp.localFaces(), patchFaceI)
 971     {
 972         label meshFaceI = pp.addressing()[patchFaceI];
 973
 974         if (addedCells[patchFaceI].size())
 975         {
 976             layerFaces_[patchFaceI].setSize(addedCells[patchFaceI].size() + 1);
 977             layerFaces_[patchFaceI][0] = meshFaceI;
 978
 979             // Get duplicated vertices on the patch face.
 980             const face& f = pp.localFaces()[patchFaceI];
 981
 982             face newFace(f.size());
 983
 984             for (label i = 0; i < addedCells[patchFaceI].size(); i++)
 985             {
 986                 forAll(f, fp)
 987                 {
 988                     if (addedPoints_[f[fp]].empty())
 989                     {
 990                         // Keep original point
 991                         newFace[fp] = meshPoints[f[fp]];
 992                     }
 993                     else
 994                     {
 995                         // Get new outside point
 996                         label offset =
 997                             addedPoints_[f[fp]].size()
 998                           - addedCells[patchFaceI].size();
 999                         newFace[fp] = addedPoints_[f[fp]][i+offset];
1000                     }
1001                 }
1002
1003
1004                 // Get new neighbour
1005                 label nei;
1006                 label patchI;
1007                 label zoneI = -1;
1008                 bool flip = false;
1009
1010
1011                 if (i == addedCells[patchFaceI].size()-1)
1012                 {
1013                     // Top layer so is patch face.
1014                     nei = -1;
1015                     patchI = patchID[patchFaceI];
1016                     zoneI = mesh_.faceZones().whichZone(meshFaceI);
1017                     if (zoneI != -1)
1018                     {
1019                         const faceZone& fz = mesh_.faceZones()[zoneI];
1020                         flip = fz.flipMap()[fz.whichFace(meshFaceI)];
1021                     }
1022                 }
1023                 else
1024                 {
1025                     // Internal face between layer i and i+1
1026                     nei = addedCells[patchFaceI][i+1];
1027                     patchI = -1;
1028                 }
1029
1030
1031                 layerFaces_[patchFaceI][i+1] = meshMod.setAction
1032                 (
1033                     polyAddFace
1034                     (
1035                         newFace,                    // face
1036                         addedCells[patchFaceI][i],  // owner
1037                         nei,                        // neighbour
1038                         -1,                         // master point
1039                         -1,                         // master edge
1040                         meshFaceI,                  // master face for addition
1041                         false,                      // flux flip
1042                         patchI,                     // patch for face
1043                         zoneI,                      // zone for face
1044                         flip                        // face zone flip
1045                     )
1046                 );
1047                 //Pout<< "Added inbetween face " << newFace
1048                 //    << " own:" << addedCells[patchFaceI][i]
1049                 //    << " nei:" << nei
1050                 //    << " patch:" << patchI
1051                 //    << endl;
1052             }
1053         }
1054     }
1055
1056     //
1057     // Modify old patch faces to be on the inside
1058     //
1059     forAll(pp, patchFaceI)
1060     {
1061         if (addedCells[patchFaceI].size())
1062         {
1063             label meshFaceI = pp.addressing()[patchFaceI];
1064
1065             meshMod.setAction
1066             (
1067                 polyModifyFace
1068                 (
1069                     pp[patchFaceI],                 // modified face
1070                     meshFaceI,                      // label of face
1071                     mesh_.faceOwner()[meshFaceI],   // owner
1072                     addedCells[patchFaceI][0],      // neighbour
1073                     false,                          // face flip
1074                     -1,                             // patch for face
1075                     true, //false,                        // remove from zone
1076                     -1, //zoneI,                          // zone for face
1077                     false                           // face flip in zone
1078                 )
1079             );
1080             //Pout<< "Modified old patch face " << meshFaceI
1081             //    << " own:" << mesh_.faceOwner()[meshFaceI]
1082             //    << " nei:" << addedCells[patchFaceI][0]
1083             //    << endl;
1084         }
1085     }
1086
1087
1088     //
1089     // Create 'side' faces, one per edge that is being extended.
1090     //
1091
1092     const labelListList& faceEdges = pp.faceEdges();
1093     const faceList& localFaces = pp.localFaces();
1094     const edgeList& edges = pp.edges();
1095
1096     // Get number of layers per edge. This is 0 if edge is not extruded;
1097     // max of connected faces otherwise.
1098     labelList edgeLayers(pp.nEdges());
1099
1100     {
1101         // Use list over mesh.nEdges() since syncTools does not yet support
1102         // partial list synchronisation.
1103         labelList meshEdgeLayers(mesh_.nEdges(), -1);
1104
1105         forAll(meshEdges, edgeI)
1106         {
1107             const edge& e = edges[edgeI];
1108
1109             label meshEdgeI = meshEdges[edgeI];
1110
1111             if ((nPointLayers[e[0]] == 0) && (nPointLayers[e[1]] == 0))
1112             {
1113                 meshEdgeLayers[meshEdgeI] = 0;
1114             }
1115             else
1116             {
1117                 const labelList& eFaces = pp.edgeFaces()[edgeI];
1118
1119                 forAll(eFaces, i)
1120                 {
1121                     meshEdgeLayers[meshEdgeI] = max
1122                     (
1123                         nFaceLayers[eFaces[i]],
1124                         meshEdgeLayers[meshEdgeI]
1125                     );
1126                 }
1127             }
1128         }
1129
1130         syncTools::syncEdgeList
1131         (
1132             mesh_,
1133             meshEdgeLayers,
1134             maxEqOp<label>(),
1135             0,                  // initial value
1136             false               // no separation
1137         );
1138
1139         forAll(meshEdges, edgeI)
1140         {
1141             edgeLayers[edgeI] = meshEdgeLayers[meshEdges[edgeI]];
1142         }
1143     }
1144
1145
1146     // Global indices engine
1147     const globalIndex globalFaces(mesh_.nFaces());
1148
1149     // Get for all pp edgeFaces a unique faceID
1150     labelListList globalEdgeFaces
1151     (
1152          calcGlobalEdgeFaces
1153          (
1154             mesh_,
1155             globalFaces,
1156             pp,
1157             meshEdges
1158         )
1159     );
1160
1161
1162     // Mark off which edges have been extruded
1163     boolList doneEdge(pp.nEdges(), false);
1164
1165
1166     // Create faces. Per face walk connected edges and find string of edges
1167     // between the same two faces and extrude string into a single face.
1168     forAll(pp, patchFaceI)
1169     {
1170         const labelList& fEdges = faceEdges[patchFaceI];
1171
1172         forAll(fEdges, fp)
1173         {
1174             // Get string of edges that needs to be extruded as a single face.
1175             // Returned as indices in fEdges.
1176             labelPair indexPair
1177             (
1178                 getEdgeString
1179                 (
1180                     pp,
1181                     globalEdgeFaces,
1182                     doneEdge,
1183                     patchFaceI,
1184                     globalFaces.toGlobal(pp.addressing()[patchFaceI])
1185                 )
1186             );
1187
1188             //Pout<< "Found unextruded edges in edges:" << fEdges
1189             //    << " start:" << indexPair[0]
1190             //    << " end:" << indexPair[1]
1191             //    << endl;
1192
1193             const label startFp = indexPair[0];
1194             const label endFp = indexPair[1];
1195
1196             if (startFp != -1)
1197             {
1198                 // Extrude edges from indexPair[0] up to indexPair[1]
1199                 // (note indexPair = indices of edges. There is one more vertex
1200                 //  than edges)
1201                 const face& f = localFaces[patchFaceI];
1202
1203                 labelList stringedVerts;
1204                 if (endFp >= startFp)
1205                 {
1206                     stringedVerts.setSize(endFp-startFp+2);
1207                 }
1208                 else
1209                 {
1210                     stringedVerts.setSize(endFp+f.size()-startFp+2);
1211                 }
1212
1213                 label fp = startFp;
1214
1215                 for (label i = 0; i < stringedVerts.size()-1; i++)
1216                 {
1217                     stringedVerts[i] = f[fp];
1218                     doneEdge[fEdges[fp]] = true;
1219                     fp = f.fcIndex(fp);
1220                 }
1221                 stringedVerts[stringedVerts.size()-1] = f[fp];
1222
1223
1224                 // Now stringedVerts contains the vertices in order of face f.
1225                 // This is consistent with the order if f becomes the owner cell
1226                 // and nbrFaceI the neighbour cell. Note that the cells get
1227                 // added in order of pp so we can just use face ordering and
1228                 // because we loop in incrementing order as well we will
1229                 // always have nbrFaceI > patchFaceI.
1230
1231                 label startEdgeI = fEdges[startFp];
1232
1233                 label meshEdgeI = meshEdges[startEdgeI];
1234
1235                 label numEdgeSideFaces = edgeLayers[startEdgeI];
1236
1237                 for (label i = 0; i < numEdgeSideFaces; i++)
1238                 {
1239                     label vEnd = stringedVerts[stringedVerts.size()-1];
1240                     label vStart = stringedVerts[0];
1241
1242                     // calculate number of points making up a face
1243                     label newFp = 2*stringedVerts.size();
1244
1245                     if (i == 0)
1246                     {
1247                         // layer 0 gets all the truncation of neighbouring
1248                         // faces with more layers.
1249                         if (addedPoints_[vEnd].size())
1250                         {
1251                             newFp +=
1252                                 addedPoints_[vEnd].size() - numEdgeSideFaces;
1253                         }
1254                         if (addedPoints_[vStart].size())
1255                         {
1256                             newFp +=
1257                                 addedPoints_[vStart].size() - numEdgeSideFaces;
1258                         }
1259                     }
1260
1261                     face newFace(newFp);
1262
1263                     newFp = 0;
1264
1265                     // For layer 0 get pp points, for all other layers get
1266                     // points of layer-1.
1267                     if (i == 0)
1268                     {
1269                         forAll(stringedVerts, stringedI)
1270                         {
1271                             label v = stringedVerts[stringedI];
1272                             addVertex(meshPoints[v], newFace, newFp);
1273                         }
1274                     }
1275                     else
1276                     {
1277                         forAll(stringedVerts, stringedI)
1278                         {
1279                             label v = stringedVerts[stringedI];
1280                             if (addedPoints_[v].size())
1281                             {
1282                                 label offset =
1283                                     addedPoints_[v].size() - numEdgeSideFaces;
1284                                 addVertex
1285                                 (
1286                                     addedPoints_[v][i+offset-1],
1287                                     newFace,
1288                                     newFp
1289                                 );
1290                             }
1291                             else
1292                             {
1293                                 addVertex(meshPoints[v], newFace, newFp);
1294                             }
1295                         }
1296                     }
1297
1298                     // add points between stringed vertices (end)
1299                     if (numEdgeSideFaces < addedPoints_[vEnd].size())
1300                     {
1301                         if (i == 0 && addedPoints_[vEnd].size())
1302                         {
1303                             label offset =
1304                                 addedPoints_[vEnd].size() - numEdgeSideFaces;
1305                             for (label ioff = 0; ioff < offset; ioff++)
1306                             {
1307                                 addVertex
1308                                 (
1309                                     addedPoints_[vEnd][ioff],
1310                                     newFace,
1311                                     newFp
1312                                 );
1313                             }
1314                         }
1315                     }
1316
1317                     forAllReverse(stringedVerts, stringedI)
1318                     {
1319                         label v = stringedVerts[stringedI];
1320                         if (addedPoints_[v].size())
1321                         {
1322                             label offset =
1323                                 addedPoints_[v].size() - numEdgeSideFaces;
1324                             addVertex
1325                             (
1326                                 addedPoints_[v][i+offset],
1327                                 newFace,
1328                                 newFp
1329                             );
1330                         }
1331                         else
1332                         {
1333                             addVertex(meshPoints[v], newFace, newFp);
1334                         }
1335                     }
1336
1337
1338                     // add points between stringed vertices (start)
1339                     if (numEdgeSideFaces < addedPoints_[vStart].size())
1340                     {
1341                         if (i == 0 && addedPoints_[vStart].size())
1342                         {
1343                             label offset =
1344                                 addedPoints_[vStart].size() - numEdgeSideFaces;
1345                             for (label ioff = offset; ioff > 0; ioff--)
1346                             {
1347                                 addVertex
1348                                 (
1349                                     addedPoints_[vStart][ioff-1],
1350                                     newFace,
1351                                     newFp
1352                                 );
1353                             }
1354                         }
1355                     }
1356
1357                     if (newFp >= 3)
1358                     {
1359                         // Add face inbetween faces patchFaceI and nbrFaceI
1360                         // (possibly -1 for external edges)
1361
1362                         newFace.setSize(newFp);
1363
1364                         label nbrFaceI = nbrFace
1365                         (
1366                             pp.edgeFaces(),
1367                             startEdgeI,
1368                             patchFaceI
1369                         );
1370
1371                         const labelList& meshFaces = mesh_.edgeFaces
1372                         (
1373                             meshEdgeI,
1374                             ef
1375                         );
1376
1377                         addSideFace
1378                         (
1379                             pp,
1380                             patchID,
1381                             addedCells,
1382                             newFace,
1383                             patchFaceI,
1384                             nbrFaceI,
1385                             startEdgeI,     // edge to inflate from
1386                             meshEdgeI,      // corresponding mesh edge
1387                             i,
1388                             numEdgeSideFaces,
1389                             meshFaces,
1390                             meshMod
1391                         );
1392                     }
1393                 }
1394             }
1395         }
1396     }
1397 }
1398
1399
1400 void Foam::addPatchCellLayer::updateMesh
1401 (
1402     const mapPolyMesh& morphMap,
1403     const labelList& faceMap,   // new to old patch faces
1404     const labelList& pointMap   // new to old patch points
1405 )
1406 {
1407     {
1408         labelListList newAddedPoints(pointMap.size());
1409
1410         forAll(newAddedPoints, newPointI)
1411         {
1412             label oldPointI = pointMap[newPointI];
1413
1414             const labelList& added = addedPoints_[oldPointI];
1415
1416             labelList& newAdded = newAddedPoints[newPointI];
1417             newAdded.setSize(added.size());
1418             label newI = 0;
1419
1420             forAll(added, i)
1421             {
1422                 label newPointI = morphMap.reversePointMap()[added[i]];
1423
1424                 if (newPointI >= 0)
1425                 {
1426                     newAdded[newI++] = newPointI;
1427                 }
1428             }
1429             newAdded.setSize(newI);
1430         }
1431         addedPoints_.transfer(newAddedPoints);
1432     }
1433
1434     {
1435         labelListList newLayerFaces(faceMap.size());
1436
1437         forAll(newLayerFaces, newFaceI)
1438         {
1439             label oldFaceI = faceMap[newFaceI];
1440
1441             const labelList& added = layerFaces_[oldFaceI];
1442
1443             labelList& newAdded = newLayerFaces[newFaceI];
1444             newAdded.setSize(added.size());
1445             label newI = 0;
1446
1447             forAll(added, i)
1448             {
1449                 label newFaceI = morphMap.reverseFaceMap()[added[i]];
1450
1451                 if (newFaceI >= 0)
1452                 {
1453                     newAdded[newI++] = newFaceI;
1454                 }
1455             }
1456             newAdded.setSize(newI);
1457         }
1458         layerFaces_.transfer(newLayerFaces);
1459     }
1460 }
1461
1462
1463 Foam::labelList Foam::addPatchCellLayer::calcMeshEdges
1464 (
1465     const primitiveMesh& mesh,
1466     const indirectPrimitivePatch& pp
1467 )
1468 {
1469     labelList meshEdges(pp.nEdges());
1470
1471     forAll(meshEdges, patchEdgeI)
1472     {
1473         const edge& e = pp.edges()[patchEdgeI];
1474
1475         label v0 = pp.meshPoints()[e[0]];
1476         label v1 = pp.meshPoints()[e[1]];
1477         meshEdges[patchEdgeI] = meshTools::findEdge
1478         (
1479             mesh.edges(),
1480             mesh.pointEdges()[v0],
1481             v0,
1482             v1
1483         );
1484     }
1485     return meshEdges;
1486 }
1487
1488
1489 // ************************************************************************* //