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