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

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of OpenFOAM.
  10
  11     OpenFOAM is free software; you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by the
  13     Free Software Foundation; either version 2 of the License, or (at your
  14     option) any later version.
  15
  16     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with OpenFOAM; if not, write to the Free Software Foundation,
  23     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
  24
  25 Description
  26     Extrude mesh from existing patch (by default outwards facing normals;
  27     optional flips faces) or from patch read from file.
  28
  29     Note: Merges close points so be careful.
  30
  31     Type of extrusion prescribed by run-time selectable model.
  32
  33 \*---------------------------------------------------------------------------*/
  34
  35 #include "argList.H"
  36 #include "Time.H"
  37 #include "dimensionedTypes.H"
  38 #include "IFstream.H"
  39 #include "faceMesh.H"
  40 #include "polyTopoChange.H"
  41 #include "polyTopoChanger.H"
  42 #include "edgeCollapser.H"
  43 #include "mathematicalConstants.H"
  44 #include "globalMeshData.H"
  45 #include "perfectInterface.H"
  46
  47 #include "extrudedMesh.H"
  48 #include "extrudeModel.H"
  49
  50 using namespace Foam;
  51
  52 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  53 // Main program:
  54
  55 int main(int argc, char *argv[])
  56 {
  57     #include "setRootCase.H"
  58     #include "createTimeExtruded.H"
  59
  60     autoPtr<extrudedMesh> meshPtr(NULL);
  61
  62     IOdictionary dict
  63     (
  64         IOobject
  65         (
  66             "extrudeProperties",
  67             runTimeExtruded.constant(),
  68             runTimeExtruded,
  69             IOobject::MUST_READ
  70         )
  71     );
  72
  73     autoPtr<extrudeModel> model(extrudeModel::New(dict));
  74
  75     const word sourceType(dict.lookup("constructFrom"));
  76
  77     autoPtr<faceMesh> fMesh;
  78
  79     if (sourceType == "patch")
  80     {
  81         fileName sourceCasePath(dict.lookup("sourceCase"));
  82         sourceCasePath.expand();
  83         fileName sourceRootDir = sourceCasePath.path();
  84         fileName sourceCaseDir = sourceCasePath.name();
  85         word patchName(dict.lookup("sourcePatch"));
  86
  87         Info<< "Extruding patch " << patchName
  88             << " on mesh " << sourceCasePath << nl
  89             << endl;
  90
  91         Time runTime
  92         (
  93             Time::controlDictName,
  94             sourceRootDir,
  95             sourceCaseDir
  96         );
  97         #include "createPolyMesh.H"
  98
  99         label patchID = mesh.boundaryMesh().findPatchID(patchName);
 100
 101         if (patchID == -1)
 102         {
 103             FatalErrorIn(args.executable())
 104                 << "Cannot find patch " << patchName
 105                 << " in the source mesh.\n"
 106                 << "Valid patch names are " << mesh.boundaryMesh().names()
 107                 << exit(FatalError);
 108         }
 109
 110         const polyPatch& pp = mesh.boundaryMesh()[patchID];
 111         fMesh.reset(new faceMesh(pp.localFaces(), pp.localPoints()));
 112
 113         {
 114             fileName surfName(runTime.path()/patchName + ".sMesh");
 115             Info<< "Writing patch as surfaceMesh to "
 116                 << surfName << nl << endl;
 117             OFstream os(surfName);
 118             os << fMesh() << nl;
 119         }
 120     }
 121     else if (sourceType == "surface")
 122     {
 123         // Read from surface
 124         fileName surfName(dict.lookup("surface"));
 125
 126         Info<< "Extruding surfaceMesh read from file " << surfName << nl
 127             << endl;
 128
 129         IFstream is(surfName);
 130
 131         fMesh.reset(new faceMesh(is));
 132
 133         Info<< "Read patch from file " << surfName << nl
 134             << endl;
 135     }
 136     else
 137     {
 138         FatalErrorIn(args.executable())
 139             << "Illegal 'constructFrom' specification. Should either be "
 140             << "patch or surface." << exit(FatalError);
 141     }
 142
 143     Switch flipNormals(dict.lookup("flipNormals"));
 144
 145     if (flipNormals)
 146     {
 147         Info<< "Flipping faces." << nl << endl;
 148
 149         faceList faces(fMesh().size());
 150         forAll(faces, i)
 151         {
 152             faces[i] = fMesh()[i].reverseFace();
 153         }
 154         fMesh.reset(new faceMesh(faces, fMesh().localPoints()));
 155     }
 156
 157
 158     Info<< "Extruding patch with :" << nl
 159             << "    points     : " << fMesh().points().size() << nl
 160             << "    faces      : " << fMesh().size() << nl
 161             << "    normals[0] : " << fMesh().faceNormals()[0]
 162             << nl
 163             << endl;
 164
 165     extrudedMesh mesh
 166     (
 167         IOobject
 168         (
 169             extrudedMesh::defaultRegion,
 170             runTimeExtruded.constant(),
 171             runTimeExtruded
 172         ),
 173         fMesh(),
 174         model()
 175     );
 176
 177
 178     const boundBox& bb = mesh.globalData().bb();
 179     const vector span = bb.span();
 180     const scalar mergeDim = 1E-4 * bb.minDim();
 181
 182     Info<< "Mesh bounding box : " << bb << nl
 183         << "        with span : " << span << nl
 184         << "Merge distance    : " << mergeDim << nl
 185         << endl;
 186
 187     const polyBoundaryMesh& patches = mesh.boundaryMesh();
 188
 189     const label origPatchID = patches.findPatchID("originalPatch");
 190     const label otherPatchID = patches.findPatchID("otherSide");
 191
 192     if (origPatchID == -1 || otherPatchID == -1)
 193     {
 194         FatalErrorIn(args.executable())
 195             << "Cannot find patch originalPatch or otherSide." << nl
 196             << "Valid patches are " << patches.names() << exit(FatalError);
 197     }
 198
 199     // Collapse edges
 200     // ~~~~~~~~~~~~~~
 201
 202     {
 203         Info<< "Collapsing edges < " << mergeDim << " ..." << nl << endl;
 204
 205         // Edge collapsing engine
 206         edgeCollapser collapser(mesh);
 207
 208         const edgeList& edges = mesh.edges();
 209         const pointField& points = mesh.points();
 210
 211         forAll(edges, edgeI)
 212         {
 213             const edge& e = edges[edgeI];
 214
 215             scalar d = e.mag(points);
 216
 217             if (d < mergeDim)
 218             {
 219                 Info<< "Merging edge " << e << " since length " << d
 220                     << " << " << mergeDim << nl;
 221
 222                 // Collapse edge to e[0]
 223                 collapser.collapseEdge(edgeI, e[0]);
 224             }
 225         }
 226
 227         // Topo change container
 228         polyTopoChange meshMod(mesh);
 229         // Put all modifications into meshMod
 230         bool anyChange = collapser.setRefinement(meshMod);
 231
 232         if (anyChange)
 233         {
 234             // Construct new mesh from polyTopoChange.
 235             autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh, false);
 236
 237             // Update fields
 238             mesh.updateMesh(map);
 239
 240             // Move mesh (if inflation used)
 241             if (map().hasMotionPoints())
 242             {
 243                 mesh.movePoints(map().preMotionPoints());
 244             }
 245         }
 246     }
 247
 248
 249     // Merging front and back patch faces
 250     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 251
 252     Switch mergeFaces(dict.lookup("mergeFaces"));
 253     if (mergeFaces)
 254     {
 255         Info<< "Assuming full 360 degree axisymmetric case;"
 256             << " stitching faces on patches "
 257             << patches[origPatchID].name() << " and "
 258             << patches[otherPatchID].name() << " together ..." << nl << endl;
 259
 260         polyTopoChanger stitcher(mesh);
 261         stitcher.setSize(1);
 262
 263         // Make list of masterPatch faces
 264         labelList isf(patches[origPatchID].size());
 265
 266         forAll (isf, i)
 267         {
 268             isf[i] = patches[origPatchID].start() + i;
 269         }
 270
 271         const word cutZoneName("originalCutFaceZone");
 272
 273         List<faceZone*> fz
 274         (
 275             1,
 276             new faceZone
 277             (
 278                 cutZoneName,
 279                 isf,
 280                 boolList(isf.size(), false),
 281                 0,
 282                 mesh.faceZones()
 283             )
 284         );
 285
 286         mesh.addZones(List<pointZone*>(0), fz, List<cellZone*>(0));
 287
 288         // Add the perfect interface mesh modifier
 289         stitcher.set
 290         (
 291             0,
 292             new perfectInterface
 293             (
 294                 "couple",
 295                 0,
 296                 stitcher,
 297                 cutZoneName,
 298                 patches[origPatchID].name(),
 299                 patches[otherPatchID].name()
 300             )
 301         );
 302
 303         // Execute all polyMeshModifiers
 304         autoPtr<mapPolyMesh> morphMap = stitcher.changeMesh(true);
 305
 306         mesh.movePoints(morphMap->preMotionPoints());
 307     }
 308
 309     if (!mesh.write())
 310     {
 311         FatalErrorIn(args.executable()) << "Failed writing mesh"
 312             << exit(FatalError);
 313     }
 314
 315     Info << "End\n" << endl;
 316
 317     return 0;
 318 }
 319
 320
 321 // ************************************************************************* //