src/conversion/meshWriter/starcd/STARCDMeshWriter.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend 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 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "STARCDMeshWriter.H"
  27
  28 #include "foamTime.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 bool Foam::meshWriters::STARCD::writeSurface
 530 (
 531     const fileName& meshName,
 532     const bool& triangulate
 533 ) const
 534 {
 535     fileName baseName(meshName);
 536
 537     if (baseName.empty())
 538     {
 539         baseName = meshWriter::defaultSurfaceName;
 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 foam 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 // ************************************************************************* //