STARCDMeshWriter.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2018 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 \*---------------------------------------------------------------------------*/
25 
26 #include "STARCDMeshWriter.H"
27 
28 #include "Time.H"
29 #include "SortableList.H"
30 #include "OFstream.H"
31 
32 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
33 
34 const char* Foam::meshWriters::STARCD::defaultBoundaryName =
35  "Default_Boundary_Region";
36 
38 {
39  { 4, 5, 2, 3, 0, 1 }, // 11 = pro-STAR hex
40  { 0, 1, 4, 5, 2, -1 }, // 12 = pro-STAR prism
41  { 5, 4, 2, 0, -1, -1 }, // 13 = pro-STAR tetra
42  { 0, 4, 3, 5, 2, -1 } // 14 = pro-STAR pyramid
43 };
44 
45 
46 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
47 
48 Foam::label Foam::meshWriters::STARCD::findDefaultBoundary() const
49 {
50  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
51 
52  label id = -1;
53 
54  // find Default_Boundary_Region if it exists
55  forAll(patches, patchi)
56  {
57  if (defaultBoundaryName == patches[patchi].name())
58  {
59  id = patchi;
60  break;
61  }
62  }
63  return id;
64 }
65 
66 
67 void Foam::meshWriters::STARCD::getCellTable()
68 {
69  // read constant/polyMesh/propertyName
70  IOList<label> ioList
71  (
72  IOobject
73  (
74  "cellTableId",
75  mesh_.time().constant(),
77  mesh_,
80  false
81  )
82  );
83 
84  bool useCellZones = false;
86 
87  // get information from constant/polyMesh/cellTableId if possible
88  if (ioList.headerOk())
89  {
90  if (ioList.size() == mesh_.nCells())
91  {
92  cellTableId_.transfer(ioList);
93 
94  if (cellTable_.empty())
95  {
96  Info<< "no cellTable information available" << endl;
97  }
98  }
99  else
100  {
102  << ioList.objectPath() << " has incorrect number of cells "
103  << " - use cellZone information"
104  << endl;
105 
106  ioList.clear();
107  useCellZones = true;
108  }
109  }
110  else
111  {
112  useCellZones = true;
113  }
114 
115 
116  if (useCellZones)
117  {
118  if (cellTable_.empty())
119  {
120  Info<< "created cellTable from cellZones" << endl;
121  cellTable_ = mesh_;
122  }
123 
124  // track if there are unzoned cells
125  label nUnzoned = mesh_.nCells();
126 
127  // get the cellZone <-> cellTable correspondence
128  Info<< "matching cellZones to cellTable" << endl;
129 
130  forAll(mesh_.cellZones(), zoneI)
131  {
132  const cellZone& cZone = mesh_.cellZones()[zoneI];
133  if (cZone.size())
134  {
135  nUnzoned -= cZone.size();
136 
137  label tableId = cellTable_.findIndex(cZone.name());
138  if (tableId < 0)
139  {
140  dictionary dict;
141 
142  dict.add("Label", cZone.name());
143  dict.add("MaterialType", "fluid");
144  tableId = cellTable_.append(dict);
145  }
146 
147  forAll(cZone, i)
148  {
149  cellTableId_[cZone[i]] = tableId;
150  }
151  }
152  }
153 
154  if (nUnzoned)
155  {
156  dictionary dict;
157 
158  dict.add("Label", "__unZonedCells__");
159  dict.add("MaterialType", "fluid");
160  label tableId = cellTable_.append(dict);
161 
162  forAll(cellTableId_, i)
163  {
164  if (cellTableId_[i] < 0)
165  {
166  cellTableId_[i] = tableId;
167  }
168  }
169  }
170  }
171 }
172 
173 
174 void Foam::meshWriters::STARCD::writeHeader(Ostream& os, const char* filetype)
175 {
176  os << "PROSTAR_" << filetype << nl
177  << 4000
178  << " " << 0
179  << " " << 0
180  << " " << 0
181  << " " << 0
182  << " " << 0
183  << " " << 0
184  << " " << 0
185  << endl;
186 }
187 
188 
189 void Foam::meshWriters::STARCD::writePoints(const fileName& prefix) const
190 {
191  OFstream os(prefix + ".vrt");
192  writeHeader(os, "VERTEX");
193 
194  // Set the precision of the points data to 10
195  os.precision(10);
196 
197  // force decimal point for Fortran input
198  os.setf(std::ios::showpoint);
199 
200  const pointField& points = mesh_.points();
201 
202  Info<< "Writing " << os.name() << " : "
203  << points.size() << " points" << endl;
204 
205  forAll(points, ptI)
206  {
207  // convert [m] -> [mm]
208  os
209  << ptI + 1 << " "
210  << scaleFactor_ * points[ptI].x() << " "
211  << scaleFactor_ * points[ptI].y() << " "
212  << scaleFactor_ * points[ptI].z() << nl;
213  }
214  os.flush();
215 
216 }
217 
218 
219 void Foam::meshWriters::STARCD::writeCells(const fileName& prefix) const
220 {
221  OFstream os(prefix + ".cel");
222  writeHeader(os, "CELL");
223 
224  // this is what we seem to need
225  // map foam cellModeller index -> star shape
226  Map<label> shapeLookupIndex;
227  shapeLookupIndex.insert(hexModel->index(), 11);
228  shapeLookupIndex.insert(prismModel->index(), 12);
229  shapeLookupIndex.insert(tetModel->index(), 13);
230  shapeLookupIndex.insert(pyrModel->index(), 14);
231 
232  const cellShapeList& shapes = mesh_.cellShapes();
233  const cellList& cells = mesh_.cells();
234  const faceList& faces = mesh_.faces();
235  const labelList& owner = mesh_.faceOwner();
236 
237  Info<< "Writing " << os.name() << " : "
238  << cells.size() << " cells" << endl;
239 
240  forAll(cells, cellId)
241  {
242  label tableId = cellTableId_[cellId];
243  label materialType = 1; // 1(fluid)
244  if (cellTable_.found(tableId))
245  {
246  const dictionary& dict = cellTable_[tableId];
247  if (dict.found("MaterialType"))
248  {
249  word matType;
250  dict.lookup("MaterialType") >> matType;
251  if (matType == "solid")
252  {
253  materialType = 2;
254  }
255 
256  }
257  }
258 
259  const cellShape& shape = shapes[cellId];
260  label mapIndex = shape.model().index();
261 
262  // a registered primitive type
263  if (shapeLookupIndex.found(mapIndex))
264  {
265  label shapeId = shapeLookupIndex[mapIndex];
266  const labelList& vrtList = shapes[cellId];
267 
268  os << cellId + 1
269  << " " << shapeId
270  << " " << vrtList.size()
271  << " " << tableId
272  << " " << materialType;
273 
274  // primitives have <= 8 vertices, but prevent overrun anyhow
275  // indent following lines for ease of reading
276  label count = 0;
277  forAll(vrtList, i)
278  {
279  if ((count % 8) == 0)
280  {
281  os << nl
282  << " " << cellId + 1;
283  }
284  os << " " << vrtList[i] + 1;
285  count++;
286  }
287  os << endl;
288 
289  }
290  else
291  {
292  label shapeId = 255; // treat as general polyhedral
293  const labelList& cFaces = cells[cellId];
294 
295  // create (beg,end) indices
296  List<label> indices(cFaces.size() + 1);
297  indices[0] = indices.size();
298 
299  label count = indices.size();
300  // determine the total number of vertices
301  forAll(cFaces, facei)
302  {
303  count += faces[cFaces[facei]].size();
304  indices[facei+1] = count;
305  }
306 
307  os << cellId + 1
308  << " " << shapeId
309  << " " << count
310  << " " << tableId
311  << " " << materialType;
312 
313  // write indices - max 8 per line
314  // indent following lines for ease of reading
315  count = 0;
316  forAll(indices, i)
317  {
318  if ((count % 8) == 0)
319  {
320  os << nl
321  << " " << cellId + 1;
322  }
323  os << " " << indices[i];
324  count++;
325  }
326 
327  // write faces - max 8 per line
328  forAll(cFaces, facei)
329  {
330  label meshFace = cFaces[facei];
331  face f;
332 
333  if (owner[meshFace] == cellId)
334  {
335  f = faces[meshFace];
336  }
337  else
338  {
339  f = faces[meshFace].reverseFace();
340  }
341 
342  forAll(f, i)
343  {
344  if ((count % 8) == 0)
345  {
346  os << nl
347  << " " << cellId + 1;
348  }
349 
350  os << " " << f[i] + 1;
351  count++;
352  }
353  }
354 
355  os << endl;
356  }
357  }
358 }
359 
360 
361 void Foam::meshWriters::STARCD::writeBoundary(const fileName& prefix) const
362 {
363  OFstream os(prefix + ".bnd");
364  writeHeader(os, "BOUNDARY");
365 
366  const cellShapeList& shapes = mesh_.cellShapes();
367  const cellList& cells = mesh_.cells();
368  const faceList& faces = mesh_.faces();
369  const labelList& owner = mesh_.faceOwner();
370  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
371 
372  // this is what we seem to need
373  // these MUST correspond to foamToStarFaceAddr
374  //
375  Map<label> faceLookupIndex;
376  faceLookupIndex.insert(hexModel->index(), 0);
377  faceLookupIndex.insert(prismModel->index(), 1);
378  faceLookupIndex.insert(tetModel->index(), 2);
379  faceLookupIndex.insert(pyrModel->index(), 3);
380 
381  Info<< "Writing " << os.name() << " : "
382  << (mesh_.nFaces() - patches[0].start()) << " boundaries" << endl;
383 
384 
385  label defaultId = findDefaultBoundary();
386 
387  //
388  // write boundary faces - skip Default_Boundary_Region entirely
389  //
390  label boundId = 0;
391  forAll(patches, patchi)
392  {
393  label regionId = patchi;
394  if (regionId == defaultId)
395  {
396  continue; // skip - already written
397  }
398  else if (defaultId == -1 || regionId < defaultId)
399  {
400  regionId++;
401  }
402 
403  label patchStart = patches[patchi].start();
404  label patchSize = patches[patchi].size();
405  word bndType = boundaryRegion_.boundaryType(patches[patchi].name());
406 
407  for
408  (
409  label facei = patchStart;
410  facei < (patchStart + patchSize);
411  ++facei
412  )
413  {
414  label cellId = owner[facei];
415  const labelList& cFaces = cells[cellId];
416  const cellShape& shape = shapes[cellId];
417  label cellFaceId = findIndex(cFaces, facei);
418 
419  // Info<< "cell " << cellId + 1 << " face " << facei
420  // << " == " << faces[facei]
421  // << " is index " << cellFaceId << " from " << cFaces;
422 
423  // Unfortunately, the order of faces returned by
424  // primitiveMesh::cells() is not necessarily the same
425  // as defined by primitiveMesh::cellShapes()
426  // Thus, for registered primitive types, do the lookup ourselves.
427  // Finally, the cellModel face number is re-mapped to the
428  // STAR-CD local face number
429 
430  label mapIndex = shape.model().index();
431 
432  // a registered primitive type
433  if (faceLookupIndex.found(mapIndex))
434  {
435  const faceList sFaces = shape.faces();
436  forAll(sFaces, sFacei)
437  {
438  if (faces[facei] == sFaces[sFacei])
439  {
440  cellFaceId = sFacei;
441  break;
442  }
443  }
444 
445  mapIndex = faceLookupIndex[mapIndex];
446  cellFaceId = foamToStarFaceAddr[mapIndex][cellFaceId];
447  }
448  // Info<< endl;
449 
450  boundId++;
451 
452  os
453  << boundId
454  << " " << cellId + 1
455  << " " << cellFaceId + 1
456  << " " << regionId
457  << " " << 0
458  << " " << bndType.c_str()
459  << endl;
460  }
461  }
462 }
463 
464 
465 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
466 
467 Foam::meshWriters::STARCD::STARCD
468 (
469  const polyMesh& mesh,
470  const scalar scaleFactor
471 )
472 :
473  meshWriter(mesh, scaleFactor)
474 {
477  getCellTable();
478 }
479 
480 
481 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
482 
484 {}
485 
486 
487 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
488 
489 void Foam::meshWriters::STARCD::rmFiles(const fileName& baseName) const
490 {
491  rm(baseName + ".vrt");
492  rm(baseName + ".cel");
493  rm(baseName + ".bnd");
494  rm(baseName + ".inp");
495 }
496 
497 
499 {
500  fileName baseName(meshName);
501 
502  if (baseName.empty())
503  {
504  baseName = meshWriter::defaultMeshName;
505 
506  if
507  (
508  mesh_.time().timeName() != "0"
509  && mesh_.time().timeName() != mesh_.time().constant()
510  )
511  {
512  baseName += "_" + mesh_.time().timeName();
513  }
514  }
515 
516  rmFiles(baseName);
517  writePoints(baseName);
518  writeCells(baseName);
519 
520  if (writeBoundary_)
521  {
522  writeBoundary(baseName);
523  }
524 
525  return true;
526 }
527 
528 
529 // ************************************************************************* //
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:424
label append(const dictionary &)
Append to the end, return index.
Definition: cellTable.C:130
write OpenFOAM meshes and/or results to another CFD format
Definition: meshWriter.H:89
static string defaultMeshName
Specify a default mesh name.
Definition: meshWriter.H:135
dictionary dict
static const label foamToStarFaceAddr[4][6]
Face addressing from OpenFOAM faces -> pro-STAR faces.
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
A class for handling file names.
Definition: fileName.H:69
virtual bool write(const fileName &meshName=fileName::null) const
Write volume mesh.
static const cellModel * tetModel
Definition: meshWriter.H:124
static const cellModel * prismModel
Definition: meshWriter.H:126
label nFaces() const
static word meshSubDir
Return the mesh sub-directory name (usually "polyMesh")
Definition: polyMesh.H:312
bool writeBoundary_
Write bnd file.
Definition: meshWriter.H:111
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:163
const cellShapeList & cellShapes() const
Return cell shapes.
boundaryRegion boundaryRegion_
boundaryRegion persistent data saved as a dictionary
Definition: meshWriter.H:114
bool empty() const
Return true if the hash table is empty.
Definition: HashTableI.H:72
List< face > faceList
Definition: faceListFwd.H:43
label nCells() const
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:256
const cellList & cells() const
static word timeName(const scalar, const int precision=precision_)
Return time name of given scalar time.
Definition: Time.C:626
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:1003
void readDict(const objectRegistry &, const word &name="cellTable", const fileName &instance="constant")
Read constant/cellTable.
Definition: cellTable.C:313
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
List< cellShape > cellShapeList
List of cellShapes and PtrList of List of cellShape.
Definition: cellShapeList.H:43
static char meshName[]
Definition: globalFoam.H:7
static const cellModel * pyrModel
Definition: meshWriter.H:125
bool found(const Key &) const
Return true if hashedEntry is found in table.
Definition: HashTable.C:113
const cellZoneMesh & cellZones() const
Return cell zone mesh.
Definition: polyMesh.H:472
const word & constant() const
Return constant name.
Definition: TimePaths.H:124
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1041
List< label > labelList
A List of labels.
Definition: labelList.H:56
virtual const faceList & faces() const
Return raw faces.
Definition: polyMesh.C:1028
scalar scaleFactor_
Scaling factor for points (eg, [m] -> [mm])
Definition: meshWriter.H:108
labelList cellTableId_
cellTable IDs for each cell
Definition: meshWriter.H:120
static const char nl
Definition: Ostream.H:265
word boundaryType(const word &name) const
Return BoundaryType corresponding to patch &#39;name&#39;.
const Time & time() const
Return time.
labelList f(nPoints)
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurrence of given element and return index,.
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
label index() const
Return index of model in the model list.
Definition: cellModelI.H:44
void readDict(const objectRegistry &, const word &name="boundaryRegion", const fileName &instance="constant")
Read constant/boundaryRegion.
label size() const
Return the number of elements in the UPtrList.
Definition: UPtrListI.H:29
void rmFiles(const fileName &baseName) const
Remove STAR-CD files for the baseName.
void setSize(const label)
Reset size of List.
Definition: List.C:281
label patchi
#define WarningInFunction
Report a warning using Foam::Warning.
virtual ~STARCD()
Destructor.
label cellId
messageStream Info
static const cellModel * hexModel
Definition: meshWriter.H:127
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
bool rm(const fileName &)
Remove a file, returning true if successful otherwise false.
Definition: POSIX.C:1008
label findIndex(const word &name) const
Return index corresponding to name.
Definition: cellTable.C:206
cellTable cellTable_
cellTable persistent data saved as a dictionary
Definition: meshWriter.H:117
void transfer(List< T > &)
Transfer the contents of the argument List into this list.
Definition: List.C:342
List< cell > cellList
list of cells
Definition: cellList.H:42
const polyMesh & mesh_
Mesh reference.
Definition: meshWriter.H:105