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