applications/utilities/mesh/conversion/kivaToFoam/readKivaGrid.H

   1 ifstream kivaFile(kivaFileName.c_str());
   2
   3 if (!kivaFile.good())
   4 {
   5     FatalErrorIn(args.executable())
   6         << "Cannot open file " << kivaFileName
   7         << exit(FatalError);
   8 }
   9
  10 Info << "Reading kiva grid from file " << kivaFileName << endl;
  11
  12
  13 char title[120];
  14 kivaFile.getline(title, 120, '\n');
  15
  16 label nPoints, nCells, nRegs;
  17 kivaFile >> nCells >> nPoints >> nRegs;
  18
  19 pointField points(nPoints);
  20 label i4;
  21 labelList idface(nPoints), fv(nPoints);
  22
  23 for(label i=0; i<nPoints; i++)
  24 {
  25     scalar ffv;
  26     kivaFile
  27         >> i4
  28         >> points[i].x() >> points[i].y() >> points[i].z()
  29         >> ffv;
  30
  31     if (kivaVersion == kiva3v)
  32     {
  33         kivaFile >> idface[i];
  34     }
  35
  36     // Convert from KIVA cgs units to SI
  37     points[i] *= 0.01;
  38
  39     fv[i] = label(ffv);
  40 }
  41
  42 labelList i1tab(nPoints), i3tab(nPoints), i8tab(nPoints), idreg(nPoints),
  43       f(nPoints), bcl(nPoints), bcf(nPoints), bcb(nPoints);
  44
  45 label nBfaces = 0;
  46
  47 for(label i=0; i<nPoints; i++)
  48 {
  49     label i1, i3, i8;
  50     scalar ff, fbcl, fbcf, fbcb;
  51
  52     kivaFile
  53         >> i1 >> i3 >> i4 >> i8
  54         >> ff >> fbcl >> fbcf >> fbcb
  55         >> idreg[i];
  56
  57     // Correct indices to start from 0
  58     i1tab[i] = i1 - 1;
  59     i3tab[i] = i3 - 1;
  60     i8tab[i] = i8 - 1;
  61
  62     // Convert scalar indices into hexLabels
  63     f[i] = label(ff);
  64     bcl[i] = label(fbcl);
  65     bcf[i] = label(fbcf);
  66     bcb[i] = label(fbcb);
  67
  68     if (bcl[i] > 0 && bcl[i] != 4) nBfaces++;
  69     if (bcf[i] > 0 && bcf[i] != 4) nBfaces++;
  70     if (bcb[i] > 0 && bcb[i] != 4) nBfaces++;
  71 }
  72
  73
  74 label mTable;
  75 kivaFile >> mTable;
  76
  77 if (mTable == 0)
  78 {
  79     FatalErrorIn(args.executable())
  80         << "mTable == 0. This is not supported."
  81            " kivaToFoam needs complete neighbour information"
  82         << exit(FatalError);
  83 }
  84
  85
  86 labelList imtab(nPoints), jmtab(nPoints), kmtab(nPoints);
  87
  88 for(label i=0; i<nPoints; i++)
  89 {
  90     label im, jm, km;
  91     kivaFile >> i4 >> im >> jm >> km;
  92
  93     // Correct indices to start from 0
  94     imtab[i] = im - 1;
  95     jmtab[i] = jm - 1;
  96     kmtab[i] = km - 1;
  97 }
  98
  99 Info << "Finished reading KIVA file" << endl;
 100
 101 cellShapeList cellShapes(nPoints);
 102
 103 labelList cellZoning(nPoints, -1);
 104
 105 const cellModel& hex = *(cellModeller::lookup("hex"));
 106 labelList hexLabels(8);
 107 label activeCells = 0;
 108
 109 // Create and set the collocated point collapse map
 110 labelList pointMap(nPoints);
 111
 112 forAll (pointMap, i)
 113 {
 114     pointMap[i] = i;
 115 }
 116
 117 // Initialise all cells to hex and search and set map for collocated points
 118 for(label i=0; i<nPoints; i++)
 119 {
 120     if (f[i] > 0.0)
 121     {
 122         hexLabels[0] = i;
 123         hexLabels[1] = i1tab[i];
 124         hexLabels[2] = i3tab[i1tab[i]];
 125         hexLabels[3] = i3tab[i];
 126         hexLabels[4] = i8tab[i];
 127         hexLabels[5] = i1tab[i8tab[i]];
 128         hexLabels[6] = i3tab[i1tab[i8tab[i]]];
 129         hexLabels[7] = i3tab[i8tab[i]];
 130
 131         cellShapes[activeCells] = cellShape(hex, hexLabels);
 132
 133         edgeList edges = cellShapes[activeCells].edges();
 134
 135         forAll (edges, ei)
 136         {
 137             if (edges[ei].mag(points) < SMALL)
 138             {
 139                 label start = pointMap[edges[ei].start()];
 140                 while (start != pointMap[start])
 141                 {
 142                     start = pointMap[start];
 143                 }
 144
 145                 label end = pointMap[edges[ei].end()];
 146                 while (end != pointMap[end])
 147                 {
 148                     end = pointMap[end];
 149                 }
 150
 151                 label minLabel = min(start, end);
 152
 153                 pointMap[start] = pointMap[end] = minLabel;
 154             }
 155         }
 156
 157         // Grab the cell zoning info
 158         cellZoning[activeCells] = idreg[i];
 159
 160         activeCells++;
 161     }
 162 }
 163
 164 // Contract list of cells to the active ones
 165 cellShapes.setSize(activeCells);
 166 cellZoning.setSize(activeCells);
 167
 168 // Map collocated points to refer to the same point and collapse cell shape
 169 // to the corresponding hex-degenerate.
 170 forAll (cellShapes, celli)
 171 {
 172     cellShape& cs = cellShapes[celli];
 173
 174     forAll (cs, i)
 175     {
 176         cs[i] = pointMap[cs[i]];
 177     }
 178
 179     cs.collapse();
 180 }
 181
 182 label bcIDs[11] = {-1, 0, 2, 4, -1, 5, -1, 6, 7, 8, 9};
 183
 184 const label nBCs = 12;
 185
 186 const word* kivaPatchTypes[nBCs] =
 187 {
 188     &wallPolyPatch::typeName,
 189     &wallPolyPatch::typeName,
 190     &wallPolyPatch::typeName,
 191     &wallPolyPatch::typeName,
 192     &symmetryPolyPatch::typeName,
 193     &wedgePolyPatch::typeName,
 194     &polyPatch::typeName,
 195     &polyPatch::typeName,
 196     &polyPatch::typeName,
 197     &polyPatch::typeName,
 198     &symmetryPolyPatch::typeName,
 199     &cyclicPolyPatch::typeName
 200 };
 201
 202 enum patchTypeNames
 203 {
 204     PISTON,
 205     VALVE,
 206     LINER,
 207     CYLINDERHEAD,
 208     AXIS,
 209     WEDGE,
 210     INFLOW,
 211     OUTFLOW,
 212     PRESIN,
 213     PRESOUT,
 214     SYMMETRYPLANE,
 215     CYCLIC
 216 };
 217
 218 const char* kivaPatchNames[nBCs] =
 219 {
 220     "piston",
 221     "valve",
 222     "liner",
 223     "cylinderHead",
 224     "axis",
 225     "wedge",
 226     "inflow",
 227     "outflow",
 228     "presin",
 229     "presout",
 230     "symmetryPlane",
 231     "cyclic"
 232 };
 233
 234
 235 List<SLList<face> > pFaces[nBCs];
 236
 237 face quadFace(4);
 238 face triFace(3);
 239
 240 for(label i=0; i<nPoints; i++)
 241 {
 242     if (f[i] > 0)
 243     {
 244         // left
 245         label bci = bcl[i];
 246         if (bci != 4)
 247         {
 248             quadFace[0] = i3tab[i];
 249             quadFace[1] = i;
 250             quadFace[2] = i8tab[i];
 251             quadFace[3] = i3tab[i8tab[i]];
 252
 253 #           include "checkPatch.H"
 254         }
 255
 256         // right
 257         bci = bcl[i1tab[i]];
 258         if (bci != 4)
 259         {
 260             quadFace[0] = i1tab[i];
 261             quadFace[1] = i3tab[i1tab[i]];
 262             quadFace[2] = i8tab[i3tab[i1tab[i]]];
 263             quadFace[3] = i8tab[i1tab[i]];
 264
 265 #           include "checkPatch.H"
 266         }
 267
 268         // front
 269         bci = bcf[i];
 270         if (bci != 4)
 271         {
 272             quadFace[0] = i;
 273             quadFace[1] = i1tab[i];
 274             quadFace[2] = i8tab[i1tab[i]];
 275             quadFace[3] = i8tab[i];
 276
 277 #           include "checkPatch.H"
 278         }
 279
 280         // derriere (back)
 281         bci = bcf[i3tab[i]];
 282         if (bci != 4)
 283         {
 284             quadFace[0] = i3tab[i];
 285             quadFace[1] = i8tab[i3tab[i]];
 286             quadFace[2] = i8tab[i1tab[i3tab[i]]];
 287             quadFace[3] = i1tab[i3tab[i]];
 288
 289 #           include "checkPatch.H"
 290         }
 291
 292         // bottom
 293         bci = bcb[i];
 294         if (bci != 4)
 295         {
 296             quadFace[0] = i;
 297             quadFace[1] = i3tab[i];
 298             quadFace[2] = i1tab[i3tab[i]];
 299             quadFace[3] = i1tab[i];
 300
 301 #           include "checkPatch.H"
 302         }
 303
 304         // top
 305         bci = bcb[i8tab[i]];
 306         if (bci != 4)
 307         {
 308             quadFace[0] = i8tab[i];
 309             quadFace[1] = i1tab[i8tab[i]];
 310             quadFace[2] = i3tab[i1tab[i8tab[i]]];
 311             quadFace[3] = i3tab[i8tab[i]];
 312
 313 #           include "checkPatch.H"
 314         }
 315     }
 316 }
 317
 318 // Tranfer liner faces that are above the minimum cylinder-head z height
 319 // into the cylinder-head region
 320 if
 321 (
 322     pFaces[LINER].size()
 323  && pFaces[LINER][0].size()
 324  && pFaces[CYLINDERHEAD].size()
 325  && pFaces[CYLINDERHEAD][0].size()
 326 )
 327 {
 328     scalar minz = GREAT;
 329
 330     for
 331     (
 332         SLList<face>::iterator iter = pFaces[CYLINDERHEAD][0].begin();
 333         iter != pFaces[CYLINDERHEAD][0].end();
 334         ++iter
 335     )
 336     {
 337         const face& pf = iter();
 338
 339         forAll (pf, pfi)
 340         {
 341             minz = min(minz, points[pf[pfi]].z());
 342         }
 343     }
 344
 345     minz += SMALL;
 346
 347     SLList<face> newLinerFaces;
 348
 349     for
 350     (
 351         SLList<face>::iterator iter = pFaces[LINER][0].begin();
 352         iter != pFaces[LINER][0].end();
 353         ++iter
 354     )
 355     {
 356         const face& pf = iter();
 357
 358         scalar minfz = GREAT;
 359         forAll (pf, pfi)
 360         {
 361             minfz = min(minfz, points[pf[pfi]].z());
 362         }
 363
 364         if (minfz > minz)
 365         {
 366             pFaces[CYLINDERHEAD][0].append(pf);
 367         }
 368         else
 369         {
 370             newLinerFaces.append(pf);
 371         }
 372     }
 373
 374     if (pFaces[LINER][0].size() != newLinerFaces.size())
 375     {
 376         Info<< "Transfered " << pFaces[LINER][0].size() - newLinerFaces.size()
 377             << " faces from liner region to cylinder head" << endl;
 378         pFaces[LINER][0] = newLinerFaces;
 379     }
 380
 381     SLList<face> newCylinderHeadFaces;
 382
 383     for
 384     (
 385         SLList<face>::iterator iter = pFaces[CYLINDERHEAD][0].begin();
 386         iter != pFaces[CYLINDERHEAD][0].end();
 387         ++iter
 388     )
 389     {
 390         const face& pf = iter();
 391
 392         scalar minfz = GREAT;
 393         forAll (pf, pfi)
 394         {
 395             minfz = min(minfz, points[pf[pfi]].z());
 396         }
 397
 398         if (minfz < zHeadMin)
 399         {
 400             pFaces[LINER][0].append(pf);
 401         }
 402         else
 403         {
 404             newCylinderHeadFaces.append(pf);
 405         }
 406     }
 407
 408     if (pFaces[CYLINDERHEAD][0].size() != newCylinderHeadFaces.size())
 409     {
 410         Info<< "Transfered faces from cylinder-head region to linder" << endl;
 411         pFaces[CYLINDERHEAD][0] = newCylinderHeadFaces;
 412     }
 413 }
 414
 415
 416 // Count the number of non-zero sized patches
 417 label nPatches = 0;
 418 for (int bci=0; bci<nBCs; bci++)
 419 {
 420     forAll (pFaces[bci], rgi)
 421     {
 422         if (pFaces[bci][rgi].size())
 423         {
 424             nPatches++;
 425         }
 426     }
 427 }
 428
 429
 430 // Sort-out the 2D/3D wedge geometry
 431 if (pFaces[WEDGE].size() && pFaces[WEDGE][0].size())
 432 {
 433     if (pFaces[WEDGE][0].size() == cellShapes.size())
 434     {
 435         // Distribute the points to be +/- 2.5deg from the x-z plane
 436
 437         scalar tanTheta = Foam::tan(2.5*mathematicalConstant::pi/180.0);
 438
 439         SLList<face>::iterator iterf = pFaces[WEDGE][0].begin();
 440         SLList<face>::iterator iterb = pFaces[WEDGE][1].begin();
 441         for
 442         (
 443             ;
 444             iterf != pFaces[WEDGE][0].end(), iterb != pFaces[WEDGE][1].end();
 445             ++iterf, ++iterb
 446         )
 447         {
 448             for (direction d=0; d<4; d++)
 449             {
 450                 points[iterf()[d]].y() = -tanTheta*points[iterf()[d]].x();
 451                 points[iterb()[d]].y() =  tanTheta*points[iterb()[d]].x();
 452             }
 453         }
 454     }
 455     else
 456     {
 457         pFaces[CYCLIC].setSize(1);
 458         pFaces[CYCLIC][0] = pFaces[WEDGE][0];
 459         for
 460         (
 461             SLList<face>::iterator iterb = pFaces[WEDGE][1].begin();
 462             iterb != pFaces[WEDGE][1].end();
 463             ++iterb
 464         )
 465         {
 466             pFaces[CYCLIC][0].append(iterb());
 467         }
 468
 469         pFaces[WEDGE].clear();
 470         nPatches--;
 471     }
 472 }
 473
 474
 475 // Build the patches
 476
 477 faceListList boundary(nPatches);
 478 wordList patchNames(nPatches);
 479 wordList patchTypes(nPatches);
 480 word defaultFacesName = "defaultFaces";
 481 word defaultFacesType = emptyPolyPatch::typeName;
 482 wordList patchPhysicalTypes(nPatches);
 483
 484 label nAddedPatches = 0;
 485
 486 for (int bci=0; bci<nBCs; bci++)
 487 {
 488     forAll (pFaces[bci], rgi)
 489     {
 490         if (pFaces[bci][rgi].size())
 491         {
 492             patchTypes[nAddedPatches] = *kivaPatchTypes[bci];
 493
 494             patchNames[nAddedPatches] = kivaPatchNames[bci];
 495
 496             if (pFaces[bci].size() > 1)
 497             {
 498                 patchNames[nAddedPatches] += name(rgi+1);
 499             }
 500
 501             boundary[nAddedPatches] = pFaces[bci][rgi];
 502             nAddedPatches++;
 503         }
 504     }
 505 }
 506
 507
 508 // Remove unused points and update cell and face labels accordingly
 509
 510 labelList pointLabels(nPoints, -1);
 511
 512 // Scan cells for used points
 513 forAll (cellShapes, celli)
 514 {
 515     forAll (cellShapes[celli], i)
 516     {
 517         pointLabels[cellShapes[celli][i]] = 1;
 518     }
 519 }
 520
 521 // Create addressing for used points and pack points array
 522 label newPointi = 0;
 523 forAll (pointLabels, pointi)
 524 {
 525     if (pointLabels[pointi] != -1)
 526     {
 527         pointLabels[pointi] = newPointi;
 528         points[newPointi++] = points[pointi];
 529     }
 530 }
 531 points.setSize(newPointi);
 532
 533 // Reset cell point labels
 534 forAll (cellShapes, celli)
 535 {
 536     cellShape& cs = cellShapes[celli];
 537
 538     forAll (cs, i)
 539     {
 540         cs[i] = pointLabels[cs[i]];
 541     }
 542 }
 543
 544 // Reset boundary-face point labels
 545 forAll (boundary, patchi)
 546 {
 547     forAll (boundary[patchi], facei)
 548     {
 549         face& f = boundary[patchi][facei];
 550
 551         forAll (f, i)
 552         {
 553             f[i] = pointLabels[f[i]];
 554         }
 555     }
 556 }
 557
 558 preservePatchTypes
 559 (
 560     runTime,
 561     runTime.constant(),
 562     polyMesh::defaultRegion,
 563     patchNames,
 564     patchTypes,
 565     defaultFacesName,
 566     defaultFacesType,
 567     patchPhysicalTypes
 568 );
 569
 570 // Build the mesh and write it out
 571 polyMesh pShapeMesh
 572 (
 573     IOobject
 574     (
 575         polyMesh::defaultRegion,
 576         runTime.constant(),
 577         runTime
 578     ),
 579     points,
 580     cellShapes,
 581     boundary,
 582     patchNames,
 583     patchTypes,
 584     defaultFacesName,
 585     defaultFacesType,
 586     patchPhysicalTypes
 587 );
 588
 589 Info << "Writing polyMesh" << endl;
 590 pShapeMesh.write();
 591
 592 fileName czPath
 593 (
 594     runTime.path()/runTime.constant()/polyMesh::defaultRegion/"cellZoning"
 595 );
 596 Info << "Writing cell zoning info to file: " << czPath << endl;
 597 OFstream cz(czPath);
 598
 599 cz << cellZoning << endl;
 600