src/conversion/meshWriter/starcd/STARCDMeshWriter.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2004-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
  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
  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",
  76             polyMesh::meshSubDir,
  77             mesh_,
  78             IOobject::READ_IF_PRESENT,
  79             IOobject::NO_WRITE,
  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
 483 Foam::meshWriters::STARCD::~STARCD()
 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
 498 bool Foam::meshWriters::STARCD::write(const fileName& meshName) const
 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
 529 // ************************************************************************* //