cfx4ToFoam.C
Go to the documentation of this file.
1/*---------------------------------------------------------------------------*\
2 ========= |
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4 \\ / O peration |
5 \\ / A nd | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2011-2016 OpenFOAM Foundation
9 Copyright (C) 2020-2021 OpenCFD Ltd.
10-------------------------------------------------------------------------------
11License
12 This file is part of OpenFOAM.
13
14 OpenFOAM is free software: you can redistribute it and/or modify it
15 under the terms of the GNU General Public License as published by
16 the Free Software Foundation, either version 3 of the License, or
17 (at your option) any later version.
18
19 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22 for more details.
23
24 You should have received a copy of the GNU General Public License
25 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26
27Application
28 cfx4ToFoam
29
30Group
31 grpMeshConversionUtilities
32
33Description
34 Convert a CFX 4 mesh to OpenFOAM format.
35
36\*---------------------------------------------------------------------------*/
37
38#include "argList.H"
39#include "Time.H"
40#include "IFstream.H"
41#include "hexBlock.H"
42#include "polyMesh.H"
43#include "wallPolyPatch.H"
44#include "symmetryPolyPatch.H"
45#include "preservePatchTypes.H"
46#include "cellShape.H"
47
48using namespace Foam;
49
50// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
51
52int main(int argc, char *argv[])
53{
54 argList::addNote
55 (
56 "Convert a CFX 4 mesh to OpenFOAM format"
57 );
58
59 argList::noParallel();
60 argList::addArgument("CFX geom file");
61 argList::addOption
62 (
63 "scale",
64 "factor",
65 "Geometry scaling factor - default is 1"
66 );
67
68 argList args(argc, argv);
69
70 if (!args.check())
71 {
73 }
74
75 const scalar scaleFactor = args.getOrDefault<scalar>("scale", 1);
76
77 #include "createTime.H"
78
79 IFstream cfxFile(args.get<fileName>(1));
80
81 // Read the cfx information using a fixed format reader.
82 // Comments in the file are in C++ style, so the stream parser will remove
83 // them with no intervention
84 label nblock, npatch, nglue, nelem, npoint;
85
86 cfxFile >> nblock >> npatch >> nglue >> nelem >> npoint;
87
88 Info<< "Reading blocks" << endl;
89
90 PtrList<hexBlock> blocks(nblock);
91
92 {
93 word blockName;
94 label nx, ny, nz;
95
96 forAll(blocks, blockI)
97 {
98 cfxFile >> blockName;
99 cfxFile >> nx >> ny >> nz;
100
101 blocks.set(blockI, new hexBlock(nx, ny, nz));
102 }
103 }
104
105 Info<< "Reading patch definitions" << endl;
106
107 wordList cfxPatchTypes(npatch);
108 wordList cfxPatchNames(npatch);
109 labelList patchMasterBlocks(npatch);
110 labelList patchDirections(npatch);
111 labelListList patchRanges(npatch);
112
113 {
114 label no, blkNo, patchLabel;
115
116 forAll(cfxPatchTypes, patchi)
117 {
118 // Grab patch type and name
119 cfxFile >> cfxPatchTypes[patchi] >> cfxPatchNames[patchi] >> no;
120
121 // Grab patch range
122 patchRanges[patchi].setSize(6);
123 labelList& curRange = patchRanges[patchi];
124
125 forAll(curRange, rI)
126 {
127 cfxFile >> curRange[rI];
128 }
129
130 // Grab patch direction and master block ID
131 // Note: direc is the direction, from the cfx manual
132 // 0 = solid (3-D patch),
133 // 1 = high i, 2 = high j, 3 = high k
134 // 4 = low i, 5 = low j, 6 = low k
135 cfxFile >> patchDirections[patchi] >> blkNo >> patchLabel;
136
137 patchMasterBlocks[patchi] = blkNo - 1;
138 }
139 }
140
141 Info<< "Reading block glueing information" << endl;
142
143 labelList glueMasterPatches(nglue, -1);
144 labelList glueSlavePatches(nglue, -1);
145
146 {
147 label masterPatch, slavePatch;
148 label dirIndex1, dirIndex2, dirIndex3, joinNumber;
149
150 for (label glueI = 0; glueI < nglue; glueI++)
151 {
152 cfxFile >> masterPatch >> slavePatch;
153 cfxFile >> dirIndex1 >> dirIndex2 >> dirIndex3 >> joinNumber;
154
155 glueMasterPatches[glueI] = masterPatch - 1;
156 glueSlavePatches[glueI] = slavePatch - 1;
157 }
158 }
159
160 Info<< "Reading block points" << endl;
161
162 forAll(blocks, blockI)
163 {
164 Info<< "block " << blockI << " is a ";
165 blocks[blockI].readPoints(cfxFile);
166 }
167
168 Info<< "Calculating block offsets" << endl;
169
170 labelList blockOffsets(nblock, -1);
171
172 blockOffsets[0] = 0;
173
174 label nMeshPoints = blocks[0].nBlockPoints();
175 label nMeshCells = blocks[0].nBlockCells();
176
177 for (label blockI = 1; blockI < nblock; blockI++)
178 {
179 nMeshPoints += blocks[blockI].nBlockPoints();
180 nMeshCells += blocks[blockI].nBlockCells();
181
182 blockOffsets[blockI] =
183 blockOffsets[blockI - 1]
184 + blocks[blockI - 1].nBlockPoints();
185 }
186
187 Info<< "Assembling patches" << endl;
188
189 faceListList rawPatches(npatch);
190
191 forAll(rawPatches, patchi)
192 {
193 const word& patchType = cfxPatchTypes[patchi];
194
195 // reject volume patches
196 if
197 (
198 patchType == "POROUS" || patchType == "SOLID"
199 || patchType == "SOLCON" || patchType == "USER3D"
200 )
201 {
202 patchMasterBlocks[patchi] = -1;
203 rawPatches[patchi].setSize(0);
204 }
205 else
206 {
207 // read and create a 2-D patch
208 rawPatches[patchi] =
209 blocks[patchMasterBlocks[patchi]].patchFaces
210 (
211 patchDirections[patchi],
212 patchRanges[patchi]
213 );
214
215 }
216 }
217
218 Info<< "Merging points ";
219
220 labelList pointMergeList(nMeshPoints, -1);
221
222 // In order to ensure robust merging, it is necessary to traverse
223 // the patch glueing list until the pointMergeList stops changing.
224 //
225
226 // For efficiency, create merge pairs in the first pass
227 labelListListList glueMergePairs(glueMasterPatches.size());
228
229 forAll(glueMasterPatches, glueI)
230 {
231 const label masterPatch = glueMasterPatches[glueI];
232 const label slavePatch = glueSlavePatches[glueI];
233
234 const label blockPlabel = patchMasterBlocks[masterPatch];
235 const label blockNlabel = patchMasterBlocks[slavePatch];
236
237 const pointField& blockPpoints = blocks[blockPlabel].points();
238 const pointField& blockNpoints = blocks[blockNlabel].points();
239
240 const faceList& blockPFaces = rawPatches[masterPatch];
241 const faceList& blockNFaces = rawPatches[slavePatch];
242
243 labelListList& curPairs = glueMergePairs[glueI];
244 curPairs.setSize(blockPFaces.size());
245
246 if (blockPFaces.size() != blockNFaces.size())
247 {
249 << "Inconsistent number of faces for glue pair "
250 << glueI << " between blocks " << blockPlabel + 1
251 << " and " << blockNlabel + 1
252 << abort(FatalError);
253 }
254
255 // Calculate sqr of the merge tolerance as 1/10th of the min
256 // sqr point to point distance on the block face. This is an
257 // N^2 algorithm, sorry but I cannot quickly come up with
258 // something better.
259
260 scalar sqrMergeTol = GREAT;
261
262 forAll(blockPFaces, blockPFaceLabel)
263 {
264 const labelList& blockPFacePoints =
265 blockPFaces[blockPFaceLabel];
266
267 forAll(blockPFacePoints, blockPFacePointi)
268 {
269 forAll(blockPFacePoints, blockPFacePointi2)
270 {
271 if (blockPFacePointi != blockPFacePointi2)
272 {
273 sqrMergeTol =
274 min
275 (
276 sqrMergeTol,
277 magSqr
278 (
279 blockPpoints
280 [blockPFacePoints[blockPFacePointi]]
281 - blockPpoints
282 [blockPFacePoints[blockPFacePointi2]]
283 )
284 );
285 }
286 }
287 }
288 }
289
290 sqrMergeTol /= 10.0;
291
292 bool found = false;
293
294 // N-squared point search over all points of all faces of
295 // master block over all point of all faces of slave block
296 forAll(blockPFaces, blockPFaceLabel)
297 {
298 const labelList& blockPFacePoints =
299 blockPFaces[blockPFaceLabel];
300
301 labelList& cp = curPairs[blockPFaceLabel];
302 cp.setSize(blockPFacePoints.size());
303
304 forAll(blockPFacePoints, blockPFacePointi)
305 {
306 found = false;
307
308 forAll(blockNFaces, blockNFaceLabel)
309 {
310 const labelList& blockNFacePoints =
311 blockNFaces[blockNFaceLabel];
312
313 forAll(blockNFacePoints, blockNFacePointi)
314 {
315 if
316 (
317 magSqr
318 (
319 blockPpoints
320 [blockPFacePoints[blockPFacePointi]]
321 - blockNpoints
322 [blockNFacePoints[blockNFacePointi]]
323 )
324 < sqrMergeTol
325 )
326 {
327 // Found a new pair
328 found = true;
329
330 cp[blockPFacePointi] =
331 blockNFacePoints[blockNFacePointi];
332
333 label PpointLabel =
334 blockPFacePoints[blockPFacePointi]
335 + blockOffsets[blockPlabel];
336
337 label NpointLabel =
338 blockNFacePoints[blockNFacePointi]
339 + blockOffsets[blockNlabel];
340
341 label minPN = min(PpointLabel, NpointLabel);
342
343 if (pointMergeList[PpointLabel] != -1)
344 {
345 minPN = min(minPN, pointMergeList[PpointLabel]);
346 }
347
348 if (pointMergeList[NpointLabel] != -1)
349 {
350 minPN = min(minPN, pointMergeList[NpointLabel]);
351 }
352
353 pointMergeList[PpointLabel]
354 = pointMergeList[NpointLabel]
355 = minPN;
356
357 break;
358 }
359 }
360 if (found) break;
361 }
362 }
363 }
364 }
365
366
367 bool changedPointMerge = false;
368 label nPasses = 0;
369
370 do
371 {
372 changedPointMerge = false;
373 nPasses++;
374
375 forAll(glueMasterPatches, glueI)
376 {
377 const label masterPatch = glueMasterPatches[glueI];
378 const label slavePatch = glueSlavePatches[glueI];
379
380 const label blockPlabel = patchMasterBlocks[masterPatch];
381 const label blockNlabel = patchMasterBlocks[slavePatch];
382
383 const faceList& blockPFaces = rawPatches[masterPatch];
384
385 const labelListList& curPairs = glueMergePairs[glueI];
386
387 forAll(blockPFaces, blockPFaceLabel)
388 {
389 const labelList& blockPFacePoints =
390 blockPFaces[blockPFaceLabel];
391
392 const labelList& cp = curPairs[blockPFaceLabel];
393
394 forAll(cp, blockPFacePointi)
395 {
396 label PpointLabel =
397 blockPFacePoints[blockPFacePointi]
398 + blockOffsets[blockPlabel];
399
400 label NpointLabel =
401 cp[blockPFacePointi]
402 + blockOffsets[blockNlabel];
403
404 if
405 (
406 pointMergeList[PpointLabel]
407 != pointMergeList[NpointLabel]
408 )
409 {
410 changedPointMerge = true;
411
412 pointMergeList[PpointLabel]
413 = pointMergeList[NpointLabel]
414 = min
415 (
416 pointMergeList[PpointLabel],
417 pointMergeList[NpointLabel]
418 );
419 }
420 }
421 }
422 }
423 Info<< "." << flush;
424 }
425 while (changedPointMerge && nPasses < 8);
426 Info<< endl;
427
428 if (changedPointMerge == true)
429 {
431 << "Point merging failed after max number of passes."
432 << abort(FatalError);
433 }
434
435
436 forAll(glueMasterPatches, glueI)
437 {
438 const label masterPatch = glueMasterPatches[glueI];
439 const label slavePatch = glueSlavePatches[glueI];
440
441 const label blockPlabel = patchMasterBlocks[masterPatch];
442 const label blockNlabel = patchMasterBlocks[slavePatch];
443
444 const faceList& blockPFaces = rawPatches[masterPatch];
445 const faceList& blockNFaces = rawPatches[slavePatch];
446
447
448 forAll(blockPFaces, blockPFaceLabel)
449 {
450 const labelList& blockPFacePoints
451 = blockPFaces[blockPFaceLabel];
452
453 forAll(blockPFacePoints, blockPFacePointi)
454 {
455 label PpointLabel =
456 blockPFacePoints[blockPFacePointi]
457 + blockOffsets[blockPlabel];
458
459 if (pointMergeList[PpointLabel] == -1)
460 {
462 << "Unable to merge point " << blockPFacePointi
463 << " of face " << blockPFaceLabel
464 << " of block " << blockPlabel
465 << abort(FatalError);
466 }
467 }
468 }
469
470 forAll(blockNFaces, blockNFaceLabel)
471 {
472 const labelList& blockNFacePoints
473 = blockNFaces[blockNFaceLabel];
474
475 forAll(blockNFacePoints, blockNFacePointi)
476 {
477 label NpointLabel =
478 blockNFacePoints[blockNFacePointi]
479 + blockOffsets[blockNlabel];
480
481 if (pointMergeList[NpointLabel] == -1)
482 {
484 << "Unable to merge point " << blockNFacePointi
485 << " of face " << blockNFaceLabel
486 << " of block " << blockNlabel
487 << abort(FatalError);
488 }
489 }
490 }
491 }
492
493
494 // sort merge list to return new point label (in new shorter list)
495 // given old point label
496 label nNewPoints = 0;
497
498 forAll(pointMergeList, pointLabel)
499 {
500 if (pointMergeList[pointLabel] > pointLabel)
501 {
503 << "ouch" << abort(FatalError);
504 }
505
506 if
507 (
508 (pointMergeList[pointLabel] == -1)
509 || pointMergeList[pointLabel] == pointLabel
510 )
511 {
512 pointMergeList[pointLabel] = nNewPoints;
513 nNewPoints++;
514 }
515 else
516 {
517 pointMergeList[pointLabel] =
518 pointMergeList[pointMergeList[pointLabel]];
519 }
520 }
521
522 nMeshPoints = nNewPoints;
523
524 Info<< "Creating points" << endl;
525
526 pointField points(nMeshPoints);
527
528 forAll(blocks, blockI)
529 {
530 const pointField& blockPoints = blocks[blockI].points();
531
532 forAll(blockPoints, blockPointLabel)
533 {
534 points
535 [
536 pointMergeList
537 [
538 blockPointLabel
539 + blockOffsets[blockI]
540 ]
541 ] = blockPoints[blockPointLabel];
542 }
543 }
544
545 // Scale the points
546 if (scaleFactor > 1.0 + SMALL || scaleFactor < 1.0 - SMALL)
547 {
548 points *= scaleFactor;
549 }
550
551 Info<< "Creating cells" << endl;
552
553 cellShapeList cellShapes(nMeshCells);
554
555 const cellModel& hex = cellModel::ref(cellModel::HEX);
556
557 label nCreatedCells = 0;
558
559 forAll(blocks, blockI)
560 {
561 labelListList curBlockCells = blocks[blockI].blockCells();
562
563 forAll(curBlockCells, blockCelli)
564 {
565 labelList cellPoints(curBlockCells[blockCelli].size());
566
567 forAll(cellPoints, pointi)
568 {
569 cellPoints[pointi] =
570 pointMergeList
571 [
572 curBlockCells[blockCelli][pointi]
573 + blockOffsets[blockI]
574 ];
575 }
576
577 cellShapes[nCreatedCells].reset(hex, cellPoints);
578
579 nCreatedCells++;
580 }
581 }
582
583 Info<< "Creating boundary patches" << endl;
584
585 faceListList boundary(npatch);
586 wordList patchNames(npatch);
587 wordList patchTypes(npatch);
588 word defaultFacesName = "defaultFaces";
589 word defaultFacesType = wallPolyPatch::typeName;
590
591 label nCreatedPatches = 0;
592
593 forAll(rawPatches, patchi)
594 {
595 if (rawPatches[patchi].size() && cfxPatchTypes[patchi] != "BLKBDY")
596 {
597 // Check if this name has been already created
598 label existingPatch = -1;
599
600 for (label oldPatchi = 0; oldPatchi < nCreatedPatches; oldPatchi++)
601 {
602 if (patchNames[oldPatchi] == cfxPatchNames[patchi])
603 {
604 existingPatch = oldPatchi;
605 break;
606 }
607 }
608
609 const faceList& curRawPatch = rawPatches[patchi];
610 label curBlock = patchMasterBlocks[patchi];
611
612 if (existingPatch >= 0)
613 {
614 Info<< "CFX patch " << patchi
615 << ", of type " << cfxPatchTypes[patchi]
616 << ", name " << cfxPatchNames[patchi]
617 << " already exists as OpenFOAM patch " << existingPatch
618 << ". Adding faces." << endl;
619
620 faceList& renumberedPatch = boundary[existingPatch];
621 label oldSize = renumberedPatch.size();
622 renumberedPatch.setSize(oldSize + curRawPatch.size());
623
624 forAll(curRawPatch, facei)
625 {
626 const face& oldFace = curRawPatch[facei];
627
628 face& newFace = renumberedPatch[oldSize + facei];
629 newFace.setSize(oldFace.size());
630
631 forAll(oldFace, pointi)
632 {
633 newFace[pointi] =
634 pointMergeList
635 [
636 oldFace[pointi]
637 + blockOffsets[curBlock]
638 ];
639 }
640 }
641 }
642 else
643 {
644 // Real patch to be created
645 faceList& renumberedPatch = boundary[nCreatedPatches];
646 renumberedPatch.setSize(curRawPatch.size());
647
648 forAll(curRawPatch, facei)
649 {
650 const face& oldFace = curRawPatch[facei];
651
652 face& newFace = renumberedPatch[facei];
653 newFace.setSize(oldFace.size());
654
655 forAll(oldFace, pointi)
656 {
657 newFace[pointi] =
658 pointMergeList
659 [
660 oldFace[pointi]
661 + blockOffsets[curBlock]
662 ];
663 }
664 }
665
666 Info<< "CFX patch " << patchi
667 << ", of type " << cfxPatchTypes[patchi]
668 << ", name " << cfxPatchNames[patchi]
669 << " converted into OpenFOAM patch " << nCreatedPatches
670 << " type ";
671
672 if (cfxPatchTypes[patchi] == "WALL")
673 {
674 Info<< "wall." << endl;
675
676 patchTypes[nCreatedPatches] = wallPolyPatch::typeName;
677 patchNames[nCreatedPatches] = cfxPatchNames[patchi];
678 nCreatedPatches++;
679 }
680 else if (cfxPatchTypes[patchi] == "SYMMET")
681 {
682 Info<< "symmetryPlane." << endl;
683
684 patchTypes[nCreatedPatches] = symmetryPolyPatch::typeName;
685 patchNames[nCreatedPatches] = cfxPatchNames[patchi];
686 nCreatedPatches++;
687 }
688 else if
689 (
690 cfxPatchTypes[patchi] == "INLET"
691 || cfxPatchTypes[patchi] == "OUTLET"
692 || cfxPatchTypes[patchi] == "PRESS"
693 || cfxPatchTypes[patchi] == "CNDBDY"
694 || cfxPatchTypes[patchi] == "USER2D"
695 )
696 {
697 Info<< "generic." << endl;
698
699 patchTypes[nCreatedPatches] = polyPatch::typeName;
700 patchNames[nCreatedPatches] = cfxPatchNames[patchi];
701 nCreatedPatches++;
702 }
703 else
704 {
706 << "Unrecognised CFX patch type "
707 << cfxPatchTypes[patchi]
708 << abort(FatalError);
709 }
710 }
711 }
712 }
713
714 boundary.setSize(nCreatedPatches);
715 patchTypes.setSize(nCreatedPatches);
716 patchNames.setSize(nCreatedPatches);
717
719
721 (
722 runTime,
723 runTime.constant(),
724 polyMesh::meshSubDir,
729 );
730
731 // Add information to dictionary
732 forAll(patchNames, patchi)
733 {
734 if (!patchDicts.set(patchi))
735 {
736 patchDicts.set(patchi, new dictionary());
737 }
738 // Add but not overwrite
739 patchDicts[patchi].add("type", patchTypes[patchi], false);
740 }
741
743 (
745 (
746 polyMesh::defaultRegion,
747 runTime.constant(),
748 runTime
749 ),
750 std::move(points),
752 boundary,
757 );
758
759 // Set the precision of the points data to 10
760 IOstream::defaultPrecision(max(10u, IOstream::defaultPrecision()));
761
762 Info<< "Writing polyMesh" << endl;
763 pShapeMesh.removeFiles();
764 pShapeMesh.write();
765
766 Info<< "End\n" << endl;
767
768 return 0;
769}
770
771
772// ************************************************************************* //
bool found
Y[inertIndex] max(0.0)
Input from file stream, using an ISstream.
Definition: IFstream.H:57
Defines the attributes of an object for which implicit objectRegistry management is supported,...
Definition: IOobject.H:170
void setSize(const label n)
Alias for resize()
Definition: List.H:218
A list of pointers to objects of type <T>, with allocation/deallocation management of the pointers....
Definition: PtrList.H:73
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
Extract command arguments and options from the supplied argc and argv parameters.
Definition: argList.H:124
T get(const label index) const
Get a value from the argument at index.
Definition: argListI.H:278
bool check(bool checkArgs=argList::argsMandatory(), bool checkOpts=true) const
Definition: argList.C:1913
T getOrDefault(const word &optName, const T &deflt) const
Get a value from the named option if present, or return default.
Definition: argListI.H:307
Maps a geometry to a set of cell primitives.
Definition: cellModel.H:73
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:126
void exit(const int errNo=1)
Exit : can be called for any error to exit program.
Definition: error.C:331
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
A class for handling file names.
Definition: fileName.H:76
Hex block definition used in the cfx converter.
Definition: hexBlock.H:53
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:81
A class for handling words, derived from Foam::string.
Definition: word.H:68
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
Namespace for OpenFOAM.
messageStream Info
Information stream (stdout output on master, null elsewhere)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33
errorManip< error > abort(error &err)
Definition: errorManip.H:144
error FatalError
Ostream & flush(Ostream &os)
Flush stream.
Definition: Ostream.H:364
preservePatchTypes
wordList patchTypes(nPatches)
wordList patchNames(nPatches)
PtrList< dictionary > patchDicts
Definition: readKivaGrid.H:532
polyMesh pShapeMesh(IOobject(polyMesh::defaultRegion, runTime.constant(), runTime), std::move(points), cellShapes, boundary, patchNames, patchDicts, defaultFacesName, defaultFacesType)
word defaultFacesName
Definition: readKivaGrid.H:455
word defaultFacesType
Definition: readKivaGrid.H:456
preservePatchTypes(runTime, runTime.constant(), polyMesh::meshSubDir, patchNames, patchDicts, defaultFacesName, defaultFacesType)
const volScalarField & cp
Foam::argList args(argc, argv)
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333