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);
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 
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.
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  &mergedCyclicPolyPatch::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<< "Transferred " << 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<< "Transferred 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;
539 (
540  runTime,
541  runTime.constant(),
542  polyMesh::meshSubDir,
543  patchNames,
544  patchDicts,
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  move(points),
569  cellShapes,
570  boundary,
571  patchNames,
572  patchDicts,
575 );
576 
577 // Un-merge any merged cyclics
579 
580 Info << "Writing polyMesh" << endl;
581 pShapeMesh.write();
582 
583 fileName czPath
584 (
585  runTime.path()/runTime.constant()/polyMesh::defaultRegion/"cellZoning"
586 );
587 Info << "Writing cell zoning info to file: " << czPath << endl;
588 OFstream cz(czPath);
589 
590 cz << cellZoning << endl;
word defaultFacesType
Definition: readKivaGrid.H:461
labelList imtab(nPoints)
List< faceList > faceListList
Definition: faceListFwd.H:45
FvWallInfoData< WallInfo, label > label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
labelList pointLabels(nPoints, -1)
error FatalError
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:306
#define forAllConstIter(Container, container, iter)
Iterate across all elements in the container object of type.
Definition: UList.H:477
#define forAllIter(Container, container, iter)
Iterate across all elements in the container object of type.
Definition: UList.H:459
uint8_t direction
Definition: direction.H:45
wordList patchTypes(nPatches)
labelList idface(nPoints)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
preservePatchTypes(runTime, runTime.constant(), polyMesh::meshSubDir, patchNames, patchDicts, defaultFacesName, defaultFacesType)
scalar degToRad(const scalar deg)
Conversion from degrees to radians.
ifstream kivaFile(kivaFileName.c_str())
labelList kmtab(nPoints)
polyMesh pShapeMesh(IOobject(polyMesh::defaultRegion, runTime.constant(), runTime), move(points), cellShapes, boundary, patchNames, patchDicts, defaultFacesName, defaultFacesType)
labelList jmtab(nPoints)
labelList i8tab(nPoints)
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
stressControl lookup("compactNormalStress") >> compactNormalStress
label nAddedPatches
Definition: readKivaGrid.H:463
labelList bcb(nPoints)
List< cellShape > cellShapeList
List of cellShapes and PtrList of List of cellShape.
Definition: cellShapeList.H:43
labelList i3tab(nPoints)
List< edge > edgeList
Definition: edgeList.H:38
label nPoints
labelList fv(nPoints)
const cellShapeList & cellShapes
face triFace(3)
wordList patchNames(nPatches)
label mTable
Definition: readKivaGrid.H:74
label i4
Definition: readKivaGrid.H:20
layerAndWeight min(const layerAndWeight &a, const layerAndWeight &b)
labelList idreg(nPoints)
List< label > labelList
A List of labels.
Definition: labelList.H:56
faceListList boundary(nPatches)
forAll(cellShapes, celli)
Definition: readKivaGrid.H:492
labelList bcl(nPoints)
labelList f(nPoints)
word defaultFacesName
Definition: readKivaGrid.H:460
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
face quadFace(4)
List< word > wordList
A List of words.
Definition: fileName.H:54
label patchi
label newPointi
Definition: readKivaGrid.H:501
pointField points(nPoints)
Info<< "Reading kiva grid from file "<< kivaFileName<< endl;char title[120];kivaFile.getline(title, 120, '\n');label nPoints, nCells, nRegs;kivaFile > nCells nPoints nRegs
Definition: readKivaGrid.H:17
messageStream Info
dimensioned< scalar > mag(const dimensioned< Type > &)
labelList bcf(nPoints)
PtrList< dictionary > patchDicts
Definition: readKivaGrid.H:537
polyMeshUnMergeCyclics(pShapeMesh)
dimensionedScalar tan(const dimensionedScalar &ds)
labelList i1tab(nPoints)
const FieldField< fvPatchField, Type > & ff(const FieldField< fvPatchField, Type > &bf)