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