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