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     &oldCyclicPolyPatch::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 // Transfer 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     forAllConstIter(SLList<face>, pFaces[CYLINDERHEAD][0], iter)
 331     {
 332         const face& pf = iter();
 333
 334         forAll(pf, pfi)
 335         {
 336             minz = min(minz, points[pf[pfi]].z());
 337         }
 338     }
 339
 340     minz += SMALL;
 341
 342     SLList<face> newLinerFaces;
 343
 344     forAllConstIter(SLList<face>, pFaces[LINER][0], iter)
 345     {
 346         const face& pf = iter();
 347
 348         scalar minfz = GREAT;
 349         forAll(pf, pfi)
 350         {
 351             minfz = min(minfz, points[pf[pfi]].z());
 352         }
 353
 354         if (minfz > minz)
 355         {
 356             pFaces[CYLINDERHEAD][0].append(pf);
 357         }
 358         else
 359         {
 360             newLinerFaces.append(pf);
 361         }
 362     }
 363
 364     if (pFaces[LINER][0].size() != newLinerFaces.size())
 365     {
 366         Info<< "Transfered " << pFaces[LINER][0].size() - newLinerFaces.size()
 367             << " faces from liner region to cylinder head" << endl;
 368         pFaces[LINER][0] = newLinerFaces;
 369     }
 370
 371     SLList<face> newCylinderHeadFaces;
 372
 373     forAllConstIter(SLList<face>, pFaces[CYLINDERHEAD][0], iter)
 374     {
 375         const face& pf = iter();
 376
 377         scalar minfz = GREAT;
 378         forAll(pf, pfi)
 379         {
 380             minfz = min(minfz, points[pf[pfi]].z());
 381         }
 382
 383         if (minfz < zHeadMin)
 384         {
 385             pFaces[LINER][0].append(pf);
 386         }
 387         else
 388         {
 389             newCylinderHeadFaces.append(pf);
 390         }
 391     }
 392
 393     if (pFaces[CYLINDERHEAD][0].size() != newCylinderHeadFaces.size())
 394     {
 395         Info<< "Transfered faces from cylinder-head region to linder" << endl;
 396         pFaces[CYLINDERHEAD][0] = newCylinderHeadFaces;
 397     }
 398 }
 399
 400
 401 // Count the number of non-zero sized patches
 402 label nPatches = 0;
 403 for (int bci=0; bci<nBCs; bci++)
 404 {
 405     forAll(pFaces[bci], rgi)
 406     {
 407         if (pFaces[bci][rgi].size())
 408         {
 409             nPatches++;
 410         }
 411     }
 412 }
 413
 414
 415 // Sort-out the 2D/3D wedge geometry
 416 if (pFaces[WEDGE].size() && pFaces[WEDGE][0].size())
 417 {
 418     if (pFaces[WEDGE][0].size() == cellShapes.size())
 419     {
 420         // Distribute the points to be +/- 2.5deg from the x-z plane
 421
 422         scalar tanTheta = Foam::tan(degToRad(2.5));
 423
 424         SLList<face>::iterator iterf = pFaces[WEDGE][0].begin();
 425         SLList<face>::iterator iterb = pFaces[WEDGE][1].begin();
 426         for
 427         (
 428             ;
 429             iterf != pFaces[WEDGE][0].end(), iterb != pFaces[WEDGE][1].end();
 430             ++iterf, ++iterb
 431         )
 432         {
 433             for (direction d=0; d<4; d++)
 434             {
 435                 points[iterf()[d]].y() = -tanTheta*points[iterf()[d]].x();
 436                 points[iterb()[d]].y() =  tanTheta*points[iterb()[d]].x();
 437             }
 438         }
 439     }
 440     else
 441     {
 442         pFaces[CYCLIC].setSize(1);
 443         pFaces[CYCLIC][0] = pFaces[WEDGE][0];
 444         forAllIter(SLList<face>, pFaces[WEDGE][1], iterb)
 445         {
 446             pFaces[CYCLIC][0].append(iterb());
 447         }
 448
 449         pFaces[WEDGE].clear();
 450         nPatches--;
 451     }
 452 }
 453
 454
 455 // Build the patches
 456
 457 faceListList boundary(nPatches);
 458 wordList patchNames(nPatches);
 459 wordList patchTypes(nPatches);
 460 word defaultFacesName = "defaultFaces";
 461 word defaultFacesType = emptyPolyPatch::typeName;
 462
 463 label nAddedPatches = 0;
 464
 465 for (int bci=0; bci<nBCs; bci++)
 466 {
 467     forAll(pFaces[bci], rgi)
 468     {
 469         if (pFaces[bci][rgi].size())
 470         {
 471             patchTypes[nAddedPatches] = *kivaPatchTypes[bci];
 472
 473             patchNames[nAddedPatches] = kivaPatchNames[bci];
 474
 475             if (pFaces[bci].size() > 1)
 476             {
 477                 patchNames[nAddedPatches] += name(rgi+1);
 478             }
 479
 480             boundary[nAddedPatches] = pFaces[bci][rgi];
 481             nAddedPatches++;
 482         }
 483     }
 484 }
 485
 486
 487 // Remove unused points and update cell and face labels accordingly
 488
 489 labelList pointLabels(nPoints, -1);
 490
 491 // Scan cells for used points
 492 forAll(cellShapes, celli)
 493 {
 494     forAll(cellShapes[celli], i)
 495     {
 496         pointLabels[cellShapes[celli][i]] = 1;
 497     }
 498 }
 499
 500 // Create addressing for used points and pack points array
 501 label newPointi = 0;
 502 forAll(pointLabels, pointi)
 503 {
 504     if (pointLabels[pointi] != -1)
 505     {
 506         pointLabels[pointi] = newPointi;
 507         points[newPointi++] = points[pointi];
 508     }
 509 }
 510 points.setSize(newPointi);
 511
 512 // Reset cell point labels
 513 forAll(cellShapes, celli)
 514 {
 515     cellShape& cs = cellShapes[celli];
 516
 517     forAll(cs, i)
 518     {
 519         cs[i] = pointLabels[cs[i]];
 520     }
 521 }
 522
 523 // Reset boundary-face point labels
 524 forAll(boundary, patchi)
 525 {
 526     forAll(boundary[patchi], facei)
 527     {
 528         face& f = boundary[patchi][facei];
 529
 530         forAll(f, i)
 531         {
 532             f[i] = pointLabels[f[i]];
 533         }
 534     }
 535 }
 536
 537 PtrList<dictionary> patchDicts;
 538 preservePatchTypes
 539 (
 540     runTime,
 541     runTime.constant(),
 542     polyMesh::meshSubDir,
 543     patchNames,
 544     patchDicts,
 545     defaultFacesName,
 546     defaultFacesType
 547 );
 548 // Add information to dictionary
 549 forAll(patchNames, patchI)
 550 {
 551     if (!patchDicts.set(patchI))
 552     {
 553         patchDicts.set(patchI, new dictionary());
 554     }
 555     // Add but not overwrite
 556     patchDicts[patchI].add("type", patchTypes[patchI], false);
 557 }
 558
 559 // Build the mesh and write it out
 560 polyMesh pShapeMesh
 561 (
 562     IOobject
 563     (
 564         polyMesh::defaultRegion,
 565         runTime.constant(),
 566         runTime
 567     ),
 568     xferMove(points),
 569     cellShapes,
 570     boundary,
 571     patchNames,
 572     patchDicts,
 573     defaultFacesName,
 574     defaultFacesType
 575 );
 576
 577 Info << "Writing polyMesh" << endl;
 578 pShapeMesh.write();
 579
 580 fileName czPath
 581 (
 582     runTime.path()/runTime.constant()/polyMesh::defaultRegion/"cellZoning"
 583 );
 584 Info << "Writing cell zoning info to file: " << czPath << endl;
 585 OFstream cz(czPath);
 586
 587 cz << cellZoning << endl;
 588