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