applications/utilities/mesh/conversion/cfx4ToFoam/cfx4ToFoam.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2009 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 Application
  26     cfx4ToFoam
  27
  28 Description
  29     Converts a CFX 4 mesh to FOAM format
  30
  31 \*---------------------------------------------------------------------------*/
  32
  33 #include "argList.H"
  34 #include "Time.H"
  35 #include "IFstream.H"
  36 #include "hexBlock.H"
  37 #include "polyMesh.H"
  38 #include "wallPolyPatch.H"
  39 #include "symmetryPolyPatch.H"
  40 #include "preservePatchTypes.H"
  41 #include "cellShape.H"
  42 #include "cellModeller.H"
  43
  44 using namespace Foam;
  45
  46 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  47 // Main program:
  48
  49 int main(int argc, char *argv[])
  50 {
  51     argList::noParallel();
  52     argList::validArgs.append("CFX geom file");
  53     argList::validOptions.insert("scale", "scale factor");
  54
  55     argList args(argc, argv);
  56
  57     if (!args.check())
  58     {
  59          FatalError.exit();
  60     }
  61
  62     scalar scaleFactor = 1.0;
  63     args.optionReadIfPresent("scale", scaleFactor);
  64
  65 #   include "createTime.H"
  66
  67     IFstream cfxFile(args.additionalArgs()[0]);
  68
  69     // Read the cfx information using a fixed format reader.
  70     // Comments in the file are in C++ style, so the stream parser will remove
  71     // them with no intervention
  72     label nblock, npatch, nglue, nelem, npoint;
  73
  74     cfxFile >> nblock >> npatch >> nglue >> nelem >> npoint;
  75
  76     Info << "Reading blocks" << endl;
  77
  78     PtrList<hexBlock> blocks(nblock);
  79
  80     {
  81         word blockName;
  82         label nx, ny, nz;
  83
  84         forAll (blocks, blockI)
  85         {
  86             cfxFile >> blockName;
  87             cfxFile >> nx >> ny >> nz;
  88
  89             blocks.set(blockI, new hexBlock(nx, ny, nz));
  90         }
  91     }
  92
  93     Info << "Reading patch definitions" << endl;
  94
  95     wordList cfxPatchTypes(npatch);
  96     wordList cfxPatchNames(npatch);
  97     labelList patchMasterBlocks(npatch);
  98     labelList patchDirections(npatch);
  99     labelListList patchRanges(npatch);
 100
 101     {
 102         label no, blkNo, patchLabel;
 103
 104         forAll (cfxPatchTypes, patchI)
 105         {
 106             // Grab patch type and name
 107             cfxFile >> cfxPatchTypes[patchI] >> cfxPatchNames[patchI] >> no;
 108
 109             // Grab patch range
 110             patchRanges[patchI].setSize(6);
 111             labelList& curRange = patchRanges[patchI];
 112
 113             forAll (curRange, rI)
 114             {
 115                 cfxFile >> curRange[rI];
 116             }
 117
 118             // Grab patch direction and master block ID
 119             // Note: direc is the direction, from the cfx manual
 120             // 0 = solid (3-D patch),
 121             // 1 = high i, 2 = high j, 3 = high k
 122             // 4 = low i, 5 = low j, 6 = low k
 123             cfxFile >> patchDirections[patchI] >> blkNo >> patchLabel;
 124
 125             patchMasterBlocks[patchI] = blkNo - 1;
 126         }
 127     }
 128
 129     Info << "Reading block glueing information" << endl;
 130
 131     labelList glueMasterPatches(nglue, -1);
 132     labelList glueSlavePatches(nglue, -1);
 133
 134     {
 135         label masterPatch, slavePatch;
 136         label dirIndex1, dirIndex2, dirIndex3, joinNumber;
 137
 138         for (label glueI = 0; glueI < nglue; glueI++)
 139         {
 140             cfxFile >> masterPatch >> slavePatch;
 141             cfxFile >> dirIndex1 >> dirIndex2 >> dirIndex3 >> joinNumber;
 142
 143             glueMasterPatches[glueI] = masterPatch - 1;
 144             glueSlavePatches[glueI] = slavePatch - 1;
 145         }
 146     }
 147
 148     Info << "Reading block points" << endl;
 149
 150     forAll (blocks, blockI)
 151     {
 152         Info << "block " << blockI << " is a ";
 153         blocks[blockI].readPoints(cfxFile);
 154     }
 155
 156     Info << "Calculating block offsets" << endl;
 157
 158     labelList blockOffsets(nblock, -1);
 159
 160     blockOffsets[0] = 0;
 161
 162     label nMeshPoints = blocks[0].nBlockPoints();
 163     label nMeshCells = blocks[0].nBlockCells();
 164
 165     for (label blockI = 1; blockI < nblock; blockI++)
 166     {
 167         nMeshPoints += blocks[blockI].nBlockPoints();
 168         nMeshCells +=  blocks[blockI].nBlockCells();
 169
 170         blockOffsets[blockI] =
 171             blockOffsets[blockI - 1]
 172           + blocks[blockI - 1].nBlockPoints();
 173     }
 174
 175     Info << "Assembling patches" << endl;
 176
 177     faceListList rawPatches(npatch);
 178
 179     forAll (rawPatches, patchI)
 180     {
 181         const word& patchType = cfxPatchTypes[patchI];
 182
 183         // reject volume patches
 184         if
 185         (
 186             patchType == "POROUS" || patchType == "SOLID"
 187          || patchType == "SOLCON" || patchType == "USER3D"
 188         )
 189         {
 190             patchMasterBlocks[patchI] = -1;
 191             rawPatches[patchI].setSize(0);
 192         }
 193         else
 194         {
 195             // read and create a 2-D patch
 196             rawPatches[patchI] =
 197                 blocks[patchMasterBlocks[patchI]].patchFaces
 198                 (
 199                     patchDirections[patchI],
 200                     patchRanges[patchI]
 201                 );
 202
 203         }
 204     }
 205
 206     Info << "Merging points ";
 207
 208     labelList pointMergeList(nMeshPoints, -1);
 209
 210     // In order to ensure robust merging, it is necessary to traverse
 211     // the patch glueing list until the pointMergeList stops changing.
 212     //
 213
 214     // For efficiency, create merge pairs in the first pass
 215     labelListListList glueMergePairs(glueMasterPatches.size());
 216
 217     forAll (glueMasterPatches, glueI)
 218     {
 219         const label masterPatch = glueMasterPatches[glueI];
 220         const label slavePatch = glueSlavePatches[glueI];
 221
 222         const label blockPlabel = patchMasterBlocks[masterPatch];
 223         const label blockNlabel = patchMasterBlocks[slavePatch];
 224
 225         const pointField& blockPpoints = blocks[blockPlabel].points();
 226         const pointField& blockNpoints = blocks[blockNlabel].points();
 227
 228         const faceList& blockPFaces = rawPatches[masterPatch];
 229         const faceList& blockNFaces = rawPatches[slavePatch];
 230
 231         labelListList& curPairs = glueMergePairs[glueI];
 232         curPairs.setSize(blockPFaces.size());
 233
 234         if (blockPFaces.size() != blockNFaces.size())
 235         {
 236             FatalErrorIn(args.executable())
 237                 << "Inconsistent number of faces for glue pair "
 238                 << glueI << " between blocks " << blockPlabel + 1
 239                 << " and " << blockNlabel + 1
 240                 << abort(FatalError);
 241         }
 242
 243         // Calculate sqr of the merge tolerance as 1/10th of the min
 244         // sqr point to point distance on the block face.  This is an
 245         // N^2 algorithm, sorry but I cannot quickly come up with
 246         // something better.
 247
 248         scalar sqrMergeTol = GREAT;
 249
 250         forAll (blockPFaces, blockPFaceLabel)
 251         {
 252             const labelList& blockPFacePoints =
 253                 blockPFaces[blockPFaceLabel];
 254
 255             forAll (blockPFacePoints, blockPFacePointI)
 256             {
 257                 forAll (blockPFacePoints, blockPFacePointI2)
 258                 {
 259                     if (blockPFacePointI != blockPFacePointI2)
 260                     {
 261                         sqrMergeTol =
 262                             min
 263                             (
 264                                 sqrMergeTol,
 265                                 magSqr
 266                                 (
 267                                     blockPpoints
 268                                         [blockPFacePoints[blockPFacePointI]]
 269                                   - blockPpoints
 270                                         [blockPFacePoints[blockPFacePointI2]]
 271                                 )
 272                             );
 273                     }
 274                 }
 275             }
 276         }
 277
 278         sqrMergeTol /= 10.0;
 279
 280         register bool found = false;
 281
 282         // N-squared point search over all points of all faces of
 283         // master block over all point of all faces of slave block
 284         forAll (blockPFaces, blockPFaceLabel)
 285         {
 286             const labelList& blockPFacePoints =
 287                 blockPFaces[blockPFaceLabel];
 288
 289             labelList& cp = curPairs[blockPFaceLabel];
 290             cp.setSize(blockPFacePoints.size());
 291
 292         forAll (blockPFacePoints, blockPFacePointI)
 293         {
 294             found = false;
 295
 296             forAll (blockNFaces, blockNFaceLabel)
 297             {
 298                 const labelList& blockNFacePoints =
 299                     blockNFaces[blockNFaceLabel];
 300
 301             forAll (blockNFacePoints, blockNFacePointI)
 302             {
 303                 if
 304                 (
 305                     magSqr
 306                     (
 307                         blockPpoints
 308                             [blockPFacePoints[blockPFacePointI]]
 309                       - blockNpoints
 310                             [blockNFacePoints[blockNFacePointI]]
 311                     )
 312                   < sqrMergeTol
 313                 )
 314                 {
 315                     // Found a new pair
 316                     found = true;
 317
 318                     cp[blockPFacePointI] =
 319                         blockNFacePoints[blockNFacePointI];
 320
 321                     label PpointLabel =
 322                         blockPFacePoints[blockPFacePointI]
 323                       + blockOffsets[blockPlabel];
 324
 325                     label NpointLabel =
 326                         blockNFacePoints[blockNFacePointI]
 327                       + blockOffsets[blockNlabel];
 328
 329                     label minPN = min(PpointLabel, NpointLabel);
 330
 331                     if (pointMergeList[PpointLabel] != -1)
 332                     {
 333                         minPN = min(minPN, pointMergeList[PpointLabel]);
 334                     }
 335
 336                     if (pointMergeList[NpointLabel] != -1)
 337                     {
 338                         minPN = min(minPN, pointMergeList[NpointLabel]);
 339                     }
 340
 341                     pointMergeList[PpointLabel]
 342                   = pointMergeList[NpointLabel]
 343                   = minPN;
 344
 345                     break;
 346                 }
 347             }
 348             if (found) break;
 349             }
 350         }
 351         }
 352     }
 353
 354
 355     register bool changedPointMerge = false;
 356     label nPasses = 0;
 357
 358     do
 359     {
 360         changedPointMerge = false;
 361         nPasses++;
 362
 363         forAll (glueMasterPatches, glueI)
 364         {
 365             const label masterPatch = glueMasterPatches[glueI];
 366             const label slavePatch = glueSlavePatches[glueI];
 367
 368             const label blockPlabel = patchMasterBlocks[masterPatch];
 369             const label blockNlabel = patchMasterBlocks[slavePatch];
 370
 371             const faceList& blockPFaces = rawPatches[masterPatch];
 372
 373             const labelListList& curPairs = glueMergePairs[glueI];
 374
 375             forAll (blockPFaces, blockPFaceLabel)
 376             {
 377                 const labelList& blockPFacePoints =
 378                     blockPFaces[blockPFaceLabel];
 379
 380                 const labelList& cp = curPairs[blockPFaceLabel];
 381
 382                 forAll (cp, blockPFacePointI)
 383                 {
 384                     label PpointLabel =
 385                         blockPFacePoints[blockPFacePointI]
 386                       + blockOffsets[blockPlabel];
 387
 388                     label NpointLabel =
 389                         cp[blockPFacePointI]
 390                       + blockOffsets[blockNlabel];
 391
 392                     if
 393                     (
 394                         pointMergeList[PpointLabel]
 395                      != pointMergeList[NpointLabel]
 396                     )
 397                     {
 398                         changedPointMerge = true;
 399
 400                         pointMergeList[PpointLabel]
 401                       = pointMergeList[NpointLabel]
 402                       = min
 403                         (
 404                             pointMergeList[PpointLabel],
 405                             pointMergeList[NpointLabel]
 406                         );
 407                     }
 408                 }
 409             }
 410         }
 411         Info << "." << flush;
 412     }
 413     while (changedPointMerge && nPasses < 8);
 414     Info << endl;
 415
 416     if (changedPointMerge == true)
 417     {
 418         FatalErrorIn(args.executable())
 419             << "Point merging failed after max number of passes."
 420             << abort(FatalError);
 421     }
 422
 423
 424     forAll (glueMasterPatches, glueI)
 425     {
 426         const label masterPatch = glueMasterPatches[glueI];
 427         const label slavePatch = glueSlavePatches[glueI];
 428
 429         const label blockPlabel = patchMasterBlocks[masterPatch];
 430         const label blockNlabel = patchMasterBlocks[slavePatch];
 431
 432         const faceList& blockPFaces = rawPatches[masterPatch];
 433         const faceList& blockNFaces = rawPatches[slavePatch];
 434
 435
 436         forAll (blockPFaces, blockPFaceLabel)
 437         {
 438             const labelList& blockPFacePoints
 439                 = blockPFaces[blockPFaceLabel];
 440
 441             forAll (blockPFacePoints, blockPFacePointI)
 442             {
 443                 label PpointLabel =
 444                     blockPFacePoints[blockPFacePointI]
 445                   + blockOffsets[blockPlabel];
 446
 447                 if (pointMergeList[PpointLabel] == -1)
 448                 {
 449                     FatalErrorIn(args.executable())
 450                         << "Unable to merge point " << blockPFacePointI
 451                         << " of face " << blockPFaceLabel
 452                         << " of block " << blockPlabel
 453                         << abort(FatalError);
 454                 }
 455             }
 456         }
 457
 458         forAll (blockNFaces, blockNFaceLabel)
 459         {
 460             const labelList& blockNFacePoints
 461                 = blockNFaces[blockNFaceLabel];
 462
 463             forAll (blockNFacePoints, blockNFacePointI)
 464             {
 465                 label NpointLabel =
 466                     blockNFacePoints[blockNFacePointI]
 467                   + blockOffsets[blockNlabel];
 468
 469                 if (pointMergeList[NpointLabel] == -1)
 470                 {
 471                     FatalErrorIn(args.executable())
 472                         << "Unable to merge point " << blockNFacePointI
 473                         << " of face " << blockNFaceLabel
 474                         << " of block " << blockNlabel
 475                         << abort(FatalError);
 476                 }
 477             }
 478         }
 479     }
 480
 481
 482     // sort merge list to return new point label (in new shorter list)
 483     // given old point label
 484     label nNewPoints = 0;
 485
 486     forAll (pointMergeList, pointLabel)
 487     {
 488         if (pointMergeList[pointLabel] > pointLabel)
 489         {
 490             FatalErrorIn(args.executable())
 491                 << "ouch" << abort(FatalError);
 492         }
 493
 494         if
 495         (
 496             (pointMergeList[pointLabel] == -1)
 497          || pointMergeList[pointLabel] == pointLabel
 498         )
 499         {
 500             pointMergeList[pointLabel] = nNewPoints;
 501             nNewPoints++;
 502         }
 503         else
 504         {
 505             pointMergeList[pointLabel] =
 506                 pointMergeList[pointMergeList[pointLabel]];
 507         }
 508     }
 509
 510     nMeshPoints = nNewPoints;
 511
 512     Info << "Creating points" << endl;
 513
 514     pointField points(nMeshPoints);
 515
 516     forAll (blocks, blockI)
 517     {
 518         const pointField& blockPoints = blocks[blockI].points();
 519
 520         forAll (blockPoints, blockPointLabel)
 521         {
 522             points
 523             [
 524                 pointMergeList
 525                 [
 526                     blockPointLabel
 527                   + blockOffsets[blockI]
 528                 ]
 529             ] = blockPoints[blockPointLabel];
 530         }
 531     }
 532
 533     // Scale the points
 534     if (scaleFactor > 1.0 + SMALL || scaleFactor < 1.0 - SMALL)
 535     {
 536         points *= scaleFactor;
 537     }
 538
 539     Info << "Creating cells" << endl;
 540
 541     cellShapeList cellShapes(nMeshCells);
 542
 543     const cellModel& hex = *(cellModeller::lookup("hex"));
 544
 545     label nCreatedCells = 0;
 546
 547     forAll (blocks, blockI)
 548     {
 549         labelListList curBlockCells = blocks[blockI].blockCells();
 550
 551         forAll (curBlockCells, blockCellI)
 552         {
 553             labelList cellPoints(curBlockCells[blockCellI].size());
 554
 555             forAll (cellPoints, pointI)
 556             {
 557                 cellPoints[pointI] =
 558                     pointMergeList
 559                     [
 560                         curBlockCells[blockCellI][pointI]
 561                       + blockOffsets[blockI]
 562                     ];
 563             }
 564
 565             cellShapes[nCreatedCells] = cellShape(hex, cellPoints);
 566
 567             nCreatedCells++;
 568         }
 569     }
 570
 571     Info << "Creating boundary patches" << endl;
 572
 573     faceListList boundary(npatch);
 574     wordList patchNames(npatch);
 575     wordList patchTypes(npatch);
 576     word defaultFacesName = "defaultFaces";
 577     word defaultFacesType = wallPolyPatch::typeName;
 578     wordList patchPhysicalTypes(npatch);
 579
 580     label nCreatedPatches = 0;
 581
 582     forAll (rawPatches, patchI)
 583     {
 584         if (rawPatches[patchI].size() && cfxPatchTypes[patchI] != "BLKBDY")
 585         {
 586             // Check if this name has been already created
 587             label existingPatch = -1;
 588
 589             for (label oldPatchI = 0; oldPatchI < nCreatedPatches; oldPatchI++)
 590             {
 591                 if (patchNames[oldPatchI] == cfxPatchNames[patchI])
 592                 {
 593                     existingPatch = oldPatchI;
 594                     break;
 595                 }
 596             }
 597
 598             const faceList& curRawPatch = rawPatches[patchI];
 599             label curBlock = patchMasterBlocks[patchI];
 600
 601             if (existingPatch >= 0)
 602             {
 603                 Info << "CFX patch " << patchI
 604                     << ", of type " << cfxPatchTypes[patchI]
 605                     << ", name " << cfxPatchNames[patchI]
 606                     << " already exists as FOAM patch " << existingPatch
 607                     << ".  Adding faces." << endl;
 608
 609                 faceList& renumberedPatch = boundary[existingPatch];
 610                 label oldSize = renumberedPatch.size();
 611                 renumberedPatch.setSize(oldSize + curRawPatch.size());
 612
 613                 forAll (curRawPatch, faceI)
 614                 {
 615                     const face& oldFace = curRawPatch[faceI];
 616
 617                     face& newFace = renumberedPatch[oldSize + faceI];
 618                     newFace.setSize(oldFace.size());
 619
 620                     forAll (oldFace, pointI)
 621                     {
 622                         newFace[pointI] =
 623                             pointMergeList
 624                             [
 625                                 oldFace[pointI]
 626                               + blockOffsets[curBlock]
 627                             ];
 628                     }
 629                 }
 630             }
 631             else
 632             {
 633                 // Real patch to be created
 634                 faceList& renumberedPatch = boundary[nCreatedPatches];
 635                 renumberedPatch.setSize(curRawPatch.size());
 636
 637                 forAll (curRawPatch, faceI)
 638                 {
 639                     const face& oldFace = curRawPatch[faceI];
 640
 641                     face& newFace = renumberedPatch[faceI];
 642                     newFace.setSize(oldFace.size());
 643
 644                     forAll (oldFace, pointI)
 645                     {
 646                         newFace[pointI] =
 647                             pointMergeList
 648                             [
 649                                 oldFace[pointI]
 650                               + blockOffsets[curBlock]
 651                             ];
 652                     }
 653                 }
 654
 655                 Info << "CFX patch " << patchI
 656                     << ", of type " << cfxPatchTypes[patchI]
 657                     << ", name " << cfxPatchNames[patchI]
 658                     << " converted into FOAM patch " << nCreatedPatches
 659                     << " type ";
 660
 661                 if (cfxPatchTypes[patchI] == "WALL")
 662                 {
 663                     Info << "wall." << endl;
 664
 665                     patchTypes[nCreatedPatches] = wallPolyPatch::typeName;
 666                     patchNames[nCreatedPatches] = cfxPatchNames[patchI];
 667                     nCreatedPatches++;
 668                 }
 669                 else if (cfxPatchTypes[patchI] == "SYMMET")
 670                 {
 671                     Info << "symmetryPlane." << endl;
 672
 673                     patchTypes[nCreatedPatches] = symmetryPolyPatch::typeName;
 674                     patchNames[nCreatedPatches] = cfxPatchNames[patchI];
 675                     nCreatedPatches++;
 676                 }
 677                 else if
 678                 (
 679                     cfxPatchTypes[patchI] == "INLET"
 680                  || cfxPatchTypes[patchI] == "OUTLET"
 681                  || cfxPatchTypes[patchI] == "PRESS"
 682                  || cfxPatchTypes[patchI] == "CNDBDY"
 683                  || cfxPatchTypes[patchI] == "USER2D"
 684                 )
 685                 {
 686                     Info << "generic." << endl;
 687
 688                     patchTypes[nCreatedPatches] = polyPatch::typeName;
 689                     patchNames[nCreatedPatches] = cfxPatchNames[patchI];
 690                     nCreatedPatches++;
 691                 }
 692                 else
 693                 {
 694                     FatalErrorIn(args.executable())
 695                         << "Unrecognised CFX patch type "
 696                         << cfxPatchTypes[patchI]
 697                         << abort(FatalError);
 698                 }
 699             }
 700         }
 701     }
 702
 703     boundary.setSize(nCreatedPatches);
 704     patchTypes.setSize(nCreatedPatches);
 705     patchNames.setSize(nCreatedPatches);
 706
 707     preservePatchTypes
 708     (
 709         runTime,
 710         runTime.constant(),
 711         polyMesh::defaultRegion,
 712         patchNames,
 713         patchTypes,
 714         defaultFacesName,
 715         defaultFacesType,
 716         patchPhysicalTypes
 717     );
 718
 719     polyMesh pShapeMesh
 720     (
 721         IOobject
 722         (
 723             polyMesh::defaultRegion,
 724             runTime.constant(),
 725             runTime
 726         ),
 727         xferMove(points),
 728         cellShapes,
 729         boundary,
 730         patchNames,
 731         patchTypes,
 732         defaultFacesName,
 733         defaultFacesType,
 734         patchPhysicalTypes
 735     );
 736
 737     // Set the precision of the points data to 10
 738     IOstream::defaultPrecision(10);
 739
 740     Info << "Writing polyMesh" << endl;
 741     pShapeMesh.write();
 742
 743     Info << "End\n" << endl;
 744
 745     return 0;
 746 }
 747
 748
 749 // ************************************************************************* //