gmshToFoam.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-2019 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 Application
25  gmshToFoam
26 
27 Description
28  Reads .msh file as written by Gmsh.
29 
30  Needs surface elements on mesh to be present and aligned with outside faces
31  of the mesh. I.e. if the mesh is hexes, the outside faces need to be
32  quads.
33 
34  Note: There is something seriously wrong with the ordering written in the
35  .msh file. Normal operation is to check the ordering and invert prisms
36  and hexes if found to be wrong way round.
37  Use the -keepOrientation to keep the raw information.
38 
39  Note: The code now uses the element (cell,face) physical region id number
40  to create cell zones and faces zones (similar to
41  fluentMeshWithInternalFaces).
42 
43  A use of the cell zone information, is for field initialization with the
44  "setFields" utility. see the classes: topoSetSource, zoneToCell.
45 
46 \*---------------------------------------------------------------------------*/
47 
48 #include "argList.H"
49 #include "polyMesh.H"
50 #include "Time.H"
51 #include "polyMesh.H"
52 #include "IFstream.H"
53 #include "cellModeller.H"
54 #include "repatchPolyTopoChanger.H"
55 #include "cellSet.H"
56 #include "faceSet.H"
57 
58 using namespace Foam;
59 
60 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
61 
62 // Element type numbers
63 static label MSHTRI = 2;
64 static label MSHQUAD = 3;
65 static label MSHTET = 4;
66 static label MSHPYR = 7;
67 static label MSHPRISM = 6;
68 static label MSHHEX = 5;
69 
70 
71 // Skips till end of section. Returns false if end of file.
72 bool skipSection(IFstream& inFile)
73 {
74  string line;
75  do
76  {
77  inFile.getLine(line);
78 
79  if (!inFile.good())
80  {
81  return false;
82  }
83  }
84  while (line.size() < 4 || line.substr(0, 4) != "$End");
85 
86  return true;
87 }
88 
89 
90 void renumber
91 (
92  const Map<label>& mshToFoam,
93  labelList& labels
94 )
95 {
96  forAll(labels, labelI)
97  {
98  labels[labelI] = mshToFoam[labels[labelI]];
99  }
100 }
101 
102 
103 // Find face in pp which uses all vertices in meshF (in mesh point labels)
104 label findFace(const primitivePatch& pp, const labelList& meshF)
105 {
106  const Map<label>& meshPointMap = pp.meshPointMap();
107 
108  // meshF[0] in pp labels.
109  if (!meshPointMap.found(meshF[0]))
110  {
111  Warning<< "Not using gmsh face " << meshF
112  << " since zero vertex is not on boundary of polyMesh" << endl;
113  return -1;
114  }
115 
116  // Find faces using first point
117  const labelList& pFaces = pp.pointFaces()[meshPointMap[meshF[0]]];
118 
119  // Go through all these faces and check if there is one which uses all of
120  // meshF vertices (in any order ;-)
121  forAll(pFaces, i)
122  {
123  label facei = pFaces[i];
124 
125  const face& f = pp[facei];
126 
127  // Count uses of vertices of meshF for f
128  label nMatched = 0;
129 
130  forAll(f, fp)
131  {
132  if (findIndex(meshF, f[fp]) != -1)
133  {
134  nMatched++;
135  }
136  }
137 
138  if (nMatched == meshF.size())
139  {
140  return facei;
141  }
142  }
143 
144  return -1;
145 }
146 
147 
148 // Same but find internal face. Expensive addressing.
149 label findInternalFace(const primitiveMesh& mesh, const labelList& meshF)
150 {
151  const labelList& pFaces = mesh.pointFaces()[meshF[0]];
152 
153  forAll(pFaces, i)
154  {
155  label facei = pFaces[i];
156 
157  const face& f = mesh.faces()[facei];
158 
159  // Count uses of vertices of meshF for f
160  label nMatched = 0;
161 
162  forAll(f, fp)
163  {
164  if (findIndex(meshF, f[fp]) != -1)
165  {
166  nMatched++;
167  }
168  }
169 
170  if (nMatched == meshF.size())
171  {
172  return facei;
173  }
174  }
175  return -1;
176 }
177 
178 
179 // Determine whether cell is inside-out by checking for any wrong-oriented
180 // face.
181 bool correctOrientation(const pointField& points, const cellShape& shape)
182 {
183  // Get centre of shape.
184  point cc(shape.centre(points));
185 
186  // Get outwards pointing faces.
187  faceList faces(shape.faces());
188 
189  forAll(faces, i)
190  {
191  const face& f = faces[i];
192 
193  const vector a(f.area(points));
194 
195  // Check if vector from any point on face to cc points outwards
196  if (((points[f[0]] - cc) & a) < 0)
197  {
198  // Incorrectly oriented
199  return false;
200  }
201  }
202 
203  return true;
204 }
205 
206 
207 void storeCellInZone
208 (
209  const label regPhys,
210  const label celli,
211  Map<label>& physToZone,
212 
213  labelList& zoneToPhys,
214  List<DynamicList<label>>& zoneCells
215 )
216 {
217  Map<label>::const_iterator zoneFnd = physToZone.find(regPhys);
218 
219  if (zoneFnd == physToZone.end())
220  {
221  // New region. Allocate zone for it.
222  label zoneI = zoneCells.size();
223  zoneCells.setSize(zoneI+1);
224  zoneToPhys.setSize(zoneI+1);
225 
226  Info<< "Mapping region " << regPhys << " to Foam cellZone "
227  << zoneI << endl;
228  physToZone.insert(regPhys, zoneI);
229 
230  zoneToPhys[zoneI] = regPhys;
231  zoneCells[zoneI].append(celli);
232  }
233  else
234  {
235  // Existing zone for region
236  zoneCells[zoneFnd()].append(celli);
237  }
238 }
239 
240 
241 // Reads mesh format
242 scalar readMeshFormat(IFstream& inFile)
243 {
244  Info<< "Starting to read mesh format at line "
245  << inFile.lineNumber()
246  << endl;
247 
248  string line;
249  inFile.getLine(line);
250  IStringStream lineStr(line);
251 
252  scalar version;
253  label asciiFlag, nBytes;
254  lineStr >> version >> asciiFlag >> nBytes;
255 
256  Info<< "Read format version " << version << " ascii " << asciiFlag << endl;
257 
258  if (asciiFlag != 0)
259  {
260  FatalIOErrorInFunction(inFile)
261  << "Can only read ascii msh files."
262  << exit(FatalIOError);
263  }
264 
265  inFile.getLine(line);
266  IStringStream tagStr(line);
267  word tag(tagStr);
268 
269  if (tag != "$EndMeshFormat")
270  {
271  FatalIOErrorInFunction(inFile)
272  << "Did not find $ENDNOD tag on line "
273  << inFile.lineNumber() << exit(FatalIOError);
274  }
275  Info<< endl;
276 
277  return version;
278 }
279 
280 
281 // Reads points and map
282 void readPoints(IFstream& inFile, pointField& points, Map<label>& mshToFoam)
283 {
284  Info<< "Starting to read points at line " << inFile.lineNumber() << endl;
285 
286  string line;
287  inFile.getLine(line);
288  IStringStream lineStr(line);
289 
290  label nVerts;
291  lineStr >> nVerts;
292 
293  Info<< "Vertices to be read:" << nVerts << endl;
294 
295  points.setSize(nVerts);
296  mshToFoam.resize(2*nVerts);
297 
298  for (label pointi = 0; pointi < nVerts; pointi++)
299  {
300  label mshLabel;
301  scalar xVal, yVal, zVal;
302 
303  string line;
304  inFile.getLine(line);
305  IStringStream lineStr(line);
306 
307  lineStr >> mshLabel >> xVal >> yVal >> zVal;
308 
309  point& pt = points[pointi];
310 
311  pt.x() = xVal;
312  pt.y() = yVal;
313  pt.z() = zVal;
314 
315  mshToFoam.insert(mshLabel, pointi);
316  }
317 
318  Info<< "Vertices read:" << mshToFoam.size() << endl;
319 
320  inFile.getLine(line);
321  IStringStream tagStr(line);
322  word tag(tagStr);
323 
324  if (tag != "$ENDNOD" && tag != "$EndNodes")
325  {
326  FatalIOErrorInFunction(inFile)
327  << "Did not find $ENDNOD tag on line "
328  << inFile.lineNumber() << exit(FatalIOError);
329  }
330  Info<< endl;
331 }
332 
333 
334 // Reads physical names
335 void readPhysNames(IFstream& inFile, Map<word>& physicalNames)
336 {
337  Info<< "Starting to read physical names at line " << inFile.lineNumber()
338  << endl;
339 
340  string line;
341  inFile.getLine(line);
342  IStringStream lineStr(line);
343 
344  label nNames;
345  lineStr >> nNames;
346 
347  Info<< "Physical names:" << nNames << endl;
348 
349  physicalNames.resize(nNames);
350 
351  for (label i = 0; i < nNames; i++)
352  {
353  label regionI;
354  string regionName;
355 
356  string line;
357  inFile.getLine(line);
358  IStringStream lineStr(line);
359  label nSpaces = lineStr.str().count(' ');
360 
361  if (nSpaces == 1)
362  {
363  lineStr >> regionI >> regionName;
364 
365  Info<< " " << regionI << '\t'
366  << string::validate<word>(regionName) << endl;
367  }
368  else if (nSpaces == 2)
369  {
370  // >= Gmsh2.4 physical types has tag in front.
371  label physType;
372  lineStr >> physType >> regionI >> regionName;
373  if (physType == 1)
374  {
375  Info<< " " << "Line " << regionI << '\t'
376  << string::validate<word>(regionName) << endl;
377  }
378  else if (physType == 2)
379  {
380  Info<< " " << "Surface " << regionI << '\t'
381  << string::validate<word>(regionName) << endl;
382  }
383  else if (physType == 3)
384  {
385  Info<< " " << "Volume " << regionI << '\t'
386  << string::validate<word>(regionName) << endl;
387  }
388  }
389 
390  physicalNames.insert(regionI, string::validate<word>(regionName));
391  }
392 
393  inFile.getLine(line);
394  IStringStream tagStr(line);
395  word tag(tagStr);
396 
397  if (tag != "$EndPhysicalNames")
398  {
399  FatalIOErrorInFunction(inFile)
400  << "Did not find $EndPhysicalNames tag on line "
401  << inFile.lineNumber() << exit(FatalIOError);
402  }
403  Info<< endl;
404 }
405 
406 
407 // Reads cells and patch faces
408 void readCells
409 (
410  const scalar versionFormat,
411  const bool keepOrientation,
412  const pointField& points,
413  const Map<label>& mshToFoam,
414  IFstream& inFile,
415  cellShapeList& cells,
416 
417  labelList& patchToPhys,
418  List<DynamicList<face>>& patchFaces,
419 
420  labelList& zoneToPhys,
421  List<DynamicList<label>>& zoneCells
422 )
423 {
424  Info<< "Starting to read cells at line " << inFile.lineNumber() << endl;
425 
426  const cellModel& hex = *(cellModeller::lookup("hex"));
427  const cellModel& prism = *(cellModeller::lookup("prism"));
428  const cellModel& pyr = *(cellModeller::lookup("pyr"));
429  const cellModel& tet = *(cellModeller::lookup("tet"));
430 
431  face triPoints(3);
432  face quadPoints(4);
433  labelList tetPoints(4);
434  labelList pyrPoints(5);
435  labelList prismPoints(6);
436  labelList hexPoints(8);
437 
438 
439  string line;
440  inFile.getLine(line);
441  IStringStream lineStr(line);
442 
443  label nElems;
444  lineStr >> nElems;
445 
446  Info<< "Cells to be read:" << nElems << endl << endl;
447 
448 
449  // Storage for all cells. Too big. Shrink later
450  cells.setSize(nElems);
451 
452  label celli = 0;
453  label nTet = 0;
454  label nPyr = 0;
455  label nPrism = 0;
456  label nHex = 0;
457 
458 
459  // From gmsh physical region to Foam patch
460  Map<label> physToPatch;
461 
462  // From gmsh physical region to Foam cellZone
463  Map<label> physToZone;
464 
465 
466  for (label elemI = 0; elemI < nElems; elemI++)
467  {
468  string line;
469  inFile.getLine(line);
470  IStringStream lineStr(line);
471 
472  label elmNumber, elmType, regPhys;
473 
474  if (versionFormat >= 2)
475  {
476  lineStr >> elmNumber >> elmType;
477 
478  label nTags;
479  lineStr>> nTags;
480 
481  if (nTags > 0)
482  {
483  // Assume the first tag is the physical surface
484  lineStr >> regPhys;
485  for (label i = 1; i < nTags; i++)
486  {
487  label dummy;
488  lineStr>> dummy;
489  }
490  }
491  }
492  else
493  {
494  label regElem, nNodes;
495  lineStr >> elmNumber >> elmType >> regPhys >> regElem >> nNodes;
496  }
497 
498  // regPhys on surface elements is region number.
499 
500  if (elmType == MSHTRI)
501  {
502  lineStr >> triPoints[0] >> triPoints[1] >> triPoints[2];
503 
504  renumber(mshToFoam, triPoints);
505 
506  Map<label>::iterator regFnd = physToPatch.find(regPhys);
507 
508  label patchi = -1;
509  if (regFnd == physToPatch.end())
510  {
511  // New region. Allocate patch for it.
512  patchi = patchFaces.size();
513 
514  patchFaces.setSize(patchi + 1);
515  patchToPhys.setSize(patchi + 1);
516 
517  Info<< "Mapping region " << regPhys << " to Foam patch "
518  << patchi << endl;
519  physToPatch.insert(regPhys, patchi);
520  patchToPhys[patchi] = regPhys;
521  }
522  else
523  {
524  // Existing patch for region
525  patchi = regFnd();
526  }
527 
528  // Add triangle to correct patchFaces.
529  patchFaces[patchi].append(triPoints);
530  }
531  else if (elmType == MSHQUAD)
532  {
533  lineStr
534  >> quadPoints[0] >> quadPoints[1] >> quadPoints[2]
535  >> quadPoints[3];
536 
537  renumber(mshToFoam, quadPoints);
538 
539  Map<label>::iterator regFnd = physToPatch.find(regPhys);
540 
541  label patchi = -1;
542  if (regFnd == physToPatch.end())
543  {
544  // New region. Allocate patch for it.
545  patchi = patchFaces.size();
546 
547  patchFaces.setSize(patchi + 1);
548  patchToPhys.setSize(patchi + 1);
549 
550  Info<< "Mapping region " << regPhys << " to Foam patch "
551  << patchi << endl;
552  physToPatch.insert(regPhys, patchi);
553  patchToPhys[patchi] = regPhys;
554  }
555  else
556  {
557  // Existing patch for region
558  patchi = regFnd();
559  }
560 
561  // Add quad to correct patchFaces.
562  patchFaces[patchi].append(quadPoints);
563  }
564  else if (elmType == MSHTET)
565  {
566  storeCellInZone
567  (
568  regPhys,
569  celli,
570  physToZone,
571  zoneToPhys,
572  zoneCells
573  );
574 
575  lineStr
576  >> tetPoints[0] >> tetPoints[1] >> tetPoints[2]
577  >> tetPoints[3];
578 
579  renumber(mshToFoam, tetPoints);
580 
581  cells[celli++] = cellShape(tet, tetPoints);
582 
583  nTet++;
584  }
585  else if (elmType == MSHPYR)
586  {
587  storeCellInZone
588  (
589  regPhys,
590  celli,
591  physToZone,
592  zoneToPhys,
593  zoneCells
594  );
595 
596  lineStr
597  >> pyrPoints[0] >> pyrPoints[1] >> pyrPoints[2]
598  >> pyrPoints[3] >> pyrPoints[4];
599 
600  renumber(mshToFoam, pyrPoints);
601 
602  cells[celli++] = cellShape(pyr, pyrPoints);
603 
604  nPyr++;
605  }
606  else if (elmType == MSHPRISM)
607  {
608  storeCellInZone
609  (
610  regPhys,
611  celli,
612  physToZone,
613  zoneToPhys,
614  zoneCells
615  );
616 
617  lineStr
618  >> prismPoints[0] >> prismPoints[1] >> prismPoints[2]
619  >> prismPoints[3] >> prismPoints[4] >> prismPoints[5];
620 
621  renumber(mshToFoam, prismPoints);
622 
623  cells[celli] = cellShape(prism, prismPoints);
624 
625  const cellShape& cell = cells[celli];
626 
627  if (!keepOrientation && !correctOrientation(points, cell))
628  {
629  Info<< "Inverting prism " << celli << endl;
630  // Reorder prism.
631  prismPoints[0] = cell[0];
632  prismPoints[1] = cell[2];
633  prismPoints[2] = cell[1];
634  prismPoints[3] = cell[3];
635  prismPoints[4] = cell[4];
636  prismPoints[5] = cell[5];
637 
638  cells[celli] = cellShape(prism, prismPoints);
639  }
640 
641  celli++;
642 
643  nPrism++;
644  }
645  else if (elmType == MSHHEX)
646  {
647  storeCellInZone
648  (
649  regPhys,
650  celli,
651  physToZone,
652  zoneToPhys,
653  zoneCells
654  );
655 
656  lineStr
657  >> hexPoints[0] >> hexPoints[1]
658  >> hexPoints[2] >> hexPoints[3]
659  >> hexPoints[4] >> hexPoints[5]
660  >> hexPoints[6] >> hexPoints[7];
661 
662  renumber(mshToFoam, hexPoints);
663 
664  cells[celli] = cellShape(hex, hexPoints);
665 
666  const cellShape& cell = cells[celli];
667 
668  if (!keepOrientation && !correctOrientation(points, cell))
669  {
670  Info<< "Inverting hex " << celli << endl;
671  // Reorder hex.
672  hexPoints[0] = cell[4];
673  hexPoints[1] = cell[5];
674  hexPoints[2] = cell[6];
675  hexPoints[3] = cell[7];
676  hexPoints[4] = cell[0];
677  hexPoints[5] = cell[1];
678  hexPoints[6] = cell[2];
679  hexPoints[7] = cell[3];
680 
681  cells[celli] = cellShape(hex, hexPoints);
682  }
683 
684  celli++;
685 
686  nHex++;
687  }
688  else
689  {
690  Info<< "Unhandled element " << elmType << " at line "
691  << inFile.lineNumber() << endl;
692  }
693  }
694 
695 
696  inFile.getLine(line);
697  IStringStream tagStr(line);
698  word tag(tagStr);
699 
700  if (tag != "$ENDELM" && tag != "$EndElements")
701  {
702  FatalIOErrorInFunction(inFile)
703  << "Did not find $ENDELM tag on line "
704  << inFile.lineNumber() << exit(FatalIOError);
705  }
706 
707 
708  cells.setSize(celli);
709 
710  forAll(patchFaces, patchi)
711  {
712  patchFaces[patchi].shrink();
713  }
714 
715 
716  Info<< "Cells:" << endl
717  << " total:" << cells.size() << endl
718  << " hex :" << nHex << endl
719  << " prism:" << nPrism << endl
720  << " pyr :" << nPyr << endl
721  << " tet :" << nTet << endl
722  << endl;
723 
724  if (cells.size() == 0)
725  {
726  FatalIOErrorInFunction(inFile)
727  << "No cells read from file " << inFile.name() << nl
728  << "Does your file specify any 3D elements (hex=" << MSHHEX
729  << ", prism=" << MSHPRISM << ", pyramid=" << MSHPYR
730  << ", tet=" << MSHTET << ")?" << nl
731  << "Perhaps you have not exported the 3D elements?"
732  << exit(FatalIOError);
733  }
734 
735  Info<< "CellZones:" << nl
736  << "Zone\tSize" << endl;
737 
738  forAll(zoneCells, zoneI)
739  {
740  zoneCells[zoneI].shrink();
741 
742  const labelList& zCells = zoneCells[zoneI];
743 
744  if (zCells.size())
745  {
746  Info<< " " << zoneI << '\t' << zCells.size() << endl;
747  }
748  }
749  Info<< endl;
750 }
751 
752 
753 
754 int main(int argc, char *argv[])
755 {
757  argList::validArgs.append(".msh file");
759  (
760  "keepOrientation",
761  "retain raw orientation for prisms/hexs"
762  );
763 
764  #include "addRegionOption.H"
765 
766  #include "setRootCase.H"
767  #include "createTime.H"
768 
770 
771  if (args.optionReadIfPresent("region", regionName))
772  {
773  Foam::Info
774  << "Creating polyMesh for region " << regionName << endl;
775  }
776  else
777  {
778  regionName = Foam::polyMesh::defaultRegion;
779  }
780 
781  const bool keepOrientation = args.optionFound("keepOrientation");
782  IFstream inFile(args[1]);
783 
784  // Storage for points
786  Map<label> mshToFoam;
787 
788  // Storage for all cells.
790 
791  // Map from patch to gmsh physical region
792  labelList patchToPhys;
793  // Storage for patch faces.
794  List<DynamicList<face>> patchFaces(0);
795 
796  // Map from cellZone to gmsh physical region
797  labelList zoneToPhys;
798  // Storage for cell zones.
799  List<DynamicList<label>> zoneCells(0);
800 
801  // Name per physical region
802  Map<word> physicalNames;
803 
804  // Version 1 or 2 format
805  scalar versionFormat = 1;
806 
807 
808  while (inFile.good())
809  {
810  string line;
811  inFile.getLine(line);
812  IStringStream lineStr(line);
813 
814  word tag(lineStr);
815 
816  if (tag == "$MeshFormat")
817  {
818  versionFormat = readMeshFormat(inFile);
819  }
820  else if (tag == "$PhysicalNames")
821  {
822  readPhysNames(inFile, physicalNames);
823  }
824  else if (tag == "$NOD" || tag == "$Nodes")
825  {
826  readPoints(inFile, points, mshToFoam);
827  }
828  else if (tag == "$ELM" || tag == "$Elements")
829  {
830  readCells
831  (
832  versionFormat,
833  keepOrientation,
834  points,
835  mshToFoam,
836  inFile,
837  cells,
838  patchToPhys,
839  patchFaces,
840  zoneToPhys,
841  zoneCells
842  );
843  }
844  else
845  {
846  Info<< "Skipping tag " << tag << " at line "
847  << inFile.lineNumber()
848  << endl;
849 
850  if (!skipSection(inFile))
851  {
852  break;
853  }
854  }
855  }
856 
857 
858  label nValidCellZones = 0;
859 
860  forAll(zoneCells, zoneI)
861  {
862  if (zoneCells[zoneI].size())
863  {
864  nValidCellZones++;
865  }
866  }
867 
868 
869  // Problem is that the orientation of the patchFaces does not have to
870  // be consistent with the outwards orientation of the mesh faces. So
871  // we have to construct the mesh in two stages:
872  // 1. define mesh with all boundary faces in one patch
873  // 2. use the read patchFaces to find the corresponding boundary face
874  // and repatch it.
875 
876 
877  // Create correct number of patches
878  // (but without any faces in it)
879  faceListList boundaryFaces(patchFaces.size());
880 
881  wordList boundaryPatchNames(boundaryFaces.size());
882 
883  forAll(boundaryPatchNames, patchi)
884  {
885  label physReg = patchToPhys[patchi];
886 
887  Map<word>::const_iterator iter = physicalNames.find(physReg);
888 
889  if (iter != physicalNames.end())
890  {
891  boundaryPatchNames[patchi] = iter();
892  }
893  else
894  {
895  boundaryPatchNames[patchi] = word("patch") + name(patchi);
896  }
897  Info<< "Patch " << patchi << " gets name "
898  << boundaryPatchNames[patchi] << endl;
899  }
900  Info<< endl;
901 
902  wordList boundaryPatchTypes(boundaryFaces.size(), polyPatch::typeName);
903  word defaultFacesName = "defaultFaces";
904  word defaultFacesType = polyPatch::typeName;
905  wordList boundaryPatchPhysicalTypes
906  (
907  boundaryFaces.size(),
908  polyPatch::typeName
909  );
910 
911  polyMesh mesh
912  (
913  IOobject
914  (
915  regionName,
916  runTime.constant(),
917  runTime
918  ),
919  move(points),
920  cells,
921  boundaryFaces,
922  boundaryPatchNames,
923  boundaryPatchTypes,
926  boundaryPatchPhysicalTypes
927  );
928 
929  repatchPolyTopoChanger repatcher(mesh);
930 
931  // Now use the patchFaces to patch up the outside faces of the mesh.
932 
933  // Get the patch for all the outside faces (= default patch added as last)
934  const polyPatch& pp = mesh.boundaryMesh().last();
935 
936  // Storage for faceZones.
937  List<DynamicList<label>> zoneFaces(patchFaces.size());
938 
939 
940  // Go through all the patchFaces and find corresponding face in pp.
941  forAll(patchFaces, patchi)
942  {
943  const DynamicList<face>& pFaces = patchFaces[patchi];
944 
945  Info<< "Finding faces of patch " << patchi << endl;
946 
947  forAll(pFaces, i)
948  {
949  const face& f = pFaces[i];
950 
951  // Find face in pp using all vertices of f.
952  label patchFacei = findFace(pp, f);
953 
954  if (patchFacei != -1)
955  {
956  label meshFacei = pp.start() + patchFacei;
957 
958  repatcher.changePatchID(meshFacei, patchi);
959  }
960  else
961  {
962  // Maybe internal face? If so add to faceZone with same index
963  // - might be useful.
964  label meshFacei = findInternalFace(mesh, f);
965 
966  if (meshFacei != -1)
967  {
968  zoneFaces[patchi].append(meshFacei);
969  }
970  else
971  {
973  << "Could not match gmsh face " << f
974  << " to any of the interior or exterior faces"
975  << " that share the same 0th point" << endl;
976  }
977  }
978  }
979  }
980  Info<< nl;
981 
982  // Face zones
983  label nValidFaceZones = 0;
984 
985  Info<< "FaceZones:" << nl
986  << "Zone\tSize" << endl;
987 
988  forAll(zoneFaces, zoneI)
989  {
990  zoneFaces[zoneI].shrink();
991 
992  const labelList& zFaces = zoneFaces[zoneI];
993 
994  if (zFaces.size())
995  {
996  nValidFaceZones++;
997 
998  Info<< " " << zoneI << '\t' << zFaces.size() << endl;
999  }
1000  }
1001  Info<< endl;
1002 
1003 
1004  // Get polyMesh to write to constant
1005 
1007 
1008  repatcher.repatch();
1009 
1010  List<cellZone*> cz;
1011  List<faceZone*> fz;
1012 
1013  // Construct and add the zones. Note that cell ordering does not change
1014  // because of repatch() and neither does internal faces so we can
1015  // use the zoneCells/zoneFaces as is.
1016 
1017  if (nValidCellZones > 0)
1018  {
1019  cz.setSize(nValidCellZones);
1020 
1021  nValidCellZones = 0;
1022 
1023  forAll(zoneCells, zoneI)
1024  {
1025  if (zoneCells[zoneI].size())
1026  {
1027  label physReg = zoneToPhys[zoneI];
1028 
1029  Map<word>::const_iterator iter = physicalNames.find(physReg);
1030 
1031  word zoneName = "cellZone_" + name(zoneI);
1032  if (iter != physicalNames.end())
1033  {
1034  zoneName = iter();
1035  }
1036 
1037  Info<< "Writing zone " << zoneI << " to cellZone "
1038  << zoneName << " and cellSet"
1039  << endl;
1040 
1041  cellSet cset(mesh, zoneName, zoneCells[zoneI]);
1042  cset.write();
1043 
1044  cz[nValidCellZones] = new cellZone
1045  (
1046  zoneName,
1047  zoneCells[zoneI],
1048  nValidCellZones,
1049  mesh.cellZones()
1050  );
1051  nValidCellZones++;
1052  }
1053  }
1054  }
1055 
1056  if (nValidFaceZones > 0)
1057  {
1058  fz.setSize(nValidFaceZones);
1059 
1060  nValidFaceZones = 0;
1061 
1062  forAll(zoneFaces, zoneI)
1063  {
1064  if (zoneFaces[zoneI].size())
1065  {
1066  label physReg = patchToPhys[zoneI];
1067 
1068  Map<word>::const_iterator iter = physicalNames.find(physReg);
1069 
1070  word zoneName = "faceZone_" + name(zoneI);
1071  if (iter != physicalNames.end())
1072  {
1073  zoneName = iter();
1074  }
1075 
1076  Info<< "Writing zone " << zoneI << " to faceZone "
1077  << zoneName << " and faceSet"
1078  << endl;
1079 
1080  faceSet fset(mesh, zoneName, zoneFaces[zoneI]);
1081  fset.write();
1082 
1083  fz[nValidFaceZones] = new faceZone
1084  (
1085  zoneName,
1086  zoneFaces[zoneI],
1087  boolList(zoneFaces[zoneI].size(), true),
1088  nValidFaceZones,
1089  mesh.faceZones()
1090  );
1091  nValidFaceZones++;
1092  }
1093  }
1094  }
1095 
1096  if (cz.size() || fz.size())
1097  {
1098  mesh.addZones(List<pointZone*>(0), fz, cz);
1099  }
1100 
1101  // Remove empty defaultFaces
1102  label defaultPatchID = mesh.boundaryMesh().findPatchID(defaultFacesName);
1103  if (mesh.boundaryMesh()[defaultPatchID].size() == 0)
1104  {
1105  List<polyPatch*> newPatchPtrList((mesh.boundaryMesh().size() - 1));
1106  label newPatchi = 0;
1107  forAll(mesh.boundaryMesh(), patchi)
1108  {
1109  if (patchi != defaultPatchID)
1110  {
1111  const polyPatch& patch = mesh.boundaryMesh()[patchi];
1112 
1113  newPatchPtrList[newPatchi] = patch.clone
1114  (
1115  mesh.boundaryMesh(),
1116  newPatchi,
1117  patch.size(),
1118  patch.start()
1119  ).ptr();
1120 
1121  newPatchi++;
1122  }
1123  }
1124  repatcher.changePatches(newPatchPtrList);
1125  }
1126 
1127  mesh.write();
1128 
1129  Info<< "End\n" << endl;
1130 
1131  return 0;
1132 }
1133 
1134 
1135 // ************************************************************************* //
ListType renumber(const labelUList &oldToNew, const ListType &)
Renumber the values (not the indices) of a list.
word defaultFacesType
Definition: readKivaGrid.H:461
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
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
static const cellModel * lookup(const word &)
Look up a model by name and return a pointer to the model or nullptr.
Definition: cellModeller.C:100
A list of face labels.
Definition: faceSet.H:48
An STL-conforming const_iterator.
Definition: HashTable.H:481
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
static iteratorEnd end()
iteratorEnd set to beyond the end of any HashTable
Definition: HashTable.H:112
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
Cell-face mesh analysis engine.
Definition: primitiveMesh.H:74
Foam::word regionName
An analytical geometric cellShape.
Definition: cellShape.H:69
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:163
engineTime & runTime
static word defaultRegion
Return the default region name.
Definition: polyMesh.H:309
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:256
bool optionFound(const word &opt) const
Return true if the named option is found.
Definition: argListI.H:108
static void noParallel()
Remove the parallel options.
Definition: argList.C:174
static SLList< string > validArgs
A list of valid (mandatory) arguments.
Definition: argList.H:153
const Cmpt & z() const
Definition: VectorI.H:87
label lineNumber() const
Return current stream line number.
Definition: IOstream.H:438
bool good() const
Return true if next operation might succeed.
Definition: IOstream.H:333
bool optionReadIfPresent(const word &opt, T &) const
Read a value from the named option if present.
Definition: argListI.H:198
label size() const
Return number of elements in table.
Definition: HashTableI.H:65
const Cmpt & y() const
Definition: VectorI.H:81
List< bool > boolList
Bool container classes.
Definition: boolList.H:50
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:142
virtual autoPtr< polyPatch > clone(const polyBoundaryMesh &bm) const
Construct and return a clone, resetting the boundary mesh.
Definition: polyPatch.H:222
A list of faces which address into the list of points.
const Map< label > & meshPointMap() const
Mesh point map. Given the global point index find its.
faceList faces() const
Faces of this cell.
Definition: cellShapeI.H:180
static const labelSphericalTensor labelI(1)
Identity labelTensor.
dynamicFvMesh & mesh
const cellShapeList & cells
A mesh which allows changes in the patch distribution of the boundary faces. The change in patching i...
const pointField & points
bool found(const Key &) const
Return true if hashedEntry is found in table.
Definition: HashTable.C:113
A class for handling words, derived from string.
Definition: word.H:59
const fileName & name() const
Return the name of the stream.
Definition: IFstream.H:116
void append(const T &)
Append an element at the end of the list.
Definition: ListI.H:177
const word & constant() const
Return constant name.
Definition: TimePaths.H:124
virtual void setTime(const Time &)
Reset the time and time-index to those of the given time.
Definition: Time.C:873
const Cmpt & x() const
Definition: VectorI.H:75
const labelListList & pointFaces() const
Return point-face addressing.
vector area(const pointField &) const
Return vector area.
Definition: face.C:552
static const char nl
Definition: Ostream.H:265
Input from file stream.
Definition: IFstream.H:81
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurrence of given element and return index,.
A subset of mesh cells.
Definition: cellZone.H:61
word defaultFacesName
Definition: readKivaGrid.H:460
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
messageStream Warning
An instant of time. Contains the time value and name.
Definition: instant.H:66
void setSize(const label)
Reset size of List.
Definition: List.C:281
label patchi
#define WarningInFunction
Report a warning using Foam::Warning.
#define FatalIOErrorInFunction(ios)
Report an error message using Foam::FatalIOError.
Definition: error.H:331
Input from memory buffer stream.
Definition: IStringStream.H:49
A collection of cell labels.
Definition: cellSet.H:48
virtual const faceList & faces() const =0
Return faces.
point centre(const pointField &) const
Centroid of the cell.
Definition: cellShapeI.H:259
label start() const
Return start label of this patch in the polyMesh face list.
Definition: polyPatch.H:303
Maps a geometry to a set of cell primitives, which enables geometric cell data to be calculated witho...
Definition: cellModel.H:64
messageStream Info
const labelListList & pointFaces() const
void resize(const label newSize)
Resize the hash table for efficiency.
Definition: HashTable.C:432
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
static void addBoolOption(const word &opt, const string &usage="")
Add to a bool option to validOptions with usage information.
Definition: argList.C:117
A subset of mesh faces organised as a primitive patch.
Definition: faceZone.H:64
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:66
Foam::argList args(argc, argv)
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:92
Namespace for OpenFOAM.
ISstream & getLine(string &)
Raw, low-level getline into a string function.
Definition: ISstreamI.H:77
IOerror FatalIOError