src/conversion/meshWriter/elmer/ElmerMeshWriter.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 "ElmerMeshWriter.H"
  27
  28 #include "foamTime.H"
  29 #include "SortableList.H"
  30 #include "OFstream.H"
  31
  32 #include "hexMatcher.H"
  33 #include "wedgeMatcher.H"
  34 #include "prismMatcher.H"
  35 #include "pyrMatcher.H"
  36 #include "tetWedgeMatcher.H"
  37 #include "tetMatcher.H"
  38
  39
  40 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  41
  42 void Foam::meshWriters::Elmer::getCellTable()
  43 {
  44     // read constant/polyMesh/propertyName
  45     IOList<label> ioList
  46     (
  47         IOobject
  48         (
  49             "cellTableId",
  50             "constant",
  51             polyMesh::meshSubDir,
  52             mesh_,
  53             IOobject::READ_IF_PRESENT,
  54             IOobject::NO_WRITE,
  55             false
  56         )
  57     );
  58
  59     bool useCellZones = false;
  60     cellTableId_.setSize(mesh_.nCells(), -1);
  61
  62     // get information from constant/polyMesh/cellTableId if possible
  63     if (ioList.headerOk())
  64     {
  65         if (ioList.size() == mesh_.nCells())
  66         {
  67             cellTableId_.transfer(ioList);
  68
  69             if (cellTable_.empty())
  70             {
  71                 Info<< "No cellTable information available" << endl;
  72             }
  73         }
  74         else
  75         {
  76             WarningIn("Elmer::getCellTable()")
  77                 << ioList.objectPath() << " has incorrect number of cells "
  78                 << " - use cellZone information"
  79                 << endl;
  80
  81             ioList.clear();
  82             useCellZones = true;
  83         }
  84     }
  85     else
  86     {
  87         useCellZones = true;
  88     }
  89
  90
  91     if (useCellZones)
  92     {
  93         if (cellTable_.empty())
  94         {
  95             Info<< "Created cellTable from cellZones" << endl;
  96             cellTable_ = mesh_;
  97         }
  98
  99         // track if there are unzoned cells
 100         label nUnzoned = mesh_.nCells();
 101
 102         // get the cellZone <-> cellTable correspondence
 103         Info<< "Matching cellZones to cellTable" << endl;
 104
 105         forAll (mesh_.cellZones(), zoneI)
 106         {
 107             const cellZone& cZone = mesh_.cellZones()[zoneI];
 108             if (cZone.size())
 109             {
 110                 nUnzoned -= cZone.size();
 111
 112                 label tableId = cellTable_.findIndex(cZone.name());
 113                 if (tableId < 0)
 114                 {
 115                     dictionary dict;
 116
 117                     dict.add("Label", cZone.name());
 118                     dict.add("MaterialType", "fluid");
 119                     tableId = cellTable_.append(dict);
 120                 }
 121
 122                 forAll (cZone, i)
 123                 {
 124                     cellTableId_[cZone[i]] = tableId;
 125                 }
 126             }
 127         }
 128
 129         if (nUnzoned)
 130         {
 131             dictionary dict;
 132
 133             dict.add("Label", "__unZonedCells__");
 134             dict.add("MaterialType", "fluid");
 135             label tableId = cellTable_.append(dict);
 136
 137             forAll (cellTableId_, i)
 138             {
 139                 if (cellTableId_[i] < 0)
 140                 {
 141                     cellTableId_[i] = tableId;
 142                 }
 143             }
 144         }
 145     }
 146 }
 147
 148
 149 Foam::label Foam::meshWriters::Elmer::getFaceType
 150 (
 151     const label nvert,
 152     const word& zname
 153 ) const
 154 {
 155     switch(nvert)
 156     {
 157         case 3:
 158             return ELMER_ETYPE_TRIA;
 159
 160         case 4:
 161             return ELMER_ETYPE_QUAD;
 162
 163         default:
 164             Info<< "Found a bad face with " << nvert
 165                 << " vertices on zone " << zname
 166                 << "." << endl;
 167             return ELMER_ETYPE_BAD;
 168     }
 169 }
 170
 171
 172 void Foam::meshWriters::Elmer::writeNames() const
 173 {
 174     OFstream os("mesh.names");
 175
 176     const cellZoneMesh& czones = mesh_.cellZones();
 177     const faceZoneMesh& fzones = mesh_.faceZones();
 178     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 179     label boundaryID = 0;
 180
 181
 182     Info<< "Writing " << os.name() << "." << endl;
 183
 184     os << "! ----- Names for bodies -----" << nl;
 185     forAll (czones, zoneI)
 186     {
 187         os  << "$ " << czones[zoneI].name() << " = " << zoneI+1 << nl;
 188     }
 189
 190     os << "! ----- Names for exterior boundaries -----" << nl;
 191     forAll(patches, patchI)
 192     {
 193         os  << "$ " << patches[patchI].name() << " = " << ++boundaryID << nl;
 194     }
 195
 196     os << "! ----- Names for interior boundaries -----" << nl;
 197     forAll(fzones, fzoneI)
 198     {
 199         if (!faceZoneExcludePattern.match(fzones[fzoneI].name()))
 200         {
 201             os  << "$ " << fzones[fzoneI].name() << " = " << ++boundaryID << nl;
 202         }
 203     }
 204 }
 205
 206
 207 bool Foam::meshWriters::Elmer::writeHeader() const
 208 {
 209     const pointField& points = mesh_.points();
 210     const cellList& cells = mesh_.cells();
 211     const faceList& faces = mesh_.faces();
 212     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 213     const faceZoneMesh& fzones = mesh_.faceZones();
 214     label extZones = 0; // Count exterior boundary zones
 215     label extFaces = 0; // Count exterior boundary faces
 216     label intZones = 0; // Used interior boundary zones
 217     label intFaces = 0; // Used interior boundary faces
 218
 219     // Count the different cell types
 220     hexMatcher hex;
 221     prismMatcher pri;
 222     pyrMatcher pyr;
 223     tetMatcher tet;
 224
 225     label nHex = 0;
 226     label nPri = 0;
 227     label nPyr = 0;
 228     label nTet = 0;
 229     label nTri = 0;
 230     label nQua = 0;
 231     label nBad = 0;
 232
 233
 234     // Count all exterior boundary zone + faces + the types
 235     forAll(patches, patchI)
 236     {
 237         const label patchStart = patches[patchI].start();
 238         const label patchEnd = patchStart + patches[patchI].size();
 239
 240         extZones++;
 241         extFaces += patches[patchI].size();
 242
 243         // Count triangles and quads
 244         for(label faceI = patchStart; faceI < patchEnd; ++faceI)
 245         {
 246             if (faces[faceI].size() == 4)
 247             {
 248                nQua++;
 249             }
 250             else
 251             {
 252                nTri++;
 253             }
 254         }
 255     }
 256
 257     // Count all interior boundary zones + faces with names not matching
 258     // the exclude pattern
 259     forAll(fzones, fzoneI)
 260     {
 261         if (faceZoneExcludePattern.match(fzones[fzoneI].name()))
 262         {
 263            continue;
 264         }
 265
 266         const faceZone& fzone = fzones[fzoneI];
 267
 268         intZones++;
 269         intFaces += fzone.size();
 270
 271         // Count triangles and quads
 272         forAll(fzone, i)
 273         {
 274             if (faces[fzone[i]].size() == 4)
 275             {
 276                 nQua++;
 277             }
 278             else
 279             {
 280                 nTri++;
 281             }
 282         }
 283     }
 284
 285     // Count the volume element types
 286     for(label cellI = 0; cellI < mesh_.nCells(); cellI++)
 287     {
 288         if (hex.isA(mesh_, cellI))
 289         {
 290             nHex++;
 291         }
 292         else if (tet.isA(mesh_, cellI))
 293         {
 294             nTet++;
 295         }
 296         else if (pyr.isA(mesh_, cellI))
 297         {
 298             nPyr++;
 299         }
 300         else if (pri.isA(mesh_, cellI))
 301         {
 302             nPri++;
 303         }
 304         else
 305         {
 306             nBad++; // Not a valid element type for Elmer
 307         }
 308     }
 309
 310     OFstream os("mesh.header");
 311
 312     os << points.size() << " "
 313        << cells.size()  << " "
 314        << extFaces+intFaces << nl;
 315
 316     Info<< "Writing " << os.name() << "." << endl;
 317
 318     // Save the number of different element types used
 319     os  << (nTet>0) + (nPri>0) + (nPyr>0) +
 320            (nHex>0) + (nQua>0) + (nTri>0) + (nBad>0) << endl;
 321
 322     // Save type and count if count > 0
 323     if (nTet > 0) os << ELMER_ETYPE_TET   << " " << nTet << endl;
 324     if (nPri > 0) os << ELMER_ETYPE_PRISM << " " << nPri << endl;
 325     if (nPyr > 0) os << ELMER_ETYPE_PYRAM << " " << nPyr << endl;
 326     if (nHex > 0) os << ELMER_ETYPE_HEX   << " " << nHex << endl;
 327     if (nQua > 0) os << ELMER_ETYPE_QUAD  << " " << nQua << endl;
 328     if (nTri > 0) os << ELMER_ETYPE_TRIA  << " " << nTri << endl;
 329     if (nBad > 0) os << ELMER_ETYPE_BAD   << " " << nBad << endl;
 330
 331     Info<< "Mesh statistics:" << nl
 332         << "    Nodes          : " << points.size()  << nl
 333         << "    Elements       : " << cells.size() << nl
 334         << "        Hexahedra  : " << nHex << nl
 335         << "        Wedges     : " << nPri << nl
 336         << "        Pyramids   : " << nPyr << nl
 337         << "        Tetrahedra : " << nTet << nl
 338         << "        Polyhedra  : " << nBad << nl
 339         << "    Regions        : " << mesh_.cellZones().size() << nl
 340         << "    Ext. boundaries: " << extZones << ", "
 341         << extFaces << " faces" << nl
 342         << "    Int. boundaries: " << intZones << ", "
 343         << intFaces << " faces" << nl
 344         << "        Tri faces  : " << nTri << nl
 345         << "        Quad faces : " << nQua << nl
 346         << endl;
 347
 348     return !nBad;
 349 }
 350
 351
 352 void Foam::meshWriters::Elmer::writeNodes() const
 353 {
 354     OFstream os("mesh.nodes");
 355
 356     // Set the precision of the points data to 10
 357     os.precision(10);
 358
 359     // Force decimal point for Fortran90 input style
 360     os.setf(std::ios::showpoint);
 361
 362     const pointField& points = mesh_.points();
 363
 364     Info<< "Writing " << os.name() << ": scale=" << scaleFactor_ << endl;
 365
 366     forAll(points, ptI)
 367     {
 368         os << ptI + 1 << " -1 "
 369            << scaleFactor_ * points[ptI].x() << " "
 370            << scaleFactor_ * points[ptI].y() << " "
 371            << scaleFactor_ * points[ptI].z() << nl;
 372     }
 373 }
 374
 375
 376 void Foam::meshWriters::Elmer::writeElements() const
 377 {
 378     OFstream os("mesh.elements");
 379
 380     // map foam cellModeller index -> Elmer element types
 381     Map<label> shapeLookupIndex;
 382     shapeLookupIndex.insert(tetModel->index(), ELMER_ETYPE_TET  );
 383     shapeLookupIndex.insert(prismModel->index(), ELMER_ETYPE_PRISM);
 384     shapeLookupIndex.insert(pyrModel->index(), ELMER_ETYPE_PYRAM);
 385     shapeLookupIndex.insert(hexModel->index(), ELMER_ETYPE_HEX  );
 386
 387     const cellShapeList& shapes = mesh_.cellShapes();
 388     const cellList& cells  = mesh_.cells();
 389
 390     Info<< "Writing " << os.name() << "." << endl;
 391
 392     forAll(cells, cellId)
 393     {
 394         const cellShape& shape = shapes[cellId];
 395         label mapIndex = shape.model().index();
 396
 397         // a registered primitive type
 398         if (shapeLookupIndex.found(mapIndex))
 399         {
 400             os << cellId+1 << " "
 401                << cellTableId_[cellId] << " "
 402                << shapeLookupIndex[mapIndex];
 403
 404             const labelList& vrtList = shapes[cellId];
 405             forAll(vrtList, i)
 406             {
 407                 os << " " << vrtList[i] + 1;
 408             }
 409             os << endl;
 410         }
 411         else
 412         {
 413             Info<< "***WARNING: polyhedron " << cellId << " ignored." << endl;
 414         }
 415     }
 416 }
 417
 418
 419 void Foam::meshWriters::Elmer::writeBoundary() const
 420 {
 421     OFstream os("mesh.boundary");
 422
 423     const faceList&  faces = mesh_.faces();
 424     const labelList& owner = mesh_.faceOwner();
 425     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 426     label bfaceID = 0;
 427     label boundaryID = 0;
 428
 429
 430     Info<< "Writing " << os.name() << "." << endl;
 431
 432     // Write all patches == exterior boundary faces
 433     forAll(patches, patchI)
 434     {
 435         const word& pname = patches[patchI].name();
 436         const label patchStart = patches[patchI].start();
 437         const label patchEnd = patchStart + patches[patchI].size();
 438
 439         boundaryID++;
 440         for(label faceI = patchStart; faceI < patchEnd; ++faceI)
 441         {
 442             const labelList& vrtList = faces[faceI];
 443
 444             os << ++bfaceID
 445                << " " << boundaryID
 446                << " " << owner[faceI] + 1  // ID of parent cell
 447                << " 0 "                    // 0 == exterior boundary face
 448                << getFaceType(vrtList.size(),pname);
 449
 450             forAll(vrtList, i)
 451             {
 452                 os << " " << vrtList[i] + 1;
 453             }
 454             os << endl;
 455         }
 456     }
 457
 458
 459     // Write face zones (== interior boundary faces) with names not matching
 460     // the exclude pattern
 461     const faceZoneMesh& fzones = mesh_.faceZones();
 462     forAll(fzones, fzoneI)
 463     {
 464         if(faceZoneExcludePattern.match(fzones[fzoneI].name()))
 465         {
 466            continue;
 467         }
 468
 469         const faceZone& fzone = fzones[fzoneI];
 470         const labelList& master = fzone.masterCells();
 471         const labelList& slave = fzone.slaveCells();
 472         const word& zname = fzone.name();
 473
 474         boundaryID++;
 475         forAll(fzone, i)
 476         {
 477             const labelList& vrtList = faces[fzone[i]];
 478
 479             os << ++bfaceID
 480                << " " << boundaryID
 481                << " " << master[i] + 1 // ID of parent cell
 482                << " " << slave[i]  + 1 // ID of neigbour cell
 483                << " " << getFaceType(vrtList.size(),zname);
 484
 485             forAll(vrtList, j)
 486             {
 487                 os << " " << vrtList[j] + 1;
 488             }
 489             os << endl;
 490         }
 491     }
 492 }
 493
 494
 495 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 496
 497 Foam::meshWriters::Elmer::Elmer
 498 (
 499     const polyMesh& mesh,
 500     const wordRe& excludePattern,
 501     const scalar scaleFactor
 502 )
 503 :
 504     meshWriter(mesh, scaleFactor),
 505     faceZoneExcludePattern(excludePattern)
 506 {
 507     boundaryRegion_.readDict(mesh_);
 508     cellTable_.readDict(mesh_);
 509     getCellTable();
 510 }
 511
 512
 513 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 514
 515 Foam::meshWriters::Elmer::~Elmer()
 516 {
 517 }
 518
 519
 520 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 521
 522 bool Foam::meshWriters::Elmer::write(const fileName& dirName) const
 523 {
 524     fileName baseName(dirName);
 525
 526     if (baseName.empty())
 527     {
 528         baseName = meshWriter::defaultMeshName;
 529
 530         if
 531         (
 532             mesh_.time().timeName() != "0"
 533          && mesh_.time().timeName() != "constant"
 534         )
 535         {
 536             baseName += "_" + mesh_.time().timeName();
 537         }
 538     }
 539
 540     // Create the mesh directory (elemerMesh), chdir into and cleanup
 541     mkDir(baseName);
 542     chDir(baseName);
 543     rm("mesh.header");
 544     rm("mesh.nodes");
 545     rm("mesh.elements");
 546     rm("mesh.boundary");
 547     rm("mesh.names");
 548
 549     bool success = writeHeader();
 550     if (success)
 551     {
 552         writeNodes();
 553         writeElements();
 554         writeBoundary();
 555         writeNames();
 556     }
 557     else
 558     {
 559         rm("mesh.header");
 560     }
 561
 562     chDir("..");
 563
 564     return success;
 565 }
 566
 567 // ************************************************************************* //