applications/utilities/mesh/conversion/plot3dToFoam/plot3dToFoam.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     Plot3d mesh (ascii/formatted format) converter.
  26
  27     Work in progress! Handles ascii multiblock (and optionally singleBlock)
  28     format.
  29     By default expects blanking. Use -noBlank if none.
  30     Use -2D \a thickness if 2D.
  31
  32     Niklas Nordin has experienced a problem with lefthandedness of the blocks.
  33     The code should detect this automatically - see hexBlock::readPoints but
  34     if this goes wrong just set the blockHandedness_ variable to 'right'
  35     always.
  36
  37 \*---------------------------------------------------------------------------*/
  38
  39 #include "argList.H"
  40 #include "Time.H"
  41 #include "IFstream.H"
  42 #include "hexBlock.H"
  43 #include "polyMesh.H"
  44 #include "wallPolyPatch.H"
  45 #include "symmetryPolyPatch.H"
  46 #include "cellShape.H"
  47 #include "cellModeller.H"
  48 #include "mergePoints.H"
  49
  50 using namespace Foam;
  51
  52 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  53 // Main program:
  54
  55 int main(int argc, char *argv[])
  56 {
  57     argList::noParallel();
  58     argList::validArgs.append("PLOT3D geom file");
  59     argList::addOption
  60     (
  61         "scale",
  62         "factor",
  63         "geometry scaling factor - default is 1"
  64     );
  65     argList::addBoolOption
  66     (
  67         "noBlank",
  68         "skip blank items"
  69     );
  70     argList::addBoolOption
  71     (
  72         "singleBlock",
  73         "input is a single block"
  74     );
  75     argList::addOption
  76     (
  77         "2D",
  78         "thickness",
  79         "use when converting a 2-D geometry"
  80     );
  81
  82     argList args(argc, argv);
  83
  84     if (!args.check())
  85     {
  86          FatalError.exit();
  87     }
  88
  89     const scalar scaleFactor = args.optionLookupOrDefault("scale", 1.0);
  90
  91     bool readBlank = !args.optionFound("noBlank");
  92     bool singleBlock = args.optionFound("singleBlock");
  93     scalar twoDThickness = -1;
  94     if (args.optionReadIfPresent("2D", twoDThickness))
  95     {
  96         Info<< "Reading 2D case by extruding points by " << twoDThickness
  97             << " in z direction." << nl << endl;
  98     }
  99
 100
 101 #   include "createTime.H"
 102
 103     IFstream plot3dFile(args[1]);
 104
 105     // Read the plot3d information using a fixed format reader.
 106     // Comments in the file are in C++ style, so the stream parser will remove
 107     // them with no intervention
 108     label nblock;
 109
 110     if (singleBlock)
 111     {
 112         nblock = 1;
 113     }
 114     else
 115     {
 116         plot3dFile >> nblock;
 117     }
 118
 119     Info<< "Reading " << nblock << " blocks" << endl;
 120
 121     PtrList<hexBlock> blocks(nblock);
 122
 123     {
 124         label nx, ny, nz;
 125
 126         forAll(blocks, blockI)
 127         {
 128             if (twoDThickness > 0)
 129             {
 130                 // Fake second set of points (done in readPoints below)
 131                 plot3dFile >> nx >> ny;
 132                 nz = 2;
 133             }
 134             else
 135             {
 136                 plot3dFile >> nx >> ny >> nz;
 137             }
 138
 139             Info<< "block " << blockI << " nx:" << nx
 140                 << " ny:" << ny << " nz:" << nz << endl;
 141
 142             blocks.set(blockI, new hexBlock(nx, ny, nz));
 143         }
 144     }
 145
 146     Info<< "Reading block points" << endl;
 147     label sumPoints(0);
 148     label nMeshCells(0);
 149
 150     forAll(blocks, blockI)
 151     {
 152         Info<< "block " << blockI << ":" << nl;
 153         blocks[blockI].readPoints(readBlank, twoDThickness, plot3dFile);
 154         sumPoints += blocks[blockI].nBlockPoints();
 155         nMeshCells += blocks[blockI].nBlockCells();
 156         Info<< nl;
 157     }
 158
 159     pointField points(sumPoints);
 160     labelList blockOffsets(blocks.size());
 161     sumPoints = 0;
 162     forAll(blocks, blockI)
 163     {
 164         const pointField& blockPoints = blocks[blockI].points();
 165         blockOffsets[blockI] = sumPoints;
 166         forAll(blockPoints, i)
 167         {
 168             points[sumPoints++] = blockPoints[i];
 169         }
 170     }
 171
 172     // From old to new master point
 173     labelList oldToNew;
 174     pointField newPoints;
 175
 176     // Merge points
 177     mergePoints
 178     (
 179         points,
 180         SMALL,
 181         false,
 182         oldToNew,
 183         newPoints
 184     );
 185
 186     Info<< "Merged points within " << SMALL << " distance. Merged from "
 187         << oldToNew.size() << " down to " << newPoints.size()
 188         << " points." << endl;
 189
 190     // Scale the points
 191     if (scaleFactor > 1.0 + SMALL || scaleFactor < 1.0 - SMALL)
 192     {
 193         newPoints *= scaleFactor;
 194     }
 195
 196     Info<< "Creating cells" << endl;
 197
 198     cellShapeList cellShapes(nMeshCells);
 199
 200     const cellModel& hex = *(cellModeller::lookup("hex"));
 201
 202     label nCreatedCells = 0;
 203
 204     forAll(blocks, blockI)
 205     {
 206         labelListList curBlockCells = blocks[blockI].blockCells();
 207
 208         forAll(curBlockCells, blockCellI)
 209         {
 210             labelList cellPoints(curBlockCells[blockCellI].size());
 211
 212             forAll(cellPoints, pointI)
 213             {
 214                 cellPoints[pointI] =
 215                     oldToNew
 216                     [
 217                         curBlockCells[blockCellI][pointI]
 218                       + blockOffsets[blockI]
 219                     ];
 220             }
 221
 222             // Do automatic collapse from hex.
 223             cellShapes[nCreatedCells] = cellShape(hex, cellPoints, true);
 224
 225             nCreatedCells++;
 226         }
 227     }
 228
 229     Info<< "Creating boundary patches" << endl;
 230
 231     faceListList boundary(0);
 232     wordList patchNames(0);
 233     wordList patchTypes(0);
 234     word defaultFacesName = "defaultFaces";
 235     word defaultFacesType = wallPolyPatch::typeName;
 236     wordList patchPhysicalTypes(0);
 237
 238     polyMesh pShapeMesh
 239     (
 240         IOobject
 241         (
 242             polyMesh::defaultRegion,
 243             runTime.constant(),
 244             runTime
 245         ),
 246         xferMove(newPoints),
 247         cellShapes,
 248         boundary,
 249         patchNames,
 250         patchTypes,
 251         defaultFacesName,
 252         defaultFacesType,
 253         patchPhysicalTypes
 254     );
 255
 256     // Set the precision of the points data to 10
 257     IOstream::defaultPrecision(10);
 258
 259     Info<< "Writing polyMesh" << endl;
 260     pShapeMesh.write();
 261
 262     Info<< "End\n" << endl;
 263
 264     return 0;
 265 }
 266
 267
 268 // ************************************************************************* //