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) 2011-2016 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  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 \*---------------------------------------------------------------------------*/
37 
38 #include "argList.H"
39 #include "polyMesh.H"
40 #include "Time.H"
41 #include "IFstream.H"
42 #include "cellModeller.H"
43 #include "cellSet.H"
44 #include "faceSet.H"
45 #include "DynamicList.H"
46 
47 #include <cassert>
48 #include "MeshedSurfaces.H"
49 
50 using namespace Foam;
51 
52 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
53 
54 namespace Foam
55 {
56  template<>
57  inline unsigned Hash<face>::operator()(const face& t, unsigned seed) const
58  {
59  return Hasher(t.cdata(),t.size()*sizeof(label), seed);
60  }
61 
62  template<>
63  inline unsigned Hash<face>::operator()(const face& t) const
64  {
65  return Hash<face>::operator()(t, 0);
66  }
67 }
68 const string SEPARATOR(" -1");
69 
70 bool isSeparator(const string& line)
71 {
72  return line.substr(0, 6) == SEPARATOR;
73 }
74 
75 
76 // Reads past -1 and reads element type
77 label readTag(IFstream& is)
78 {
79  string tag;
80  do
81  {
82  if (!is.good())
83  {
84  return -1;
85  }
86 
87  string line;
88 
89  is.getLine(line);
90 
91  if (line.size() < 6)
92  {
93  return -1;
94  }
95 
96  tag = line.substr(0, 6);
97 
98  } while (tag == SEPARATOR);
99 
100  return readLabel(IStringStream(tag)());
101 }
102 
103 
104 // Reads and prints header
105 void readHeader(IFstream& is)
106 {
107  string line;
108 
109  while (is.good())
110  {
111  is.getLine(line);
112 
113  if (isSeparator(line))
114  {
115  break;
116  }
117  else
118  {
119  Info<< line << endl;
120  }
121  }
122 }
123 
124 
125 // Skip
126 void skipSection(IFstream& is)
127 {
128  Info<< "Skipping section at line " << is.lineNumber() << '.' << endl;
129 
130  string line;
131 
132  while (is.good())
133  {
134  is.getLine(line);
135 
136  if (isSeparator(line))
137  {
138  break;
139  }
140  }
141 }
142 
143 
144 scalar readUnvScalar(const string& unvString)
145 {
146  string s(unvString);
147 
148  s.replaceAll("d", "E");
149  s.replaceAll("D", "E");
150 
151  return readScalar(IStringStream(s)());
152 }
153 
154 
155 // Reads unit section
156 void readUnits
157 (
158  IFstream& is,
159  scalar& lengthScale,
160  scalar& forceScale,
161  scalar& tempScale,
162  scalar& tempOffset
163 )
164 {
165  Info<< "Starting reading units at line " << is.lineNumber() << '.' << endl;
166 
167  string line;
168  is.getLine(line);
169 
170  label l = readLabel(IStringStream(line.substr(0, 10))());
171  Info<< "l:" << l << endl;
172 
173  string units(line.substr(10, 20));
174  Info<< "units:" << units << endl;
175 
176  label unitType = readLabel(IStringStream(line.substr(30, 10))());
177  Info<< "unitType:" << unitType << endl;
178 
179  // Read lengthscales
180  is.getLine(line);
181 
182  lengthScale = readUnvScalar(line.substr(0, 25));
183  forceScale = readUnvScalar(line.substr(25, 25));
184  tempScale = readUnvScalar(line.substr(50, 25));
185 
186  is.getLine(line);
187  tempOffset = readUnvScalar(line.substr(0, 25));
188 
189  Info<< "Unit factors:" << nl
190  << " Length scale : " << lengthScale << nl
191  << " Force scale : " << forceScale << nl
192  << " Temperature scale : " << tempScale << nl
193  << " Temperature offset : " << tempOffset << nl
194  << endl;
195 }
196 
197 
198 // Reads points section. Read region as well?
199 void readPoints
200 (
201  IFstream& is,
202  DynamicList<point>& points, // coordinates
203  DynamicList<label>& unvPointID // unv index
204 )
205 {
206  Info<< "Starting reading points at line " << is.lineNumber() << '.' << endl;
207 
208  static bool hasWarned = false;
209 
210  while (true)
211  {
212  string line;
213  is.getLine(line);
214 
215  label pointi = readLabel(IStringStream(line.substr(0, 10))());
216 
217  if (pointi == -1)
218  {
219  break;
220  }
221  else if (pointi != points.size()+1 && !hasWarned)
222  {
223  hasWarned = true;
224 
226  (
227  is
228  ) << "Points not in order starting at point " << pointi
229  //<< " at line " << is.lineNumber()
230  //<< abort(FatalError);
231  << endl;
232  }
233 
234  point pt;
235  is.getLine(line);
236  pt[0] = readUnvScalar(line.substr(0, 25));
237  pt[1] = readUnvScalar(line.substr(25, 25));
238  pt[2] = readUnvScalar(line.substr(50, 25));
239 
240  unvPointID.append(pointi);
241  points.append(pt);
242  }
243 
244  points.shrink();
245  unvPointID.shrink();
246 
247  Info<< "Read " << points.size() << " points." << endl;
248 }
249 
250 void addAndExtend
251 (
252  DynamicList<label>& indizes,
253  label celli,
254  label val
255 )
256 {
257  if (indizes.size() < (celli+1))
258  {
259  indizes.setSize(celli+1,-1);
260  }
261  indizes[celli] = val;
262 }
263 
264 // Reads cells section. Read region as well? Not handled yet but should just
265 // be a matter of reading corresponding to boundaryFaces the correct property
266 // and sorting it later on.
267 void readCells
268 (
269  IFstream& is,
270  DynamicList<cellShape>& cellVerts,
271  DynamicList<label>& cellMaterial,
272  DynamicList<label>& boundaryFaceIndices,
273  DynamicList<face>& boundaryFaces,
274  DynamicList<label>& cellCorrespondence,
275  DynamicList<label>& unvPointID // unv index
276 )
277 {
278  Info<< "Starting reading cells at line " << is.lineNumber() << '.' << endl;
279 
280  // Invert point numbering.
281  label maxUnvPoint = 0;
282  forAll(unvPointID, pointi)
283  {
284  maxUnvPoint = max(maxUnvPoint, unvPointID[pointi]);
285  }
286  labelList unvToFoam(invert(maxUnvPoint+1, unvPointID));
287 
288 
289  const cellModel& hex = *(cellModeller::lookup("hex"));
290  const cellModel& prism = *(cellModeller::lookup("prism"));
291  const cellModel& tet = *(cellModeller::lookup("tet"));
292 
293  labelHashSet skippedElements;
294 
295  labelHashSet foundFeType;
296 
297  while (true)
298  {
299  string line;
300  is.getLine(line);
301 
302  if (isSeparator(line))
303  {
304  break;
305  }
306 
307  label celli, feID, physProp, matProp, colour, nNodes;
308 
309  IStringStream lineStr(line);
310  lineStr
311  >> celli >> feID >> physProp >> matProp >> colour >> nNodes;
312 
313  if (foundFeType.insert(feID))
314  {
315  Info<< "First occurrence of element type " << feID
316  << " for cell " << celli << " at line "
317  << is.lineNumber() << endl;
318  }
319 
320  if (feID == 11)
321  {
322  // Rod. Skip.
323  is.getLine(line);
324  is.getLine(line);
325  }
326  else if (feID == 171)
327  {
328  // Rod. Skip.
329  is.getLine(line);
330  }
331  else if (feID == 41 || feID == 91)
332  {
333  // Triangle. Save - used for patching later on.
334  is.getLine(line);
335 
336  face cVerts(3);
337  IStringStream lineStr(line);
338  lineStr
339  >> cVerts[0] >> cVerts[1] >> cVerts[2];
340  boundaryFaces.append(cVerts);
341  boundaryFaceIndices.append(celli);
342  }
343  else if (feID == 44 || feID == 94)
344  {
345  // Quad. Save - used for patching later on.
346  is.getLine(line);
347 
348  face cVerts(4);
349  IStringStream lineStr(line);
350  lineStr
351  >> cVerts[0] >> cVerts[1] >> cVerts[2] >> cVerts[3];
352  boundaryFaces.append(cVerts);
353  boundaryFaceIndices.append(celli);
354  }
355  else if (feID == 111)
356  {
357  // Tet.
358  is.getLine(line);
359 
360  labelList cVerts(4);
361  IStringStream lineStr(line);
362  lineStr
363  >> cVerts[0] >> cVerts[1] >> cVerts[2] >> cVerts[3];
364 
365  cellVerts.append(cellShape(tet, cVerts, true));
366  cellMaterial.append(physProp);
367  addAndExtend(cellCorrespondence,celli,cellMaterial.size()-1);
368 
369  if (cellVerts.last().size() != cVerts.size())
370  {
371  Info<< "Line:" << is.lineNumber()
372  << " element:" << celli
373  << " type:" << feID
374  << " collapsed from " << cVerts << nl
375  << " to:" << cellVerts.last()
376  << endl;
377  }
378  }
379  else if (feID == 112)
380  {
381  // Wedge.
382  is.getLine(line);
383 
384  labelList cVerts(6);
385  IStringStream lineStr(line);
386  lineStr
387  >> cVerts[0] >> cVerts[1] >> cVerts[2]
388  >> cVerts[3] >> cVerts[4] >> cVerts[5];
389 
390  cellVerts.append(cellShape(prism, cVerts, true));
391  cellMaterial.append(physProp);
392  addAndExtend(cellCorrespondence,celli,cellMaterial.size()-1);
393 
394  if (cellVerts.last().size() != cVerts.size())
395  {
396  Info<< "Line:" << is.lineNumber()
397  << " element:" << celli
398  << " type:" << feID
399  << " collapsed from " << cVerts << nl
400  << " to:" << cellVerts.last()
401  << endl;
402  }
403  }
404  else if (feID == 115)
405  {
406  // Hex.
407  is.getLine(line);
408 
409  labelList cVerts(8);
410  IStringStream lineStr(line);
411  lineStr
412  >> cVerts[0] >> cVerts[1] >> cVerts[2] >> cVerts[3]
413  >> cVerts[4] >> cVerts[5] >> cVerts[6] >> cVerts[7];
414 
415  cellVerts.append(cellShape(hex, cVerts, true));
416  cellMaterial.append(physProp);
417  addAndExtend(cellCorrespondence,celli,cellMaterial.size()-1);
418 
419  if (cellVerts.last().size() != cVerts.size())
420  {
421  Info<< "Line:" << is.lineNumber()
422  << " element:" << celli
423  << " type:" << feID
424  << " collapsed from " << cVerts << nl
425  << " to:" << cellVerts.last()
426  << endl;
427  }
428  }
429  else if (feID == 118)
430  {
431  // Parabolic Tet
432  is.getLine(line);
433 
434  labelList cVerts(4);
435  label dummy;
436  {
437  IStringStream lineStr(line);
438  lineStr
439  >> cVerts[0] >> dummy >> cVerts[1] >> dummy >> cVerts[2];
440  }
441  is.getLine(line);
442  {
443  IStringStream lineStr(line);
444  lineStr >> dummy>> cVerts[3];
445  }
446 
447  cellVerts.append(cellShape(tet, cVerts, true));
448  cellMaterial.append(physProp);
449  addAndExtend(cellCorrespondence,celli,cellMaterial.size()-1);
450 
451  if (cellVerts.last().size() != cVerts.size())
452  {
453  Info<< "Line:" << is.lineNumber()
454  << " element:" << celli
455  << " type:" << feID
456  << " collapsed from " << cVerts << nl
457  << " to:" << cellVerts.last()
458  << endl;
459  }
460  }
461  else
462  {
463  if (skippedElements.insert(feID))
464  {
466  << "Cell type " << feID << " not supported" << endl;
467  }
468  is.getLine(line); // Do nothing
469  }
470  }
471 
472  cellVerts.shrink();
473  cellMaterial.shrink();
474  boundaryFaces.shrink();
475  boundaryFaceIndices.shrink();
476  cellCorrespondence.shrink();
477 
478  Info<< "Read " << cellVerts.size() << " cells"
479  << " and " << boundaryFaces.size() << " boundary faces." << endl;
480 }
481 
482 
483 void readSets
484 (
485  IFstream& is,
486  DynamicList<word>& patchNames,
487  DynamicList<labelList>& patchFaceIndices
488 )
489 {
490  Info<< "Starting reading patches at line " << is.lineNumber() << '.'
491  << endl;
492 
493  while (true)
494  {
495  string line;
496  is.getLine(line);
497 
498  if (isSeparator(line))
499  {
500  break;
501  }
502 
503  IStringStream lineStr(line);
504  label group, constraintSet, restraintSet, loadSet, dofSet,
505  tempSet, contactSet, nFaces;
506  lineStr
507  >> group >> constraintSet >> restraintSet >> loadSet
508  >> dofSet >> tempSet >> contactSet >> nFaces;
509 
510  is.getLine(line);
511  word groupName = string::validate<word>(line);
512 
513  Info<< "For group " << group
514  << " named " << groupName
515  << " trying to read " << nFaces << " patch face indices."
516  << endl;
517 
518  labelList groupIndices(nFaces);
519  label groupType = -1;
520  nFaces = 0;
521 
522  while (nFaces < groupIndices.size())
523  {
524  is.getLine(line);
525  IStringStream lineStr(line);
526 
527  // Read one (for last face) or two entries from line.
528  label nRead = 2;
529  if (nFaces == groupIndices.size()-1)
530  {
531  nRead = 1;
532  }
533 
534  for (label i = 0; i < nRead; i++)
535  {
536  label tag, nodeLeaf, component;
537 
538  lineStr >> groupType >> tag >> nodeLeaf >> component;
539 
540  groupIndices[nFaces++] = tag;
541  }
542  }
543 
544 
545  // Store
546  if (groupType == 8)
547  {
548  patchNames.append(groupName);
549  patchFaceIndices.append(groupIndices);
550  }
551  else
552  {
554  << "When reading patches expect entity type code 8"
555  << nl << " Skipping group code " << groupType
556  << endl;
557  }
558  }
559 
560  patchNames.shrink();
561  patchFaceIndices.shrink();
562 }
563 
564 
565 
566 // Read dof set (757)
567 void readDOFS
568 (
569  IFstream& is,
570  DynamicList<word>& patchNames,
571  DynamicList<labelList>& dofVertices
572 )
573 {
574  Info<< "Starting reading contraints at line " << is.lineNumber() << '.'
575  << endl;
576 
577  string line;
578  is.getLine(line);
579  label group;
580  {
581  IStringStream lineStr(line);
582  lineStr >> group;
583  }
584 
585  is.getLine(line);
586  {
587  IStringStream lineStr(line);
588  patchNames.append(lineStr);
589  }
590 
591  Info<< "For DOF set " << group
592  << " named " << patchNames.last()
593  << " trying to read vertex indices."
594  << endl;
595 
596  DynamicList<label> vertices;
597  while (true)
598  {
599  string line;
600  is.getLine(line);
601 
602  if (isSeparator(line))
603  {
604  break;
605  }
606 
607  IStringStream lineStr(line);
608  label pointi;
609  lineStr >> pointi;
610 
611  vertices.append(pointi);
612  }
613 
614  Info<< "For DOF set " << group
615  << " named " << patchNames.last()
616  << " read " << vertices.size() << " vertex indices." << endl;
617 
618  dofVertices.append(vertices.shrink());
619 
620  patchNames.shrink();
621  dofVertices.shrink();
622 }
623 
624 
625 // Returns -1 or group that all of the vertices of f are in,
626 label findPatch(const List<labelHashSet>& dofGroups, const face& f)
627 {
628  forAll(dofGroups, patchi)
629  {
630  if (dofGroups[patchi].found(f[0]))
631  {
632  bool allInGroup = true;
633 
634  // Check rest of face
635  for (label fp = 1; fp < f.size(); fp++)
636  {
637  if (!dofGroups[patchi].found(f[fp]))
638  {
639  allInGroup = false;
640  break;
641  }
642  }
643 
644  if (allInGroup)
645  {
646  return patchi;
647  }
648  }
649  }
650  return -1;
651 }
652 
653 
654 
655 int main(int argc, char *argv[])
656 {
658  argList::validArgs.append(".unv file");
660  (
661  "dump",
662  "dump boundary faces as boundaryFaces.obj (for debugging)"
663  );
664 
665  #include "setRootCase.H"
666  #include "createTime.H"
667 
668  const fileName ideasName = args[1];
669  IFstream inFile(ideasName);
670 
671  if (!inFile.good())
672  {
674  << "Cannot open file " << ideasName
675  << exit(FatalError);
676  }
677 
678 
679  // Switch on additional debug info
680  const bool verbose = false; //true;
681 
682  // Unit scale factors
683  scalar lengthScale = 1;
684  scalar forceScale = 1;
685  scalar tempScale = 1;
686  scalar tempOffset = 0;
687 
688  // Points
690  // Original unv point label
691  DynamicList<label> unvPointID;
692 
693  // Cells
694  DynamicList<cellShape> cellVerts;
695  DynamicList<label> cellMat;
696  DynamicList<label> cellCorrespondence;
697 
698  // Boundary faces
699  DynamicList<label> boundaryFaceIndices;
700  DynamicList<face> boundaryFaces;
701 
702  // Patch names and patchFace indices.
704  DynamicList<labelList> patchFaceIndices;
705  DynamicList<labelList> dofVertIndices;
706 
707 
708  while (inFile.good())
709  {
710  label tag = readTag(inFile);
711 
712  if (tag == -1)
713  {
714  break;
715  }
716 
717  Info<< "Processing tag:" << tag << endl;
718 
719  switch (tag)
720  {
721  case 151:
722  readHeader(inFile);
723  break;
724 
725  case 164:
726  readUnits
727  (
728  inFile,
729  lengthScale,
730  forceScale,
731  tempScale,
732  tempOffset
733  );
734  break;
735 
736  case 2411:
737  readPoints(inFile, points, unvPointID);
738  break;
739 
740  case 2412:
741  readCells
742  (
743  inFile,
744  cellVerts,
745  cellMat,
746  boundaryFaceIndices,
747  boundaryFaces,
748  cellCorrespondence,
749  unvPointID
750  );
751  break;
752 
753  case 2452:
754  case 2467:
755  readSets
756  (
757  inFile,
758  patchNames,
759  patchFaceIndices
760  );
761  break;
762 
763  case 757:
764  readDOFS
765  (
766  inFile,
767  patchNames,
768  dofVertIndices
769  );
770  break;
771 
772  default:
773  Info<< "Skipping tag " << tag << " on line "
774  << inFile.lineNumber() << endl;
775  skipSection(inFile);
776  break;
777  }
778  Info<< endl;
779  }
780 
781 
782  // Invert point numbering.
783  label maxUnvPoint = 0;
784  forAll(unvPointID, pointi)
785  {
786  maxUnvPoint = max(maxUnvPoint, unvPointID[pointi]);
787  }
788  labelList unvToFoam(invert(maxUnvPoint+1, unvPointID));
789 
790 
791  // Renumber vertex numbers on cells
792 
793  forAll(cellVerts, celli)
794  {
795  labelList foamVerts
796  (
797  renumber
798  (
799  unvToFoam,
800  static_cast<labelList&>(cellVerts[celli])
801  )
802  );
803 
804  if (findIndex(foamVerts, -1) != -1)
805  {
807  << "Cell " << celli
808  << " unv vertices " << cellVerts[celli]
809  << " has some undefined vertices " << foamVerts
810  << abort(FatalError);
811  }
812 
813  // Bit nasty: replace vertex list.
814  cellVerts[celli].transfer(foamVerts);
815  }
816 
817  // Renumber vertex numbers on boundaryFaces
818 
819  forAll(boundaryFaces, bFacei)
820  {
821  labelList foamVerts(renumber(unvToFoam, boundaryFaces[bFacei]));
822 
823  if (findIndex(foamVerts, -1) != -1)
824  {
826  << "Boundary face " << bFacei
827  << " unv vertices " << boundaryFaces[bFacei]
828  << " has some undefined vertices " << foamVerts
829  << abort(FatalError);
830  }
831 
832  // Bit nasty: replace vertex list.
833  boundaryFaces[bFacei].transfer(foamVerts);
834  }
835 
836 
837 
838  // Patchify = create patchFaceVerts for use in cellShape construction
839  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
840 
841  // Works in one of two modes. Either has read boundaryFaces which
842  // just need to be sorted according to patch. Or has read point constraint
843  // sets (dofVertIndices).
844 
845  List<faceList> patchFaceVerts;
846 
847  labelList own(boundaryFaces.size(), -1);
848  labelList nei(boundaryFaces.size(), -1);
849  HashTable<label, label> faceToCell[2];
850 
851  {
852  HashTable<label, face, Hash<face>> faceToFaceID(boundaryFaces.size());
853  forAll(boundaryFaces, facei)
854  {
855  SortableList<label> sortedVerts(boundaryFaces[facei]);
856  faceToFaceID.insert(face(sortedVerts), facei);
857  }
858 
859  forAll(cellVerts, celli)
860  {
861  faceList faces = cellVerts[celli].faces();
862  forAll(faces, i)
863  {
864  SortableList<label> sortedVerts(faces[i]);
865  HashTable<label, face, Hash<face>>::const_iterator fnd =
866  faceToFaceID.find(face(sortedVerts));
867 
868  if (fnd != faceToFaceID.end())
869  {
870  label facei = fnd();
871  int stat = face::compare(faces[i], boundaryFaces[facei]);
872 
873  if (stat == 1)
874  {
875  // Same orientation. Cell is owner.
876  own[facei] = celli;
877  }
878  else if (stat == -1)
879  {
880  // Opposite orientation. Cell is neighbour.
881  nei[facei] = celli;
882  }
883  }
884  }
885  }
886 
887  label nReverse = 0;
888  forAll(own, facei)
889  {
890  if (own[facei] == -1 && nei[facei] != -1)
891  {
892  // Boundary face with incorrect orientation
893  boundaryFaces[facei] = boundaryFaces[facei].reverseFace();
894  Swap(own[facei], nei[facei]);
895  nReverse++;
896  }
897  }
898  if (nReverse > 0)
899  {
900  Info << "Found " << nReverse << " reversed boundary faces out of "
901  << boundaryFaces.size() << endl;
902  }
903 
904 
905  label cnt = 0;
906  forAll(own, facei)
907  {
908  if (own[facei] != -1 && nei[facei] != -1)
909  {
910  faceToCell[1].insert(facei, own[facei]);
911  faceToCell[0].insert(facei, nei[facei]);
912  cnt++;
913  }
914  }
915 
916  if (cnt > 0)
917  {
918  Info << "Of " << boundaryFaces.size() << " so-called"
919  << " boundary faces " << cnt << " belong to two cells "
920  << "and are therefore internal" << endl;
921  }
922  }
923 
924  HashTable<labelList,word> cellZones;
925  HashTable<labelList,word> faceZones;
926  List<bool> isAPatch(patchNames.size(),true);
927 
928  if (dofVertIndices.size())
929  {
930  // Use the vertex constraints to patch. Is of course bit dodgy since
931  // face goes on patch if all its vertices are on a constraint.
932  // Note: very inefficient since goes through all faces (including
933  // internal ones) twice. Should do a construct without patches
934  // and then repatchify.
935 
936  Info<< "Using " << dofVertIndices.size()
937  << " DOF sets to detect boundary faces."<< endl;
938 
939  // Renumber vertex numbers on contraints
940  forAll(dofVertIndices, patchi)
941  {
942  inplaceRenumber(unvToFoam, dofVertIndices[patchi]);
943  }
944 
945 
946  // Build labelHashSet of points per dofGroup/patch
947  List<labelHashSet> dofGroups(dofVertIndices.size());
948 
949  forAll(dofVertIndices, patchi)
950  {
951  const labelList& foamVerts = dofVertIndices[patchi];
952 
953  forAll(foamVerts, i)
954  {
955  dofGroups[patchi].insert(foamVerts[i]);
956  }
957  }
958 
959  List<DynamicList<face>> dynPatchFaces(dofVertIndices.size());
960 
961  forAll(cellVerts, celli)
962  {
963  const cellShape& shape = cellVerts[celli];
964 
965  const faceList shapeFaces(shape.faces());
966 
967  forAll(shapeFaces, i)
968  {
969  label patchi = findPatch(dofGroups, shapeFaces[i]);
970 
971  if (patchi != -1)
972  {
973  dynPatchFaces[patchi].append(shapeFaces[i]);
974  }
975  }
976  }
977 
978  // Transfer
979  patchFaceVerts.setSize(dynPatchFaces.size());
980 
981  forAll(dynPatchFaces, patchi)
982  {
983  patchFaceVerts[patchi].transfer(dynPatchFaces[patchi]);
984  }
985  }
986  else
987  {
988  // Use the boundary faces.
989 
990  // Construct the patch faces by sorting the boundaryFaces according to
991  // patch.
992  patchFaceVerts.setSize(patchFaceIndices.size());
993 
994  Info<< "Sorting boundary faces according to group (patch)" << endl;
995 
996  // make sure that no face is used twice on the boundary
997  // (possible for boundary-only faceZones)
998  labelHashSet alreadyOnBoundary;
999 
1000  // Construct map from boundaryFaceIndices
1001  Map<label> boundaryFaceToIndex(boundaryFaceIndices.size());
1002 
1003  forAll(boundaryFaceIndices, i)
1004  {
1005  boundaryFaceToIndex.insert(boundaryFaceIndices[i], i);
1006  }
1007 
1008  forAll(patchFaceVerts, patchi)
1009  {
1010  Info << patchi << ": " << patchNames[patchi] << " is " << flush;
1011 
1012  faceList& patchFaces = patchFaceVerts[patchi];
1013  const labelList& faceIndices = patchFaceIndices[patchi];
1014 
1015  patchFaces.setSize(faceIndices.size());
1016 
1017  bool duplicateFaces = false;
1018 
1019  label cnt = 0;
1020  forAll(patchFaces, i)
1021  {
1022  if (boundaryFaceToIndex.found(faceIndices[i]))
1023  {
1024  label bFacei = boundaryFaceToIndex[faceIndices[i]];
1025 
1026  if (own[bFacei] != -1 && nei[bFacei] == -1)
1027  {
1028  patchFaces[cnt] = boundaryFaces[bFacei];
1029  cnt++;
1030  if (alreadyOnBoundary.found(bFacei))
1031  {
1032  duplicateFaces = true;
1033  }
1034  }
1035  }
1036  }
1037 
1038  if (cnt != patchFaces.size() || duplicateFaces)
1039  {
1040  isAPatch[patchi] = false;
1041 
1042  if (verbose)
1043  {
1044  if (cnt != patchFaces.size())
1045  {
1047  << "For patch " << patchi << " there were "
1048  << patchFaces.size()-cnt
1049  << " faces not used because they seem"
1050  << " to be internal. "
1051  << "This seems to be a face or a cell-zone"
1052  << endl;
1053  }
1054  else
1055  {
1057  << "Patch "
1058  << patchi << " has faces that are already "
1059  << " in use on other boundary-patches,"
1060  << " Assuming faceZoneset." << endl;
1061  }
1062  }
1063 
1064  patchFaces.setSize(0); // Assume that this is no patch at all
1065 
1066  if (cellCorrespondence[faceIndices[0]] >= 0)
1067  {
1068  Info << "cellZone" << endl;
1069  labelList theCells(faceIndices.size());
1070  forAll(faceIndices, i)
1071  {
1072  if (cellCorrespondence[faceIndices[0]] < 0)
1073  {
1075  << "The face index " << faceIndices[i]
1076  << " was not found amongst the cells."
1077  << " This kills the theory that "
1078  << patchNames[patchi] << " is a cell zone"
1079  << endl
1080  << abort(FatalError);
1081  }
1082  theCells[i] = cellCorrespondence[faceIndices[i]];
1083  }
1084  cellZones.insert(patchNames[patchi], theCells);
1085  }
1086  else
1087  {
1088  Info << "faceZone" << endl;
1089  labelList theFaces(faceIndices.size());
1090  forAll(faceIndices, i)
1091  {
1092  theFaces[i] = boundaryFaceToIndex[faceIndices[i]];
1093  }
1094  faceZones.insert(patchNames[patchi],theFaces);
1095  }
1096  }
1097  else
1098  {
1099  Info << "patch" << endl;
1100 
1101  forAll(patchFaces, i)
1102  {
1103  label bFacei = boundaryFaceToIndex[faceIndices[i]];
1104  alreadyOnBoundary.insert(bFacei);
1105  }
1106  }
1107  }
1108  }
1109 
1110  pointField polyPoints;
1111  polyPoints.transfer(points);
1112 
1113  // Length scaling factor
1114  polyPoints /= lengthScale;
1115 
1116 
1117  // For debugging: dump boundary faces as OBJ surface mesh
1118  if (args.optionFound("dump"))
1119  {
1120  Info<< "Writing boundary faces to OBJ file boundaryFaces.obj"
1121  << nl << endl;
1122 
1123  // Create globally numbered surface
1124  meshedSurface rawSurface
1125  (
1126  xferCopy(polyPoints),
1127  xferCopyTo<faceList>(boundaryFaces)
1128  );
1129 
1130  // Write locally numbered surface
1132  (
1133  xferCopy(rawSurface.localPoints()),
1134  xferCopy(rawSurface.localFaces())
1135  ).write(runTime.path()/"boundaryFaces.obj");
1136  }
1137 
1138 
1139  Info<< "\nConstructing mesh with non-default patches of size:" << nl;
1140  DynamicList<word> usedPatchNames;
1141  DynamicList<faceList> usedPatchFaceVerts;
1142 
1143  forAll(patchNames, patchi)
1144  {
1145  if (isAPatch[patchi])
1146  {
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 // ************************************************************************* //
ListType renumber(const labelUList &oldToNew, const ListType &)
Renumber the values (not the indices) of a list.
const char *const group
Group name for atomic constants.
point centre(const pointField &) const
Centre point of face.
Definition: face.C:488
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
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 NULL.
Definition: cellModeller.C:88
A class for handling file names.
Definition: fileName.H:69
void inplaceRenumber(const labelUList &oldToNew, ListType &)
Inplace renumber the values of a list.
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:106
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
static const char tab
Definition: Ostream.H:261
error FatalError
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
A list that is sorted upon construction or when explicitly requested with the sort() method...
Definition: List.H:68
An analytical geometric cellShape.
Definition: cellShape.H:69
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:76
bool good() const
Return true if next operation might succeed.
Definition: IOstream.H:333
const dimensionedScalar c2
Second radiation constant: default SI units: [m.K].
static word defaultRegion
Return the default region name.
Definition: polyMesh.H:306
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:253
static void noParallel()
Remove the parallel options.
Definition: argList.C:146
static SLList< string > validArgs
A list of valid (mandatory) arguments.
Definition: argList.H:154
Xfer< T > xferCopy(const T &)
Construct by copying the contents of the arg.
label size() const
Return number of elements in table.
Definition: HashTableI.H:65
const vectorField & faceCentres() const
bool insert(const Key &key)
Insert a new entry.
Definition: HashSet.H:116
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:979
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:138
unsigned operator()(const PrimitiveType &p, unsigned seed) const
Definition: Hash.H:61
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.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
labelList invert(const label len, const labelUList &)
Invert one-to-one map. Unmapped elements will be -1.
Definition: ListOps.C:37
void setSize(const label)
Alter the addressed list size.
Definition: DynamicListI.H:163
dynamicFvMesh & mesh
void Swap(T &a, T &b)
Definition: Swap.H:43
const pointField & points
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects...
Definition: DynamicList.H:56
bool optionFound(const word &opt) const
Return true if the named option is found.
Definition: argListI.H:108
Xfer< T > xferMove(T &)
Construct by transferring the contents of the arg.
A class for handling words, derived from string.
Definition: word.H:59
MeshedSurface< face > meshedSurface
static int compare(const face &, const face &)
Compare faces.
Definition: face.C:298
const T * cdata() const
Return a const pointer to the first data element,.
Definition: UListI.H:142
wordList patchNames(nPatches)
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
Definition: DynamicListI.H:292
bool readScalar(const char *buf, doubleScalar &s)
Read whole of buf as a scalar. Return true if succesful.
Definition: doubleScalar.H:63
An STL-conforming hash table.
Definition: HashTable.H:61
errorManip< error > abort(error &err)
Definition: errorManip.H:131
label readLabel(Istream &is)
Definition: label.H:64
List< Key > toc() const
Return the table of contents.
Definition: HashTable.C:198
static void write(const fileName &, const MeshedSurface< Face > &)
Write to file.
bool found(const Key &) const
Return true if hashedEntry is found in table.
Definition: HashTable.C:109
DynamicList< T, SizeInc, SizeMult, SizeDiv > & shrink()
Shrink the allocated space to the number of elements used.
Definition: DynamicListI.H:240
static const char nl
Definition: Ostream.H:262
Input from file stream.
Definition: IFstream.H:81
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurence of given element and return index,.
A subset of mesh cells.
Definition: cellZone.H:61
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
List< word > wordList
A List of words.
Definition: fileName.H:54
void setSize(const label)
Reset size of List.
Definition: List.C:295
faceList faces() const
Faces of this cell.
Definition: cellShapeI.H:162
const cellZoneMesh & cellZones() const
Return cell zone mesh.
Definition: polyMesh.H:469
unsigned Hasher(const void *data, size_t len, unsigned seed=0)
Bob Jenkins&#39;s 96-bit mixer hashing function (lookup3)
Definition: Hasher.C:476
label patchi
Ostream & flush(Ostream &os)
Flush stream.
Definition: Ostream.H:245
virtual const labelList & faceNeighbour() const
Return face neighbour.
Definition: polyMesh.C:1023
#define WarningInFunction
Report a warning using Foam::Warning.
Input from memory buffer stream.
Definition: IStringStream.H:49
label lineNumber() const
Return current stream line number.
Definition: IOstream.H:438
Maps a geometry to a set of cell primitives, which enables geometric cell data to be calculated witho...
Definition: cellModel.H:64
const dimensionedScalar c1
First radiation constant: default SI units: [W/m2].
messageStream Info
dimensioned< scalar > mag(const dimensioned< Type > &)
const faceZoneMesh & faceZones() const
Return face zone mesh.
Definition: polyMesh.H:463
void addZones(const List< pointZone * > &pz, const List< faceZone * > &fz, const List< cellZone * > &cz)
Add mesh zones.
Definition: polyMesh.C:922
#define IOWarningInFunction(ios)
Report an IO warning using Foam::Warning.
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:83
A subset of mesh faces organised as a primitive patch.
Definition: faceZone.H:64
A class for managing temporary objects.
Definition: PtrList.H:54
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1017
T & last()
Return the last element of the list.
Definition: UListI.H:128
Foam::argList args(argc, argv)
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:91
void transfer(List< T > &)
Transfer contents of the argument List into this.
Definition: DynamicListI.H:259
void component(FieldField< Field, typename FieldField< Field, Type >::cmptType > &sf, const FieldField< Field, Type > &f, const direction d)
bool found
void transfer(List< T > &)
Transfer the contents of the argument List into this list.
Definition: List.C:365
virtual bool write() const
Write mesh using IO settings from time.
Definition: fvMesh.C:870
Namespace for OpenFOAM.
IOporosityModelList pZones(mesh)
ISstream & getLine(string &)
Raw, low-level getline into a string function.
Definition: ISstreamI.H:77