FreeFOAM The Cross-Platform CFD Toolkit
ccm26ToFoam.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  ccm26ToFoam
26 
27 Description
28  Reads CCM files as written by Prostar/ccm using ccm 2.6 (not 2.4)
29 
30  - does polyhedral mesh
31  - does not handle 'interfaces' (couples)
32  - does not handle cyclics. Use createPatch to recreate these
33  - does not do data
34  - does patch names only if they are in the problem description
35 
36 Usage
37 
38  - ccm26ToFoam [OPTIONS] <ccm26 file>
39 
40  @param <ccm26 file> \n
41  @todo Detailed description of argument.
42 
43  @param -case <dir>\n
44  Case directory.
45 
46  @param -help \n
47  Display help message.
48 
49  @param -doc \n
50  Display Doxygen API documentation page for this application.
51 
52  @param -srcDoc \n
53  Display Doxygen source documentation page for this application.
54 
55 \*---------------------------------------------------------------------------*/
56 
57 #include <OpenFOAM/ListOps.H>
58 #include <OpenFOAM/argList.H>
59 #include <OpenFOAM/Time.H>
60 #include <finiteVolume/fvMesh.H>
61 #include <finiteVolume/volFields.H>
64 #include <OpenFOAM/wallPolyPatch.H>
65 #include <OpenFOAM/SortableList.H>
66 #include <meshTools/cellSet.H>
67 
68 #include <libccmio/ccmio.h>
69 #include <vector>
70 
71 using namespace Foam;
72 
73 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
74 
75 static char const kDefaultState[] = "default";
76 static int const kVertOffset = 2;
77 
78 
79 // Determine upper-triangular order for internal faces:
80 labelList getInternalFaceOrder
81 (
82  const cellList& cells,
83  const labelList& owner,
84  const labelList& neighbour
85 )
86 {
87  labelList oldToNew(owner.size(), -1);
88 
89  // First unassigned face
90  label newFaceI = 0;
91 
92  forAll(cells, cellI)
93  {
94  const labelList& cFaces = cells[cellI];
95 
96  SortableList<label> nbr(cFaces.size());
97 
98  forAll(cFaces, i)
99  {
100  label faceI = cFaces[i];
101 
102  label nbrCellI = neighbour[faceI];
103 
104  if (nbrCellI != -1)
105  {
106  // Internal face. Get cell on other side.
107  if (nbrCellI == cellI)
108  {
109  nbrCellI = owner[faceI];
110  }
111 
112  if (cellI < nbrCellI)
113  {
114  // CellI is master
115  nbr[i] = nbrCellI;
116  }
117  else
118  {
119  // nbrCell is master. Let it handle this face.
120  nbr[i] = -1;
121  }
122  }
123  else
124  {
125  // External face. Do later.
126  nbr[i] = -1;
127  }
128  }
129 
130  nbr.sort();
131 
132  forAll(nbr, i)
133  {
134  if (nbr[i] != -1)
135  {
136  oldToNew[cFaces[nbr.indices()[i]]] = newFaceI++;
137  }
138  }
139  }
140 
141  // Keep boundary faces in same order.
142  for (label faceI = newFaceI; faceI < owner.size(); faceI++)
143  {
144  oldToNew[faceI] = faceI;
145  }
146 
147  return oldToNew;
148 }
149 
150 
151 void storeCellInZone
152 (
153  const label cellI,
154  const label cellType,
155  Map<label>& typeToZone,
156  List<DynamicList<label> >& zoneCells
157 )
158 {
159  if (cellType >= 0)
160  {
161  Map<label>::iterator zoneFnd = typeToZone.find(cellType);
162 
163  if (zoneFnd == typeToZone.end())
164  {
165  // New type. Allocate zone for it.
166  zoneCells.setSize(zoneCells.size() + 1);
167 
168  label zoneI = zoneCells.size()-1;
169 
170  Info<< "Mapping type " << cellType << " to Foam cellZone "
171  << zoneI << endl;
172  typeToZone.insert(cellType, zoneI);
173 
174  zoneCells[zoneI].append(cellI);
175  }
176  else
177  {
178  // Existing zone for type
179  zoneCells[zoneFnd()].append(cellI);
180  }
181  }
182 }
183 
184 
185 void CheckError(CCMIOError const &err, const Foam::string& str)
186 {
187  if (err != kCCMIONoErr)
188  {
189  FatalErrorIn("CheckError")
190  << str << exit(FatalError);
191  }
192 }
193 
194 
195 void ReadVertices
196 (
197  CCMIOError &err,
198  CCMIOID &vertices,
199  labelList& foamPointMap,
200  pointField& foamPoints
201 )
202 {
203 
204  // Read the vertices. This involves reading both the vertex data and
205  // the map, which maps the index into the data array with the ID number.
206  // As we process the vertices we need to be sure to scale them by the
207  // appropriate scaling factor. The offset is just to show you can read
208  // any chunk. Normally this would be in a for loop.
209 
210  CCMIOSize nVertices;
211  CCMIOEntitySize(&err, vertices, &nVertices, NULL);
212 
213  List<int> mapData(nVertices, 0);
214  List<float> verts(3*nVertices, 0);
215 
216  int offset = 0;
217  int offsetPlusSize = offset+nVertices;
218 
219  int dims = 1;
220  float scale;
221  CCMIOID mapID;
222  CCMIOReadVerticesf(&err, vertices, &dims, &scale, &mapID, verts.begin(),
223  offset, offsetPlusSize);
224  CCMIOReadMap(&err, mapID, mapData.begin(), offset, offsetPlusSize);
225 
226  //CCMIOSize size;
227  //CCMIOEntityDescription(&err, vertices, &size, NULL);
228  //char *desc = new char[size + 1];
229  //CCMIOEntityDescription(&err, vertices, NULL, desc);
230  //Pout<< "label: '" << desc << "'" << endl;
231  //delete [] desc;
232 
233  // Convert to foamPoints
234  foamPoints.setSize(nVertices);
235  foamPoints = vector::zero;
236  foamPointMap.setSize(nVertices);
237 
238  forAll(foamPointMap, i)
239  {
240  foamPointMap[i] = mapData[i];
241  for (direction cmpt = 0; cmpt < dims; cmpt++)
242  {
243  foamPoints[i][cmpt] = verts[dims*i + cmpt]*scale;
244  }
245  }
246 }
247 
248 
249 void ReadProblem
250 (
251  CCMIOError &err,
252  CCMIOID& problem,
253 
254  const Map<label>& prostarToFoamPatch,
255  Map<word>& foamCellTypeNames,
256  wordList& foamPatchTypes,
257  wordList& foamPatchNames
258 )
259 {
260  // ... walk through each cell type and print it...
261  CCMIOID next;
262  int i = 0;
263  while
264  (
265  CCMIONextEntity(NULL, problem, kCCMIOCellType, &i, &next)
266  == kCCMIONoErr
267  )
268  {
269  char *name;
270  int size;
271  int cellType;
272 
273  // ... if it has a material type. (Note that we do not pass in
274  // an array to get the name because we do not know how long the
275  // string is yet. Many parameters to CCMIO functions that
276  // return
277  // data can be NULL if that data is not needed.)
278  if
279  (
280  CCMIOReadOptstr(NULL, next, "MaterialType", &size, NULL)
281  == kCCMIONoErr
282  )
283  {
284  name = new char[size + 1];
285  CCMIOReadOptstr(&err, next, "MaterialType", &size, name);
286  CCMIOGetEntityIndex(&err, next, &cellType);
287 
288  foamCellTypeNames.insert(cellType, name);
289  Pout<< "Celltype:" << cellType << " name:" << name << endl;
290 
291  delete [] name;
292  }
293  }
294 
295  // ... walk through each region description and print it...
296 
297 
298  CCMIOID boundary;
299  label regionI = 0;
300  int k = 0;
301  while
302  (
303  CCMIONextEntity(NULL, problem, kCCMIOBoundaryRegion, &k, &boundary)
304  == kCCMIONoErr
305  )
306  {
307  // Index of foam patch
308  label foamPatchI = -1;
309 
310  // Read prostar id
311 
312  int prostarI = -1;
313  if
314  (
315  CCMIOReadOpti(NULL, boundary, "ProstarRegionNumber", &prostarI)
316  == kCCMIONoErr
317  )
318  {
319  Pout<< "For region:" << regionI
320  << " found ProstarRegionNumber:" << prostarI << endl;
321  }
322  else
323  {
324  prostarI = regionI;
325 
326  Pout<< "For region:" << regionI
327  << "did not find ProstarRegionNumber entry. Assuming "
328  << prostarI << endl;
329  }
330 
331 
332  if (prostarToFoamPatch.found(prostarI))
333  {
334  foamPatchI = prostarToFoamPatch[prostarI];
335 
336  // Read boundary type
337 
338  int size;
339  if
340  (
341  CCMIOReadOptstr(NULL, boundary, "BoundaryType", &size, NULL)
342  == kCCMIONoErr
343  )
344  {
345  char* s = new char[size + 1];
346  CCMIOReadOptstr(NULL, boundary, "BoundaryType", &size, s);
347  s[size] = '\0';
348  foamPatchTypes[foamPatchI] = string::validate<word>(string(s));
349  delete [] s;
350  }
351 
352 
353  //foamPatchMap.append(prostarI);
354 
355 
356  // Read boundary name:
357  // - from BoundaryName field (Prostar)
358  // - from 'Label' field (ccm+?)
359  // - make up one from prostar id.
360 
361  if
362  (
363  CCMIOReadOptstr(NULL, boundary, "BoundaryName", &size, NULL)
364  == kCCMIONoErr
365  )
366  {
367  char* name = new char[size + 1];
368  CCMIOReadOptstr(NULL, boundary, "BoundaryName", &size, name);
369  name[size] = '\0';
370  foamPatchNames[foamPatchI] = string::validate<word>(string(name));
371  delete [] name;
372  }
373  else if
374  (
375  CCMIOReadOptstr(NULL, boundary, "Label", &size, NULL)
376  == kCCMIONoErr
377  )
378  {
379  char* name = new char[size + 1];
380  CCMIOReadOptstr(NULL, boundary, "Label", &size, name);
381  name[size] = '\0';
382  foamPatchNames[foamPatchI] = string::validate<word>(string(name));
383  delete [] name;
384  }
385  else
386  {
387  foamPatchNames[foamPatchI] =
388  foamPatchTypes[foamPatchI]
389  + Foam::name(foamPatchI);
390  Pout<< "Made up name:" << foamPatchNames[foamPatchI]
391  << endl;
392  }
393 
394  Pout<< "Read patch:" << foamPatchI
395  << " name:" << foamPatchNames[foamPatchI]
396  << " foamPatchTypes:" << foamPatchTypes[foamPatchI]
397  << endl;
398  }
399 
400  regionI++;
401  }
402 }
403 
404 
405 void ReadCells
406 (
407  CCMIOError &err,
408  CCMIOID& topology,
409  labelList& foamCellMap,
410  labelList& foamCellType,
411  Map<label>& prostarToFoamPatch,
412  DynamicList<label>& foamPatchSizes,
413  DynamicList<label>& foamPatchStarts,
414  labelList& foamFaceMap,
415  labelList& foamOwner,
416  labelList& foamNeighbour,
417  faceList& foamFaces
418 )
419 {
420  // Read the cells.
421  // ~~~~~~~~~~~~~~~
422 
423  // Store cell IDs so that we know what cells are in
424  // this post data.
425  CCMIOID id;
426  CCMIOGetEntity(&err, topology, kCCMIOCells, 0, &id);
427  CCMIOSize nCells;
428  CCMIOEntitySize(&err, id, &nCells, NULL);
429 
430  std::vector<int> mapData(nCells);
431  std::vector<int> cellType(nCells);
432 
433  CCMIOID mapID;
434  CCMIOReadCells(&err, id, &mapID, &cellType[0], 0, nCells);
435  CCMIOReadMap(&err, mapID, &mapData[0], 0, nCells);
436  CheckError(err, "Error reading cells");
437 
438  foamCellMap.setSize(nCells);
439  foamCellType.setSize(nCells);
440  forAll(foamCellMap, i)
441  {
442  foamCellMap[i] = mapData[i];
443  foamCellType[i] = cellType[i];
444  }
445 
446 
447  // Read the internal faces.
448  // ~~~~~~~~~~~~~~~~~~~~~~~~
449 
450  CCMIOGetEntity(&err, topology, kCCMIOInternalFaces, 0, &id);
451  CCMIOSize nInternalFaces;
452  CCMIOEntitySize(&err, id, &nInternalFaces, NULL);
453  Pout<< "nInternalFaces:" << nInternalFaces << endl;
454 
455  // Determine patch sizes before reading internal faces
456  label foamNFaces = nInternalFaces;
457  int index = 0;
458  while
459  (
460  CCMIONextEntity(NULL, topology, kCCMIOBoundaryFaces, &index, &id)
461  == kCCMIONoErr
462  )
463  {
464  CCMIOSize size;
465  CCMIOEntitySize(&err, id, &size, NULL);
466 
467  Pout<< "Read kCCMIOBoundaryFaces entry with " << size
468  << " faces." << endl;
469 
470  foamPatchStarts.append(foamNFaces);
471  foamPatchSizes.append(size);
472  foamNFaces += size;
473  }
474  foamPatchStarts.shrink();
475  foamPatchSizes.shrink();
476 
477  Pout<< "patchSizes:" << foamPatchSizes << endl;
478  Pout<< "patchStarts:" << foamPatchStarts << endl;
479  Pout<< "nFaces:" << foamNFaces << endl;
480 
481  mapData.resize(nInternalFaces);
482  CCMIOGetEntity(&err, topology, kCCMIOInternalFaces, 0, &id);
483  CCMIOSize size;
484  CCMIOReadFaces(&err, id, kCCMIOInternalFaces, NULL, &size, NULL,
485  kCCMIOStart, kCCMIOEnd);
486  std::vector<int> faces(size);
487  CCMIOReadFaces(&err, id, kCCMIOInternalFaces, &mapID, NULL, &faces[0],
488  kCCMIOStart, kCCMIOEnd);
489  std::vector<int> faceCells(2*nInternalFaces);
490  CCMIOReadFaceCells(&err, id, kCCMIOInternalFaces, &faceCells[0],
491  kCCMIOStart, kCCMIOEnd);
492  CCMIOReadMap(&err, mapID, &mapData[0], kCCMIOStart, kCCMIOEnd);
493  CheckError(err, "Error reading internal faces");
494 
495  // Copy into Foam lists
496  foamFaceMap.setSize(foamNFaces);
497  foamFaces.setSize(foamNFaces);
498  foamOwner.setSize(foamNFaces);
499  foamNeighbour.setSize(foamNFaces);
500 
501  unsigned int pos = 0;
502 
503  for (unsigned int faceI = 0; faceI < nInternalFaces; faceI++)
504  {
505  foamFaceMap[faceI] = mapData[faceI];
506  foamOwner[faceI] = faceCells[2*faceI];
507  foamNeighbour[faceI] = faceCells[2*faceI+1];
508  face& f = foamFaces[faceI];
509 
510  f.setSize(faces[pos++]);
511  forAll(f, fp)
512  {
513  f[fp] = faces[pos++];
514  }
515  }
516 
517  // Read the boundary faces.
518  // ~~~~~~~~~~~~~~~~~~~~~~~~
519 
520  index = 0;
521  label regionI = 0;
522  while
523  (
524  CCMIONextEntity(NULL, topology, kCCMIOBoundaryFaces, &index, &id)
525  == kCCMIONoErr
526  )
527  {
528  CCMIOSize nFaces;
529  CCMIOEntitySize(&err, id, &nFaces, NULL);
530 
531  mapData.resize(nFaces);
532  faceCells.resize(nFaces);
533  CCMIOReadFaces(&err, id, kCCMIOBoundaryFaces, NULL, &size, NULL,
534  kCCMIOStart, kCCMIOEnd);
535  faces.resize(size);
536  CCMIOReadFaces(&err, id, kCCMIOBoundaryFaces, &mapID, NULL, &faces[0],
537  kCCMIOStart, kCCMIOEnd);
538  CCMIOReadFaceCells(&err, id, kCCMIOBoundaryFaces, &faceCells[0],
539  kCCMIOStart, kCCMIOEnd);
540  CCMIOReadMap(&err, mapID, &mapData[0], kCCMIOStart, kCCMIOEnd);
541  CheckError(err, "Error reading boundary faces");
542 
543  // Read prostar id
544  int prostarI;
545  if
546  (
547  CCMIOReadOpti(NULL, id, "ProstarRegionNumber", &prostarI)
548  == kCCMIONoErr
549  )
550  {
551  Pout<< "For region:" << regionI
552  << " found ProstarRegionNumber:" << prostarI << endl;
553  }
554  else
555  {
556  prostarI = regionI;
557 
558  Pout<< "For region:" << regionI
559  << " did not find ProstarRegionNumber entry. Assuming "
560  << prostarI << endl;
561  }
562  prostarToFoamPatch.insert(prostarI, regionI);
563 
564 
565  Pout<< "region:" << regionI
566  << " ProstarRegionNumber:" << prostarI
567  << " foamPatchStart:"
568  << foamPatchStarts[regionI]
569  << " size:"
570  << foamPatchSizes[regionI]
571  << endl;
572 
573  // Copy into Foam list.
574  unsigned int pos = 0;
575 
576  for (unsigned int i = 0; i < nFaces; i++)
577  {
578  label foamFaceI = foamPatchStarts[regionI] + i;
579 
580  foamFaceMap[foamFaceI] = mapData[i];
581  foamOwner[foamFaceI] = faceCells[i];
582  foamNeighbour[foamFaceI] = -1;
583  face& f = foamFaces[foamFaceI];
584 
585  f.setSize(faces[pos++]);
586  forAll(f, fp)
587  {
588  f[fp] = faces[pos++];
589  }
590  }
591  Pout<< endl;
592 
593  regionI++;
594  }
595 }
596 
597 
598 // Main program:
599 
600 int main(int argc, char *argv[])
601 {
603  argList::validArgs.append("ccm26 file");
604 
605 # include <OpenFOAM/setRootCase.H>
606 # include <OpenFOAM/createTime.H>
607 
608 
609  // Foam mesh data
610  // ~~~~~~~~~~~~~~
611 
612  // Coordinates
613  pointField foamPoints;
614  labelList foamPointMap;
615 
616  // Cell type
617  labelList foamCellType;
618  labelList foamCellMap;
619 
620  // Patching info
621  Map<label> prostarToFoamPatch;
622  DynamicList<label> foamPatchSizes;
623  DynamicList<label> foamPatchStarts;
624  // Face connectivity
625  labelList foamFaceMap;
626  labelList foamOwner;
627  labelList foamNeighbour;
628  faceList foamFaces;
629 
630  // Celltypes (not the names of the zones; just the material type)
631  // and patch type names
632  Map<word> foamCellTypeNames;
633  wordList foamPatchTypes;
634  wordList foamPatchNames;
635 
636  {
637  fileName ccmFile(args.additionalArgs()[0]);
638 
639  if (!isFile(ccmFile))
640  {
642  << "Cannot read file " << ccmFile
643  << exit(FatalError);
644  }
645 
646  word ccmExt = ccmFile.ext();
647 
648  if (ccmExt != "ccm" && ccmExt != "ccmg")
649  {
651  << "Illegal extension " << ccmExt << " for file " << ccmFile
652  << nl << "Allowed extensions are '.ccm', '.ccmg'"
653  << exit(FatalError);
654  }
655 
656  // Open the file. Because we did not initialize 'err' we need to pass
657  // in NULL (which always means kCCMIONoErr) and then assign the return
658  // value to 'err'.).
659  CCMIOID root;
660  CCMIOError err = CCMIOOpenFile(NULL, ccmFile.c_str(), kCCMIORead, &root);
661 
662  // We are going to assume that we have a state with a known name.
663  // We could instead use CCMIONextEntity() to walk through all the
664  // states in the file and present the list to the user for selection.
665  CCMIOID state;
666  int stateI = 0;
667  CCMIONextEntity(&err, root, kCCMIOState, &stateI, &state);
668  CheckError(err, "Error opening state");
669 
670  unsigned int size;
671  CCMIOEntityDescription(&err, state, &size, NULL);
672  char *desc = new char[size + 1];
673  CCMIOEntityDescription(&err, state, NULL, desc);
674  Pout<< "Reading state '" << kDefaultState << "' (" << desc << ")"
675  << endl;
676  delete [] desc;
677 
678  // Find the first processor (i has previously been initialized to 0) and
679  // read the mesh and solution information.
680  int i = 0;
681  CCMIOID processor;
682  CCMIONextEntity(&err, state, kCCMIOProcessor, &i, &processor);
683  CCMIOID solution, vertices, topology;
684  CCMIOReadProcessor
685  (
686  &err,
687  processor,
688  &vertices,
689  &topology,
690  NULL,
691  &solution
692  );
693 
694  if (err != kCCMIONoErr)
695  {
696  // Maybe no solution; try again
697  err = kCCMIONoErr;
698  CCMIOReadProcessor
699  (
700  &err,
701  processor,
702  &vertices,
703  &topology,
704  NULL,
705  NULL
706  );
707  if (err != kCCMIONoErr)
708  {
710  << "Could not read the file."
711  << exit(FatalError);
712  }
713  }
714 
715  ReadVertices(err, vertices, foamPointMap, foamPoints);
716 
717  Pout<< "nPoints:" << foamPoints.size() << endl
718  << "bounding box:" << boundBox(foamPoints) << endl
719  << endl;
720 
721  ReadCells
722  (
723  err,
724  topology,
725  foamCellMap,
726  foamCellType,
727  prostarToFoamPatch,
728  foamPatchSizes,
729  foamPatchStarts,
730  foamFaceMap,
731  foamOwner,
732  foamNeighbour,
733  foamFaces
734  );
735 
736  Pout<< "nCells:" << foamCellMap.size() << endl
737  << "nFaces:" << foamOwner.size() << endl
738  << "nPatches:" << foamPatchStarts.size() << endl
739  << "nInternalFaces:" << foamPatchStarts[0] << endl
740  << endl;
741 
742  // Create some default patch names/types. These will be overwritten
743  // by any problem desciption (if it is there)
744  foamPatchTypes.setSize(foamPatchStarts.size());
745  foamPatchNames.setSize(foamPatchStarts.size());
746 
747  forAll(foamPatchNames, i)
748  {
749  foamPatchNames[i] = word("patch") + Foam::name(i);
750  foamPatchTypes[i] = "patch";
751  }
752 
753  // Get problem description
754 
755  CCMIOID problem;
756  int problemI = 0;
757  CCMIONextEntity
758  (
759  &err,
760  root,
761  kCCMIOProblemDescription,
762  &problemI,
763  &problem
764  );
765  CheckError(err, "Error stepping to first problem description");
766 
767  if (CCMIOIsValidEntity(problem)) // if we have a problem description
768  {
769  ReadProblem
770  (
771  err,
772  problem,
773  prostarToFoamPatch,
774 
775  foamCellTypeNames,
776  foamPatchTypes,
777  foamPatchNames
778  );
779  }
780 
781 
782  CCMIOCloseFile(&err, vertices);
783  CCMIOCloseFile(&err, topology);
784  CCMIOCloseFile(&err, solution);
785  CCMIOCloseFile(&err, root);
786  }
787 
788 
789  Pout<< "foamPatchNames:" << foamPatchNames << endl;
790 
791 
792  Pout<< "foamOwner : min:" << min(foamOwner)
793  << " max:" << max(foamOwner)
794  << nl
795  << "foamNeighbour : min:" << min(foamNeighbour)
796  << " max:" << max(foamNeighbour)
797  << nl
798  << "foamCellType : min:" << min(foamCellType)
799  << " max:" << max(foamCellType)
800  << nl << endl;
801 
802 
803 
804  // We now have extracted all info from CCMIO:
805  // - coordinates (points)
806  // - face to point addressing (faces)
807  // - face to cell addressing (owner, neighbour)
808  // - cell based data (cellId)
809 
810 
811  // Renumber vertex labels to Foam point labels
812  {
813  label maxCCMPointI = max(foamPointMap);
814  labelList toFoamPoints(invert(maxCCMPointI+1, foamPointMap));
815 
816  forAll(foamFaces, faceI)
817  {
818  inplaceRenumber(toFoamPoints, foamFaces[faceI]);
819  }
820  }
821 
822  // Renumber cell labels
823  {
824  label maxCCMCellI = max(foamCellMap);
825  labelList toFoamCells(invert(maxCCMCellI+1, foamCellMap));
826 
827  inplaceRenumber(toFoamCells, foamOwner);
828  inplaceRenumber(toFoamCells, foamNeighbour);
829  }
830 
831 
832  //
833  // Basic mesh info complete. Now convert to Foam convention:
834  // - owner < neighbour
835  // - face vertices such that normal points away from owner
836  // - order faces: upper-triangular for internal faces; boundary faces after
837  // internal faces
838  //
839 
840  // Set owner/neighbour so owner < neighbour
841  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
842 
843  forAll(foamNeighbour, faceI)
844  {
845  label nbr = foamNeighbour[faceI];
846  label own = foamOwner[faceI];
847 
848  if (nbr >= foamCellType.size() || own >= foamCellType.size())
849  {
851  << "face:" << faceI
852  << " nbr:" << nbr
853  << " own:" << own
854  << " nCells:" << foamCellType.size()
855  << exit(FatalError);
856  }
857 
858  if (nbr >= 0)
859  {
860  if (nbr < own)
861  {
862  foamOwner[faceI] = foamNeighbour[faceI];
863  foamNeighbour[faceI] = own;
864  foamFaces[faceI] = foamFaces[faceI].reverseFace();
865  }
866  }
867 
868 
869  // And check the face
870  const face& f = foamFaces[faceI];
871 
872  forAll(f, fp)
873  {
874  if (f[fp] < 0 || f[fp] >= foamPoints.size())
875  {
876  FatalErrorIn(args.executable()) << "face:" << faceI
877  << " f:" << f << abort(FatalError);
878  }
879  }
880  }
881 
882 
883  // Do upper-triangular ordering
884  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
885 
886  {
887  // Create cells (inverse of face-to-cell addressing)
888  cellList foamCells;
889  primitiveMesh::calcCells
890  (
891  foamCells,
892  foamOwner,
893  foamNeighbour,
894  foamCellType.size()
895  );
896 
897  // Determine face order for upper-triangular ordering
898  labelList oldToNew
899  (
900  getInternalFaceOrder
901  (
902  foamCells,
903  foamOwner,
904  foamNeighbour
905  )
906  );
907 
908  // Reorder faces accordingly
909  forAll(foamCells, cellI)
910  {
911  inplaceRenumber(oldToNew, foamCells[cellI]);
912  }
913 
914  // Reorder faces.
915  inplaceReorder(oldToNew, foamFaces);
916  inplaceReorder(oldToNew, foamOwner);
917  inplaceReorder(oldToNew, foamNeighbour);
918  }
919 
920 
921  // Construct fvMesh
922  // ~~~~~~~~~~~~~~~~
923 
924  fvMesh mesh
925  (
926  IOobject
927  (
929  runTime.constant(),
930  runTime
931  ),
932  xferMove<pointField>(foamPoints),
933  xferMove<faceList>(foamFaces),
934  xferCopy<labelList>(foamOwner),
935  xferMove<labelList>(foamNeighbour)
936  );
937 
938  // Create patches. Use patch types to determine what Foam types to generate.
939  List<polyPatch*> newPatches(foamPatchNames.size());
940 
941  label meshFaceI = foamPatchStarts[0];
942 
943  forAll(newPatches, patchI)
944  {
945  const word& patchName = foamPatchNames[patchI];
946  const word& patchType = foamPatchTypes[patchI];
947 
948  Pout<< "Patch:" << patchName << " start at:" << meshFaceI
949  << " size:" << foamPatchSizes[patchI]
950  << " end at:" << meshFaceI+foamPatchSizes[patchI]
951  << endl;
952 
953  if (patchType == "wall")
954  {
955  newPatches[patchI] =
956  new wallPolyPatch
957  (
958  patchName,
959  foamPatchSizes[patchI],
960  meshFaceI,
961  patchI,
963  );
964  }
965  else if (patchType == "symmetryplane")
966  {
967  newPatches[patchI] =
969  (
970  patchName,
971  foamPatchSizes[patchI],
972  meshFaceI,
973  patchI,
975  );
976  }
977  else if (patchType == "empty")
978  {
979  // Note: not ccm name, introduced by us above.
980  newPatches[patchI] =
981  new emptyPolyPatch
982  (
983  patchName,
984  foamPatchSizes[patchI],
985  meshFaceI,
986  patchI,
988  );
989  }
990  else
991  {
992  // All other ccm types become straight polyPatch:
993  // 'inlet', 'outlet', 'pressured'.
994  newPatches[patchI] =
995  new polyPatch
996  (
997  patchName,
998  foamPatchSizes[patchI],
999  meshFaceI,
1000  patchI,
1001  mesh.boundaryMesh()
1002  );
1003  }
1004 
1005  meshFaceI += foamPatchSizes[patchI];
1006  }
1007 
1008  if (meshFaceI != foamOwner.size())
1009  {
1011  << "meshFaceI:" << meshFaceI
1012  << " nFaces:" << foamOwner.size()
1013  << abort(FatalError);
1014  }
1015  mesh.addFvPatches(newPatches);
1016 
1017 
1018 
1019  // Construct sets and zones from types
1020  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1021 
1022  label maxType = max(foamCellType);
1023  label minType = min(foamCellType);
1024 
1025  if (maxType > minType)
1026  {
1027  // From foamCellType physical region to Foam cellZone
1028  Map<label> typeToZone;
1029  // Storage for cell zones.
1030  List<DynamicList<label> > zoneCells(0);
1031 
1032  forAll(foamCellType, cellI)
1033  {
1034  storeCellInZone
1035  (
1036  cellI,
1037  foamCellType[cellI],
1038  typeToZone,
1039  zoneCells
1040  );
1041  }
1042 
1043  mesh.cellZones().clear();
1044  mesh.cellZones().setSize(typeToZone.size());
1045 
1046  label nValidCellZones = 0;
1047 
1048  forAllConstIter(Map<label>, typeToZone, iter)
1049  {
1050  label type = iter.key();
1051  label zoneI = iter();
1052 
1053  zoneCells[zoneI].shrink();
1054 
1055  word zoneName = "cellZone_" + name(type);
1056 
1057  Info<< "Writing zone " << type
1058  << " size " << zoneCells[zoneI].size()
1059  << " to cellZone "
1060  << zoneName << " and cellSet " << zoneName
1061  << endl;
1062 
1063  cellSet cset(mesh, zoneName, zoneCells[zoneI]);
1064  cset.write();
1065 
1066  mesh.cellZones().set
1067  (
1068  nValidCellZones,
1069  new cellZone
1070  (
1071  zoneName,
1072  zoneCells[zoneI],
1073  nValidCellZones,
1074  mesh.cellZones()
1075  )
1076  );
1077  nValidCellZones++;
1078  }
1080  }
1081 
1082 
1083  Info<< "Writing mesh to " << mesh.objectRegistry::objectPath()
1084  << "..." << nl << endl;
1085 
1086 
1087  // Construct field with calculated bc to hold Star cell Id.
1088  volScalarField cellIdField
1089  (
1090  IOobject
1091  (
1092  "cellId",
1093  runTime.timeName(),
1094  mesh,
1097  ),
1098  mesh,
1099  dimensionedScalar("cellId", dimless, 0.0)
1100  );
1101 
1102  forAll(foamCellMap, cellI)
1103  {
1104  cellIdField[cellI] = foamCellMap[cellI];
1105  }
1106 
1107  // Construct field with calculated bc to hold cell type.
1108  volScalarField cellTypeField
1109  (
1110  IOobject
1111  (
1112  "cellType",
1113  runTime.timeName(),
1114  mesh,
1117  ),
1118  mesh,
1119  dimensionedScalar("cellType", dimless, 0.0)
1120  );
1121 
1122  forAll(foamCellType, cellI)
1123  {
1124  cellTypeField[cellI] = foamCellType[cellI];
1125  }
1126 
1127  Info<< "Writing cellIds as volScalarField to " << cellIdField.objectPath()
1128  << "..." << nl << endl;
1129  mesh.write();
1130 
1131  Info<< "End\n" << endl;
1132 
1133  return 0;
1134 }
1135 
1136 
1137 // ************************ vim: set sw=4 sts=4 et: ************************ //