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