extrudeMesh.C 29.4 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
6
7
8
    \\  /    A nd           | Copyright (C) 2015-2019 OpenCFD Ltd.
     \\/     M anipulation  |
-------------------------------------------------------------------------------
                            | Copyright (C) 2011-2016 OpenFOAM Foundation
9
10
11
12
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

13
14
15
16
    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.
17
18
19
20
21
22
23

    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
24
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
25

26
27
28
Application
    extrudeMesh

29
30
31
Group
    grpMeshGenerationUtilities

32
Description
mattijs's avatar
mattijs committed
33
    Extrude mesh from existing patch (by default outwards facing normals;
34
35
    optional flips faces) or from patch read from file.

36
    Note: Merges close points so be careful.
37

38
    Type of extrusion prescribed by run-time selectable model.
39
40
41
42
43
44
45
46
47

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

#include "argList.H"
#include "Time.H"
#include "polyTopoChange.H"
#include "polyTopoChanger.H"
#include "edgeCollapser.H"
#include "perfectInterface.H"
48
49
50
#include "addPatchCellLayer.H"
#include "fvMesh.H"
#include "MeshedSurfaces.H"
51
#include "globalIndex.H"
52
#include "cellSet.H"
53
#include "fvMeshTools.H"
54
55

#include "extrudedMesh.H"
56
#include "extrudeModel.H"
57

58
59
#include "wedge.H"
#include "wedgePolyPatch.H"
60
#include "planeExtrusion.H"
61
#include "emptyPolyPatch.H"
62
#include "processorMeshes.H"
63
#include "hexRef8Data.H"
64

65
66
67
using namespace Foam;

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
68
69
70
71
72
73
74
75

enum ExtrudeMode
{
    MESH,
    PATCH,
    SURFACE
};

76
static const Enum<ExtrudeMode> ExtrudeModeNames
77
{
78
79
80
81
    { ExtrudeMode::MESH, "mesh" },
    { ExtrudeMode::PATCH, "patch" },
    { ExtrudeMode::SURFACE, "surface" },
};
82
83
84
85


label findPatchID(const polyBoundaryMesh& patches, const word& name)
{
86
    const label patchID = patches.findPatchID(name);
87
88
89

    if (patchID == -1)
    {
90
        FatalErrorInFunction
91
92
93
94
95
96
97
98
99
            << "Cannot find patch " << name
            << " in the source mesh.\n"
            << "Valid patch names are " << patches.names()
            << exit(FatalError);
    }
    return patchID;
}


mattijs's avatar
mattijs committed
100
labelList patchFaces(const polyBoundaryMesh& patches, const wordList& names)
101
{
mattijs's avatar
mattijs committed
102
    label n = 0;
103

mattijs's avatar
mattijs committed
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
    forAll(names, i)
    {
        const polyPatch& pp = patches[findPatchID(patches, names[i])];

        n += pp.size();
    }
    labelList faceLabels(n);
    n = 0;
    forAll(names, i)
    {
        const polyPatch& pp = patches[findPatchID(patches, names[i])];

        forAll(pp, j)
        {
            faceLabels[n++] = pp.start()+j;
        }
    }

    return faceLabels;
123
124
125
126
127
128
129
130
131
132
133
134
}


void updateFaceLabels(const mapPolyMesh& map, labelList& faceLabels)
{
    const labelList& reverseMap = map.reverseFaceMap();

    labelList newFaceLabels(faceLabels.size());
    label newI = 0;

    forAll(faceLabels, i)
    {
135
        label oldFacei = faceLabels[i];
136

137
        if (reverseMap[oldFacei] >= 0)
138
        {
139
            newFaceLabels[newI++] = reverseMap[oldFacei];
140
141
142
143
144
145
146
        }
    }
    newFaceLabels.setSize(newI);
    faceLabels.transfer(newFaceLabels);
}


147
148
149
150
151
152
153
154
void updateCellSet(const mapPolyMesh& map, labelHashSet& cellLabels)
{
    const labelList& reverseMap = map.reverseCellMap();

    labelHashSet newCellLabels(2*cellLabels.size());

    forAll(cellLabels, i)
    {
155
        label oldCelli = cellLabels[i];
156

157
        if (reverseMap[oldCelli] >= 0)
158
        {
159
            newCellLabels.insert(reverseMap[oldCelli]);
160
161
162
163
        }
    }
    cellLabels.transfer(newCellLabels);
}
164

165

166
167
168
169
170
171
172
173
174
175
template<class PatchType>
void changeFrontBackPatches
(
    polyMesh& mesh,
    const word& frontPatchName,
    const word& backPatchName
)
{
    const polyBoundaryMesh& patches = mesh.boundaryMesh();

176
177
    label frontPatchi = findPatchID(patches, frontPatchName);
    label backPatchi = findPatchID(patches, backPatchName);
178
179
180

    DynamicList<polyPatch*> newPatches(patches.size());

181
    forAll(patches, patchi)
182
    {
183
        const polyPatch& pp(patches[patchi]);
184

185
        if (patchi == frontPatchi || patchi == backPatchi)
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
        {
            newPatches.append
            (
                new PatchType
                (
                    pp.name(),
                    pp.size(),
                    pp.start(),
                    pp.index(),
                    mesh.boundaryMesh(),
                    PatchType::typeName
                )
            );
        }
        else
        {
            newPatches.append(pp.clone(mesh.boundaryMesh()).ptr());
        }
    }

    // Edit patches
    mesh.removeBoundary();
    mesh.addPatches(newPatches, true);
}


212
213
int main(int argc, char *argv[])
{
214
215
216
217
218
    argList::addNote
    (
        "Extrude mesh from existing patch."
    );

219
    #include "addRegionOption.H"
mattijs's avatar
mattijs committed
220
    argList::addOption("dict", "file", "Use alternative extrudeMeshDict");
mattijs's avatar
mattijs committed
221
    #include "setRootCase.H"
222
    #include "createTimeExtruded.H"
223

224
225
226
    // Get optional regionName
    word regionName;
    word regionDir;
227
    if (args.readIfPresent("region", regionName))
228
229
230
231
232
233
234
235
236
237
238
239
    {
        regionDir = regionName;
        Info<< "Create mesh " << regionName << " for time = "
            << runTimeExtruded.timeName() << nl << endl;
    }
    else
    {
        regionName = fvMesh::defaultRegion;
        Info<< "Create mesh for time = "
            << runTimeExtruded.timeName() << nl << endl;
    }

mattijs's avatar
mattijs committed
240
    const IOdictionary dict
241
    (
mattijs's avatar
mattijs committed
242
        IOobject::selectIO
243
        (
mattijs's avatar
mattijs committed
244
245
246
247
248
249
250
251
252
            IOobject
            (
                "extrudeMeshDict",
                runTimeExtruded.system(),
                runTimeExtruded,
                IOobject::MUST_READ_IF_MODIFIED,
                IOobject::NO_WRITE
            ),
            args.opt<fileName>("dict", "")
253
254
255
        )
    );

256
    // Point generator
mattijs's avatar
mattijs committed
257
258
    autoPtr<extrudeModel> model(extrudeModel::New(dict));

259
    // Whether to flip normals
Mark OLESEN's avatar
Mark OLESEN committed
260
    const bool flipNormals(dict.get<bool>("flipNormals"));
261
262

    // What to extrude
Mark OLESEN's avatar
Mark OLESEN committed
263
    const ExtrudeMode mode = ExtrudeModeNames.get("constructFrom", dict);
264

265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
    // Any merging of small edges
    const scalar mergeTol(dict.lookupOrDefault<scalar>("mergeTol", 1e-4));

    Info<< "Extruding from " << ExtrudeModeNames[mode]
        << " using model " << model().type() << endl;
    if (flipNormals)
    {
        Info<< "Flipping normals before extruding" << endl;
    }
    if (mergeTol > 0)
    {
        Info<< "Collapsing edges < " << mergeTol << " of bounding box" << endl;
    }
    else
    {
        Info<< "Not collapsing any edges after extrusion" << endl;
    }
    Info<< endl;

mattijs's avatar
mattijs committed
284

285
286
287
    // Generated mesh (one of either)
    autoPtr<fvMesh> meshFromMesh;
    autoPtr<polyMesh> meshFromSurface;
mattijs's avatar
mattijs committed
288

289
290
291
292
293
294
    // Faces on front and back for stitching (in case of mergeFaces)
    word frontPatchName;
    labelList frontPatchFaces;
    word backPatchName;
    labelList backPatchFaces;

295
296
297
    // Optional added cells (get written to cellSet)
    labelHashSet addedCellsSet;

298
299
300
    // Optional refinement data
    autoPtr<hexRef8Data> refDataPtr;

301
    if (mode == PATCH || mode == MESH)
302
    {
303
        if (flipNormals && mode == MESH)
mattijs's avatar
mattijs committed
304
        {
305
            FatalErrorInFunction
306
                << "Flipping normals not supported for extrusions from mesh."
mattijs's avatar
mattijs committed
307
308
309
                << exit(FatalError);
        }

mattijs's avatar
mattijs committed
310
        fileName sourceCasePath(dict.lookup("sourceCase"));
mattijs's avatar
mattijs committed
311
        sourceCasePath.expand();
312
313
        fileName sourceRootDir = sourceCasePath.path();
        fileName sourceCaseDir = sourceCasePath.name();
314
315
316
        if (Pstream::parRun())
        {
            sourceCaseDir =
317
                sourceCaseDir/("processor" + Foam::name(Pstream::myProcNo()));
318
        }
mattijs's avatar
mattijs committed
319
        wordList sourcePatches;
320
        dict.readEntry("sourcePatches", sourcePatches);
mattijs's avatar
mattijs committed
321
322
323
324
325

        if (sourcePatches.size() == 1)
        {
            frontPatchName = sourcePatches[0];
        }
326

mattijs's avatar
mattijs committed
327
        Info<< "Extruding patches " << sourcePatches
328
            << " on mesh " << sourceCasePath << nl
329
330
331
332
333
            << endl;

        Time runTime
        (
            Time::controlDictName,
334
335
            sourceRootDir,
            sourceCaseDir
336
        );
337

338
        #include "createNamedMesh.H"
339
340

        const polyBoundaryMesh& patches = mesh.boundaryMesh();
341
342


343
344
345
346
        // Extrusion engine. Either adding to existing mesh or
        // creating separate mesh.
        addPatchCellLayer layerExtrude(mesh, (mode == MESH));

347
        const labelList meshFaces(patchFaces(patches, sourcePatches));
348
349
350
351
352
353
354
355

        if (mode == PATCH && flipNormals)
        {
            // Cheat. Flip patch faces in mesh. This invalidates the
            // mesh (open cells) but does produce the correct extrusion.
            polyTopoChange meshMod(mesh);
            forAll(meshFaces, i)
            {
356
                label meshFacei = meshFaces[i];
357

358
359
                label patchi = patches.whichPatch(meshFacei);
                label own = mesh.faceOwner()[meshFacei];
360
                label nei = -1;
361
                if (patchi == -1)
362
                {
363
                    nei = mesh.faceNeighbour()[meshFacei];
364
365
                }

366
                label zoneI = mesh.faceZones().whichZone(meshFacei);
367
368
369
                bool zoneFlip = false;
                if (zoneI != -1)
                {
370
                    label index = mesh.faceZones()[zoneI].whichFace(meshFacei);
371
372
373
374
375
                    zoneFlip = mesh.faceZones()[zoneI].flipMap()[index];
                }

                meshMod.modifyFace
                (
376
377
                    mesh.faces()[meshFacei].reverseFace(),  // modified face
                    meshFacei,                      // label of face
378
379
380
                    own,                            // owner
                    nei,                            // neighbour
                    true,                           // face flip
381
                    patchi,                         // patch for face
382
383
384
385
386
387
388
389
390
                    zoneI,                          // zone for face
                    zoneFlip                        // face flip in zone
                );
            }

            // Change the mesh. No inflation.
            autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh, false);

            // Update fields
391
            mesh.updateMesh(map());
392
393
394
395
396
397
398
399
400

            // Move mesh (since morphing does not do this)
            if (map().hasMotionPoints())
            {
                mesh.movePoints(map().preMotionPoints());
            }
        }


401
402
403
404
405
        indirectPrimitivePatch extrudePatch
        (
            IndirectList<face>
            (
                mesh.faces(),
406
                meshFaces
407
408
409
410
            ),
            mesh.points()
        );

411
412
413
        // Determine extrudePatch normal
        pointField extrudePatchPointNormals
        (
414
            PatchTools::pointNormals(mesh, extrudePatch)
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
        );


        // Precalculate mesh edges for pp.edges.
        const labelList meshEdges
        (
            extrudePatch.meshEdges
            (
                mesh.edges(),
                mesh.pointEdges()
            )
        );

        // Global face indices engine
        const globalIndex globalFaces(mesh.nFaces());

        // Determine extrudePatch.edgeFaces in global numbering (so across
        // coupled patches)
        labelListList edgeGlobalFaces
        (
            addPatchCellLayer::globalEdgeFaces
            (
                mesh,
                globalFaces,
                extrudePatch
            )
        );


        // Determine what patches boundary edges need to get extruded into.
        // This might actually cause edge-connected processors to become
        // face-connected so might need to introduce new processor boundaries.
        // Calculates:
        //  - per pp.edge the patch to extrude into
        //  - any additional processor boundaries (patchToNbrProc = map from
        //    new patchID to neighbour processor)
        //  - number of new patches (nPatches)

453
454
455
456
        labelList edgePatchID;
        labelList edgeZoneID;
        boolList edgeFlip;
        labelList inflateFaceID;
457
458
459
        label nPatches;
        Map<label> nbrProcToPatch;
        Map<label> patchToNbrProc;
460
        addPatchCellLayer::calcExtrudeInfo
461
        (
462
463
            true,   // for internal edges get zone info from any face

464
465
466
467
468
            mesh,
            globalFaces,
            edgeGlobalFaces,
            extrudePatch,

469
            edgePatchID,
470
471
            nPatches,
            nbrProcToPatch,
472
473
474
475
476
            patchToNbrProc,

            edgeZoneID,
            edgeFlip,
            inflateFaceID
477
478
479
480
481
482
483
484
485
486
487
488
489
        );


        // Add any patches.

        label nAdded = nPatches - mesh.boundaryMesh().size();
        reduce(nAdded, sumOp<label>());

        Info<< "Adding overall " << nAdded << " processor patches." << endl;

        if (nAdded > 0)
        {
            DynamicList<polyPatch*> newPatches(nPatches);
490
            forAll(mesh.boundaryMesh(), patchi)
491
492
493
            {
                newPatches.append
                (
494
                    mesh.boundaryMesh()[patchi].clone
495
496
497
498
499
500
501
                    (
                        mesh.boundaryMesh()
                    ).ptr()
                );
            }
            for
            (
502
503
504
                label patchi = mesh.boundaryMesh().size();
                patchi < nPatches;
                patchi++
505
506
            )
            {
507
                label nbrProci = patchToNbrProc[patchi];
508

509
                Pout<< "Adding patch " << patchi
510
                    << " between " << Pstream::myProcNo()
511
                    << " and " << nbrProci
512
513
514
515
516
517
518
519
                    << endl;

                newPatches.append
                (
                    new processorPolyPatch
                    (
                        0,                  // size
                        mesh.nFaces(),      // start
520
                        patchi,             // index
521
522
                        mesh.boundaryMesh(),// polyBoundaryMesh
                        Pstream::myProcNo(),// myProcNo
523
                        nbrProci            // neighbProcNo
524
525
526
527
528
529
530
531
532
533
534
                    )
                );
            }

            // Add patches. Do no parallel updates.
            mesh.removeFvBoundary();
            mesh.addFvPatches(newPatches, true);
        }



535
        // Only used for addPatchCellLayer into new mesh
536
        labelList exposedPatchID;
537
        if (mode == PATCH)
538
        {
539
            dict.readEntry("exposedPatchName", backPatchName);
540
            exposedPatchID.setSize
541
542
543
544
            (
                extrudePatch.size(),
                findPatchID(patches, backPatchName)
            );
545
546
        }

547
        // Determine points and extrusion
548
549
        pointField layer0Points(extrudePatch.nPoints());
        pointField displacement(extrudePatch.nPoints());
550
        forAll(displacement, pointi)
551
        {
552
            const vector& patchNormal = extrudePatchPointNormals[pointi];
553
554

            // layer0 point
555
            layer0Points[pointi] = model()
556
            (
557
                extrudePatch.localPoints()[pointi],
558
559
560
561
562
563
                patchNormal,
                0
            );
            // layerN point
            point extrudePt = model()
            (
564
                extrudePatch.localPoints()[pointi],
565
566
567
                patchNormal,
                model().nLayers()
            );
568
            displacement[pointi] = extrudePt - layer0Points[pointi];
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
        }


        // Check if wedge (has layer0 different from original patch points)
        // If so move the mesh to starting position.
        if (gMax(mag(layer0Points-extrudePatch.localPoints())) > SMALL)
        {
            Info<< "Moving mesh to layer0 points since differ from original"
                << " points - this can happen for wedge extrusions." << nl
                << endl;

            pointField newPoints(mesh.points());
            forAll(extrudePatch.meshPoints(), i)
            {
                newPoints[extrudePatch.meshPoints()[i]] = layer0Points[i];
            }
            mesh.movePoints(newPoints);
        }


        // Layers per face
        labelList nFaceLayers(extrudePatch.size(), model().nLayers());
591

592
593
        // Layers per point
        labelList nPointLayers(extrudePatch.nPoints(), model().nLayers());
594

595
        // Displacement for first layer
596
        vectorField firstLayerDisp(displacement*model().sumThickness(1));
597

598
599
600
        // Expansion ratio not used.
        scalarField ratio(extrudePatch.nPoints(), 1.0);

601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626

        // Load any refinement data
        if (mode == MESH)
        {
            // Tricky: register hexRef8 onto the database
            // since the mesh does not yet exist. It should not be registered
            // onto the input mesh.
            refDataPtr.reset
            (
                new hexRef8Data
                (
                    IOobject
                    (
                        "dummy",
                        mesh.facesInstance(),
                        polyMesh::meshSubDir,
                        runTimeExtruded,        //mesh,
                        IOobject::READ_IF_PRESENT,
                        IOobject::NO_WRITE,
                        false
                    )
                )
            );
        }


627
628
629
630
631
632
633
634
635
636
637
        // Topo change container. Either copy an existing mesh or start
        // with empty storage (number of patches only needed for checking)
        autoPtr<polyTopoChange> meshMod
        (
            (
                mode == MESH
              ? new polyTopoChange(mesh)
              : new polyTopoChange(patches.size())
            )
        );

638
639
        layerExtrude.setRefinement
        (
640
641
642
            globalFaces,
            edgeGlobalFaces,

643
644
            ratio,              // expansion ratio
            extrudePatch,       // patch faces to extrude
645
646
647
648
649
650

            edgePatchID,        // if boundary edge: patch for extruded face
            edgeZoneID,         // optional zone for extruded face
            edgeFlip,
            inflateFaceID,      // mesh face that zone/patch info is from

651
            exposedPatchID,     // if new mesh: patches for exposed faces
652
653
654
655
656
657
658
            nFaceLayers,
            nPointLayers,
            firstLayerDisp,
            meshMod()
        );

        // Reset points according to extrusion model
659
        forAll(layerExtrude.addedPoints(), pointi)
660
        {
661
662
            const labelList& pPoints = layerExtrude.addedPoints()[pointi];
            forAll(pPoints, pPointi)
663
            {
664
                label meshPointi = pPoints[pPointi];
665
666
667
668
669

                point modelPt
                (
                    model()
                    (
670
671
672
                        extrudePatch.localPoints()[pointi],
                        extrudePatchPointNormals[pointi],
                        pPointi+1       // layer
673
674
675
676
677
678
                    )
                );

                const_cast<DynamicList<point>&>
                (
                    meshMod().points()
679
                )[meshPointi] = modelPt;
680
681
682
683
684
685
686
687
            }
        }

        // Store faces on front and exposed patch (if mode=patch there are
        // only added faces so cannot used map to old faces)
        const labelListList& layerFaces = layerExtrude.layerFaces();
        backPatchFaces.setSize(layerFaces.size());
        frontPatchFaces.setSize(layerFaces.size());
688
        forAll(backPatchFaces, patchFacei)
689
        {
690
691
            backPatchFaces[patchFacei]  = layerFaces[patchFacei].first();
            frontPatchFaces[patchFacei] = layerFaces[patchFacei].last();
692
693
        }

694

695
        // Create dummy fvSchemes, fvSolution
696
        fvMeshTools::createDummyFvMeshFiles(mesh, regionDir, true);
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713

        // Create actual mesh from polyTopoChange container
        autoPtr<mapPolyMesh> map = meshMod().makeMesh
        (
            meshFromMesh,
            IOobject
            (
                regionName,
                runTimeExtruded.constant(),
                runTimeExtruded,
                IOobject::NO_READ,
                IOobject::AUTO_WRITE,
                false
            ),
            mesh
        );

714
715
716
717
718
719
720
        layerExtrude.updateMesh
        (
            map(),
            identity(extrudePatch.size()),
            identity(extrudePatch.nPoints())
        );

721
722
723
724
725
726
727
728
729
730
731
        // Calculate face labels for front and back.
        frontPatchFaces = renumber
        (
            map().reverseFaceMap(),
            frontPatchFaces
        );
        backPatchFaces = renumber
        (
            map().reverseFaceMap(),
            backPatchFaces
        );
732

733
734
735
736
737
738
        // Update
        if (refDataPtr.valid())
        {
            refDataPtr().updateMesh(map());
        }

739
740
741
742
743
744
745
        // Store added cells
        if (mode == MESH)
        {
            const labelListList addedCells
            (
                layerExtrude.addedCells
                (
746
                    *meshFromMesh,
747
748
749
                    layerExtrude.layerFaces()
                )
            );
750
            forAll(addedCells, facei)
751
            {
752
                const labelList& aCells = addedCells[facei];
753
                addedCellsSet.insert(aCells);
754
755
            }
        }
756
    }
757
    else
758
759
    {
        // Read from surface
mattijs's avatar
mattijs committed
760
        fileName surfName(dict.lookup("surface"));
761
        surfName.expand();
762
763
764
765

        Info<< "Extruding surfaceMesh read from file " << surfName << nl
            << endl;

766
        MeshedSurface<face> fMesh(surfName);
767

768
769
770
        if (flipNormals)
        {
            Info<< "Flipping faces." << nl << endl;
771
            faceList& faces = const_cast<faceList&>(fMesh.surfFaces());
772
773
774
775
776
            forAll(faces, i)
            {
                faces[i] = fMesh[i].reverseFace();
            }
        }
777

778
779
780
781
782
783
        Info<< "Extruding surface with :" << nl
                << "    points     : " << fMesh.points().size() << nl
                << "    faces      : " << fMesh.size() << nl
                << "    normals[0] : " << fMesh.faceNormals()[0]
                << nl
                << endl;
mattijs's avatar
mattijs committed
784

785
786
787
788
789
790
791
792
793
794
795
796
797
798
        meshFromSurface.reset
        (
            new extrudedMesh
            (
                IOobject
                (
                    extrudedMesh::defaultRegion,
                    runTimeExtruded.constant(),
                    runTimeExtruded
                ),
                fMesh,
                model()
            )
        );
mattijs's avatar
mattijs committed
799
800


801
802
803
804
805
        // Get the faces on front and back
        frontPatchName = "originalPatch";
        frontPatchFaces = patchFaces
        (
            meshFromSurface().boundaryMesh(),
mattijs's avatar
mattijs committed
806
            wordList(1, frontPatchName)
807
808
809
810
811
        );
        backPatchName = "otherSide";
        backPatchFaces = patchFaces
        (
            meshFromSurface().boundaryMesh(),
mattijs's avatar
mattijs committed
812
            wordList(1, backPatchName)
813
        );
mattijs's avatar
mattijs committed
814
815
    }

816

817
    polyMesh& mesh =
mattijs's avatar
mattijs committed
818
    (
819
820
821
        meshFromMesh.valid()
      ? meshFromMesh()
      : meshFromSurface()
mattijs's avatar
mattijs committed
822
    );
823
824


825
    const boundBox& bb = mesh.bounds();
826
    const vector span = bb.span();
827
    const scalar mergeDim = mergeTol * bb.minDim();
828

mattijs's avatar
mattijs committed
829
830
831
    Info<< "Mesh bounding box : " << bb << nl
        << "        with span : " << span << nl
        << "Merge distance    : " << mergeDim << nl
832
833
834
835
836
837
        << endl;


    // Collapse edges
    // ~~~~~~~~~~~~~~

838
    if (mergeDim > 0)
839
    {
840
        Info<< "Collapsing edges < " << mergeDim << " ..." << nl << endl;
841
842
843
844
845
846
847

        // Edge collapsing engine
        edgeCollapser collapser(mesh);

        const edgeList& edges = mesh.edges();
        const pointField& points = mesh.points();

848
        bitSet collapseEdge(mesh.nEdges());
849
850
        Map<point> collapsePointToLocation(mesh.nPoints());

851
852
853
854
855
856
857
858
        forAll(edges, edgeI)
        {
            const edge& e = edges[edgeI];

            scalar d = e.mag(points);

            if (d < mergeDim)
            {
859
                Info<< "Merging edge " << e << " since length " << d
860
861
                    << " << " << mergeDim << nl;

862
                collapseEdge.set(edgeI);
863
                collapsePointToLocation.set(e[1], points[e[0]]);
864
865
866
            }
        }

867
868
        List<pointEdgeCollapse> allPointInfo;
        const globalIndex globalPoints(mesh.nPoints());
869
        labelList pointPriority(mesh.nPoints(), Zero);
870
871
872
873
874
875
876
877
878
879

        collapser.consistentCollapse
        (
            globalPoints,
            pointPriority,
            collapsePointToLocation,
            collapseEdge,
            allPointInfo
        );

880
881
        // Topo change container
        polyTopoChange meshMod(mesh);
882

883
        // Put all modifications into meshMod
884
        bool anyChange = collapser.setRefinement(allPointInfo, meshMod);
885
886
887
888
889
890
891

        if (anyChange)
        {
            // Construct new mesh from polyTopoChange.
            autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh, false);

            // Update fields
892
            mesh.updateMesh(map());
893

894
895
896
            // Update stored data
            updateFaceLabels(map(), frontPatchFaces);
            updateFaceLabels(map(), backPatchFaces);
897
            updateCellSet(map(), addedCellsSet);
898

899
900
901
902
903
            if (refDataPtr.valid())
            {
                refDataPtr().updateMesh(map());
            }

904
905
906
907
908
909
910
911
912
            // Move mesh (if inflation used)
            if (map().hasMotionPoints())
            {
                mesh.movePoints(map().preMotionPoints());
            }
        }
    }


913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
    // Change the front and back patch types as required
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    word frontBackType(word::null);

    if (isType<extrudeModels::wedge>(model()))
    {
        changeFrontBackPatches<wedgePolyPatch>
        (
            mesh,
            frontPatchName,
            backPatchName
        );
    }
    else if (isType<extrudeModels::plane>(model()))
    {
        changeFrontBackPatches<emptyPolyPatch>
        (
            mesh,
            frontPatchName,
            backPatchName
        );
    }


938
939
940
    // Merging front and back patch faces
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Mark OLESEN's avatar
Mark OLESEN committed
941
    if (dict.get<bool>("mergeFaces"))
942
    {
943
944
        if (mode == MESH)
        {
945
            FatalErrorInFunction
946
947
948
949
950
                << "Cannot stitch front and back of extrusion since"
                << " in 'mesh' mode (extrusion appended to mesh)."
                << exit(FatalError);
        }

951
952
        Info<< "Assuming full 360 degree axisymmetric case;"
            << " stitching faces on patches "
953
954
            << frontPatchName << " and "
            << backPatchName << " together ..." << nl << endl;
955

956
        if (frontPatchFaces.size() != backPatchFaces.size())
957
        {
958
            FatalErrorInFunction
959
960
961
962
                << "Differing number of faces on front ("
                << frontPatchFaces.size() << ") and back ("
                << backPatchFaces.size() << ")"
                << exit(FatalError);
963
964
        }

965
966
967
968

        polyTopoChanger stitcher(mesh);
        stitcher.setSize(1);

969
970
971
972
973
974
975
976
        const word cutZoneName("originalCutFaceZone");

        List<faceZone*> fz
        (
            1,
            new faceZone
            (
                cutZoneName,
977
                frontPatchFaces,
Mark OLESEN's avatar
Mark OLESEN committed
978
                false, // none are flipped
979
980
981
982
983
984
985
986
                0,
                mesh.faceZones()
            )
        );

        mesh.addZones(List<pointZone*>(0), fz, List<cellZone*>(0));

        // Add the perfect interface mesh modifier
987
        perfectInterface perfectStitcher
988
        (
989
            "couple",
990
            0,
991
992
993
994
995
996
997
998
999
1000
            stitcher,
            cutZoneName,
            word::null,         // dummy patch name
            word::null          // dummy patch name
        );

        // Topo change container
        polyTopoChange meshMod(mesh);

        perfectStitcher.setRefinement
For faster browsing, not all history is shown. View entire blame