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) 1991-2008 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 the
  13     Free Software Foundation; either version 2 of the License, or (at your
  14     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, write to the Free Software Foundation,
  23     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
  24
  25 \*---------------------------------------------------------------------------*/
  26
  27 #include "STARCDMeshWriter.H"
  28
  29 #include "Time.H"
  30 #include "SortableList.H"
  31 #include "OFstream.H"
  32
  33 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  34
  35 const char* Foam::meshWriters::STARCD::defaultBoundaryName =
  36     "Default_Boundary_Region";
  37
  38 const Foam::label Foam::meshWriters::STARCD::foamToStarFaceAddr[4][6] =
  39 {
  40     { 4, 5, 2, 3, 0, 1 },     // 11 = pro-STAR hex
  41     { 0, 1, 4, 5, 2, -1 },    // 12 = pro-STAR prism
  42     { 5, 4, 2, 0, -1, -1 },   // 13 = pro-STAR tetra
  43     { 0, 4, 3, 5, 2, -1 }     // 14 = pro-STAR pyramid
  44 };
  45
  46
  47 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  48
  49 Foam::label Foam::meshWriters::STARCD::findDefaultBoundary() const
  50 {
  51     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
  52
  53     label id = -1;
  54
  55     // find Default_Boundary_Region if it exists
  56     forAll(patches, patchI)
  57     {
  58         if (defaultBoundaryName == patches[patchI].name())
  59         {
  60             id = patchI;
  61             break;
  62         }
  63     }
  64     return id;
  65 }
  66
  67
  68 void Foam::meshWriters::STARCD::getCellTable()
  69 {
  70     // read constant/polyMesh/propertyName
  71     IOList<label> ioList
  72     (
  73         IOobject
  74         (
  75             "cellTableId",
  76             "constant",
  77             polyMesh::meshSubDir,
  78             mesh_,
  79             IOobject::READ_IF_PRESENT,
  80             IOobject::NO_WRITE,
  81             false
  82         )
  83     );
  84
  85     bool useCellZones = false;
  86     cellTableId_.setSize(mesh_.nCells(), -1);
  87
  88     // get information from constant/polyMesh/cellTableId if possible
  89     if (ioList.headerOk())
  90     {
  91         if (ioList.size() == mesh_.nCells())
  92         {
  93             cellTableId_.transfer(ioList);
  94
  95             if (!cellTable_.size())
  96             {
  97                 Info<< "no cellTable information available" << endl;
  98             }
  99         }
 100         else
 101         {
 102             WarningIn("STARCD::getCellTable()")
 103                 << ioList.objectPath() << " has incorrect number of cells "
 104                 << " - use cellZone information"
 105                 << endl;
 106
 107             ioList.clear();
 108             useCellZones = true;
 109         }
 110     }
 111     else
 112     {
 113         useCellZones = true;
 114     }
 115
 116
 117     if (useCellZones)
 118     {
 119         if (!cellTable_.size())
 120         {
 121             Info<< "created cellTable from cellZones" << endl;
 122             cellTable_ = mesh_;
 123         }
 124
 125         // track if there are unzoned cells
 126         label nUnzoned = mesh_.nCells();
 127
 128         // get the cellZone <-> cellTable correspondence
 129         Info<< "matching cellZones to cellTable" << endl;
 130
 131         forAll (mesh_.cellZones(), zoneI)
 132         {
 133             const cellZone& cZone = mesh_.cellZones()[zoneI];
 134             if (cZone.size())
 135             {
 136                 nUnzoned -= cZone.size();
 137
 138                 label tableId = cellTable_.findIndex(cZone.name());
 139                 if (tableId < 0)
 140                 {
 141                     dictionary dict;
 142
 143                     dict.add("Label", cZone.name());
 144                     dict.add("MaterialType", "fluid");
 145                     tableId = cellTable_.append(dict);
 146                 }
 147
 148                 forAll (cZone, i)
 149                 {
 150                     cellTableId_[cZone[i]] = tableId;
 151                 }
 152             }
 153         }
 154
 155         if (nUnzoned)
 156         {
 157             dictionary dict;
 158
 159             dict.add("Label", "__unZonedCells__");
 160             dict.add("MaterialType", "fluid");
 161             label tableId = cellTable_.append(dict);
 162
 163             forAll (cellTableId_, i)
 164             {
 165                 if (cellTableId_[i] < 0)
 166                 {
 167                     cellTableId_[i] = tableId;
 168                 }
 169             }
 170         }
 171     }
 172 }
 173
 174
 175 void Foam::meshWriters::STARCD::writeHeader(Ostream& os, const char* filetype)
 176 {
 177     os  << "PROSTAR_" << filetype << nl
 178         << 4000
 179         << " " << 0
 180         << " " << 0
 181         << " " << 0
 182         << " " << 0
 183         << " " << 0
 184         << " " << 0
 185         << " " << 0
 186         << endl;
 187 }
 188
 189
 190 void Foam::meshWriters::STARCD::writePoints(const fileName& prefix) const
 191 {
 192     OFstream os(prefix + ".vrt");
 193     writeHeader(os, "VERTEX");
 194
 195     // Set the precision of the points data to 10
 196     os.precision(10);
 197
 198     // force decimal point for Fortran input
 199     os.setf(std::ios::showpoint);
 200
 201     const pointField& points = mesh_.points();
 202
 203     Info<< "Writing " << os.name() << " : "
 204         << points.size() << " points" << endl;
 205
 206     forAll(points, ptI)
 207     {
 208         // convert [m] -> [mm]
 209         os
 210             << ptI + 1 << " "
 211             << scaleFactor_ * points[ptI].x() << " "
 212             << scaleFactor_ * points[ptI].y() << " "
 213             << scaleFactor_ * points[ptI].z() << nl;
 214     }
 215     os.flush();
 216
 217 }
 218
 219
 220 void Foam::meshWriters::STARCD::writeCells(const fileName& prefix) const
 221 {
 222     OFstream os(prefix + ".cel");
 223     writeHeader(os, "CELL");
 224
 225     // this is what we seem to need
 226     // map foam cellModeller index -> star shape
 227     Map<label> shapeLookupIndex;
 228     shapeLookupIndex.insert(hexModel->index(), 11);
 229     shapeLookupIndex.insert(prismModel->index(), 12);
 230     shapeLookupIndex.insert(tetModel->index(), 13);
 231     shapeLookupIndex.insert(pyrModel->index(), 14);
 232
 233     const cellShapeList& shapes = mesh_.cellShapes();
 234     const cellList& cells  = mesh_.cells();
 235     const faceList& faces  = mesh_.faces();
 236     const labelList& owner = mesh_.faceOwner();
 237
 238     Info<< "Writing " << os.name() << " : "
 239         << cells.size() << " cells" << endl;
 240
 241     forAll(cells, cellId)
 242     {
 243         label tableId = cellTableId_[cellId];
 244         label materialType  = 1;        // 1(fluid)
 245         if (cellTable_.found(tableId))
 246         {
 247             const dictionary& dict = cellTable_[tableId];
 248             if (dict.found("MaterialType"))
 249             {
 250                 word matType;
 251                 dict.lookup("MaterialType") >> matType;
 252                 if (matType == "solid")
 253                 {
 254                     materialType = 2;
 255                 }
 256
 257             }
 258         }
 259
 260         const cellShape& shape = shapes[cellId];
 261         label mapIndex = shape.model().index();
 262
 263         // a registered primitive type
 264         if (shapeLookupIndex.found(mapIndex))
 265         {
 266             label shapeId = shapeLookupIndex[mapIndex];
 267             const labelList& vrtList = shapes[cellId];
 268
 269             os  << cellId + 1
 270                 << " " << shapeId
 271                 << " " << vrtList.size()
 272                 << " " << tableId
 273                 << " " << materialType;
 274
 275             // primitives have <= 8 vertices, but prevent overrun anyhow
 276             // indent following lines for ease of reading
 277             label count = 0;
 278             forAll(vrtList, i)
 279             {
 280                 if ((count % 8) == 0)
 281                 {
 282                     os  << nl
 283                         << "  " << cellId + 1;
 284                 }
 285                 os << " " << vrtList[i] + 1;
 286                 count++;
 287             }
 288             os << endl;
 289
 290         }
 291         else
 292         {
 293             label shapeId = 255;        // treat as general polyhedral
 294             const labelList& cFaces  = cells[cellId];
 295
 296             // create (beg,end) indices
 297             List<label> indices(cFaces.size() + 1);
 298             indices[0] = indices.size();
 299
 300             label count = indices.size();
 301             // determine the total number of vertices
 302             forAll(cFaces, faceI)
 303             {
 304                 count += faces[cFaces[faceI]].size();
 305                 indices[faceI+1] = count;
 306             }
 307
 308             os  << cellId + 1
 309                 << " " << shapeId
 310                 << " " << count
 311                 << " " << tableId
 312                 << " " << materialType;
 313
 314             // write indices - max 8 per line
 315             // indent following lines for ease of reading
 316             count = 0;
 317             forAll(indices, i)
 318             {
 319                 if ((count % 8) == 0)
 320                 {
 321                     os  << nl
 322                         << "  " << cellId + 1;
 323                 }
 324                 os << " " << indices[i];
 325                 count++;
 326             }
 327
 328             // write faces - max 8 per line
 329             forAll(cFaces, faceI)
 330             {
 331                 label meshFace = cFaces[faceI];
 332                 face f;
 333
 334                 if (owner[meshFace] == cellId)
 335                 {
 336                     f = faces[meshFace];
 337                 }
 338                 else
 339                 {
 340                     f = faces[meshFace].reverseFace();
 341                 }
 342
 343                 forAll(f, i)
 344                 {
 345                     if ((count % 8) == 0)
 346                     {
 347                         os  << nl
 348                             << "  " << cellId + 1;
 349                     }
 350
 351                     os << " " << f[i] + 1;
 352                     count++;
 353                 }
 354             }
 355
 356             os << endl;
 357         }
 358     }
 359 }
 360
 361
 362 void Foam::meshWriters::STARCD::writeBoundary(const fileName& prefix) const
 363 {
 364     OFstream os(prefix + ".bnd");
 365     writeHeader(os, "BOUNDARY");
 366
 367     const cellShapeList& shapes = mesh_.cellShapes();
 368     const cellList& cells  = mesh_.cells();
 369     const faceList& faces  = mesh_.faces();
 370     const labelList& owner = mesh_.faceOwner();
 371     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 372
 373     // this is what we seem to need
 374     // these MUST correspond to foamToStarFaceAddr
 375     //
 376     Map<label> faceLookupIndex;
 377     faceLookupIndex.insert(hexModel->index(), 0);
 378     faceLookupIndex.insert(prismModel->index(), 1);
 379     faceLookupIndex.insert(tetModel->index(), 2);
 380     faceLookupIndex.insert(pyrModel->index(), 3);
 381
 382     Info<< "Writing " << os.name() << " : "
 383         << (mesh_.nFaces() - patches[0].start()) << " boundaries" << endl;
 384
 385
 386     label defaultId = findDefaultBoundary();
 387
 388     //
 389     // write boundary faces - skip Default_Boundary_Region entirely
 390     //
 391     label boundId = 0;
 392     forAll(patches, patchI)
 393     {
 394         label regionId = patchI;
 395         if (regionId == defaultId)
 396         {
 397             continue;   // skip - already written
 398         }
 399         else if (defaultId == -1 || regionId < defaultId)
 400         {
 401             regionId++;
 402         }
 403
 404         label patchStart = patches[patchI].start();
 405         label patchSize  = patches[patchI].size();
 406         word  bndType = boundaryRegion_.boundaryType(patches[patchI].name());
 407
 408         for
 409         (
 410             label faceI = patchStart;
 411             faceI < (patchStart + patchSize);
 412             ++faceI
 413         )
 414         {
 415             label cellId = owner[faceI];
 416             const labelList& cFaces  = cells[cellId];
 417             const cellShape& shape = shapes[cellId];
 418             label cellFaceId = findIndex(cFaces, faceI);
 419
 420             //      Info<< "cell " << cellId + 1 << " face " << faceI
 421             //          << " == " << faces[faceI]
 422             //          << " is index " << cellFaceId << " from " << cFaces;
 423
 424             // Unfortunately, the order of faces returned by
 425             //   primitiveMesh::cells() is not necessarily the same
 426             //   as defined by primitiveMesh::cellShapes()
 427             // Thus, for registered primitive types, do the lookup ourselves.
 428             // Finally, the cellModel face number is re-mapped to the
 429             // STAR-CD local face number
 430
 431             label mapIndex = shape.model().index();
 432
 433             // a registered primitive type
 434             if (faceLookupIndex.found(mapIndex))
 435             {
 436                 const faceList sFaces = shape.faces();
 437                 forAll(sFaces, sFaceI)
 438                 {
 439                     if (faces[faceI] == sFaces[sFaceI])
 440                     {
 441                         cellFaceId = sFaceI;
 442                         break;
 443                     }
 444                 }
 445
 446                 mapIndex = faceLookupIndex[mapIndex];
 447                 cellFaceId = foamToStarFaceAddr[mapIndex][cellFaceId];
 448             }
 449             // Info<< endl;
 450
 451             boundId++;
 452
 453             os
 454                 << boundId
 455                 << " " << cellId + 1
 456                 << " " << cellFaceId + 1
 457                 << " " << regionId
 458                 << " " << 0
 459                 << " " << bndType.c_str()
 460                 << endl;
 461         }
 462     }
 463 }
 464
 465
 466 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 467
 468 Foam::meshWriters::STARCD::STARCD
 469 (
 470     const polyMesh& mesh,
 471     const scalar scaleFactor
 472 )
 473 :
 474     meshWriter(mesh, scaleFactor)
 475 {
 476     boundaryRegion_.readDict(mesh_);
 477     cellTable_.readDict(mesh_);
 478     getCellTable();
 479 }
 480
 481
 482 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 483
 484 Foam::meshWriters::STARCD::~STARCD()
 485 {}
 486
 487
 488 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 489
 490 void Foam::meshWriters::STARCD::rmFiles(const fileName& baseName) const
 491 {
 492     rm(baseName + ".vrt");
 493     rm(baseName + ".cel");
 494     rm(baseName + ".bnd");
 495     rm(baseName + ".inp");
 496 }
 497
 498
 499 bool Foam::meshWriters::STARCD::write(const fileName& meshName) const
 500 {
 501     fileName baseName(meshName);
 502
 503     if (!baseName.size())
 504     {
 505         baseName = meshWriter::defaultMeshName;
 506
 507         if
 508         (
 509             mesh_.time().timeName() != "0"
 510          && mesh_.time().timeName() != "constant"
 511         )
 512         {
 513             baseName += "_" + mesh_.time().timeName();
 514         }
 515     }
 516
 517     rmFiles(baseName);
 518     writePoints(baseName);
 519     writeCells(baseName);
 520
 521     if (writeBoundary_)
 522     {
 523         writeBoundary(baseName);
 524     }
 525
 526     return true;
 527 }
 528
 529
 530 bool Foam::meshWriters::STARCD::writeSurface
 531 (
 532     const fileName& meshName,
 533     const bool& triangulate
 534 ) const
 535 {
 536     fileName baseName(meshName);
 537
 538     if (!baseName.size())
 539     {
 540         baseName = meshWriter::defaultSurfaceName;
 541
 542         if
 543         (
 544             mesh_.time().timeName() != "0"
 545          && mesh_.time().timeName() != "constant"
 546         )
 547         {
 548             baseName += "_" + mesh_.time().timeName();
 549         }
 550     }
 551
 552     rmFiles(baseName);
 553
 554     OFstream celFile(baseName + ".cel");
 555     writeHeader(celFile, "CELL");
 556
 557     Info << "Writing " << celFile.name() << endl;
 558
 559     // mesh and patch info
 560     const pointField& points = mesh_.points();
 561     const labelList& owner = mesh_.faceOwner();
 562     const faceList& meshFaces = mesh_.faces();
 563     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 564
 565     label shapeId = 3;  // shell/baffle element
 566     label typeId  = 4;  // 4(shell)
 567
 568     // remember which points need to be written
 569     labelHashSet pointHash;
 570
 571     // write boundary faces as normal STAR-CD mesh
 572     if (triangulate)
 573     {
 574         // cell Id has no particular meaning - just increment
 575         // use the cellTable id from the patch Number
 576         label cellId = 0;
 577
 578         forAll(patches, patchI)
 579         {
 580             label patchStart = patches[patchI].start();
 581             label patchSize  = patches[patchI].size();
 582
 583             label ctableId = patchI + 1;
 584
 585             for
 586             (
 587                 label faceI = patchStart;
 588                 faceI < (patchStart + patchSize);
 589                 ++faceI
 590             )
 591             {
 592                 const face& f = meshFaces[faceI];
 593
 594                 label nTri = f.nTriangles(points);
 595                 faceList triFaces;
 596
 597                 // triangulate polygons, but not quads
 598                 if (nTri <= 2)
 599                 {
 600                     triFaces.setSize(1);
 601                     triFaces[0] = f;
 602                 }
 603                 else
 604                 {
 605                     triFaces.setSize(nTri);
 606                     nTri = 0;
 607                     f.triangles(points, nTri, triFaces);
 608                 }
 609
 610                 forAll(triFaces, faceI)
 611                 {
 612                     const labelList& vrtList = triFaces[faceI];
 613
 614                     celFile
 615                         << cellId + 1 << " "
 616                         << shapeId << " "
 617                         << vrtList.size() << " "
 618                         << ctableId << " "
 619                         << typeId;
 620
 621                     // must be 3 (triangle) but could be quad
 622                     label count = 0;
 623                     forAll(vrtList, i)
 624                     {
 625                         if ((count % 8) == 0)
 626                         {
 627                             celFile
 628                                 << nl
 629                                 << "  " << cellId + 1;
 630                         }
 631                         // remember which points we'll need to write
 632                         pointHash.insert(vrtList[i]);
 633                         celFile << " " << vrtList[i] + 1;
 634                         count++;
 635                     }
 636                     celFile << endl;
 637
 638                     cellId++;
 639                 }
 640             }
 641         }
 642     }
 643     else
 644     {
 645         // cell Id is the OpenFOAM face Id
 646         // use the cellTable id from the face owner
 647         // - allows separation of parts
 648         forAll(patches, patchI)
 649         {
 650             label patchStart = patches[patchI].start();
 651             label patchSize  = patches[patchI].size();
 652
 653             for
 654             (
 655                 label faceI = patchStart;
 656                 faceI < (patchStart + patchSize);
 657                 ++faceI
 658             )
 659             {
 660                 const labelList& vrtList = meshFaces[faceI];
 661                 label cellId = faceI;
 662
 663                 celFile
 664                     << cellId + 1 << " "
 665                     << shapeId << " "
 666                     << vrtList.size() << " "
 667                     << cellTableId_[owner[faceI]] << " "
 668                     << typeId;
 669
 670                 // likely <= 8 vertices, but prevent overrun anyhow
 671                 label count = 0;
 672                 forAll(vrtList, i)
 673                 {
 674                     if ((count % 8) == 0)
 675                     {
 676                         celFile
 677                             << nl
 678                             << "  " << cellId + 1;
 679                     }
 680                     // remember which points we'll need to write
 681                     pointHash.insert(vrtList[i]);
 682                     celFile << " " << vrtList[i] + 1;
 683                     count++;
 684                 }
 685                 celFile << endl;
 686             }
 687         }
 688     }
 689
 690     OFstream vrtFile(baseName + ".vrt");
 691     writeHeader(vrtFile, "VERTEX");
 692
 693     vrtFile.precision(10);
 694     vrtFile.setf(std::ios::showpoint);  // force decimal point for Fortran
 695
 696     Info << "Writing " << vrtFile.name() << endl;
 697
 698     // build sorted table of contents
 699     SortableList<label> toc(pointHash.size());
 700     {
 701         label i = 0;
 702         forAllConstIter(labelHashSet, pointHash, iter)
 703         {
 704             toc[i++] = iter.key();
 705         }
 706     }
 707     toc.sort();
 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 // ************************************************************************* //