FreeFOAM The Cross-Platform CFD Toolkit
STARCDMeshWriter.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 \*---------------------------------------------------------------------------*/
25 
27 
28 #include <OpenFOAM/Time.H>
29 #include <OpenFOAM/SortableList.H>
30 #include <OpenFOAM/OFstream.H>
31 
32 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
33 
34 const char* Foam::meshWriters::STARCD::defaultBoundaryName =
35  "Default_Boundary_Region";
36 
37 const Foam::label Foam::meshWriters::STARCD::foamToStarFaceAddr[4][6] =
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  "constant",
77  mesh_,
80  false
81  )
82  );
83 
84  bool useCellZones = false;
85  cellTableId_.setSize(mesh_.nCells(), -1);
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  {
101  WarningIn("STARCD::getCellTable()")
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 {
475  boundaryRegion_.readDict(mesh_);
476  cellTable_.readDict(mesh_);
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() != "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 
530 (
531  const fileName& meshName,
532  const bool& triangulate
533 ) const
534 {
535  fileName baseName(meshName);
536 
537  if (baseName.empty())
538  {
540 
541  if
542  (
543  mesh_.time().timeName() != "0"
544  && mesh_.time().timeName() != "constant"
545  )
546  {
547  baseName += "_" + mesh_.time().timeName();
548  }
549  }
550 
551  rmFiles(baseName);
552 
553  OFstream celFile(baseName + ".cel");
554  writeHeader(celFile, "CELL");
555 
556  Info << "Writing " << celFile.name() << endl;
557 
558  // mesh and patch info
559  const pointField& points = mesh_.points();
560  const labelList& owner = mesh_.faceOwner();
561  const faceList& meshFaces = mesh_.faces();
562  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
563 
564  label shapeId = 3; // shell/baffle element
565  label typeId = 4; // 4(shell)
566 
567  // remember which points need to be written
568  labelHashSet pointHash;
569 
570  // write boundary faces as normal STAR-CD mesh
571  if (triangulate)
572  {
573  // cell Id has no particular meaning - just increment
574  // use the cellTable id from the patch Number
575  label cellId = 0;
576 
577  forAll(patches, patchI)
578  {
579  label patchStart = patches[patchI].start();
580  label patchSize = patches[patchI].size();
581 
582  label ctableId = patchI + 1;
583 
584  for
585  (
586  label faceI = patchStart;
587  faceI < (patchStart + patchSize);
588  ++faceI
589  )
590  {
591  const face& f = meshFaces[faceI];
592 
593  label nTri = f.nTriangles(points);
594  faceList triFaces;
595 
596  // triangulate polygons, but not quads
597  if (nTri <= 2)
598  {
599  triFaces.setSize(1);
600  triFaces[0] = f;
601  }
602  else
603  {
604  triFaces.setSize(nTri);
605  nTri = 0;
606  f.triangles(points, nTri, triFaces);
607  }
608 
609  forAll(triFaces, faceI)
610  {
611  const labelList& vrtList = triFaces[faceI];
612 
613  celFile
614  << cellId + 1 << " "
615  << shapeId << " "
616  << vrtList.size() << " "
617  << ctableId << " "
618  << typeId;
619 
620  // must be 3 (triangle) but could be quad
621  label count = 0;
622  forAll(vrtList, i)
623  {
624  if ((count % 8) == 0)
625  {
626  celFile
627  << nl
628  << " " << cellId + 1;
629  }
630  // remember which points we'll need to write
631  pointHash.insert(vrtList[i]);
632  celFile << " " << vrtList[i] + 1;
633  count++;
634  }
635  celFile << endl;
636 
637  cellId++;
638  }
639  }
640  }
641  }
642  else
643  {
644  // cell Id is the OpenFOAM face Id
645  // use the cellTable id from the face owner
646  // - allows separation of parts
647  forAll(patches, patchI)
648  {
649  label patchStart = patches[patchI].start();
650  label patchSize = patches[patchI].size();
651 
652  for
653  (
654  label faceI = patchStart;
655  faceI < (patchStart + patchSize);
656  ++faceI
657  )
658  {
659  const labelList& vrtList = meshFaces[faceI];
660  label cellId = faceI;
661 
662  celFile
663  << cellId + 1 << " "
664  << shapeId << " "
665  << vrtList.size() << " "
666  << cellTableId_[owner[faceI]] << " "
667  << typeId;
668 
669  // likely <= 8 vertices, but prevent overrun anyhow
670  label count = 0;
671  forAll(vrtList, i)
672  {
673  if ((count % 8) == 0)
674  {
675  celFile
676  << nl
677  << " " << cellId + 1;
678  }
679  // remember which points we'll need to write
680  pointHash.insert(vrtList[i]);
681  celFile << " " << vrtList[i] + 1;
682  count++;
683  }
684  celFile << endl;
685  }
686  }
687  }
688 
689  OFstream vrtFile(baseName + ".vrt");
690  writeHeader(vrtFile, "VERTEX");
691 
692  vrtFile.precision(10);
693  vrtFile.setf(std::ios::showpoint); // force decimal point for Fortran
694 
695  Info << "Writing " << vrtFile.name() << endl;
696 
697  // build sorted table of contents
698  SortableList<label> toc(pointHash.size());
699  {
700  label i = 0;
701  forAllConstIter(labelHashSet, pointHash, iter)
702  {
703  toc[i++] = iter.key();
704  }
705  }
706  toc.sort();
707  toc.shrink();
708  pointHash.clear();
709 
710  // write points in sorted order
711  forAll(toc, i)
712  {
713  label vrtId = toc[i];
714  vrtFile
715  << vrtId + 1
716  << " " << scaleFactor_ * points[vrtId].x()
717  << " " << scaleFactor_ * points[vrtId].y()
718  << " " << scaleFactor_ * points[vrtId].z()
719  << endl;
720  }
721 
722  return true;
723 }
724 
725 
726 // ************************ vim: set sw=4 sts=4 et: ************************ //