readKivaGrid.H
Go to the documentation of this file.
1std::ifstream kivaFile(kivaFileName);
2
3if (!kivaFile.good())
4{
6 << "Cannot open file " << kivaFileName
7 << exit(FatalError);
8}
9
10Info << "Reading kiva grid from file " << kivaFileName << endl;
11
12
13char title[120];
14kivaFile.getline(title, 120, '\n');
15
16label nPoints, nCells, nRegs;
17kivaFile >> nCells >> nPoints >> nRegs;
18
19pointField points(nPoints);
20label i4;
21labelList idface(nPoints), fv(nPoints);
22
23for (label i=0; i<nPoints; i++)
24{
25 scalar ffv;
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
45label nBfaces = 0;
46
47for (label i=0; i<nPoints; i++)
48{
49 label i1, i3, i8;
50 scalar ff, fbcl, fbcf, fbcb;
51
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
74label mTable;
76
77if (mTable == 0)
78{
80 << "mTable == 0. This is not supported."
81 " kivaToFoam needs complete neighbour information"
82 << exit(FatalError);
83}
84
85
87
88for (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
99Info << "Finished reading KIVA file" << endl;
100
102
103labelList cellZoning(nPoints, -1);
104
105const cellModel& hex = cellModel::ref(cellModel::HEX);
106labelList hexLabels(8);
107label activeCells = 0;
108
109// Create and set the collocated point collapse map
110labelList pointMap(nPoints);
111
112forAll(pointMap, i)
113{
114 pointMap[i] = i;
115}
116
117// Initialise all cells to hex and search and set map for collocated points
118for (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].reset(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
165cellShapes.setSize(activeCells);
166cellZoning.setSize(activeCells);
167
168// Map collocated points to refer to the same point and collapse cell shape
169// to the corresponding hex-degenerate.
170forAll(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
182label bcIDs[11] = {-1, 0, 2, 4, -1, 5, -1, 6, 7, 8, 9};
183
184const label nBCs = 12;
185
186const 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
202enum 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
218const 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
235List<SLList<face>> pFaces[nBCs];
236
237face quadFace(4);
238face triFace(3);
239
240for (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
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 (const face& pf : pFaces[CYLINDERHEAD][0])
331 {
332 forAll(pf, pfi)
333 {
334 minz = min(minz, points[pf[pfi]].z());
335 }
336 }
337
338 minz += SMALL;
339
340 SLList<face> newLinerFaces;
341
342 for (const face& pf : pFaces[LINER][0])
343 {
344 scalar minfz = GREAT;
345 forAll(pf, pfi)
346 {
347 minfz = min(minfz, points[pf[pfi]].z());
348 }
349
350 if (minfz > minz)
351 {
352 pFaces[CYLINDERHEAD][0].append(pf);
353 }
354 else
355 {
356 newLinerFaces.append(pf);
357 }
358 }
359
360 if (pFaces[LINER][0].size() != newLinerFaces.size())
361 {
362 Info<< "Transferred " << pFaces[LINER][0].size() - newLinerFaces.size()
363 << " faces from liner region to cylinder head" << endl;
364 pFaces[LINER][0] = newLinerFaces;
365 }
366
367 SLList<face> newCylinderHeadFaces;
368
369 for (const face& pf : pFaces[CYLINDERHEAD][0])
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
396label nPatches = 0;
397for (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
410if (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 const scalar tanTheta = Foam::tan(degToRad(2.5));
417
418 auto iterf = pFaces[WEDGE][0].begin();
419 auto iterb = pFaces[WEDGE][1].begin();
420
421 for
422 (
423 ;
424 iterf.good() && iterb.good();
425 ++iterf, ++iterb
426 )
427 {
428 for (direction d=0; d<4; d++)
429 {
430 points[iterf()[d]].y() = -tanTheta*points[iterf()[d]].x();
431 points[iterb()[d]].y() = tanTheta*points[iterb()[d]].x();
432 }
433 }
434 }
435 else
436 {
437 pFaces[CYCLIC].setSize(1);
438 pFaces[CYCLIC][0] = pFaces[WEDGE][0];
439 for (const face& pf : pFaces[WEDGE][1])
440 {
441 pFaces[CYCLIC][0].append(pf);
442 }
443
444 pFaces[WEDGE].clear();
445 nPatches--;
446 }
447}
448
449
450// Build the patches
451
452faceListList boundary(nPatches);
455word defaultFacesName = "defaultFaces";
456word defaultFacesType = emptyPolyPatch::typeName;
457
459
460for (int bci=0; bci<nBCs; bci++)
461{
462 forAll(pFaces[bci], rgi)
463 {
464 if (pFaces[bci][rgi].size())
465 {
466 patchTypes[nAddedPatches] = *kivaPatchTypes[bci];
467
468 patchNames[nAddedPatches] = kivaPatchNames[bci];
469
470 if (pFaces[bci].size() > 1)
471 {
472 patchNames[nAddedPatches] += name(rgi+1);
473 }
474
475 boundary[nAddedPatches] = pFaces[bci][rgi];
477 }
478 }
479}
480
481
482// Remove unused points and update cell and face labels accordingly
483
484labelList pointLabels(nPoints, -1);
485
486// Scan cells for used points
488{
489 forAll(cellShapes[celli], i)
490 {
491 pointLabels[cellShapes[celli][i]] = 1;
492 }
493}
494
495// Create addressing for used points and pack points array
496label newPointi = 0;
498{
499 if (pointLabels[pointi] != -1)
500 {
501 pointLabels[pointi] = newPointi;
502 points[newPointi++] = points[pointi];
503 }
504}
506
507// Reset cell point labels
508forAll(cellShapes, celli)
509{
510 cellShape& cs = cellShapes[celli];
511
512 forAll(cs, i)
513 {
514 cs[i] = pointLabels[cs[i]];
515 }
516}
517
518// Reset boundary-face point labels
520{
521 forAll(boundary[patchi], facei)
522 {
523 face& f = boundary[patchi][facei];
524
525 forAll(f, i)
526 {
527 f[i] = pointLabels[f[i]];
528 }
529 }
530}
531
532PtrList<dictionary> patchDicts;
534(
535 runTime,
536 runTime.constant(),
537 polyMesh::meshSubDir,
542);
543// Add information to dictionary
545{
546 if (!patchDicts.set(patchi))
547 {
548 patchDicts.set(patchi, new dictionary());
549 }
550 // Add but not overwrite
551 patchDicts[patchi].add("type", patchTypes[patchi], false);
552}
553
554// Build the mesh and write it out
556(
557 IOobject
558 (
559 polyMesh::defaultRegion,
560 runTime.constant(),
561 runTime
562 ),
563 std::move(points),
565 boundary,
570);
571
572// Set the precision of the points data to 10
573IOstream::defaultPrecision(max(10u, IOstream::defaultPrecision()));
574
575Info << "Writing polyMesh" << endl;
576pShapeMesh.removeFiles();
577pShapeMesh.write();
578
579fileName czPath
580(
581 runTime.path()/runTime.constant()/polyMesh::defaultRegion/"cellZoning"
582);
583Info << "Writing cell zoning info to file: " << czPath << endl;
584OFstream cz(czPath);
585
586cz << cellZoning << endl;
Y[inertIndex] max(0.0)
const cellModel & hex
cellShapeList cellShapes
faceListList boundary
engineTime & runTime
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
const pointField & points
label nPoints
List< cellShape > cellShapeList
List of cellShapes and PtrList of List of cellShape.
Definition: cellShapeList.H:45
dimensionedScalar tan(const dimensionedScalar &ds)
List< label > labelList
A List of labels.
Definition: List.H:66
messageStream Info
Information stream (stdout output on master, null elsewhere)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33
List< edge > edgeList
A List of edges.
Definition: edgeList.H:63
constexpr auto end(C &c) -> decltype(c.end())
Return iterator to the end of the container c.
Definition: stdFoam.H:158
wordList patchTypes(nPatches)
label nBfaces
Definition: readKivaGrid.H:45
labelList i3tab(nPoints)
labelList kmtab(nPoints)
labelList idreg(nPoints)
Info<< "Finished reading KIVA file"<< endl;cellShapeList cellShapes(nPoints);labelList cellZoning(nPoints, -1);const cellModel &hex=cellModel::ref(cellModel::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].reset(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, &oldCyclicPolyPatch::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:235
wordList patchNames(nPatches)
labelList bcl(nPoints)
label i4
Definition: readKivaGrid.H:20
labelList i1tab(nPoints)
PtrList< dictionary > patchDicts
Definition: readKivaGrid.H:532
label newPointi
Definition: readKivaGrid.H:496
face quadFace(4)
labelList imtab(nPoints)
labelList jmtab(nPoints)
polyMesh pShapeMesh(IOobject(polyMesh::defaultRegion, runTime.constant(), runTime), std::move(points), cellShapes, boundary, patchNames, patchDicts, defaultFacesName, defaultFacesType)
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:458
word defaultFacesName
Definition: readKivaGrid.H:455
labelList i8tab(nPoints)
labelList bcf(nPoints)
word defaultFacesType
Definition: readKivaGrid.H:456
preservePatchTypes(runTime, runTime.constant(), polyMesh::meshSubDir, patchNames, patchDicts, defaultFacesName, defaultFacesType)
labelList bcb(nPoints)
labelList fv(nPoints)
label mTable
Definition: readKivaGrid.H:74
std::ifstream kivaFile(kivaFileName)
label nPatches
Definition: readKivaGrid.H:396
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333