FreeFOAM The Cross-Platform CFD Toolkit
readKivaGrid.H
Go to the documentation of this file.
1 ifstream kivaFile(kivaFileName.c_str());
2 
3 if (!kivaFile.good())
4 {
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 {
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.
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 
203 {
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;
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
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 
559 (
560  runTime,
561  runTime.constant(),
562  polyMesh::defaultRegion,
563  patchNames,
564  patchTypes,
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  xferMove(points),
580  cellShapes,
581  boundary,
582  patchNames,
583  patchTypes,
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 
601 
602 // ************************ vim: set sw=4 sts=4 et: ************************ //