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