1 /*---------------------------------------------------------------------------*\
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
5 \\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
7 -------------------------------------------------------------------------------
9 This file is part of OpenFOAM.
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.
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
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
25 \*---------------------------------------------------------------------------*/
27 #include "STARCDMeshReader.H"
28 #include <OpenFOAM/cyclicPolyPatch.H>
29 #include <OpenFOAM/emptyPolyPatch.H>
30 #include <OpenFOAM/wallPolyPatch.H>
31 #include <OpenFOAM/symmetryPolyPatch.H>
32 #include <OpenFOAM/cellModeller.H>
33 #include <OpenFOAM/ListOps.H>
34 #include <OpenFOAM/IFstream.H>
35 #include <OpenFOAM/IOMap.H>
37 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
39 const char* Foam::meshReaders::STARCD::defaultBoundaryName =
40 "Default_Boundary_Region";
42 const char* Foam::meshReaders::STARCD::defaultSolidBoundaryName =
43 "Default_Boundary_Solid";
45 bool Foam::meshReaders::STARCD::keepSolids = false;
47 const int Foam::meshReaders::STARCD::starToFoamFaceAddr[4][6] =
49 { 4, 5, 2, 3, 0, 1 }, // 11 = pro-STAR hex
50 { 0, 1, 4, -1, 2, 3 }, // 12 = pro-STAR prism
51 { 3, -1, 2, -1, 1, 0 }, // 13 = pro-STAR tetra
52 { 0, -1, 4, 2, 1, 3 } // 14 = pro-STAR pyramid
56 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
58 void Foam::meshReaders::STARCD::readToNewline(IFstream& is)
65 while ((is) && ch != '\n');
69 bool Foam::meshReaders::STARCD::readHeader(IFstream& is, word fileSignature)
73 FatalErrorIn("meshReaders::STARCD::readHeader()")
74 << "cannot read " << fileSignature << " " << is.name()
84 // skip the rest of the line
87 // add other checks ...
88 if (header != fileSignature)
90 Info<< "header mismatch " << fileSignature << " " << is.name()
98 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
100 void Foam::meshReaders::STARCD::readAux(const objectRegistry& registry)
102 boundaryRegion_.readDict(registry);
103 cellTable_.readDict(registry);
107 // read in the points from the .vrt file
109 /*---------------------------------------------------------------------------*\
111 PROSTAR_VERTEX [newline]
114 <version> 0 0 0 0 0 0 0 [newline]
117 <vertexId> <x> <y> <z> [newline]
119 \*---------------------------------------------------------------------------*/
120 void Foam::meshReaders::STARCD::readPoints
122 const fileName& inputName,
123 const scalar scaleFactor
126 const word fileSignature = "PROSTAR_VERTEX";
127 label nPoints = 0, maxId = 0;
130 // get # points and maximum vertex label
132 IFstream is(inputName);
133 readHeader(is, fileSignature);
138 while ((is >> lineLabel).good())
141 maxId = max(maxId, lineLabel);
146 Info<< "Number of points = " << nPoints << endl;
148 // set sizes and reset to invalid values
150 points_.setSize(nPoints);
151 mapToFoamPointId_.setSize(maxId+1);
153 //- original Point number for a given vertex
154 // might need again in the future
155 //// labelList origPointId(nPoints);
156 //// origPointId = -1;
158 mapToFoamPointId_ = -1;
161 // construct pointList and conversion table
162 // from Star vertex numbers to Foam point labels
165 IFstream is(inputName);
166 readHeader(is, fileSignature);
171 while ((is >> lineLabel).good())
173 is >> points_[pointI].x()
174 >> points_[pointI].y()
175 >> points_[pointI].z();
177 // might need again in the future
178 //// origPointId[pointI] = lineLabel;
179 mapToFoamPointId_[lineLabel] = pointI;
183 if (nPoints > pointI)
186 points_.setSize(nPoints);
187 // might need again in the future
188 //// origPointId.setSize(nPoints);
191 if (scaleFactor > 1.0 + SMALL || scaleFactor < 1.0 - SMALL)
193 points_ *= scaleFactor;
198 FatalErrorIn("meshReaders::STARCD::readPoints()")
199 << "no points in file " << inputName
200 << abort(FatalError);
206 // read in the cells from the .cel file
208 /*---------------------------------------------------------------------------*\
210 PROSTAR_CELL [newline]
213 <version> 0 0 0 0 0 0 0 [newline]
216 <cellId> <shapeId> <nLabels> <cellTableId> <typeId> [newline]
217 <cellId> <int1> .. <int8>
218 <cellId> <int9> .. <int16>
238 For primitive cell shapes, the number of vertices will never exceed 8 (hexa)
239 and corresponds to <nLabels>.
240 For polyhedral, <nLabels> includess an index table comprising beg/end pairs
243 Strictly speaking, we only need the cellModeller for adding boundaries.
244 \*---------------------------------------------------------------------------*/
246 void Foam::meshReaders::STARCD::readCells(const fileName& inputName)
248 const word fileSignature = "PROSTAR_CELL";
249 label nFluids = 0, nSolids = 0, nBaffles = 0, nShells = 0;
252 bool unknownVertices = false;
256 // count nFluids, nSolids, nBaffle, nShell and maxId
257 // also see if polyhedral cells were used
259 IFstream is(inputName);
260 readHeader(is, fileSignature);
262 label lineLabel, shapeId, nLabels, cellTableId, typeId;
264 while ((is >> lineLabel).good())
266 label starCellId = lineLabel;
272 // skip the rest of the line
275 // max 8 indices per line
282 if (typeId == starcdFluidType)
285 maxId = max(maxId, starCellId);
287 if (!cellTable_.found(cellTableId))
289 cellTable_.setName(cellTableId);
290 cellTable_.setMaterial(cellTableId, "fluid");
293 else if (typeId == starcdSolidType)
298 maxId = max(maxId, starCellId);
301 if (!cellTable_.found(cellTableId))
303 cellTable_.setName(cellTableId);
304 cellTable_.setMaterial(cellTableId, "solid");
308 else if (typeId == starcdBaffleType)
310 // baffles have no cellTable entry
312 maxId = max(maxId, starCellId);
314 else if (typeId == starcdShellType)
317 if (!cellTable_.found(cellTableId))
319 cellTable_.setName(cellTableId);
320 cellTable_.setMaterial(cellTableId, "shell");
327 Info<< "Number of fluids = " << nFluids << nl
328 << "Number of baffles = " << nBaffles << nl;
331 Info<< "Number of solids = " << nSolids << nl;
335 Info<< "Ignored solids = " << nSolids << nl;
337 Info<< "Ignored shells = " << nShells << endl;
343 nCells = nFluids + nSolids;
350 cellFaces_.setSize(nCells);
351 cellShapes_.setSize(nCells);
352 cellTableId_.setSize(nCells);
354 // information for the interfaces
355 baffleFaces_.setSize(nBaffles);
357 // extra space for baffles
358 origCellId_.setSize(nCells + nBaffles);
359 mapToFoamCellId_.setSize(maxId+1);
360 mapToFoamCellId_ = -1;
363 // avoid undefined shapes for polyhedra
364 cellShape genericShape(*unknownModel, labelList(0));
367 // construct cellFaces_ and possibly cellShapes_
370 FatalErrorIn("meshReaders::STARCD::readCells()")
371 << "no cells in file " << inputName
372 << abort(FatalError);
376 IFstream is(inputName);
377 readHeader(is, fileSignature);
379 labelList starLabels(64);
380 label lineLabel, shapeId, nLabels, cellTableId, typeId;
385 while ((is >> lineLabel).good())
387 label starCellId = lineLabel;
393 if (nLabels > starLabels.size())
395 starLabels.setSize(nLabels);
399 // read indices - max 8 per line
400 for (label i = 0; i < nLabels; ++i)
410 if (typeId == starcdSolidType && !keepSolids)
415 // determine the foam cell shape
416 const cellModel* curModelPtr = NULL;
422 curModelPtr = hexModel;
425 curModelPtr = prismModel;
428 curModelPtr = tetModel;
431 curModelPtr = pyrModel;
437 // primitive cell - use shapes
439 // convert orig vertex Id to point label
441 for (label i=0; i < nLabels; ++i)
443 label pointId = mapToFoamPointId_[starLabels[i]];
446 Info<< "Cells inconsistent with vertex file. "
447 << "Star vertex " << starLabels[i]
448 << " does not exist" << endl;
450 unknownVertices = true;
452 starLabels[i] = pointId;
460 // record original cell number and lookup
461 origCellId_[cellI] = starCellId;
462 mapToFoamCellId_[starCellId] = cellI;
464 cellTableId_[cellI] = cellTableId;
465 cellShapes_[cellI] = cellShape
468 SubList<label>(starLabels, nLabels)
471 cellFaces_[cellI] = cellShapes_[cellI].faces();
474 else if (shapeId == starcdPoly)
477 label nFaces = starLabels[0] - 1;
479 // convert orig vertex id to point label
480 // start with offset (skip the index table)
482 for (label i=starLabels[0]; i < nLabels; ++i)
484 label pointId = mapToFoamPointId_[starLabels[i]];
487 Info<< "Cells inconsistent with vertex file. "
488 << "Star vertex " << starLabels[i]
489 << " does not exist" << endl;
491 unknownVertices = true;
493 starLabels[i] = pointId;
501 // traverse beg/end indices
502 faceList faces(nFaces);
504 for (label i=0; i < nFaces; ++i)
506 label beg = starLabels[i];
507 label n = starLabels[i+1] - beg;
511 SubList<label>(starLabels, n, beg)
525 Info<< "star cell " << starCellId << " has "
527 << " empty faces - could cause boundary "
528 << "addressing problems"
532 faces.setSize(nFaces);
537 FatalErrorIn("meshReaders::STARCD::readCells()")
538 << "star cell " << starCellId << " has " << nFaces
539 << abort(FatalError);
542 // record original cell number and lookup
543 origCellId_[cellI] = starCellId;
544 mapToFoamCellId_[starCellId] = cellI;
546 cellTableId_[cellI] = cellTableId;
547 cellShapes_[cellI] = genericShape;
548 cellFaces_[cellI] = faces;
551 else if (typeId == starcdBaffleType)
555 // convert orig vertex id to point label
557 for (label i=0; i < nLabels; ++i)
559 label pointId = mapToFoamPointId_[starLabels[i]];
562 Info<< "Baffles inconsistent with vertex file. "
563 << "Star vertex " << starLabels[i]
564 << " does not exist" << endl;
566 unknownVertices = true;
568 starLabels[i] = pointId;
579 SubList<label>(starLabels, nLabels)
587 baffleFaces_[baffleI] = f;
588 // insert lookup addressing in normal list
589 mapToFoamCellId_[starCellId] = nCells + baffleI;
590 origCellId_[nCells + baffleI] = starCellId;
596 baffleFaces_.setSize(baffleI);
601 FatalErrorIn("meshReaders::STARCD::readCells()")
602 << "cells with unknown vertices"
603 << abort(FatalError);
609 Info<< "CELLS READ" << endl;
613 mapToFoamPointId_.clear();
617 // read in the boundaries from the .bnd file
619 /*---------------------------------------------------------------------------*\
621 PROSTAR_BOUNDARY [newline]
624 <version> 0 0 0 0 0 0 0 [newline]
627 <boundId> <cellId> <cellFace> <regionId> 0 <boundaryType> [newline]
629 where boundaryType is truncated to 4 characters from one of the following:
635 \*---------------------------------------------------------------------------*/
637 void Foam::meshReaders::STARCD::readBoundary(const fileName& inputName)
639 const word fileSignature = "PROSTAR_BOUNDARY";
640 label nPatches = 0, nFaces = 0, nBafflePatches = 0, maxId = 0;
641 label lineLabel, starCellId, cellFaceId, starRegion, configNumber;
644 labelList mapToFoamPatchId(1000, -1);
645 labelList nPatchFaces(1000, 0);
646 labelList origRegion(1000, 0);
647 patchTypes_.setSize(1000);
649 // this is what we seem to need
650 // these MUST correspond to starToFoamFaceAddr
652 Map<label> faceLookupIndex;
654 faceLookupIndex.insert(hexModel->index(), 0);
655 faceLookupIndex.insert(prismModel->index(), 1);
656 faceLookupIndex.insert(tetModel->index(), 2);
657 faceLookupIndex.insert(pyrModel->index(), 3);
661 // no. of faces (nFaces), no. of patches (nPatches)
662 // and for each of these patches the number of faces
663 // (nPatchFaces[patchLabel])
665 // and a conversion table from Star regions to (Foam) patchLabels
667 // additionally note the no. of baffle patches (nBafflePatches)
668 // so that we sort these to the end of the patch list
669 // - this makes it easier to transfer them to an adjacent patch if reqd
671 IFstream is(inputName);
675 readHeader(is, fileSignature);
677 while ((is >> lineLabel).good())
686 // Build translation table to convert star patch to foam patch
687 label patchLabel = mapToFoamPatchId[starRegion];
688 if (patchLabel == -1)
690 patchLabel = nPatches;
691 mapToFoamPatchId[starRegion] = patchLabel;
692 origRegion[patchLabel] = starRegion;
693 patchTypes_[patchLabel] = patchType;
695 maxId = max(maxId, starRegion);
697 if (patchType == "BAFF") // should actually be case-insensitive
704 nPatchFaces[patchLabel]++;
709 Info<< "No boundary faces in file " << inputName << endl;
714 Info<< "Could not read boundary file " << inputName << endl;
718 // keep empty patch region in reserve
720 Info<< "Number of patches = " << nPatches
721 << " (including extra for missing)" << endl;
724 origRegion.setSize(nPatches);
725 patchTypes_.setSize(nPatches);
726 patchNames_.setSize(nPatches);
727 nPatchFaces.setSize(nPatches);
729 // add our empty patch
730 origRegion[nPatches-1] = 0;
731 nPatchFaces[nPatches-1] = 0;
732 patchTypes_[nPatches-1] = "none";
735 // - use 'Label' entry from "constant/boundaryRegion" dictionary
736 forAll(patchTypes_, patchI)
738 bool foundName = false, foundType = false;
740 Map<dictionary>::const_iterator
741 iter = boundaryRegion_.find(origRegion[patchI]);
745 iter != boundaryRegion_.end()
748 foundType = iter().readIfPresent
754 foundName = iter().readIfPresent
761 // consistent names, in long form and in lowercase
765 forAllIter(string, patchTypes_[patchI], i)
770 if (patchTypes_[patchI] == "symp")
772 patchTypes_[patchI] = "symplane";
774 else if (patchTypes_[patchI] == "cycl")
776 patchTypes_[patchI] = "cyclic";
778 else if (patchTypes_[patchI] == "baff")
780 patchTypes_[patchI] = "baffle";
782 else if (patchTypes_[patchI] == "moni")
784 patchTypes_[patchI] = "monitoring";
788 // create a name if needed
791 patchNames_[patchI] =
792 patchTypes_[patchI] + "_" + name(origRegion[patchI]);
796 // enforce name "Default_Boundary_Region"
797 patchNames_[nPatches-1] = defaultBoundaryName;
799 // sort according to ascending region numbers, but leave
800 // Default_Boundary_Region as the final patch
802 labelList sortedIndices;
803 sortedOrder(SubList<label>(origRegion, nPatches-1), sortedIndices);
805 labelList oldToNew = identity(nPatches);
806 forAll(sortedIndices, i)
808 oldToNew[sortedIndices[i]] = i;
811 inplaceReorder(oldToNew, origRegion);
812 inplaceReorder(oldToNew, patchTypes_);
813 inplaceReorder(oldToNew, patchNames_);
814 inplaceReorder(oldToNew, nPatchFaces);
817 // re-sort to have baffles near the end
821 labelList oldToNew = identity(nPatches);
823 label baffleIndex = (nPatches-1 - nBafflePatches);
825 for (label i=0; i < oldToNew.size()-1; ++i)
827 if (patchTypes_[i] == "baffle")
829 oldToNew[i] = baffleIndex++;
833 oldToNew[i] = newIndex++;
837 inplaceReorder(oldToNew, origRegion);
838 inplaceReorder(oldToNew, patchTypes_);
839 inplaceReorder(oldToNew, patchNames_);
840 inplaceReorder(oldToNew, nPatchFaces);
843 mapToFoamPatchId.setSize(maxId+1, -1);
844 forAll(origRegion, patchI)
846 mapToFoamPatchId[origRegion[patchI]] = patchI;
849 boundaryIds_.setSize(nPatches);
850 forAll(boundaryIds_, patchI)
852 boundaryIds_[patchI].setSize(nPatchFaces[patchI]);
853 nPatchFaces[patchI] = 0;
859 if (nPatches > 1 && mapToFoamCellId_.size() > 1)
861 IFstream is(inputName);
862 readHeader(is, fileSignature);
864 while ((is >> lineLabel).good())
873 label patchI = mapToFoamPatchId[starRegion];
875 // zero-based indexing
880 // convert to FOAM cell number
881 if (starCellId < mapToFoamCellId_.size())
883 cellId = mapToFoamCellId_[starCellId];
889 << "Boundaries inconsistent with cell file. "
890 << "Star cell " << starCellId << " does not exist"
895 // restrict lookup to volume cells (no baffles)
896 if (cellId < cellShapes_.size())
898 label index = cellShapes_[cellId].model().index();
899 if (faceLookupIndex.found(index))
901 index = faceLookupIndex[index];
902 cellFaceId = starToFoamFaceAddr[index][cellFaceId];
907 // we currently use cellId >= nCells to tag baffles,
908 // we can also use a negative face number
912 boundaryIds_[patchI][nPatchFaces[patchI]] =
913 cellFaceIdentifier(cellId, cellFaceId);
915 #ifdef DEBUG_BOUNDARY
916 Info<< "bnd " << cellId << " " << cellFaceId << endl;
918 // increment counter of faces in current patch
919 nPatchFaces[patchI]++;
924 // retain original information in patchPhysicalTypes_ - overwrite latter
925 patchPhysicalTypes_.setSize(patchTypes_.size());
928 forAll(boundaryIds_, patchI)
930 // resize - avoid invalid boundaries
931 if (nPatchFaces[patchI] < boundaryIds_[patchI].size())
933 boundaryIds_[patchI].setSize(nPatchFaces[patchI]);
936 word origType = patchTypes_[patchI];
937 patchPhysicalTypes_[patchI] = origType;
939 if (origType == "symplane")
941 patchTypes_[patchI] = symmetryPolyPatch::typeName;
942 patchPhysicalTypes_[patchI] = patchTypes_[patchI];
944 else if (origType == "wall")
946 patchTypes_[patchI] = wallPolyPatch::typeName;
947 patchPhysicalTypes_[patchI] = patchTypes_[patchI];
949 else if (origType == "cyclic")
951 // incorrect. should be cyclicPatch but this
952 // requires info on connected faces.
953 patchTypes_[patchI] = cyclicPolyPatch::typeName;
954 patchPhysicalTypes_[patchI] = patchTypes_[patchI];
956 else if (origType == "baffle")
958 // incorrect. tag the patch until we get proper support.
959 // set physical type to a canonical "baffle"
960 patchTypes_[patchI] = emptyPolyPatch::typeName;
961 patchPhysicalTypes_[patchI] = "baffle";
965 patchTypes_[patchI] = polyPatch::typeName;
968 Info<< "patch " << patchI
969 << " (region " << origRegion[patchI]
970 << ": " << origType << ") type: '" << patchTypes_[patchI]
971 << "' name: " << patchNames_[patchI] << endl;
975 mapToFoamCellId_.clear();
981 // remove unused points
983 void Foam::meshReaders::STARCD::cullPoints()
985 label nPoints = points_.size();
986 labelList oldToNew(nPoints, -1);
988 // loop through cell faces and note which points are being used
989 forAll(cellFaces_, cellI)
991 const faceList& faces = cellFaces_[cellI];
994 const labelList& labels = faces[i];
997 oldToNew[labels[j]]++;
1002 // the new ordering and the count of unused points
1006 if (oldToNew[i] >= 0)
1008 oldToNew[i] = pointI++;
1012 // report unused points
1013 if (nPoints > pointI)
1015 Info<< "Unused points = " << (nPoints - pointI) << endl;
1018 // adjust points and truncate
1019 inplaceReorder(oldToNew, points_);
1020 points_.setSize(nPoints);
1022 // adjust cellFaces - with mesh shapes this might be faster
1023 forAll(cellFaces_, cellI)
1025 faceList& faces = cellFaces_[cellI];
1028 inplaceRenumber(oldToNew, faces[i]);
1033 forAll(baffleFaces_, faceI)
1035 inplaceRenumber(oldToNew, baffleFaces_[faceI]);
1041 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
1043 bool Foam::meshReaders::STARCD::readGeometry(const scalar scaleFactor)
1045 // Info<< "called meshReaders::STARCD::readGeometry" << endl;
1047 readPoints(geometryFile_ + ".vrt", scaleFactor);
1048 readCells(geometryFile_ + ".cel");
1050 readBoundary(geometryFile_ + ".bnd");
1056 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
1058 Foam::meshReaders::STARCD::STARCD
1060 const fileName& prefix,
1061 const objectRegistry& registry,
1062 const scalar scaleFactor
1065 meshReader(prefix, scaleFactor),
1067 mapToFoamPointId_(0),
1074 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
1076 Foam::meshReaders::STARCD::~STARCD()
1080 // ************************ vim: set sw=4 sts=4 et: ************************ //