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