FreeFOAM The Cross-Platform CFD Toolkit
ideasUnvToFoam.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-2011 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  ideasUnvToFoam
26 
27 Description
28  I-Deas unv format mesh conversion.
29 
30  Uses either
31  - DOF sets (757) or
32  - face groups (2452(Cubit), 2467)
33  to do patching.
34  Works without but then puts all faces in defaultFaces patch.
35 
36 Usage
37 
38  - ideasUnvToFoam [OPTIONS] <.unv file>
39 
40  @param <.unv file> \n
41  @todo Detailed description of argument.
42 
43  @param -case <dir>\n
44  Case directory.
45 
46  @param -help \n
47  Display help message.
48 
49  @param -doc \n
50  Display Doxygen API documentation page for this application.
51 
52  @param -srcDoc \n
53  Display Doxygen source documentation page for this application.
54 
55 \*---------------------------------------------------------------------------*/
56 
57 #include <OpenFOAM/argList.H>
58 #include <OpenFOAM/polyMesh.H>
59 #include <OpenFOAM/Time.H>
60 #include <OpenFOAM/IFstream.H>
61 #include <OpenFOAM/cellModeller.H>
62 #include <meshTools/cellSet.H>
63 #include <meshTools/faceSet.H>
64 #include <OpenFOAM/DynamicList.H>
65 #include <triSurface/triSurface.H>
66 #include <OpenFOAM/SortableList.H>
67 
68 #include <cassert>
69 
70 using namespace Foam;
71 
72 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
73 
74 namespace Foam
75 {
76  template<>
77  inline unsigned Hash<face>::operator()(const face& t, unsigned seed) const
78  {
79  return Hasher(t.cdata(),t.size()*sizeof(label), seed);
80  }
81 
82  template<>
83  inline unsigned Hash<face>::operator()(const face& t) const
84  {
85  return Hash<face>::operator()(t, 0);
86  }
87 }
88 const string SEPARATOR(" -1");
89 
90 bool isSeparator(const string& line)
91 {
92  return line.substr(0, 6) == SEPARATOR;
93 }
94 
95 
96 // Reads past -1 and reads element type
97 label readTag(IFstream& is)
98 {
99  string tag;
100  do
101  {
102  if (!is.good())
103  {
104  return -1;
105  }
106 
107  string line;
108 
109  is.getLine(line);
110 
111  if (line.size() < 6)
112  {
113  return -1;
114  }
115 
116  tag = line.substr(0, 6);
117 
118  } while (tag == SEPARATOR);
119 
120  return readLabel(IStringStream(tag)());
121 }
122 
123 
124 // Reads and prints header
125 void readHeader(IFstream& is)
126 {
127  string line;
128 
129  while (is.good())
130  {
131  is.getLine(line);
132 
133  if (isSeparator(line))
134  {
135  break;
136  }
137  else
138  {
139  Sout<< line << endl;
140  }
141  }
142 }
143 
144 
145 // Skip
146 void skipSection(IFstream& is)
147 {
148  Sout<< "Skipping section at line " << is.lineNumber() << '.' << endl;
149 
150  string line;
151 
152  while (is.good())
153  {
154  is.getLine(line);
155 
156  if (isSeparator(line))
157  {
158  break;
159  }
160  }
161 }
162 
163 
164 scalar readUnvScalar(const string& unvString)
165 {
166  string s(unvString);
167 
168  s.replaceAll("d", "E");
169  s.replaceAll("D", "E");
170 
171  return readScalar(IStringStream(s)());
172 }
173 
174 
175 // Reads unit section
176 void readUnits
177 (
178  IFstream& is,
179  scalar& lengthScale,
180  scalar& forceScale,
181  scalar& tempScale,
182  scalar& tempOffset
183 )
184 {
185  Sout<< "Starting reading units at line " << is.lineNumber() << '.' << endl;
186 
187  string line;
188  is.getLine(line);
189 
190  label l = readLabel(IStringStream(line.substr(0, 10))());
191  Sout<< "l:" << l << endl;
192 
193  string units(line.substr(10, 20));
194  Sout<< "units:" << units << endl;
195 
196  label unitType = readLabel(IStringStream(line.substr(30, 10))());
197  Sout<< "unitType:" << unitType << endl;
198 
199  // Read lengthscales
200  is.getLine(line);
201 
202  lengthScale = readUnvScalar(line.substr(0, 25));
203  forceScale = readUnvScalar(line.substr(25, 25));
204  tempScale = readUnvScalar(line.substr(50, 25));
205 
206  is.getLine(line);
207  tempOffset = readUnvScalar(line.substr(0, 25));
208 
209  Sout<< "Unit factors:" << nl
210  << " Length scale : " << lengthScale << nl
211  << " Force scale : " << forceScale << nl
212  << " Temperature scale : " << tempScale << nl
213  << " Temperature offset : " << tempOffset << nl
214  << endl;
215 }
216 
217 
218 // Reads points section. Read region as well?
219 void readPoints
220 (
221  IFstream& is,
222  DynamicList<point>& points, // coordinates
223  DynamicList<label>& unvPointID // unv index
224 )
225 {
226  Sout<< "Starting reading points at line " << is.lineNumber() << '.' << endl;
227 
228  static bool hasWarned = false;
229 
230  while (true)
231  {
232  string line;
233  is.getLine(line);
234 
235  label pointI = readLabel(IStringStream(line.substr(0, 10))());
236 
237  if (pointI == -1)
238  {
239  break;
240  }
241  else if (pointI != points.size()+1 && !hasWarned)
242  {
243  hasWarned = true;
244 
246  (
247  "readPoints(IFstream&, label&, DynamicList<point>"
248  ", DynamicList<label>&)",
249  is
250  ) << "Points not in order starting at point " << pointI
251  //<< " at line " << is.lineNumber()
252  //<< abort(FatalError);
253  << endl;
254  }
255 
256  point pt;
257  is.getLine(line);
258  pt[0] = readUnvScalar(line.substr(0, 25));
259  pt[1] = readUnvScalar(line.substr(25, 25));
260  pt[2] = readUnvScalar(line.substr(50, 25));
261 
262  unvPointID.append(pointI);
263  points.append(pt);
264  }
265 
266  points.shrink();
267  unvPointID.shrink();
268 
269  Sout<< "Read " << points.size() << " points." << endl;
270 }
271 
272 void addAndExtend
273 (
274  DynamicList<label>& indizes,
275  label cellI,
276  label val
277 )
278 {
279  if (indizes.size() < (cellI+1))
280  {
281  indizes.setSize(cellI+1,-1);
282  }
283  indizes[cellI] = val;
284 }
285 
286 // Reads cells section. Read region as well? Not handled yet but should just
287 // be a matter of reading corresponding to boundaryFaces the correct property
288 // and sorting it later on.
289 void readCells
290 (
291  IFstream& is,
292  DynamicList<cellShape>& cellVerts,
293  DynamicList<label>& cellMaterial,
294  DynamicList<label>& boundaryFaceIndices,
295  DynamicList<face>& boundaryFaces,
296  DynamicList<label>& cellCorrespondence,
297  DynamicList<label>& unvPointID // unv index
298 )
299 {
300  Sout<< "Starting reading cells at line " << is.lineNumber() << '.' << endl;
301 
302  // Invert point numbering.
303  label maxUnvPoint = 0;
304  forAll(unvPointID, pointI)
305  {
306  maxUnvPoint = max(maxUnvPoint, unvPointID[pointI]);
307  }
308  labelList unvToFoam(invert(maxUnvPoint+1, unvPointID));
309 
310 
311  const cellModel& hex = *(cellModeller::lookup("hex"));
312  const cellModel& prism = *(cellModeller::lookup("prism"));
313  const cellModel& tet = *(cellModeller::lookup("tet"));
314 
315  labelHashSet skippedElements;
316 
317  labelHashSet foundFeType;
318 
319  while (true)
320  {
321  string line;
322  is.getLine(line);
323 
324  if (isSeparator(line))
325  {
326  break;
327  }
328 
329  IStringStream lineStr(line);
330  label cellI, feID, physProp, matProp, colour, nNodes;
331  lineStr >> cellI >> feID >> physProp >> matProp >> colour >> nNodes;
332 
333  if (foundFeType.insert(feID))
334  {
335  Info<< "First occurrence of element type " << feID
336  << " for cell " << cellI << " at line "
337  << is.lineNumber() << endl;
338  }
339 
340  if (feID == 11)
341  {
342  // Rod. Skip.
343  is.getLine(line);
344  is.getLine(line);
345  }
346  else if (feID == 171)
347  {
348  // Rod. Skip.
349  is.getLine(line);
350  }
351  else if (feID == 41 || feID == 91)
352  {
353  // Triangle. Save - used for patching later on.
354  is.getLine(line);
355 
356  face cVerts(3);
357  IStringStream lineStr(line);
358  lineStr >> cVerts[0] >> cVerts[1] >> cVerts[2];
359  boundaryFaces.append(cVerts);
360  boundaryFaceIndices.append(cellI);
361  }
362  else if (feID == 44 || feID == 94)
363  {
364  // Quad. Save - used for patching later on.
365  is.getLine(line);
366 
367  face cVerts(4);
368  IStringStream lineStr(line);
369  lineStr >> cVerts[0] >> cVerts[1] >> cVerts[2] >> cVerts[3];
370  boundaryFaces.append(cVerts);
371  boundaryFaceIndices.append(cellI);
372  }
373  else if (feID == 111)
374  {
375  // Tet.
376  is.getLine(line);
377 
378  labelList cVerts(4);
379  IStringStream lineStr(line);
380  lineStr >> cVerts[0] >> cVerts[1] >> cVerts[2] >> cVerts[3];
381 
382  cellVerts.append(cellShape(tet, cVerts, true));
383  cellMaterial.append(physProp);
384  addAndExtend(cellCorrespondence,cellI,cellMaterial.size()-1);
385 
386  if (cellVerts[cellVerts.size()-1].size() != cVerts.size())
387  {
388  Pout<< "Line:" << is.lineNumber()
389  << " element:" << cellI
390  << " type:" << feID
391  << " collapsed from " << cVerts << nl
392  << " to:" << cellVerts[cellVerts.size()-1]
393  << endl;
394  }
395  }
396  else if (feID == 112)
397  {
398  // Wedge.
399  is.getLine(line);
400 
401  labelList cVerts(6);
402  IStringStream lineStr(line);
403  lineStr >> cVerts[0] >> cVerts[1] >> cVerts[2] >> cVerts[3]
404  >> cVerts[4] >> cVerts[5];
405 
406  cellVerts.append(cellShape(prism, cVerts, true));
407  cellMaterial.append(physProp);
408  addAndExtend(cellCorrespondence,cellI,cellMaterial.size()-1);
409 
410  if (cellVerts[cellVerts.size()-1].size() != cVerts.size())
411  {
412  Pout<< "Line:" << is.lineNumber()
413  << " element:" << cellI
414  << " type:" << feID
415  << " collapsed from " << cVerts << nl
416  << " to:" << cellVerts[cellVerts.size()-1]
417  << endl;
418  }
419  }
420  else if (feID == 115)
421  {
422  // Hex.
423  is.getLine(line);
424 
425  labelList cVerts(8);
426  IStringStream lineStr(line);
427  lineStr
428  >> cVerts[0] >> cVerts[1] >> cVerts[2] >> cVerts[3]
429  >> cVerts[4] >> cVerts[5] >> cVerts[6] >> cVerts[7];
430 
431  cellVerts.append(cellShape(hex, cVerts, true));
432  cellMaterial.append(physProp);
433  addAndExtend(cellCorrespondence,cellI,cellMaterial.size()-1);
434 
435  if (cellVerts[cellVerts.size()-1].size() != cVerts.size())
436  {
437  Pout<< "Line:" << is.lineNumber()
438  << " element:" << cellI
439  << " type:" << feID
440  << " collapsed from " << cVerts << nl
441  << " to:" << cellVerts[cellVerts.size()-1]
442  << endl;
443  }
444  }
445  else if (feID == 118)
446  {
447  // Parabolic Tet
448 
449  is.getLine(line);
450 
451  labelList cVerts(4);
452  label dummy;
453  {
454  IStringStream lineStr(line);
455  lineStr
456  >> cVerts[0] >> dummy >> cVerts[1] >> dummy
457  >> cVerts[2];
458  }
459  is.getLine(line);
460  {
461  IStringStream lineStr(line);
462  lineStr >> dummy>> cVerts[3];
463  }
464 
465  cellVerts.append(cellShape(tet, cVerts, true));
466  cellMaterial.append(physProp);
467  addAndExtend(cellCorrespondence,cellI,cellMaterial.size()-1);
468 
469  if (cellVerts[cellVerts.size()-1].size() != cVerts.size())
470  {
471  Pout<< "Line:" << is.lineNumber()
472  << " element:" << cellI
473  << " type:" << feID
474  << " collapsed from " << cVerts << nl
475  << " to:" << cellVerts[cellVerts.size()-1]
476  << endl;
477  }
478  }
479  else
480  {
481  if (skippedElements.insert(feID))
482  {
483  IOWarningIn("readCells(IFstream&, label&)", is)
484  << "Cell type " << feID << " not supported" << endl;
485  }
486  is.getLine(line); //Do nothing
487  }
488  }
489 
490  cellVerts.shrink();
491  cellMaterial.shrink();
492  boundaryFaces.shrink();
493  boundaryFaceIndices.shrink();
494  cellCorrespondence.shrink();
495 
496  Sout<< "Read " << cellVerts.size() << " cells"
497  << " and " << boundaryFaces.size() << " boundary faces." << endl;
498 }
499 
500 
501 void readSets
502 (
503  IFstream& is,
505  DynamicList<labelList>& patchFaceIndices
506 )
507 {
508  Sout<< "Starting reading patches at line " << is.lineNumber() << '.'
509  << endl;
510 
511  while(true)
512  {
513  string line;
514  is.getLine(line);
515 
516  if (isSeparator(line))
517  {
518  break;
519  }
520 
521  IStringStream lineStr(line);
522  label group, constraintSet, restraintSet, loadSet, dofSet,
523  tempSet, contactSet, nFaces;
524  lineStr
525  >> group >> constraintSet >> restraintSet >> loadSet
526  >> dofSet >> tempSet >> contactSet >> nFaces;
527 
528  is.getLine(line);
529  word groupName = string::validate<word>(line);
530 
531  Info<< "For group " << group
532  << " named " << groupName
533  << " trying to read " << nFaces << " patch face indices."
534  << endl;
535 
536  labelList groupIndices(nFaces);
537  label groupType = -1;
538  nFaces = 0;
539 
540  while (nFaces < groupIndices.size())
541  {
542  is.getLine(line);
543  IStringStream lineStr(line);
544 
545  // Read one (for last face) or two entries from line.
546  label nRead = 2;
547  if (nFaces == groupIndices.size()-1)
548  {
549  nRead = 1;
550  }
551 
552  for (label i = 0; i < nRead; i++)
553  {
554  label tag, nodeLeaf, component;
555 
556  lineStr >> groupType >> tag >> nodeLeaf >> component;
557 
558  groupIndices[nFaces++] = tag;
559  }
560  }
561 
562 
563  // Store
564  if (groupType == 8)
565  {
566  patchNames.append(groupName);
567  patchFaceIndices.append(groupIndices);
568  }
569  else
570  {
571  IOWarningIn("readSets(..)", is)
572  << "When reading patches expect entity type code 8"
573  << nl << " Skipping group code " << groupType
574  << endl;
575  }
576  }
577 
578  patchNames.shrink();
579  patchFaceIndices.shrink();
580 }
581 
582 
583 
584 // Read dof set (757)
585 void readDOFS
586 (
587  IFstream& is,
588  DynamicList<word>& patchNames,
589  DynamicList<labelList>& dofVertices
590 )
591 {
592  Sout<< "Starting reading constraints at line " << is.lineNumber() << '.'
593  << endl;
594 
595  string line;
596  is.getLine(line);
597  label group;
598  {
599  IStringStream lineStr(line);
600  lineStr >> group;
601  }
602 
603  is.getLine(line);
604  {
605  IStringStream lineStr(line);
606  patchNames.append(lineStr);
607  }
608 
609  Info<< "For DOF set " << group
610  << " named " << patchNames[patchNames.size()-1]
611  << " trying to read vertex indices."
612  << endl;
613 
614  DynamicList<label> vertices;
615  while (true)
616  {
617  string line;
618  is.getLine(line);
619 
620  if (isSeparator(line))
621  {
622  break;
623  }
624 
625  IStringStream lineStr(line);
626  label pointI;
627  lineStr >> pointI;
628 
629  vertices.append(pointI);
630  }
631 
632  Info<< "For DOF set " << group
633  << " named " << patchNames[patchNames.size()-1]
634  << " read " << vertices.size() << " vertex indices." << endl;
635 
636  dofVertices.append(vertices.shrink());
637 
638  patchNames.shrink();
639  dofVertices.shrink();
640 }
641 
642 
643 // Returns -1 or group that all of the vertices of f are in,
644 label findPatch(const List<labelHashSet>& dofGroups, const face& f)
645 {
646  forAll(dofGroups, patchI)
647  {
648  if (dofGroups[patchI].found(f[0]))
649  {
650  bool allInGroup = true;
651 
652  // Check rest of face
653  for (label fp = 1; fp < f.size(); fp++)
654  {
655  if (!dofGroups[patchI].found(f[fp]))
656  {
657  allInGroup = false;
658  break;
659  }
660  }
661 
662  if (allInGroup)
663  {
664  return patchI;
665  }
666  }
667  }
668  return -1;
669 }
670 
671 // Additional output while we're unsure
672 const bool verbose=false;
673 
674 // Main program:
675 
676 int main(int argc, char *argv[])
677 {
679  argList::validArgs.append(".unv file");
680  argList::validOptions.insert("dump", "");
681 
682 # include <OpenFOAM/setRootCase.H>
683 # include <OpenFOAM/createTime.H>
684 
685  fileName ideasName(args.additionalArgs()[0]);
686 
687  IFstream inFile(ideasName.c_str());
688 
689  if (!inFile.good())
690  {
692  << "Cannot open file " << ideasName
693  << exit(FatalError);
694  }
695 
696 
697  // Unit scale factors
698  scalar lengthScale = 1;
699  scalar forceScale = 1;
700  scalar tempScale = 1;
701  scalar tempOffset = 0;
702 
703  // Points
705  // Original unv point label
706  DynamicList<label> unvPointID;
707 
708  // Cells
709  DynamicList<cellShape> cellVerts;
710  DynamicList<label> cellMat;
711  DynamicList<label> cellCorrespondence;
712 
713  // Boundary faces
714  DynamicList<label> boundaryFaceIndices;
715  DynamicList<face> boundaryFaces;
716 
717  // Patch names and patchFace indices.
719  DynamicList<labelList> patchFaceIndices;
720  DynamicList<labelList> dofVertIndices;
721 
722 
723  while (inFile.good())
724  {
725  label tag = readTag(inFile);
726 
727  if (tag == -1)
728  {
729  break;
730  }
731 
732  Sout<< "Processing tag:" << tag << endl;
733 
734  switch (tag)
735  {
736  case 151:
737  readHeader(inFile);
738  break;
739 
740  case 164:
741  readUnits
742  (
743  inFile,
744  lengthScale,
745  forceScale,
746  tempScale,
747  tempOffset
748  );
749  break;
750 
751  case 2411:
752  readPoints(inFile, points, unvPointID);
753  break;
754 
755  case 2412:
756  readCells
757  (
758  inFile,
759  cellVerts,
760  cellMat,
761  boundaryFaceIndices,
762  boundaryFaces,
763  cellCorrespondence,
764  unvPointID
765  );
766  break;
767 
768  case 2452:
769  case 2467:
770  readSets
771  (
772  inFile,
773  patchNames,
774  patchFaceIndices
775  );
776  break;
777 
778  case 757:
779  readDOFS
780  (
781  inFile,
782  patchNames,
783  dofVertIndices
784  );
785  break;
786 
787  default:
788  Sout<< "Skipping tag " << tag << " on line "
789  << inFile.lineNumber() << endl;
790  skipSection(inFile);
791  break;
792  }
793  Sout<< endl;
794  }
795 
796 
797  // Invert point numbering.
798  label maxUnvPoint = 0;
799  forAll(unvPointID, pointI)
800  {
801  maxUnvPoint = max(maxUnvPoint, unvPointID[pointI]);
802  }
803  labelList unvToFoam(invert(maxUnvPoint+1, unvPointID));
804 
805 
806  // Renumber vertex numbers on cells
807 
808  forAll(cellVerts, cellI)
809  {
810  labelList foamVerts
811  (
812  renumber
813  (
814  unvToFoam,
815  static_cast<labelList&>(cellVerts[cellI])
816  )
817  );
818 
819  if (findIndex(foamVerts, -1) != -1)
820  {
822  << "Cell " << cellI
823  << " unv vertices " << cellVerts[cellI]
824  << " has some undefined vertices " << foamVerts
825  << abort(FatalError);
826  }
827 
828  // Bit nasty: replace vertex list.
829  cellVerts[cellI].transfer(foamVerts);
830  }
831 
832  // Renumber vertex numbers on boundaryFaces
833 
834  forAll(boundaryFaces, bFaceI)
835  {
836  labelList foamVerts(renumber(unvToFoam, boundaryFaces[bFaceI]));
837 
838  if (findIndex(foamVerts, -1) != -1)
839  {
841  << "Boundary face " << bFaceI
842  << " unv vertices " << boundaryFaces[bFaceI]
843  << " has some undefined vertices " << foamVerts
844  << abort(FatalError);
845  }
846 
847  // Bit nasty: replace vertex list.
848  boundaryFaces[bFaceI].transfer(foamVerts);
849  }
850 
851 
852 
853  // Patchify = create patchFaceVerts for use in cellShape construction
854  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
855 
856  // Works in one of two modes. Either has read boundaryFaces which
857  // just need to be sorted according to patch. Or has read point constraint
858  // sets (dofVertIndices).
859 
860  List<faceList> patchFaceVerts;
861 
862  labelList nrFaceCells(boundaryFaces.size(),0);
864 
865  {
866  HashTable<label, face, Hash<face> > faceToFaceID(boundaryFaces.size());
867  forAll(boundaryFaces,faceI)
868  {
869  SortableList<label> foo(boundaryFaces[faceI]);
870  face theFace(foo);
871  faceToFaceID.insert(theFace,faceI);
872  }
873 
874  forAll(cellVerts,cellI)
875  {
876  faceList faces = cellVerts[cellI].faces();
877  forAll(faces,i)
878  {
879  SortableList<label> foo(faces[i]);
880  face theFace(foo);
881  if (faceToFaceID.found(theFace))
882  {
883  label faceI = faceToFaceID[theFace];
884  if (nrFaceCells[faceI] < 2)
885  {
886  faceToCell[nrFaceCells[faceI]].insert(faceI,cellI);
887  }
888  nrFaceCells[faceI]++;
889  }
890  }
891  }
892 
893  label cnt = 0;
894  forAll(nrFaceCells, faceI)
895  {
896  assert(nrFaceCells[faceI] == 1 || nrFaceCells[faceI] == 2);
897  if (nrFaceCells[faceI]>1)
898  {
899  cnt++;
900  }
901  }
902 
903  if (cnt>0)
904  {
905  Info << "Of " << boundaryFaces.size() << " so-called"
906  << " boundary faces " << cnt << " belong to two cells "
907  << "and are therefore internal" << endl;
908  }
909  }
910 
911  HashTable<labelList,word> cellZones;
912  HashTable<labelList,word> faceZones;
913  List<bool> isAPatch(patchNames.size(),true);
914 
915  if (dofVertIndices.size())
916  {
917  // Use the vertex constraints to patch. Is of course bit dodgy since
918  // face goes on patch if all its vertices are on a constraint.
919  // Note: very inefficient since goes through all faces (including
920  // internal ones) twice. Should do a construct without patches
921  // and then repatchify.
922 
923  Info<< "Using " << dofVertIndices.size()
924  << " DOF sets to detect boundary faces."<< endl;
925 
926  // Renumber vertex numbers on constraints
927  forAll(dofVertIndices, patchI)
928  {
929  inplaceRenumber(unvToFoam, dofVertIndices[patchI]);
930  }
931 
932 
933  // Build labelHashSet of points per dofGroup/patch
934  List<labelHashSet> dofGroups(dofVertIndices.size());
935 
936  forAll(dofVertIndices, patchI)
937  {
938  const labelList& foamVerts = dofVertIndices[patchI];
939 
940  forAll(foamVerts, i)
941  {
942  dofGroups[patchI].insert(foamVerts[i]);
943  }
944  }
945 
946  List<DynamicList<face> > dynPatchFaces(dofVertIndices.size());
947 
948  forAll(cellVerts, cellI)
949  {
950  const cellShape& shape = cellVerts[cellI];
951 
952  const faceList shapeFaces(shape.faces());
953 
954  forAll(shapeFaces, i)
955  {
956  label patchI = findPatch(dofGroups, shapeFaces[i]);
957 
958  if (patchI != -1)
959  {
960  dynPatchFaces[patchI].append(shapeFaces[i]);
961  }
962  }
963  }
964 
965  // Transfer
966  patchFaceVerts.setSize(dynPatchFaces.size());
967 
968  forAll(dynPatchFaces, patchI)
969  {
970  patchFaceVerts[patchI].transfer(dynPatchFaces[patchI]);
971  }
972  }
973  else
974  {
975  // Use the boundary faces.
976 
977  // Construct the patch faces by sorting the boundaryFaces according to
978  // patch.
979  patchFaceVerts.setSize(patchFaceIndices.size());
980 
981  Info<< "Sorting boundary faces according to group (patch)" << endl;
982 
983  // make sure that no face is used twice on the boundary
984  // (possible for boundary-only faceZones)
985  labelHashSet alreadyOnBoundary;
986 
987  // Construct map from boundaryFaceIndices
988  Map<label> boundaryFaceToIndex(boundaryFaceIndices.size());
989 
990  forAll(boundaryFaceIndices, i)
991  {
992  boundaryFaceToIndex.insert(boundaryFaceIndices[i], i);
993  }
994 
995  forAll(patchFaceVerts, patchI)
996  {
997  Info << patchI << ": " << patchNames[patchI] << " is " << flush;
998 
999  faceList& patchFaces = patchFaceVerts[patchI];
1000  const labelList& faceIndices = patchFaceIndices[patchI];
1001 
1002  patchFaces.setSize(faceIndices.size());
1003 
1004  bool duplicateFaces = false;
1005 
1006  label cnt = 0;
1007  forAll(patchFaces, i)
1008  {
1009  if (boundaryFaceToIndex.found(faceIndices[i]))
1010  {
1011  label bFaceI = boundaryFaceToIndex[faceIndices[i]];
1012  if (nrFaceCells[bFaceI] == 1)
1013  {
1014  patchFaces[cnt] = boundaryFaces[bFaceI];
1015  cnt++;
1016  if (alreadyOnBoundary.found(bFaceI))
1017  {
1018  duplicateFaces = true;
1019  }
1020  }
1021  }
1022  }
1023 
1024  if (cnt!=patchFaces.size() || duplicateFaces)
1025  {
1026  isAPatch[patchI] = false;
1027 
1028  if (verbose)
1029  {
1030  if (cnt != patchFaces.size())
1031  {
1033  << "For patch " << patchI << " there were "
1034  << patchFaces.size()-cnt
1035  << " faces not used because they seem"
1036  << " to be internal. "
1037  << "This seems to be a face or a cell-zone"
1038  << endl;
1039  }
1040  else
1041  {
1043  << "Patch "
1044  << patchI << " has faces that are already "
1045  << " in use on other boundary-patches,"
1046  << " Assuming faceZoneset." << endl;
1047  }
1048  }
1049 
1050  patchFaces.setSize(0); // Assume that this is no patch at all
1051 
1052  if (cellCorrespondence[faceIndices[0]] >= 0)
1053  {
1054  Info << "cellZone" << endl;
1055  labelList theCells(faceIndices.size());
1056  forAll(faceIndices, i)
1057  {
1058  if (cellCorrespondence[faceIndices[0]] < 0)
1059  {
1061  << "The face index " << faceIndices[i]
1062  << " was not found amongst the cells."
1063  << " This kills the theory that "
1064  << patchNames[patchI] << " is a cell zone"
1065  << endl
1066  << abort(FatalError);
1067  }
1068  theCells[i] = cellCorrespondence[faceIndices[i]];
1069  }
1070  cellZones.insert(patchNames[patchI], theCells);
1071  }
1072  else
1073  {
1074  Info << "faceZone" << endl;
1075  labelList theFaces(faceIndices.size());
1076  forAll(faceIndices, i)
1077  {
1078  theFaces[i] = boundaryFaceToIndex[faceIndices[i]];
1079  }
1080  faceZones.insert(patchNames[patchI],theFaces);
1081  }
1082  }
1083  else
1084  {
1085  Info << "patch" << endl;
1086 
1087  forAll(patchFaces, i)
1088  {
1089  label bFaceI = boundaryFaceToIndex[faceIndices[i]];
1090  alreadyOnBoundary.insert(bFaceI);
1091  }
1092  }
1093  }
1094  }
1095 
1096  pointField polyPoints;
1097  polyPoints.transfer(points);
1098 
1099  // Length scaling factor
1100  polyPoints /= lengthScale;
1101 
1102 
1103  // For debugging: dump boundary faces as triSurface
1104  if (args.optionFound("dump"))
1105  {
1106  DynamicList<labelledTri> triangles(boundaryFaces.size());
1107 
1108  forAll(boundaryFaces, i)
1109  {
1110  const face& f = boundaryFaces[i];
1111 
1112  faceList triFaces(f.nTriangles(polyPoints));
1113  label nTri = 0;
1114  f.triangles(polyPoints, nTri, triFaces);
1115 
1116  forAll(triFaces, triFaceI)
1117  {
1118  const face& f = triFaces[triFaceI];
1119  triangles.append(labelledTri(f[0], f[1], f[2], 0));
1120  }
1121  }
1122 
1123  // Create globally numbered tri surface
1124  triSurface rawSurface(triangles.shrink(), polyPoints);
1125 
1126  // Create locally numbered tri surface
1127  triSurface surface
1128  (
1129  rawSurface.localFaces(),
1130  rawSurface.localPoints()
1131  );
1132 
1133  Info<< "Writing boundary faces to STL file boundaryFaces.stl"
1134  << nl << endl;
1135 
1136  surface.write(runTime.path()/"boundaryFaces.stl");
1137  }
1138 
1139 
1140  Info<< "\nConstructing mesh with non-default patches of size:" << nl;
1141  DynamicList<word> usedPatchNames;
1142  DynamicList<faceList> usedPatchFaceVerts;
1143 
1144  forAll(patchNames, patchI)
1145  {
1146  if (isAPatch[patchI]) {
1147  Info<< " " << patchNames[patchI] << '\t'
1148  << patchFaceVerts[patchI].size() << nl;
1149  usedPatchNames.append(patchNames[patchI]);
1150  usedPatchFaceVerts.append(patchFaceVerts[patchI]);
1151  }
1152  }
1153  usedPatchNames.shrink();
1154  usedPatchFaceVerts.shrink();
1155 
1156  Info<< endl;
1157 
1158 
1159 
1160  // Construct mesh
1161  polyMesh mesh
1162  (
1163  IOobject
1164  (
1166  runTime.constant(),
1167  runTime
1168  ),
1169  xferMove(polyPoints),
1170  cellVerts,
1171  usedPatchFaceVerts, // boundaryFaces,
1172  usedPatchNames, // boundaryPatchNames,
1173  wordList(patchNames.size(), polyPatch::typeName), // boundaryPatchTypes,
1174  "defaultFaces", // defaultFacesName
1175  polyPatch::typeName, // defaultFacesType,
1176  wordList(0) // boundaryPatchPhysicalTypes
1177  );
1178 
1179 
1180  if (faceZones.size()>0 || cellZones.size()>0)
1181  {
1182  Info << "Adding cell and face zones" << endl;
1183 
1185  List<faceZone*> fZones(faceZones.size());
1186  List<cellZone*> cZones(cellZones.size());
1187 
1188  if (cellZones.size()>0)
1189  {
1190  forAll(cellZones.toc(),cnt)
1191  {
1192  word name = cellZones.toc()[cnt];
1193  Info<< " Cell Zone " << name << " " << tab
1194  << cellZones[name].size() << endl;
1195 
1196  cZones[cnt] = new cellZone
1197  (
1198  name,
1199  cellZones[name],
1200  cnt,
1201  mesh.cellZones()
1202  );
1203  }
1204  }
1205  if (faceZones.size() > 0)
1206  {
1207  const labelList& own = mesh.faceOwner();
1208  const labelList& nei = mesh.faceNeighbour();
1209  const pointField& centers = mesh.faceCentres();
1210  const pointField& points = mesh.points();
1211 
1212  forAll(faceZones.toc(), cnt)
1213  {
1214  word name = faceZones.toc()[cnt];
1215  const labelList& oldIndizes = faceZones[name];
1216  labelList indizes(oldIndizes.size());
1217 
1218  Info<< " Face Zone " << name << " " << tab
1219  << oldIndizes.size() << endl;
1220 
1221  forAll(indizes, i)
1222  {
1223  const label old = oldIndizes[i];
1224  label noveau = -1;
1225  label c1 = -1, c2 = -1;
1226  if (faceToCell[0].found(old))
1227  {
1228  c1 = faceToCell[0][old];
1229  }
1230  if (faceToCell[1].found(old))
1231  {
1232  c2 = faceToCell[1][old];
1233  }
1234  if (c1<c2)
1235  {
1236  label tmp = c1;
1237  c1 = c2;
1238  c2 = tmp;
1239  }
1240  if (c2 == -1)
1241  {
1242  // Boundary face is part of the faceZone
1243  forAll(own, j)
1244  {
1245  if (own[j] == c1)
1246  {
1247  const face& f = boundaryFaces[old];
1248  if (mag(centers[j]-f.centre(points))<SMALL)
1249  {
1250  noveau = j;
1251  break;
1252  }
1253  }
1254  }
1255  }
1256  else
1257  {
1258  forAll(nei, j)
1259  {
1260  if
1261  (
1262  (c1 == own[j] && c2 == nei[j])
1263  || (c2 == own[j] && c1 == nei[j])
1264  )
1265  {
1266  noveau = j;
1267  break;
1268  }
1269  }
1270  }
1271  assert(noveau>-1);
1272  indizes[i] = noveau;
1273  }
1274  fZones[cnt] = new faceZone
1275  (
1276  faceZones.toc()[cnt],
1277  indizes,
1278  boolList(indizes.size(),false),
1279  cnt,
1280  mesh.faceZones()
1281  );
1282  }
1283  }
1284  mesh.addZones(pZones, fZones, cZones);
1285 
1286  Info << endl;
1287  }
1288 
1289  mesh.write();
1290 
1291  Info<< "End\n" << endl;
1292 
1293  return 0;
1294 }
1295 
1296 
1297 // ************************ vim: set sw=4 sts=4 et: ************************ //