createShellMesh.C 28.3 KB
Newer Older
1 2 3 4
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
OpenFOAM bot's avatar
OpenFOAM bot committed
5
    \\  /    A nd           | www.openfoam.com
6
     \\/     M anipulation  |
OpenFOAM bot's avatar
OpenFOAM bot committed
7
-------------------------------------------------------------------------------
OpenFOAM bot's avatar
OpenFOAM bot committed
8 9
    Copyright (C) 2011-2016 OpenFOAM Foundation
    Copyright (C) 2019 OpenCFD Ltd.
10 11 12 13
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

14 15 16 17
    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.
18 19 20 21 22 23 24

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
25
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
26 27 28 29 30 31 32 33 34 35 36

\*---------------------------------------------------------------------------*/

#include "createShellMesh.H"
#include "polyTopoChange.H"
#include "meshTools.H"
#include "mapPolyMesh.H"
#include "polyAddPoint.H"
#include "polyAddFace.H"
#include "polyModifyFace.H"
#include "polyAddCell.H"
37 38 39 40 41 42
#include "labelPair.H"
#include "indirectPrimitivePatch.H"
#include "mapDistribute.H"
#include "globalMeshData.H"
#include "PatchTools.H"
#include "globalIndex.H"
43 44 45

// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

46 47
namespace Foam
{
48 49
defineTypeNameAndDebug(createShellMesh, 0);

50 51 52 53 54 55 56 57 58 59 60 61
template<>
class minEqOp<labelPair>
{
public:
    void operator()(labelPair& x, const labelPair& y) const
    {
        x[0] = min(x[0], y[0]);
        x[1] = min(x[1], y[1]);
    }
};
}

62 63 64

// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //

65 66 67 68 69 70 71
// Synchronise edges
void Foam::createShellMesh::syncEdges
(
    const globalMeshData& globalData,

    const labelList& patchEdges,
    const labelList& coupledEdges,
72
    const bitSet& sameEdgeOrientation,
73
    const bool syncNonCollocated,
74

75
    bitSet& isChangedEdge,
76 77 78 79 80
    DynamicList<label>& changedEdges,
    labelPairList& allEdgeData
)
{
    const mapDistribute& map = globalData.globalEdgeSlavesMap();
81
    const bitSet& cppOrientation = globalData.globalEdgeOrientation();
82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117

    // Convert patch-edge data into cpp-edge data
    labelPairList cppEdgeData
    (
        map.constructSize(),
        labelPair(labelMax, labelMax)
    );

    forAll(patchEdges, i)
    {
        label patchEdgeI = patchEdges[i];
        label coupledEdgeI = coupledEdges[i];

        if (isChangedEdge[patchEdgeI])
        {
            const labelPair& data = allEdgeData[patchEdgeI];

            // Patch-edge data needs to be converted into coupled-edge data
            // (optionally flipped) and consistent in orientation with
            // other coupled edge (optionally flipped)
            if (sameEdgeOrientation[i] == cppOrientation[coupledEdgeI])
            {
                cppEdgeData[coupledEdgeI] = data;
            }
            else
            {
                cppEdgeData[coupledEdgeI] = labelPair(data[1], data[0]);
            }
        }
    }

    // Synchronise
    globalData.syncData
    (
        cppEdgeData,
        globalData.globalEdgeSlaves(),
118 119 120 121 122
        (
            syncNonCollocated
          ? globalData.globalEdgeTransformedSlaves()    // transformed elems
          : labelListList(globalData.globalEdgeSlaves().size()) //no transformed
        ),
123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148
        map,
        minEqOp<labelPair>()
    );

    // Back from cpp-edge to patch-edge data
    forAll(patchEdges, i)
    {
        label patchEdgeI = patchEdges[i];
        label coupledEdgeI = coupledEdges[i];

        if (cppEdgeData[coupledEdgeI] != labelPair(labelMax, labelMax))
        {
            const labelPair& data = cppEdgeData[coupledEdgeI];

            if (sameEdgeOrientation[i] == cppOrientation[coupledEdgeI])
            {
                allEdgeData[patchEdgeI] = data;
            }
            else
            {
                allEdgeData[patchEdgeI] = labelPair(data[1], data[0]);
            }

            if (!isChangedEdge[patchEdgeI])
            {
                changedEdges.append(patchEdgeI);
149
                isChangedEdge.set(patchEdgeI);
150 151 152 153 154 155
            }
        }
    }
}


156 157
void Foam::createShellMesh::calcPointRegions
(
158
    const globalMeshData& globalData,
159
    const primitiveFacePatch& patch,
160
    const bitSet& nonManifoldEdge,
161
    const bool syncNonCollocated,
162 163 164 165

    faceList& pointGlobalRegions,
    faceList& pointLocalRegions,
    labelList& localToGlobalRegion
166 167
)
{
168 169 170 171 172
    const indirectPrimitivePatch& cpp = globalData.coupledPatch();

    // Calculate correspondence between patch and globalData.coupledPatch.
    labelList patchEdges;
    labelList coupledEdges;
173
    bitSet sameEdgeOrientation;
174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190
    PatchTools::matchEdges
    (
        cpp,
        patch,

        coupledEdges,
        patchEdges,
        sameEdgeOrientation
    );


    // Initial unique regions
    // ~~~~~~~~~~~~~~~~~~~~~~
    // These get merged later on across connected edges.

    // 1. Count
    label nMaxRegions = 0;
191
    forAll(patch.localFaces(), facei)
192
    {
193
        const face& f = patch.localFaces()[facei];
194
        nMaxRegions += f.size();
195 196
    }

197 198 199
    const globalIndex globalRegions(nMaxRegions);

    // 2. Assign unique regions
200 201
    label nRegions = 0;

202
    pointGlobalRegions.setSize(patch.size());
203
    forAll(pointGlobalRegions, facei)
204
    {
205 206
        const face& f = patch.localFaces()[facei];
        labelList& pRegions = pointGlobalRegions[facei];
207 208 209 210 211 212
        pRegions.setSize(f.size());
        forAll(pRegions, fp)
        {
            pRegions[fp] = globalRegions.toGlobal(nRegions++);
        }
    }
213

214 215 216

    DynamicList<label> changedEdges(patch.nEdges());
    labelPairList allEdgeData(patch.nEdges(), labelPair(labelMax, labelMax));
217
    bitSet isChangedEdge(patch.nEdges());
218 219 220 221 222 223 224 225 226 227 228


    // Fill initial seed
    // ~~~~~~~~~~~~~~~~~

    forAll(patch.edgeFaces(), edgeI)
    {
        if (!nonManifoldEdge[edgeI])
        {
            // Take over value from one face only.
            const edge& e = patch.edges()[edgeI];
229 230
            label facei = patch.edgeFaces()[edgeI][0];
            const face& f = patch.localFaces()[facei];
231

232 233
            label fp0 = f.find(e[0]);
            label fp1 = f.find(e[1]);
234 235
            allEdgeData[edgeI] = labelPair
            (
236 237
                pointGlobalRegions[facei][fp0],
                pointGlobalRegions[facei][fp1]
238 239 240 241
            );
            if (!isChangedEdge[edgeI])
            {
                changedEdges.append(edgeI);
242
                isChangedEdge.set(edgeI);
243 244 245 246 247 248 249 250 251 252 253 254
            }
        }
    }


    syncEdges
    (
        globalData,

        patchEdges,
        coupledEdges,
        sameEdgeOrientation,
255
        syncNonCollocated,
256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272

        isChangedEdge,
        changedEdges,
        allEdgeData
    );


    // Edge-Face-Edge walk across patch
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    // Across edge minimum regions win

    while (true)
    {
        // From edge to face
        // ~~~~~~~~~~~~~~~~~

        DynamicList<label> changedFaces(patch.size());
273
        bitSet isChangedFace(patch.size());
274 275 276 277 278 279 280 281 282 283 284

        forAll(changedEdges, changedI)
        {
            label edgeI = changedEdges[changedI];
            const labelPair& edgeData = allEdgeData[edgeI];

            const edge& e = patch.edges()[edgeI];
            const labelList& eFaces = patch.edgeFaces()[edgeI];

            forAll(eFaces, i)
            {
285 286
                label facei = eFaces[i];
                const face& f = patch.localFaces()[facei];
287 288

                // Combine edgeData with face data
289
                label fp0 = f.find(e[0]);
290
                if (pointGlobalRegions[facei][fp0] > edgeData[0])
291
                {
292 293
                    pointGlobalRegions[facei][fp0] = edgeData[0];
                    if (!isChangedFace[facei])
294
                    {
295
                        isChangedFace.set(facei);
296
                        changedFaces.append(facei);
297 298 299
                    }
                }

300
                label fp1 = f.find(e[1]);
301
                if (pointGlobalRegions[facei][fp1] > edgeData[1])
302
                {
303 304
                    pointGlobalRegions[facei][fp1] = edgeData[1];
                    if (!isChangedFace[facei])
305
                    {
306
                        isChangedFace.set(facei);
307
                        changedFaces.append(facei);
308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327
                    }
                }
            }
        }


        label nChangedFaces = returnReduce(changedFaces.size(), sumOp<label>());
        if (nChangedFaces == 0)
        {
            break;
        }


        // From face to edge
        // ~~~~~~~~~~~~~~~~~

        isChangedEdge = false;
        changedEdges.clear();

        forAll(changedFaces, i)
328
        {
329 330 331
            label facei = changedFaces[i];
            const face& f = patch.localFaces()[facei];
            const labelList& fEdges = patch.faceEdges()[facei];
332 333

            forAll(fEdges, fp)
334
            {
335 336 337
                label edgeI = fEdges[fp];

                if (!nonManifoldEdge[edgeI])
338
                {
339
                    const edge& e = patch.edges()[edgeI];
340
                    label fp0 = f.find(e[0]);
341
                    label region0 = pointGlobalRegions[facei][fp0];
342
                    label fp1 = f.find(e[1]);
343
                    label region1 = pointGlobalRegions[facei][fp1];
344

345 346 347 348 349
                    if
                    (
                        (allEdgeData[edgeI][0] > region0)
                     || (allEdgeData[edgeI][1] > region1)
                    )
350
                    {
351 352
                        allEdgeData[edgeI] = labelPair(region0, region1);
                        if (!isChangedEdge[edgeI])
mattijs's avatar
mattijs committed
353
                        {
354
                            changedEdges.append(edgeI);
355
                            isChangedEdge.set(edgeI);
356 357 358 359 360 361
                        }
                    }
                }
            }
        }

362 363 364
        syncEdges
        (
            globalData,
365

366 367 368
            patchEdges,
            coupledEdges,
            sameEdgeOrientation,
369
            syncNonCollocated,
370

371 372 373 374 375 376 377 378
            isChangedEdge,
            changedEdges,
            allEdgeData
        );


        label nChangedEdges = returnReduce(changedEdges.size(), sumOp<label>());
        if (nChangedEdges == 0)
379
        {
380
            break;
381 382 383 384
        }
    }


385 386 387 388 389 390 391 392

    // Assign local regions
    // ~~~~~~~~~~~~~~~~~~~~

    // Calculate addressing from global region back to local region
    pointLocalRegions.setSize(patch.size());
    Map<label> globalToLocalRegion(globalRegions.localSize()/4);
    DynamicList<label> dynLocalToGlobalRegion(globalToLocalRegion.size());
393
    forAll(patch.localFaces(), facei)
394
    {
395 396
        const face& f = patch.localFaces()[facei];
        face& pRegions = pointLocalRegions[facei];
397 398
        pRegions.setSize(f.size());
        forAll(f, fp)
399
        {
400
            label globalRegionI = pointGlobalRegions[facei][fp];
401

402
            const auto fnd = globalToLocalRegion.cfind(globalRegionI);
403

404
            if (fnd.found())
405 406 407 408 409 410 411 412 413 414 415
            {
                // Already encountered this global region. Assign same local one
                pRegions[fp] = fnd();
            }
            else
            {
                // Region not yet seen. Create new one
                label localRegionI = globalToLocalRegion.size();
                pRegions[fp] = localRegionI;
                globalToLocalRegion.insert(globalRegionI, localRegionI);
                dynLocalToGlobalRegion.append(globalRegionI);
416 417 418
            }
        }
    }
419
    localToGlobalRegion.transfer(dynLocalToGlobalRegion);
420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437
}


// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

Foam::createShellMesh::createShellMesh
(
    const primitiveFacePatch& patch,
    const faceList& pointRegions,
    const labelList& regionPoints
)
:
    patch_(patch),
    pointRegions_(pointRegions),
    regionPoints_(regionPoints)
{
    if (pointRegions_.size() != patch_.size())
    {
438
        FatalErrorInFunction
439 440 441 442 443 444 445 446 447 448 449
            << "nFaces:" << patch_.size()
            << " pointRegions:" << pointRegions.size()
            << exit(FatalError);
    }
}


// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

void Foam::createShellMesh::setRefinement
(
mattijs's avatar
mattijs committed
450 451 452
    const pointField& firstLayerDisp,
    const scalar expansionRatio,
    const label nLayers,
453 454
    const labelList& topPatchID,
    const labelList& bottomPatchID,
455 456 457 458
    const labelListList& extrudeEdgePatches,
    polyTopoChange& meshMod
)
{
mattijs's avatar
mattijs committed
459
    if (firstLayerDisp.size() != regionPoints_.size())
460
    {
461
        FatalErrorInFunction
462
            << "nRegions:" << regionPoints_.size()
mattijs's avatar
mattijs committed
463
            << " firstLayerDisp:" << firstLayerDisp.size()
464 465 466 467 468
            << exit(FatalError);
    }

    if
    (
469 470
        topPatchID.size() != patch_.size()
     && bottomPatchID.size() != patch_.size()
471 472
    )
    {
473
        FatalErrorInFunction
474
            << "nFaces:" << patch_.size()
475 476
            << " topPatchID:" << topPatchID.size()
            << " bottomPatchID:" << bottomPatchID.size()
477 478 479 480 481
            << exit(FatalError);
    }

    if (extrudeEdgePatches.size() != patch_.nEdges())
    {
482
        FatalErrorInFunction
483 484 485 486 487 488 489 490
            << "nEdges:" << patch_.nEdges()
            << " extrudeEdgePatches:" << extrudeEdgePatches.size()
            << exit(FatalError);
    }



    // From cell to patch (trivial)
mattijs's avatar
mattijs committed
491
    DynamicList<label> cellToFaceMap(nLayers*patch_.size());
492
    // From face to patch+turning index
mattijs's avatar
mattijs committed
493 494 495 496
    DynamicList<label> faceToFaceMap
    (
        (nLayers+1)*(patch_.size()+patch_.nEdges())
    );
497
    // From face to patch edge index
mattijs's avatar
mattijs committed
498
    DynamicList<label> faceToEdgeMap(nLayers*(patch_.nEdges()+patch_.nEdges()));
499
    // From point to patch point index
mattijs's avatar
mattijs committed
500
    DynamicList<label> pointToPointMap((nLayers+1)*patch_.nPoints());
501 502 503 504 505


    // Introduce new cell for every face
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

mattijs's avatar
mattijs committed
506
    labelList addedCells(nLayers*patch_.size());
507
    forAll(patch_, facei)
508
    {
mattijs's avatar
mattijs committed
509 510
        for (label layerI = 0; layerI < nLayers; layerI++)
        {
511
            addedCells[nLayers*facei+layerI] = meshMod.addCell
mattijs's avatar
mattijs committed
512 513 514 515 516 517 518
            (
                -1,                     // masterPointID
                -1,                     // masterEdgeID
                -1,                     // masterFaceID
                cellToFaceMap.size(),   // masterCellID
                -1                      // zoneID
            );
519
            cellToFaceMap.append(facei);
mattijs's avatar
mattijs committed
520
        }
521 522
    }

523

524 525 526 527
    // Introduce original points
    // ~~~~~~~~~~~~~~~~~~~~~~~~~

    // Original point numbers in local point ordering so no need to store.
528
    forAll(patch_.localPoints(), pointi)
529
    {
530
        //label addedPointi =
531 532
        meshMod.addPoint
        (
533
            patch_.localPoints()[pointi],   // point
mattijs's avatar
mattijs committed
534
            pointToPointMap.size(),         // masterPointID
535 536 537
            -1,                             // zoneID
            true                            // inCell
        );
538
        pointToPointMap.append(pointi);
539

540 541 542
        //Pout<< "Added bottom point " << addedPointi
        //    << " at " << patch_.localPoints()[pointi]
        //    << "  from point " << pointi
543
        //    << endl;
544 545
    }

546

547 548 549
    // Introduce new points (one for every region)
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

mattijs's avatar
mattijs committed
550
    labelList addedPoints(nLayers*regionPoints_.size());
551 552
    forAll(regionPoints_, regionI)
    {
553
        label pointi = regionPoints_[regionI];
554

555
        point pt = patch_.localPoints()[pointi];
mattijs's avatar
mattijs committed
556 557 558 559
        point disp = firstLayerDisp[regionI];
        for (label layerI = 0; layerI < nLayers; layerI++)
        {
            pt += disp;
560

mattijs's avatar
mattijs committed
561 562 563 564 565 566 567
            addedPoints[nLayers*regionI+layerI] = meshMod.addPoint
            (
                pt,                     // point
                pointToPointMap.size(), // masterPointID - used only addressing
                -1,                     // zoneID
                true                    // inCell
            );
568
            pointToPointMap.append(pointi);
mattijs's avatar
mattijs committed
569 570 571

            disp *= expansionRatio;
        }
572 573 574
    }


575
    // Add face on bottom side
576
    forAll(patch_.localFaces(), facei)
577 578 579
    {
        meshMod.addFace
        (
580 581
            patch_.localFaces()[facei].reverseFace(),// vertices
            addedCells[nLayers*facei],  // own
582 583 584
            -1,                         // nei
            -1,                         // masterPointID
            -1,                         // masterEdgeID
585
            faceToFaceMap.size(),       // masterFaceID : current facei
586
            true,                       // flipFaceFlux
587
            bottomPatchID[facei],       // patchID
588 589 590
            -1,                         // zoneID
            false                       // zoneFlip
        );
591
        faceToFaceMap.append(-facei-1); // points to flipped original face
592 593
        faceToEdgeMap.append(-1);

594
        //const face newF(patch_.localFaces()[facei].reverseFace());
595
        //Pout<< "Added bottom face "
mattijs's avatar
mattijs committed
596 597
        //    << newF
        //    << " coords:" << UIndirectList<point>(meshMod.points(), newF)
598 599 600
        //    << " own " << addedCells[facei]
        //    << " patch:" << bottomPatchID[facei]
        //    << "  at " << patch_.faceCentres()[facei]
601 602 603
        //    << endl;
    }

mattijs's avatar
mattijs committed
604
    // Add inbetween faces and face on top
605
    forAll(patch_.localFaces(), facei)
606 607
    {
        // Get face in original ordering
608
        const face& f = patch_.localFaces()[facei];
609 610

        face newF(f.size());
mattijs's avatar
mattijs committed
611 612

        for (label layerI = 0; layerI < nLayers; layerI++)
613
        {
mattijs's avatar
mattijs committed
614 615 616
            // Pick up point based on region and layer
            forAll(f, fp)
            {
617
                label region = pointRegions_[facei][fp];
mattijs's avatar
mattijs committed
618 619
                newF[fp] = addedPoints[region*nLayers+layerI];
            }
620

621
            label own = addedCells[facei*nLayers+layerI];
mattijs's avatar
mattijs committed
622
            label nei;
623
            label patchi;
mattijs's avatar
mattijs committed
624 625 626
            if (layerI == nLayers-1)
            {
                nei = -1;
627
                patchi = topPatchID[facei];
mattijs's avatar
mattijs committed
628 629 630
            }
            else
            {
631
                nei = addedCells[facei*nLayers+layerI+1];
632
                patchi = -1;
mattijs's avatar
mattijs committed
633
            }
634

mattijs's avatar
mattijs committed
635 636 637 638 639 640 641
            meshMod.addFace
            (
                newF,                       // vertices
                own,                        // own
                nei,                        // nei
                -1,                         // masterPointID
                -1,                         // masterEdgeID
642
                faceToFaceMap.size(),       // masterFaceID : current facei
mattijs's avatar
mattijs committed
643
                false,                      // flipFaceFlux
644
                patchi,                     // patchID
mattijs's avatar
mattijs committed
645 646 647
                -1,                         // zoneID
                false                       // zoneFlip
            );
648
            faceToFaceMap.append(facei+1);  // unflipped
mattijs's avatar
mattijs committed
649 650 651 652 653 654 655
            faceToEdgeMap.append(-1);

            //Pout<< "Added inbetween face " << newF
            //    << " coords:" << UIndirectList<point>(meshMod.points(), newF)
            //    << " at layer " << layerI
            //    << " own " << own
            //    << " nei " << nei
656
            //    << "  at " << patch_.faceCentres()[facei]
mattijs's avatar
mattijs committed
657 658
            //    << endl;
        }
659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674
    }


    // Add side faces
    // ~~~~~~~~~~~~~~
    // Note that we loop over edges multiple times so for edges with
    // two cyclic faces they get added in two passes (for correct ordering)

    // Pass1. Internal edges and first face of other edges
    forAll(extrudeEdgePatches, edgeI)
    {
        const labelList& eFaces = patch_.edgeFaces()[edgeI];
        const labelList& ePatches = extrudeEdgePatches[edgeI];

        if (ePatches.size() == 0)
        {
mattijs's avatar
mattijs committed
675
            // Internal face
676 677
            if (eFaces.size() != 2)
            {
678
                FatalErrorInFunction
679 680 681 682 683 684 685 686 687
                    << "edge:" << edgeI
                    << " not internal but does not have side-patches defined."
                    << exit(FatalError);
            }
        }
        else
        {
            if (eFaces.size() != ePatches.size())
            {
688
                FatalErrorInFunction
689
                    << "external/feature edge:" << edgeI
690
                    << " has " << eFaces.size() << " connected extruded faces "
691 692 693
                    << " but only " << ePatches.size()
                    << " boundary faces defined." << exit(FatalError);
            }
mattijs's avatar
mattijs committed
694 695 696
        }


697

mattijs's avatar
mattijs committed
698 699 700
        // Make face pointing in to eFaces[0] so out of new master face
        const face& f = patch_.localFaces()[eFaces[0]];
        const edge& e = patch_.edges()[edgeI];
701

702
        label fp0 = f.find(e[0]);
mattijs's avatar
mattijs committed
703 704 705 706 707 708 709
        label fp1 = f.fcIndex(fp0);

        if (f[fp1] != e[1])
        {
            fp1 = fp0;
            fp0 = f.rcIndex(fp1);
        }
710

mattijs's avatar
mattijs committed
711
        face newF(4);
712

mattijs's avatar
mattijs committed
713 714 715 716
        for (label layerI = 0; layerI < nLayers; layerI++)
        {
            label region0 = pointRegions_[eFaces[0]][fp0];
            label region1 = pointRegions_[eFaces[0]][fp1];
717

718
            // Pick up points with correct normal
mattijs's avatar
mattijs committed
719
            if (layerI == 0)
720
            {
mattijs's avatar
mattijs committed
721 722 723 724 725 726 727 728 729 730 731
                newF[0] = f[fp0];
                newF[1] = f[fp1];
                newF[2] = addedPoints[nLayers*region1+layerI];
                newF[3] = addedPoints[nLayers*region0+layerI];
            }
            else
            {
                newF[0] = addedPoints[nLayers*region0+layerI-1];
                newF[1] = addedPoints[nLayers*region1+layerI-1];
                newF[2] = addedPoints[nLayers*region1+layerI];
                newF[3] = addedPoints[nLayers*region0+layerI];
732 733
            }

734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749
            // Optionally rotate so e[0] is always 0th vertex. Note that
            // this normally is automatically done by coupled face ordering
            // but with NOORDERING we have to do it ourselves.
            if (f[fp0] != e[0])
            {
                // rotate one back to get newF[1] (originating from e[0])
                // into newF[0]
                label v0 = newF[0];
                for (label i = 0; i < newF.size()-1; i++)
                {
                    newF[i] = newF[newF.fcIndex(i)];
                }
                newF.last() = v0;
            }


750 751
            label minCelli = addedCells[nLayers*eFaces[0]+layerI];
            label maxCelli;
752
            label patchi;
mattijs's avatar
mattijs committed
753 754
            if (ePatches.size() == 0)
            {
755 756
                maxCelli = addedCells[nLayers*eFaces[1]+layerI];
                if (minCelli > maxCelli)
mattijs's avatar
mattijs committed
757 758
                {
                    // Swap
759
                    Swap(minCelli, maxCelli);
mattijs's avatar
mattijs committed
760 761
                    newF = newF.reverseFace();
                }
762
                patchi = -1;
mattijs's avatar
mattijs committed
763 764 765
            }
            else
            {
766
                maxCelli = -1;
767
                patchi = ePatches[0];
mattijs's avatar
mattijs committed
768
            }
769

mattijs's avatar
mattijs committed
770 771 772 773 774 775 776 777 778
            //{
            //    Pout<< "Adding from face:" << patch_.faceCentres()[eFaces[0]]
            //        << " from edge:"
            //        << patch_.localPoints()[f[fp0]]
            //        << patch_.localPoints()[f[fp1]]
            //        << " at layer:" << layerI
            //        << " with new points:" << newF
            //        << " locations:"
            //        << UIndirectList<point>(meshMod.points(), newF)
779 780
            //        << " own:" << minCelli
            //        << " nei:" << maxCelli
mattijs's avatar
mattijs committed
781 782
            //        << endl;
            //}
783 784 785 786 787 788


            // newF already outwards pointing.
            meshMod.addFace
            (
                newF,                   // vertices
789 790
                minCelli,               // own
                maxCelli,               // nei
791 792 793 794
                -1,                     // masterPointID
                -1,                     // masterEdgeID
                faceToFaceMap.size(),   // masterFaceID
                false,                  // flipFaceFlux
795
                patchi,                 // patchID
796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814
                -1,                     // zoneID
                false                   // zoneFlip
            );
            faceToFaceMap.append(0);
            faceToEdgeMap.append(edgeI);
        }
    }

    // Pass2. Other faces of boundary edges
    forAll(extrudeEdgePatches, edgeI)
    {
        const labelList& eFaces = patch_.edgeFaces()[edgeI];
        const labelList& ePatches = extrudeEdgePatches[edgeI];

        if (ePatches.size() >= 2)
        {
            for (label i = 1; i < ePatches.size(); i++)
            {
                // Extrude eFaces[i]
815
                label minFacei = eFaces[i];
816 817

                // Make face pointing in to eFaces[0] so out of new master face
818
                const face& f = patch_.localFaces()[minFacei];
819 820

                const edge& e = patch_.edges()[edgeI];
821
                label fp0 = f.find(e[0]);
822 823 824 825 826 827 828 829 830
                label fp1 = f.fcIndex(fp0);

                if (f[fp1] != e[1])
                {
                    fp1 = fp0;
                    fp0 = f.rcIndex(fp1);
                }

                face newF(4);
mattijs's avatar
mattijs committed
831 832
                for (label layerI = 0; layerI < nLayers; layerI++)
                {
833 834
                    label region0 = pointRegions_[minFacei][fp0];
                    label region1 = pointRegions_[minFacei][fp1];
mattijs's avatar
mattijs committed
835 836 837 838 839 840 841 842 843 844 845 846 847 848 849 850

                    if (layerI == 0)
                    {
                        newF[0] = f[fp0];
                        newF[1] = f[fp1];
                        newF[2] = addedPoints[nLayers*region1+layerI];
                        newF[3] = addedPoints[nLayers*region0+layerI];
                    }
                    else
                    {
                        newF[0] = addedPoints[nLayers*region0+layerI-1];
                        newF[1] = addedPoints[nLayers*region1+layerI-1];
                        newF[2] = addedPoints[nLayers*region1+layerI];
                        newF[3] = addedPoints[nLayers*region0+layerI];
                    }

851 852 853 854 855 856 857 858 859 860 861 862 863 864 865

                    // Optionally rotate so e[0] is always 0th vertex. Note that
                    // this normally is automatically done by coupled face
                    // ordering but with NOORDERING we have to do it ourselves.
                    if (f[fp0] != e[0])
                    {
                        // rotate one back to get newF[1] (originating
                        // from e[0]) into newF[0].
                        label v0 = newF[0];
                        for (label i = 0; i < newF.size()-1; i++)
                        {
                            newF[i] = newF[newF.fcIndex(i)];
                        }
                        newF.last() = v0;
                    }
mattijs's avatar
mattijs committed
866 867 868
                    ////if (ePatches.size() == 0)
                    //{
                    //    Pout<< "Adding from MULTI face:"
869
                    //        << patch_.faceCentres()[minFacei]
mattijs's avatar
mattijs committed
870 871 872 873 874 875 876 877 878 879 880 881 882 883
                    //        << " from edge:"
                    //        << patch_.localPoints()[f[fp0]]
                    //        << patch_.localPoints()[f[fp1]]
                    //        << " at layer:" << layerI
                    //        << " with new points:" << newF
                    //        << " locations:"
                    //        << UIndirectList<point>(meshMod.points(), newF)
                    //        << endl;
                    //}

                    // newF already outwards pointing.
                    meshMod.addFace
                    (
                        newF,                   // vertices
884
                        addedCells[nLayers*minFacei+layerI],   // own
mattijs's avatar
mattijs committed
885 886 887 888 889 890 891 892 893 894 895 896
                        -1,                     // nei
                        -1,                     // masterPointID
                        -1,                     // masterEdgeID
                        faceToFaceMap.size(),   // masterFaceID
                        false,                  // flipFaceFlux
                        ePatches[i],            // patchID
                        -1,                     // zoneID
                        false                   // zoneFlip
                    );
                    faceToFaceMap.append(0);
                    faceToEdgeMap.append(edgeI);
                }
897 898 899 900 901 902 903 904 905 906 907 908 909 910 911 912 913 914 915 916 917 918
            }
        }
    }


    cellToFaceMap_.transfer(cellToFaceMap);
    faceToFaceMap_.transfer(faceToFaceMap);
    faceToEdgeMap_.transfer(faceToEdgeMap);
    pointToPointMap_.transfer(pointToPointMap);
}


void Foam::createShellMesh::updateMesh(const mapPolyMesh& map)
{
    inplaceReorder(map.reverseCellMap(), cellToFaceMap_);
    inplaceReorder(map.reverseFaceMap(), faceToFaceMap_);
    inplaceReorder(map.reverseFaceMap(), faceToEdgeMap_);
    inplaceReorder(map.reversePointMap(), pointToPointMap_);
}


// ************************************************************************* //