applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMesh.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
   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
  13     the Free Software Foundation, either version 3 of the License, or
  14     (at your 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, see <http://www.gnu.org/licenses/>.
  23
  24 Description
  25     Extrude mesh from existing patch (by default outwards facing normals;
  26     optional flips faces) or from patch read from file.
  27
  28     Note: Merges close points so be careful.
  29
  30     Type of extrusion prescribed by run-time selectable model.
  31
  32 \*---------------------------------------------------------------------------*/
  33
  34 #include "argList.H"
  35 #include "Time.H"
  36 #include "dimensionedTypes.H"
  37 #include "IFstream.H"
  38 #include "polyTopoChange.H"
  39 #include "polyTopoChanger.H"
  40 #include "edgeCollapser.H"
  41 #include "globalMeshData.H"
  42 #include "perfectInterface.H"
  43 #include "addPatchCellLayer.H"
  44 #include "fvMesh.H"
  45 #include "MeshedSurfaces.H"
  46 #include "globalIndex.H"
  47
  48 #include "extrudedMesh.H"
  49 #include "extrudeModel.H"
  50
  51 using namespace Foam;
  52
  53 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  54
  55 enum ExtrudeMode
  56 {
  57     MESH,
  58     PATCH,
  59     SURFACE
  60 };
  61
  62 namespace Foam
  63 {
  64     template<>
  65     const char* NamedEnum<ExtrudeMode, 3>::names[] =
  66     {
  67         "mesh",
  68         "patch",
  69         "surface"
  70     };
  71 }
  72
  73 static const NamedEnum<ExtrudeMode, 3> ExtrudeModeNames;
  74
  75
  76 void createDummyFvMeshFiles(const polyMesh& mesh, const word& regionName)
  77 {
  78     // Create dummy system/fv*
  79     {
  80         IOobject io
  81         (
  82             "fvSchemes",
  83             mesh.time().system(),
  84             regionName,
  85             mesh,
  86             IOobject::NO_READ,
  87             IOobject::NO_WRITE,
  88             false
  89         );
  90
  91         Info<< "Testing:" << io.objectPath() << endl;
  92
  93         if (!io.headerOk())
  94         {
  95             Info<< "Writing dummy " << regionName/io.name() << endl;
  96             dictionary dummyDict;
  97             dictionary divDict;
  98             dummyDict.add("divSchemes", divDict);
  99             dictionary gradDict;
 100             dummyDict.add("gradSchemes", gradDict);
 101             dictionary laplDict;
 102             dummyDict.add("laplacianSchemes", laplDict);
 103
 104             IOdictionary(io, dummyDict).regIOobject::write();
 105         }
 106     }
 107     {
 108         IOobject io
 109         (
 110             "fvSolution",
 111             mesh.time().system(),
 112             regionName,
 113             mesh,
 114             IOobject::NO_READ,
 115             IOobject::NO_WRITE,
 116             false
 117         );
 118
 119         if (!io.headerOk())
 120         {
 121             Info<< "Writing dummy " << regionName/io.name() << endl;
 122             dictionary dummyDict;
 123             IOdictionary(io, dummyDict).regIOobject::write();
 124         }
 125     }
 126 }
 127
 128
 129 label findPatchID(const polyBoundaryMesh& patches, const word& name)
 130 {
 131     const label patchID = patches.findPatchID(name);
 132
 133     if (patchID == -1)
 134     {
 135         FatalErrorIn("findPatchID(const polyBoundaryMesh&, const word&)")
 136             << "Cannot find patch " << name
 137             << " in the source mesh.\n"
 138             << "Valid patch names are " << patches.names()
 139             << exit(FatalError);
 140     }
 141     return patchID;
 142 }
 143
 144
 145 labelList patchFaces(const polyBoundaryMesh& patches, const wordList& names)
 146 {
 147     label n = 0;
 148
 149     forAll(names, i)
 150     {
 151         const polyPatch& pp = patches[findPatchID(patches, names[i])];
 152
 153         n += pp.size();
 154     }
 155     labelList faceLabels(n);
 156     n = 0;
 157     forAll(names, i)
 158     {
 159         const polyPatch& pp = patches[findPatchID(patches, names[i])];
 160
 161         forAll(pp, j)
 162         {
 163             faceLabels[n++] = pp.start()+j;
 164         }
 165     }
 166
 167     return faceLabels;
 168 }
 169
 170
 171 void updateFaceLabels(const mapPolyMesh& map, labelList& faceLabels)
 172 {
 173     const labelList& reverseMap = map.reverseFaceMap();
 174
 175     labelList newFaceLabels(faceLabels.size());
 176     label newI = 0;
 177
 178     forAll(faceLabels, i)
 179     {
 180         label oldFaceI = faceLabels[i];
 181
 182         if (reverseMap[oldFaceI] >= 0)
 183         {
 184             newFaceLabels[newI++] = reverseMap[oldFaceI];
 185         }
 186     }
 187     newFaceLabels.setSize(newI);
 188     faceLabels.transfer(newFaceLabels);
 189 }
 190
 191
 192
 193 // Main program:
 194
 195 int main(int argc, char *argv[])
 196 {
 197     #include "addRegionOption.H"
 198     #include "setRootCase.H"
 199     #include "createTimeExtruded.H"
 200
 201     // Get optional regionName
 202     word regionName;
 203     word regionDir;
 204     if (args.optionReadIfPresent("region", regionName))
 205     {
 206         regionDir = regionName;
 207         Info<< "Create mesh " << regionName << " for time = "
 208             << runTimeExtruded.timeName() << nl << endl;
 209     }
 210     else
 211     {
 212         regionName = fvMesh::defaultRegion;
 213         Info<< "Create mesh for time = "
 214             << runTimeExtruded.timeName() << nl << endl;
 215     }
 216
 217
 218     IOdictionary dict
 219     (
 220         IOobject
 221         (
 222             "extrudeMeshDict",
 223             runTimeExtruded.system(),
 224             runTimeExtruded,
 225             IOobject::MUST_READ_IF_MODIFIED
 226         )
 227     );
 228
 229     // Point generator
 230     autoPtr<extrudeModel> model(extrudeModel::New(dict));
 231
 232     // Whether to flip normals
 233     const Switch flipNormals(dict.lookup("flipNormals"));
 234
 235     // What to extrude
 236     const ExtrudeMode mode = ExtrudeModeNames.read
 237     (
 238         dict.lookup("constructFrom")
 239     );
 240
 241
 242     // Generated mesh (one of either)
 243     autoPtr<fvMesh> meshFromMesh;
 244     autoPtr<polyMesh> meshFromSurface;
 245
 246     // Faces on front and back for stitching (in case of mergeFaces)
 247     word frontPatchName;
 248     labelList frontPatchFaces;
 249     word backPatchName;
 250     labelList backPatchFaces;
 251
 252     if (mode == PATCH || mode == MESH)
 253     {
 254         if (flipNormals)
 255         {
 256             FatalErrorIn(args.executable())
 257                 << "Flipping normals not supported for extrusions from patch."
 258                 << exit(FatalError);
 259         }
 260
 261         fileName sourceCasePath(dict.lookup("sourceCase"));
 262         sourceCasePath.expand();
 263         fileName sourceRootDir = sourceCasePath.path();
 264         fileName sourceCaseDir = sourceCasePath.name();
 265         if (Pstream::parRun())
 266         {
 267             sourceCaseDir =
 268                 sourceCaseDir
 269                /"processor" + Foam::name(Pstream::myProcNo());
 270         }
 271         wordList sourcePatches;
 272         dict.lookup("sourcePatches") >> sourcePatches;
 273
 274         if (sourcePatches.size() == 1)
 275         {
 276             frontPatchName = sourcePatches[0];
 277         }
 278
 279         Info<< "Extruding patches " << sourcePatches
 280             << " on mesh " << sourceCasePath << nl
 281             << endl;
 282
 283         Time runTime
 284         (
 285             Time::controlDictName,
 286             sourceRootDir,
 287             sourceCaseDir
 288         );
 289
 290         #include "createMesh.H"
 291
 292         const polyBoundaryMesh& patches = mesh.boundaryMesh();
 293
 294
 295         // Extrusion engine. Either adding to existing mesh or
 296         // creating separate mesh.
 297         addPatchCellLayer layerExtrude(mesh, (mode == MESH));
 298
 299         const labelList meshFaces(patchFaces(patches, sourcePatches));
 300         indirectPrimitivePatch extrudePatch
 301         (
 302             IndirectList<face>
 303             (
 304                 mesh.faces(),
 305                 meshFaces
 306             ),
 307             mesh.points()
 308         );
 309
 310         // Determine extrudePatch normal
 311         pointField extrudePatchPointNormals
 312         (
 313             PatchTools::pointNormals    //calcNormals
 314             (
 315                 mesh,
 316                 extrudePatch,
 317                 meshFaces
 318             )
 319         );
 320
 321
 322         // Precalculate mesh edges for pp.edges.
 323         const labelList meshEdges
 324         (
 325             extrudePatch.meshEdges
 326             (
 327                 mesh.edges(),
 328                 mesh.pointEdges()
 329             )
 330         );
 331
 332         // Global face indices engine
 333         const globalIndex globalFaces(mesh.nFaces());
 334
 335         // Determine extrudePatch.edgeFaces in global numbering (so across
 336         // coupled patches)
 337         labelListList edgeGlobalFaces
 338         (
 339             addPatchCellLayer::globalEdgeFaces
 340             (
 341                 mesh,
 342                 globalFaces,
 343                 extrudePatch
 344             )
 345         );
 346
 347
 348         // Determine what patches boundary edges need to get extruded into.
 349         // This might actually cause edge-connected processors to become
 350         // face-connected so might need to introduce new processor boundaries.
 351         // Calculates:
 352         //  - per pp.edge the patch to extrude into
 353         //  - any additional processor boundaries (patchToNbrProc = map from
 354         //    new patchID to neighbour processor)
 355         //  - number of new patches (nPatches)
 356
 357         labelList sidePatchID;
 358         label nPatches;
 359         Map<label> nbrProcToPatch;
 360         Map<label> patchToNbrProc;
 361         addPatchCellLayer::calcSidePatch
 362         (
 363             mesh,
 364             globalFaces,
 365             edgeGlobalFaces,
 366             extrudePatch,
 367
 368             sidePatchID,
 369             nPatches,
 370             nbrProcToPatch,
 371             patchToNbrProc
 372         );
 373
 374
 375         // Add any patches.
 376
 377         label nAdded = nPatches - mesh.boundaryMesh().size();
 378         reduce(nAdded, sumOp<label>());
 379
 380         Info<< "Adding overall " << nAdded << " processor patches." << endl;
 381
 382         if (nAdded > 0)
 383         {
 384             DynamicList<polyPatch*> newPatches(nPatches);
 385             forAll(mesh.boundaryMesh(), patchI)
 386             {
 387                 newPatches.append
 388                 (
 389                     mesh.boundaryMesh()[patchI].clone
 390                     (
 391                         mesh.boundaryMesh()
 392                     ).ptr()
 393                 );
 394             }
 395             for
 396             (
 397                 label patchI = mesh.boundaryMesh().size();
 398                 patchI < nPatches;
 399                 patchI++
 400             )
 401             {
 402                 label nbrProcI = patchToNbrProc[patchI];
 403                 word name =
 404                         "procBoundary"
 405                       + Foam::name(Pstream::myProcNo())
 406                       + "to"
 407                       + Foam::name(nbrProcI);
 408
 409                 Pout<< "Adding patch " << patchI
 410                     << " name:" << name
 411                     << " between " << Pstream::myProcNo()
 412                     << " and " << nbrProcI
 413                     << endl;
 414
 415
 416                 newPatches.append
 417                 (
 418                     new processorPolyPatch
 419                     (
 420                         name,
 421                         0,                  // size
 422                         mesh.nFaces(),      // start
 423                         patchI,             // index
 424                         mesh.boundaryMesh(),// polyBoundaryMesh
 425                         Pstream::myProcNo(),// myProcNo
 426                         nbrProcI            // neighbProcNo
 427                     )
 428                 );
 429             }
 430
 431             // Add patches. Do no parallel updates.
 432             mesh.removeFvBoundary();
 433             mesh.addFvPatches(newPatches, true);
 434         }
 435
 436
 437
 438         // Only used for addPatchCellLayer into new mesh
 439         labelList exposedPatchID;
 440         if (mode == PATCH)
 441         {
 442             dict.lookup("exposedPatchName") >> backPatchName;
 443             exposedPatchID.setSize
 444             (
 445                 extrudePatch.size(),
 446                 findPatchID(patches, backPatchName)
 447             );
 448         }
 449
 450         // Determine points and extrusion
 451         pointField layer0Points(extrudePatch.nPoints());
 452         pointField displacement(extrudePatch.nPoints());
 453         forAll(displacement, pointI)
 454         {
 455             const vector& patchNormal = extrudePatchPointNormals[pointI];
 456
 457             // layer0 point
 458             layer0Points[pointI] = model()
 459             (
 460                 extrudePatch.localPoints()[pointI],
 461                 patchNormal,
 462                 0
 463             );
 464             // layerN point
 465             point extrudePt = model()
 466             (
 467                 extrudePatch.localPoints()[pointI],
 468                 patchNormal,
 469                 model().nLayers()
 470             );
 471             displacement[pointI] = extrudePt - layer0Points[pointI];
 472         }
 473
 474         if (flipNormals)
 475         {
 476             Info<< "Flipping faces." << nl << endl;
 477             displacement = -displacement;
 478         }
 479
 480         // Check if wedge (has layer0 different from original patch points)
 481         // If so move the mesh to starting position.
 482         if (gMax(mag(layer0Points-extrudePatch.localPoints())) > SMALL)
 483         {
 484             Info<< "Moving mesh to layer0 points since differ from original"
 485                 << " points - this can happen for wedge extrusions." << nl
 486                 << endl;
 487
 488             pointField newPoints(mesh.points());
 489             forAll(extrudePatch.meshPoints(), i)
 490             {
 491                 newPoints[extrudePatch.meshPoints()[i]] = layer0Points[i];
 492             }
 493             mesh.movePoints(newPoints);
 494         }
 495
 496
 497         // Layers per face
 498         labelList nFaceLayers(extrudePatch.size(), model().nLayers());
 499         // Layers per point
 500         labelList nPointLayers(extrudePatch.nPoints(), model().nLayers());
 501         // Displacement for first layer
 502         vectorField firstLayerDisp(displacement*model().sumThickness(1));
 503
 504         // Expansion ratio not used.
 505         scalarField ratio(extrudePatch.nPoints(), 1.0);
 506
 507         // Topo change container. Either copy an existing mesh or start
 508         // with empty storage (number of patches only needed for checking)
 509         autoPtr<polyTopoChange> meshMod
 510         (
 511             (
 512                 mode == MESH
 513               ? new polyTopoChange(mesh)
 514               : new polyTopoChange(patches.size())
 515             )
 516         );
 517
 518         layerExtrude.setRefinement
 519         (
 520             globalFaces,
 521             edgeGlobalFaces,
 522
 523             ratio,              // expansion ratio
 524             extrudePatch,       // patch faces to extrude
 525             sidePatchID,        // if boundary edge: patch to use
 526             exposedPatchID,     // if new mesh: patches for exposed faces
 527             nFaceLayers,
 528             nPointLayers,
 529             firstLayerDisp,
 530             meshMod()
 531         );
 532
 533         // Reset points according to extrusion model
 534         forAll(layerExtrude.addedPoints(), pointI)
 535         {
 536             const labelList& pPoints = layerExtrude.addedPoints()[pointI];
 537             forAll(pPoints, pPointI)
 538             {
 539                 label meshPointI = pPoints[pPointI];
 540
 541                 point modelPt
 542                 (
 543                     model()
 544                     (
 545                         extrudePatch.localPoints()[pointI],
 546                         extrudePatchPointNormals[pointI],
 547                         pPointI+1       // layer
 548                     )
 549                 );
 550
 551                 const_cast<DynamicList<point>&>
 552                 (
 553                     meshMod().points()
 554                 )[meshPointI] = modelPt;
 555             }
 556         }
 557
 558         // Store faces on front and exposed patch (if mode=patch there are
 559         // only added faces so cannot used map to old faces)
 560         const labelListList& layerFaces = layerExtrude.layerFaces();
 561         backPatchFaces.setSize(layerFaces.size());
 562         frontPatchFaces.setSize(layerFaces.size());
 563         forAll(backPatchFaces, i)
 564         {
 565             backPatchFaces[i]  = layerFaces[i].first();
 566             frontPatchFaces[i] = layerFaces[i].last();
 567         }
 568
 569         // Create dummy fvSchemes, fvSolution
 570         createDummyFvMeshFiles(mesh, regionDir);
 571
 572         // Create actual mesh from polyTopoChange container
 573         autoPtr<mapPolyMesh> map = meshMod().makeMesh
 574         (
 575             meshFromMesh,
 576             IOobject
 577             (
 578                 regionName,
 579                 runTimeExtruded.constant(),
 580                 runTimeExtruded,
 581                 IOobject::NO_READ,
 582                 IOobject::AUTO_WRITE,
 583                 false
 584             ),
 585             mesh
 586         );
 587
 588         // Calculate face labels for front and back.
 589         frontPatchFaces = renumber
 590         (
 591             map().reverseFaceMap(),
 592             frontPatchFaces
 593         );
 594         backPatchFaces = renumber
 595         (
 596             map().reverseFaceMap(),
 597             backPatchFaces
 598         );
 599     }
 600     else
 601     {
 602         // Read from surface
 603         fileName surfName(dict.lookup("surface"));
 604
 605         Info<< "Extruding surfaceMesh read from file " << surfName << nl
 606             << endl;
 607
 608         MeshedSurface<face> fMesh(surfName);
 609
 610         if (flipNormals)
 611         {
 612             Info<< "Flipping faces." << nl << endl;
 613             faceList& faces = const_cast<faceList&>(fMesh.faces());
 614             forAll(faces, i)
 615             {
 616                 faces[i] = fMesh[i].reverseFace();
 617             }
 618         }
 619
 620         Info<< "Extruding surface with :" << nl
 621                 << "    points     : " << fMesh.points().size() << nl
 622                 << "    faces      : " << fMesh.size() << nl
 623                 << "    normals[0] : " << fMesh.faceNormals()[0]
 624                 << nl
 625                 << endl;
 626
 627         meshFromSurface.reset
 628         (
 629             new extrudedMesh
 630             (
 631                 IOobject
 632                 (
 633                     extrudedMesh::defaultRegion,
 634                     runTimeExtruded.constant(),
 635                     runTimeExtruded
 636                 ),
 637                 fMesh,
 638                 model()
 639             )
 640         );
 641
 642
 643         // Get the faces on front and back
 644         frontPatchName = "originalPatch";
 645         frontPatchFaces = patchFaces
 646         (
 647             meshFromSurface().boundaryMesh(),
 648             wordList(1, frontPatchName)
 649         );
 650         backPatchName = "otherSide";
 651         backPatchFaces = patchFaces
 652         (
 653             meshFromSurface().boundaryMesh(),
 654             wordList(1, backPatchName)
 655         );
 656     }
 657
 658
 659     polyMesh& mesh =
 660     (
 661         meshFromMesh.valid()
 662       ? meshFromMesh()
 663       : meshFromSurface()
 664     );
 665
 666
 667     const boundBox& bb = mesh.bounds();
 668     const vector span = bb.span();
 669     const scalar mergeDim = 1E-4 * bb.minDim();
 670
 671     Info<< "Mesh bounding box : " << bb << nl
 672         << "        with span : " << span << nl
 673         << "Merge distance    : " << mergeDim << nl
 674         << endl;
 675
 676
 677     // Collapse edges
 678     // ~~~~~~~~~~~~~~
 679
 680     {
 681         Info<< "Collapsing edges < " << mergeDim << " ..." << nl << endl;
 682
 683         // Edge collapsing engine
 684         edgeCollapser collapser(mesh);
 685
 686         const edgeList& edges = mesh.edges();
 687         const pointField& points = mesh.points();
 688
 689         forAll(edges, edgeI)
 690         {
 691             const edge& e = edges[edgeI];
 692
 693             scalar d = e.mag(points);
 694
 695             if (d < mergeDim)
 696             {
 697                 Info<< "Merging edge " << e << " since length " << d
 698                     << " << " << mergeDim << nl;
 699
 700                 // Collapse edge to e[0]
 701                 collapser.collapseEdge(edgeI, e[0]);
 702             }
 703         }
 704
 705         // Topo change container
 706         polyTopoChange meshMod(mesh);
 707         // Put all modifications into meshMod
 708         bool anyChange = collapser.setRefinement(meshMod);
 709
 710         if (anyChange)
 711         {
 712             // Construct new mesh from polyTopoChange.
 713             autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh, false);
 714
 715             // Update fields
 716             mesh.updateMesh(map);
 717
 718             // Update stored data
 719             updateFaceLabels(map(), frontPatchFaces);
 720             updateFaceLabels(map(), backPatchFaces);
 721
 722             // Move mesh (if inflation used)
 723             if (map().hasMotionPoints())
 724             {
 725                 mesh.movePoints(map().preMotionPoints());
 726             }
 727         }
 728     }
 729
 730
 731     // Merging front and back patch faces
 732     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 733
 734     Switch mergeFaces(dict.lookup("mergeFaces"));
 735     if (mergeFaces)
 736     {
 737         if (mode == MESH)
 738         {
 739             FatalErrorIn(args.executable())
 740                 << "Cannot stitch front and back of extrusion since"
 741                 << " in 'mesh' mode (extrusion appended to mesh)."
 742                 << exit(FatalError);
 743         }
 744
 745         Info<< "Assuming full 360 degree axisymmetric case;"
 746             << " stitching faces on patches "
 747             << frontPatchName << " and "
 748             << backPatchName << " together ..." << nl << endl;
 749
 750         if (frontPatchFaces.size() != backPatchFaces.size())
 751         {
 752             FatalErrorIn(args.executable())
 753                 << "Differing number of faces on front ("
 754                 << frontPatchFaces.size() << ") and back ("
 755                 << backPatchFaces.size() << ")"
 756                 << exit(FatalError);
 757         }
 758
 759
 760
 761         polyTopoChanger stitcher(mesh);
 762         stitcher.setSize(1);
 763
 764         const word cutZoneName("originalCutFaceZone");
 765
 766         List<faceZone*> fz
 767         (
 768             1,
 769             new faceZone
 770             (
 771                 cutZoneName,
 772                 frontPatchFaces,
 773                 boolList(frontPatchFaces.size(), false),
 774                 0,
 775                 mesh.faceZones()
 776             )
 777         );
 778
 779         mesh.addZones(List<pointZone*>(0), fz, List<cellZone*>(0));
 780
 781         // Add the perfect interface mesh modifier
 782         perfectInterface perfectStitcher
 783         (
 784             "couple",
 785             0,
 786             stitcher,
 787             cutZoneName,
 788             word::null,         // dummy patch name
 789             word::null          // dummy patch name
 790         );
 791
 792         // Topo change container
 793         polyTopoChange meshMod(mesh);
 794
 795         perfectStitcher.setRefinement
 796         (
 797             indirectPrimitivePatch
 798             (
 799                 IndirectList<face>
 800                 (
 801                     mesh.faces(),
 802                     frontPatchFaces
 803                 ),
 804                 mesh.points()
 805             ),
 806             indirectPrimitivePatch
 807             (
 808                 IndirectList<face>
 809                 (
 810                     mesh.faces(),
 811                     backPatchFaces
 812                 ),
 813                 mesh.points()
 814             ),
 815             meshMod
 816         );
 817
 818         // Construct new mesh from polyTopoChange.
 819         autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh, false);
 820
 821         // Update fields
 822         mesh.updateMesh(map);
 823
 824         // Move mesh (if inflation used)
 825         if (map().hasMotionPoints())
 826         {
 827             mesh.movePoints(map().preMotionPoints());
 828         }
 829     }
 830
 831     mesh.setInstance(runTimeExtruded.constant());
 832     Info<< "Writing mesh to " << mesh.objectPath() << nl << endl;
 833
 834     if (!mesh.write())
 835     {
 836         FatalErrorIn(args.executable()) << "Failed writing mesh"
 837             << exit(FatalError);
 838     }
 839
 840     Info<< "End\n" << endl;
 841
 842     return 0;
 843 }
 844
 845
 846 // ************************************************************************* //