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 | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2011-2017 OpenFOAM Foundation
9 Copyright (C) 2021 OpenCFD Ltd.
10-------------------------------------------------------------------------------
11License
12 This file is part of OpenFOAM.
13
14 OpenFOAM is free software: you can redistribute it and/or modify it
15 under the terms of the GNU General Public License as published by
16 the Free Software Foundation, either version 3 of the License, or
17 (at your option) any later version.
18
19 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22 for more details.
23
24 You should have received a copy of the GNU General Public License
25 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26
27Application
28 ideasUnvToFoam
29
30Group
31 grpMeshConversionUtilities
32
33Description
34 I-Deas unv format mesh conversion.
35
36 Uses either
37 - DOF sets (757) or
38 - face groups (2452(Cubit), 2467)
39 to do patching.
40 Works without but then puts all faces in defaultFaces patch.
41
42\*---------------------------------------------------------------------------*/
43
44#include "argList.H"
45#include "polyMesh.H"
46#include "Time.H"
47#include "IFstream.H"
48#include "cellModel.H"
49#include "cellSet.H"
50#include "faceSet.H"
51#include "DynamicList.H"
52
53#include <cassert>
54#include "MeshedSurfaces.H"
55
56using namespace Foam;
57
58// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
59
60const string SEPARATOR(" -1");
61
62bool isSeparator(const std::string& line)
63{
64 return line.substr(0, 6) == SEPARATOR;
65}
66
67
68// Reads past -1 and reads element type
69label readTag(IFstream& is)
70{
71 string tag;
72 do
73 {
74 if (!is.good())
75 {
76 return -1;
77 }
78
79 string line;
80
81 is.getLine(line);
82
83 if (line.size() < 6)
84 {
85 return -1;
86 }
87
88 tag = line.substr(0, 6);
89
90 } while (tag == SEPARATOR);
91
92 return readLabel(tag);
93}
94
95
96// Reads and prints header
97void readHeader(IFstream& is)
98{
99 string line;
100
101 while (is.good())
102 {
103 is.getLine(line);
104
105 if (isSeparator(line))
106 {
107 break;
108 }
109 else
110 {
111 Info<< line << endl;
112 }
113 }
114}
115
116
117// Skip
118void skipSection(IFstream& is)
119{
120 Info<< "Skipping section at line " << is.lineNumber() << '.' << endl;
121
122 string line;
123
124 while (is.good())
125 {
126 is.getLine(line);
127
128 if (isSeparator(line))
129 {
130 break;
131 }
132 }
133}
134
135
136scalar readUnvScalar(const std::string& unvString)
137{
138 string s(unvString);
139
140 s.replaceAll("d", "E");
141 s.replaceAll("D", "E");
142
143 return readScalar(s);
144}
145
146
147// Reads unit section
148void readUnits
149(
150 IFstream& is,
151 scalar& lengthScale,
152 scalar& forceScale,
153 scalar& tempScale,
154 scalar& tempOffset
155)
156{
157 Info<< "Starting reading units at line " << is.lineNumber() << '.' << endl;
158
159 string line;
160 is.getLine(line);
161
162 label l = readLabel(line.substr(0, 10));
163 Info<< "l:" << l << endl;
164
165 string units(line.substr(10, 20));
166 Info<< "units:" << units << endl;
167
168 label unitType = readLabel(line.substr(30, 10));
169 Info<< "unitType:" << unitType << endl;
170
171 // Read lengthscales
172 is.getLine(line);
173
174 lengthScale = readUnvScalar(line.substr(0, 25));
175 forceScale = readUnvScalar(line.substr(25, 25));
176 tempScale = readUnvScalar(line.substr(50, 25));
177
178 is.getLine(line);
179 tempOffset = readUnvScalar(line.substr(0, 25));
180
181 Info<< "Unit factors:" << nl
182 << " Length scale : " << lengthScale << nl
183 << " Force scale : " << forceScale << nl
184 << " Temperature scale : " << tempScale << nl
185 << " Temperature offset : " << tempOffset << nl
186 << endl;
187}
188
189
190// Reads points section. Read region as well?
191void readPoints
192(
193 IFstream& is,
194 DynamicList<point>& points, // coordinates
195 DynamicList<label>& unvPointID // unv index
196)
197{
198 Info<< "Starting reading points at line " << is.lineNumber() << '.' << endl;
199
200 static bool hasWarned = false;
201
202 while (true)
203 {
204 string line;
205 is.getLine(line);
206
207 label pointi = readLabel(line.substr(0, 10));
208
209 if (pointi == -1)
210 {
211 break;
212 }
213 else if (pointi != points.size()+1 && !hasWarned)
214 {
215 hasWarned = true;
216
218 << "Points not in order starting at point " << pointi
219 << endl;
220 }
221
222 point pt;
223 is.getLine(line);
224 pt[0] = readUnvScalar(line.substr(0, 25));
225 pt[1] = readUnvScalar(line.substr(25, 25));
226 pt[2] = readUnvScalar(line.substr(50, 25));
227
228 unvPointID.append(pointi);
229 points.append(pt);
230 }
231
232 points.shrink();
233 unvPointID.shrink();
234
235 Info<< "Read " << points.size() << " points." << endl;
236}
237
238void addAndExtend
239(
240 DynamicList<label>& indizes,
241 label celli,
242 label val
243)
244{
245 if (indizes.size() < (celli+1))
246 {
247 indizes.setSize(celli+1,-1);
248 }
249 indizes[celli] = val;
250}
251
252// Reads cells section. Read region as well? Not handled yet but should just
253// be a matter of reading corresponding to boundaryFaces the correct property
254// and sorting it later on.
255void readCells
256(
257 IFstream& is,
258 DynamicList<cellShape>& cellVerts,
259 DynamicList<label>& cellMaterial,
260 DynamicList<label>& boundaryFaceIndices,
261 DynamicList<face>& boundaryFaces,
262 DynamicList<label>& cellCorrespondence,
263 DynamicList<label>& unvPointID // unv index
264)
265{
266 Info<< "Starting reading cells at line " << is.lineNumber() << '.' << endl;
267
268 // Invert point numbering.
269 label maxUnvPoint = 0;
270 forAll(unvPointID, pointi)
271 {
272 maxUnvPoint = max(maxUnvPoint, unvPointID[pointi]);
273 }
274 labelList unvToFoam(invert(maxUnvPoint+1, unvPointID));
275
276
277 const cellModel& hex = cellModel::ref(cellModel::HEX);
278 const cellModel& prism = cellModel::ref(cellModel::PRISM);
279 const cellModel& tet = cellModel::ref(cellModel::TET);
280
281 labelHashSet skippedElements;
282
283 labelHashSet foundFeType;
284
285 while (true)
286 {
287 string line;
288 is.getLine(line);
289
290 if (isSeparator(line))
291 {
292 break;
293 }
294
295 label celli, feID, physProp, matProp, colour, nNodes;
296
297 IStringStream lineStr(line);
298 lineStr
299 >> celli >> feID >> physProp >> matProp >> colour >> nNodes;
300
301 if (foundFeType.insert(feID))
302 {
303 Info<< "First occurrence of element type " << feID
304 << " for cell " << celli << " at line "
305 << is.lineNumber() << endl;
306 }
307
308 if (feID == 11)
309 {
310 // Rod. Skip.
311 is.getLine(line);
312 is.getLine(line);
313 }
314 else if (feID == 171)
315 {
316 // Rod. Skip.
317 is.getLine(line);
318 }
319 else if (feID == 41 || feID == 91)
320 {
321 // Triangle. Save - used for patching later on.
322 is.getLine(line);
323
324 face cVerts(3);
325 IStringStream lineStr(line);
326 lineStr
327 >> cVerts[0] >> cVerts[1] >> cVerts[2];
328 boundaryFaces.append(cVerts);
329 boundaryFaceIndices.append(celli);
330 }
331 else if (feID == 44 || feID == 94)
332 {
333 // Quad. Save - used for patching later on.
334 is.getLine(line);
335
336 face cVerts(4);
337 IStringStream lineStr(line);
338 lineStr
339 >> cVerts[0] >> cVerts[1] >> cVerts[2] >> cVerts[3];
340 boundaryFaces.append(cVerts);
341 boundaryFaceIndices.append(celli);
342 }
343 else if (feID == 111)
344 {
345 // Tet.
346 is.getLine(line);
347
348 labelList cVerts(4);
349 IStringStream lineStr(line);
350 lineStr
351 >> cVerts[0] >> cVerts[1] >> cVerts[2] >> cVerts[3];
352
353 cellVerts.append(cellShape(tet, cVerts, true));
354 cellMaterial.append(physProp);
355 addAndExtend(cellCorrespondence,celli,cellMaterial.size()-1);
356
357 if (cellVerts.last().size() != cVerts.size())
358 {
359 Info<< "Line:" << is.lineNumber()
360 << " element:" << celli
361 << " type:" << feID
362 << " collapsed from " << cVerts << nl
363 << " to:" << cellVerts.last()
364 << endl;
365 }
366 }
367 else if (feID == 112)
368 {
369 // Wedge.
370 is.getLine(line);
371
372 labelList cVerts(6);
373 IStringStream lineStr(line);
374 lineStr
375 >> cVerts[0] >> cVerts[1] >> cVerts[2]
376 >> cVerts[3] >> cVerts[4] >> cVerts[5];
377
378 cellVerts.append(cellShape(prism, cVerts, true));
379 cellMaterial.append(physProp);
380 addAndExtend(cellCorrespondence,celli,cellMaterial.size()-1);
381
382 if (cellVerts.last().size() != cVerts.size())
383 {
384 Info<< "Line:" << is.lineNumber()
385 << " element:" << celli
386 << " type:" << feID
387 << " collapsed from " << cVerts << nl
388 << " to:" << cellVerts.last()
389 << endl;
390 }
391 }
392 else if (feID == 115)
393 {
394 // Hex.
395 is.getLine(line);
396
397 labelList cVerts(8);
398 IStringStream lineStr(line);
399 lineStr
400 >> cVerts[0] >> cVerts[1] >> cVerts[2] >> cVerts[3]
401 >> cVerts[4] >> cVerts[5] >> cVerts[6] >> cVerts[7];
402
403 cellVerts.append(cellShape(hex, cVerts, true));
404 cellMaterial.append(physProp);
405 addAndExtend(cellCorrespondence,celli,cellMaterial.size()-1);
406
407 if (cellVerts.last().size() != cVerts.size())
408 {
409 Info<< "Line:" << is.lineNumber()
410 << " element:" << celli
411 << " type:" << feID
412 << " collapsed from " << cVerts << nl
413 << " to:" << cellVerts.last()
414 << endl;
415 }
416 }
417 else if (feID == 118)
418 {
419 // Parabolic Tet
420 is.getLine(line);
421
422 labelList cVerts(4);
423 label dummy;
424 {
425 IStringStream lineStr(line);
426 lineStr
427 >> cVerts[0] >> dummy >> cVerts[1] >> dummy >> cVerts[2];
428 }
429 is.getLine(line);
430 {
431 IStringStream lineStr(line);
432 lineStr >> dummy>> cVerts[3];
433 }
434
435 cellVerts.append(cellShape(tet, cVerts, true));
436 cellMaterial.append(physProp);
437 addAndExtend(cellCorrespondence,celli,cellMaterial.size()-1);
438
439 if (cellVerts.last().size() != cVerts.size())
440 {
441 Info<< "Line:" << is.lineNumber()
442 << " element:" << celli
443 << " type:" << feID
444 << " collapsed from " << cVerts << nl
445 << " to:" << cellVerts.last()
446 << endl;
447 }
448 }
449 else
450 {
451 if (skippedElements.insert(feID))
452 {
454 << "Cell type " << feID << " not supported" << endl;
455 }
456
457 is.getLine(line);
458 }
459 }
460
461 cellVerts.shrink();
462 cellMaterial.shrink();
463 boundaryFaces.shrink();
464 boundaryFaceIndices.shrink();
465 cellCorrespondence.shrink();
466
467 Info<< "Read " << cellVerts.size() << " cells"
468 << " and " << boundaryFaces.size() << " boundary faces." << endl;
469
470 if (!cellVerts.size())
471 {
473 << "There are no cells in the mesh." << endl
474 << " Note: 2D meshes are not supported."<< endl;
475 }
476}
477
478
479void readSets
480(
481 IFstream& is,
483 DynamicList<labelList>& patchFaceIndices
484)
485{
486 Info<< "Starting reading patches at line " << is.lineNumber() << '.'
487 << endl;
488
489 while (true)
490 {
491 string line;
492 is.getLine(line);
493
494 if (isSeparator(line))
495 {
496 break;
497 }
498
499 IStringStream lineStr(line);
500 label group, constraintSet, restraintSet, loadSet, dofSet,
501 tempSet, contactSet, nFaces;
502 lineStr
503 >> group >> constraintSet >> restraintSet >> loadSet
504 >> dofSet >> tempSet >> contactSet >> nFaces;
505
506 is.getLine(line);
507 const word groupName = word::validate(line);
508
509 Info<< "For group " << group
510 << " named " << groupName
511 << " trying to read " << nFaces << " patch face indices."
512 << endl;
513
514 labelList groupIndices(nFaces);
515 label groupType = -1;
516 nFaces = 0;
517
518 while (nFaces < groupIndices.size())
519 {
520 is.getLine(line);
521 IStringStream lineStr(line);
522
523 // Read one (for last face) or two entries from line.
524 label nRead = 2;
525 if (nFaces == groupIndices.size()-1)
526 {
527 nRead = 1;
528 }
529
530 for (label i = 0; i < nRead; i++)
531 {
532 label tag, nodeLeaf, component;
533
534 lineStr >> groupType >> tag >> nodeLeaf >> component;
535
536 groupIndices[nFaces++] = tag;
537 }
538 }
539
540
541 // Store
542 if (groupType == 8)
543 {
544 patchNames.append(groupName);
545 patchFaceIndices.append(groupIndices);
546 }
547 else
548 {
550 << "When reading patches expect entity type code 8"
551 << nl << " Skipping group code " << groupType
552 << endl;
553 }
554 }
555
556 patchNames.shrink();
557 patchFaceIndices.shrink();
558}
559
560
561
562// Read dof set (757)
563void readDOFS
564(
565 IFstream& is,
567 DynamicList<labelList>& dofVertices
568)
569{
570 Info<< "Starting reading constraints at line " << is.lineNumber() << '.'
571 << endl;
572
573 string line;
574 is.getLine(line);
575 label group;
576 {
577 IStringStream lineStr(line);
578 lineStr >> group;
579 }
580
581 is.getLine(line);
582 {
583 IStringStream lineStr(line);
584 word pName(lineStr);
585 patchNames.append(pName);
586 }
587
588 Info<< "For DOF set " << group
589 << " named " << patchNames.last()
590 << " trying to read vertex indices."
591 << endl;
592
594 while (true)
595 {
596 string line;
597 is.getLine(line);
598
599 if (isSeparator(line))
600 {
601 break;
602 }
603
604 IStringStream lineStr(line);
605 label pointi;
606 lineStr >> pointi;
607
608 vertices.append(pointi);
609 }
610
611 Info<< "For DOF set " << group
612 << " named " << patchNames.last()
613 << " read " << vertices.size() << " vertex indices." << endl;
614
615 dofVertices.append(vertices.shrink());
616
617 patchNames.shrink();
618 dofVertices.shrink();
619}
620
621
622// Returns -1 or group that all of the vertices of f are in,
623label findPatch(const List<labelHashSet>& dofGroups, const face& f)
624{
625 forAll(dofGroups, patchi)
626 {
627 if (dofGroups[patchi].found(f[0]))
628 {
629 bool allInGroup = true;
630
631 // Check rest of face
632 for (label fp = 1; fp < f.size(); fp++)
633 {
634 if (!dofGroups[patchi].found(f[fp]))
635 {
636 allInGroup = false;
637 break;
638 }
639 }
640
641 if (allInGroup)
642 {
643 return patchi;
644 }
645 }
646 }
647 return -1;
648}
649
650
651
652int main(int argc, char *argv[])
653{
654 argList::addNote
655 (
656 "Convert I-Deas unv format to OpenFOAM"
657 );
658 argList::noParallel();
659 argList::addArgument(".unv file");
660 argList::addBoolOption
661 (
662 "dump",
663 "Dump boundary faces as boundaryFaces.obj (for debugging)"
664 );
665
666 #include "setRootCase.H"
667 #include "createTime.H"
668
669 const auto ideasName = args.get<fileName>(1);
670 IFstream inFile(ideasName);
671
672 if (!inFile.good())
673 {
675 << "Cannot open file " << ideasName
676 << exit(FatalError);
677 }
678
679
680 // Switch on additional debug info
681 const bool verbose = false; //true;
682
683 // Unit scale factors
684 scalar lengthScale = 1;
685 scalar forceScale = 1;
686 scalar tempScale = 1;
687 scalar tempOffset = 0;
688
689 // Points
691 // Original unv point label
692 DynamicList<label> unvPointID;
693
694 // Cells
695 DynamicList<cellShape> cellVerts;
696 DynamicList<label> cellMat;
697 DynamicList<label> cellCorrespondence;
698
699 // Boundary faces
700 DynamicList<label> boundaryFaceIndices;
701 DynamicList<face> boundaryFaces;
702
703 // Patch names and patchFace indices.
705 DynamicList<labelList> patchFaceIndices;
706 DynamicList<labelList> dofVertIndices;
707
708
709 while (inFile.good())
710 {
711 label tag = readTag(inFile);
712
713 if (tag == -1)
714 {
715 break;
716 }
717
718 Info<< "Processing tag:" << tag << endl;
719
720 switch (tag)
721 {
722 case 151:
723 readHeader(inFile);
724 break;
725
726 case 164:
727 readUnits
728 (
729 inFile,
730 lengthScale,
731 forceScale,
732 tempScale,
733 tempOffset
734 );
735 break;
736
737 case 2411:
738 readPoints(inFile, points, unvPointID);
739 break;
740
741 case 2412:
742 readCells
743 (
744 inFile,
745 cellVerts,
746 cellMat,
747 boundaryFaceIndices,
748 boundaryFaces,
749 cellCorrespondence,
750 unvPointID
751 );
752 break;
753
754 case 2452:
755 case 2467:
756 readSets
757 (
758 inFile,
760 patchFaceIndices
761 );
762 break;
763
764 case 757:
765 readDOFS
766 (
767 inFile,
769 dofVertIndices
770 );
771 break;
772
773 default:
774 Info<< "Skipping tag " << tag << " on line "
775 << inFile.lineNumber() << endl;
776 skipSection(inFile);
777 break;
778 }
779 Info<< endl;
780 }
781
782
783 // Invert point numbering.
784 label maxUnvPoint = 0;
785 forAll(unvPointID, pointi)
786 {
787 maxUnvPoint = max(maxUnvPoint, unvPointID[pointi]);
788 }
789 labelList unvToFoam(invert(maxUnvPoint+1, unvPointID));
790
791
792 // Renumber vertex numbers on cells
793
794 forAll(cellVerts, celli)
795 {
796 labelList foamVerts
797 (
798 renumber
799 (
800 unvToFoam,
801 static_cast<labelList&>(cellVerts[celli])
802 )
803 );
804
805 if (foamVerts.found(-1))
806 {
808 << "Cell " << celli
809 << " unv vertices " << cellVerts[celli]
810 << " has some undefined vertices " << foamVerts
811 << abort(FatalError);
812 }
813
814 // Bit nasty: replace vertex list.
815 cellVerts[celli].transfer(foamVerts);
816 }
817
818 // Renumber vertex numbers on boundaryFaces
819
820 forAll(boundaryFaces, bFacei)
821 {
822 labelList foamVerts(renumber(unvToFoam, boundaryFaces[bFacei]));
823
824 if (foamVerts.found(-1))
825 {
827 << "Boundary face " << bFacei
828 << " unv vertices " << boundaryFaces[bFacei]
829 << " has some undefined vertices " << foamVerts
830 << abort(FatalError);
831 }
832
833 // Bit nasty: replace vertex list.
834 boundaryFaces[bFacei].transfer(foamVerts);
835 }
836
837
838
839 // Patchify = create patchFaceVerts for use in cellShape construction
840 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
841
842 // Works in one of two modes. Either has read boundaryFaces which
843 // just need to be sorted according to patch. Or has read point constraint
844 // sets (dofVertIndices).
845
846 List<faceList> patchFaceVerts;
847
848 labelList own(boundaryFaces.size(), -1);
849 labelList nei(boundaryFaces.size(), -1);
851
852 {
853 // Can use face::symmHasher or use sorted indices instead
854 // - choose the latter in case UNV has anything odd
855 HashTable<label, face> faceToFaceID(2*boundaryFaces.size());
856
857 forAll(boundaryFaces, bfacei)
858 {
859 face sortedVerts(boundaryFaces[bfacei]);
860 Foam::sort(sortedVerts);
861 faceToFaceID.insert(sortedVerts, bfacei);
862 }
863
864 forAll(cellVerts, celli)
865 {
866 const cellShape& shape = cellVerts[celli];
867
868 for (const face& f : shape.faces())
869 {
870 face sortedVerts(f);
871 Foam::sort(sortedVerts);
872 const label bfacei = faceToFaceID.lookup(sortedVerts, -1);
873 if (bfacei != -1)
874 {
875 const int cmp = face::compare(f, boundaryFaces[bfacei]);
876
877 if (cmp == 1)
878 {
879 // Same orientation. Cell is owner.
880 own[bfacei] = celli;
881 }
882 else if (cmp == -1)
883 {
884 // Opposite orientation. Cell is neighbour.
885 nei[bfacei] = celli;
886 }
887 }
888 }
889 }
890
891 label nReverse = 0;
892 forAll(own, facei)
893 {
894 if (own[facei] == -1 && nei[facei] != -1)
895 {
896 // Boundary face with incorrect orientation
897 boundaryFaces[facei].flip();
898 std::swap(own[facei], nei[facei]);
899 nReverse++;
900 }
901 }
902 if (nReverse > 0)
903 {
904 Info << "Found " << nReverse << " reversed boundary faces out of "
905 << boundaryFaces.size() << endl;
906 }
907
908
909 label cnt = 0;
910 forAll(own, facei)
911 {
912 if (own[facei] != -1 && nei[facei] != -1)
913 {
914 faceToCell[1].insert(facei, own[facei]);
915 faceToCell[0].insert(facei, nei[facei]);
916 cnt++;
917 }
918 }
919
920 if (cnt > 0)
921 {
922 Info << "Of " << boundaryFaces.size() << " so-called"
923 << " boundary faces " << cnt << " belong to two cells "
924 << "and are therefore internal" << endl;
925 }
926 }
927
928 HashTable<labelList> cellZones;
929 HashTable<labelList> faceZones;
930 List<bool> isAPatch(patchNames.size(),true);
931
932 if (dofVertIndices.size())
933 {
934 // Use the vertex constraints to patch. Is of course bit dodgy since
935 // face goes on patch if all its vertices are on a constraint.
936 // Note: very inefficient since goes through all faces (including
937 // internal ones) twice. Should do a construct without patches
938 // and then repatchify.
939
940 Info<< "Using " << dofVertIndices.size()
941 << " DOF sets to detect boundary faces."<< endl;
942
943 // Renumber vertex numbers on constraints
944 forAll(dofVertIndices, patchi)
945 {
946 inplaceRenumber(unvToFoam, dofVertIndices[patchi]);
947 }
948
949
950 // Build labelHashSet of points per dofGroup/patch
951 List<labelHashSet> dofGroups(dofVertIndices.size());
952
953 forAll(dofVertIndices, patchi)
954 {
955 const labelList& foamVerts = dofVertIndices[patchi];
956 dofGroups[patchi].insert(foamVerts);
957 }
958
959 List<DynamicList<face>> dynPatchFaces(dofVertIndices.size());
960
961 forAll(cellVerts, celli)
962 {
963 const cellShape& shape = cellVerts[celli];
964
965 for (const face& f : shape.faces())
966 {
967 label patchi = findPatch(dofGroups, f);
968
969 if (patchi != -1)
970 {
971 dynPatchFaces[patchi].append(f);
972 }
973 }
974 }
975
976 // Transfer
977 patchFaceVerts.setSize(dynPatchFaces.size());
978
979 forAll(dynPatchFaces, patchi)
980 {
981 patchFaceVerts[patchi].transfer(dynPatchFaces[patchi]);
982 }
983 }
984 else
985 {
986 // Use the boundary faces.
987
988 // Construct the patch faces by sorting the boundaryFaces according to
989 // patch.
990 patchFaceVerts.setSize(patchFaceIndices.size());
991
992 Info<< "Sorting boundary faces according to group (patch)" << endl;
993
994 // make sure that no face is used twice on the boundary
995 // (possible for boundary-only faceZones)
996 labelHashSet alreadyOnBoundary;
997
998 // Construct map from boundaryFaceIndices
999 Map<label> boundaryFaceToIndex(boundaryFaceIndices.size());
1000
1001 forAll(boundaryFaceIndices, i)
1002 {
1003 boundaryFaceToIndex.insert(boundaryFaceIndices[i], i);
1004 }
1005
1006 forAll(patchFaceVerts, patchi)
1007 {
1008 Info << patchi << ": " << patchNames[patchi] << " is " << flush;
1009
1010 faceList& patchFaces = patchFaceVerts[patchi];
1011 const labelList& faceIndices = patchFaceIndices[patchi];
1012
1013 patchFaces.setSize(faceIndices.size());
1014
1015 bool duplicateFaces = false;
1016
1017 label cnt = 0;
1018 forAll(patchFaces, i)
1019 {
1020 if (boundaryFaceToIndex.found(faceIndices[i]))
1021 {
1022 label bFacei = boundaryFaceToIndex[faceIndices[i]];
1023
1024 if (own[bFacei] != -1 && nei[bFacei] == -1)
1025 {
1026 patchFaces[cnt] = boundaryFaces[bFacei];
1027 cnt++;
1028 if (alreadyOnBoundary.found(bFacei))
1029 {
1030 duplicateFaces = true;
1031 }
1032 }
1033 }
1034 }
1035
1036 if (cnt != patchFaces.size() || duplicateFaces)
1037 {
1038 isAPatch[patchi] = false;
1039
1040 if (verbose)
1041 {
1042 if (cnt != patchFaces.size())
1043 {
1045 << "For patch " << patchi << " there were "
1046 << patchFaces.size()-cnt
1047 << " faces not used because they seem"
1048 << " to be internal. "
1049 << "This seems to be a face or a cell-zone"
1050 << endl;
1051 }
1052 else
1053 {
1055 << "Patch "
1056 << patchi << " has faces that are already "
1057 << " in use on other boundary-patches,"
1058 << " Assuming faceZoneset." << endl;
1059 }
1060 }
1061
1062 patchFaces.setSize(0); // Assume that this is no patch at all
1063
1064 if (cellCorrespondence[faceIndices[0]] >= 0)
1065 {
1066 Info << "cellZone" << endl;
1067 labelList theCells(faceIndices.size());
1068 forAll(faceIndices, i)
1069 {
1070 if (cellCorrespondence[faceIndices[0]] < 0)
1071 {
1073 << "The face index " << faceIndices[i]
1074 << " was not found amongst the cells."
1075 << " This kills the theory that "
1076 << patchNames[patchi] << " is a cell zone"
1077 << endl
1078 << abort(FatalError);
1079 }
1080 theCells[i] = cellCorrespondence[faceIndices[i]];
1081 }
1082 cellZones.insert(patchNames[patchi], theCells);
1083 }
1084 else
1085 {
1086 Info << "faceZone" << endl;
1087 labelList theFaces(faceIndices.size());
1088 forAll(faceIndices, i)
1089 {
1090 theFaces[i] = boundaryFaceToIndex[faceIndices[i]];
1091 }
1092 faceZones.insert(patchNames[patchi],theFaces);
1093 }
1094 }
1095 else
1096 {
1097 Info << "patch" << endl;
1098
1099 forAll(patchFaces, i)
1100 {
1101 label bFacei = boundaryFaceToIndex[faceIndices[i]];
1102 alreadyOnBoundary.insert(bFacei);
1103 }
1104 }
1105 }
1106 }
1107
1108 pointField polyPoints;
1109 polyPoints.transfer(points);
1110
1111 // Length scaling factor
1112 polyPoints /= lengthScale;
1113
1114
1115 // For debugging: dump boundary faces as OBJ surface mesh
1116 if (args.found("dump"))
1117 {
1118 Info<< "Writing boundary faces to OBJ file boundaryFaces.obj"
1119 << nl << endl;
1120
1121 // Create globally numbered surface
1122 meshedSurface rawSurface(polyPoints, boundaryFaces);
1123
1124 // Write locally numbered surface
1126 (
1127 rawSurface.localPoints(),
1128 rawSurface.localFaces()
1129 ).write(runTime.path()/"boundaryFaces.obj");
1130 }
1131
1132
1133 Info<< "\nConstructing mesh with non-default patches of size:" << nl;
1134 DynamicList<word> usedPatchNames;
1135 DynamicList<faceList> usedPatchFaceVerts;
1136
1137 forAll(patchNames, patchi)
1138 {
1139 if (isAPatch[patchi])
1140 {
1141 Info<< " " << patchNames[patchi] << '\t'
1142 << patchFaceVerts[patchi].size() << nl;
1143 usedPatchNames.append(patchNames[patchi]);
1144 usedPatchFaceVerts.append(patchFaceVerts[patchi]);
1145 }
1146 }
1147 usedPatchNames.shrink();
1148 usedPatchFaceVerts.shrink();
1149
1150 Info<< endl;
1151
1152 // Construct mesh
1154 (
1155 IOobject
1156 (
1157 polyMesh::defaultRegion,
1158 runTime.constant(),
1159 runTime
1160 ),
1161 std::move(polyPoints),
1162 cellVerts,
1163 usedPatchFaceVerts, // boundaryFaces,
1164 usedPatchNames, // boundaryPatchNames,
1165 wordList(patchNames.size(), polyPatch::typeName), // boundaryPatchTypes,
1166 "defaultFaces", // defaultFacesName
1167 polyPatch::typeName, // defaultFacesType,
1168 wordList() // boundaryPatchPhysicalTypes
1169 );
1170
1171 // Remove files now, to ensure all mesh files written are consistent.
1172 mesh.removeFiles();
1173
1174 if (faceZones.size() || cellZones.size())
1175 {
1176 Info << "Adding cell and face zones" << endl;
1177
1179 List<faceZone*> fZones(faceZones.size());
1180 List<cellZone*> cZones(cellZones.size());
1181
1182 if (cellZones.size())
1183 {
1184 forAll(cellZones.toc(), cnt)
1185 {
1186 word name = cellZones.toc()[cnt];
1187 Info<< " Cell Zone " << name << " " << tab
1188 << cellZones[name].size() << endl;
1189
1190 cZones[cnt] = new cellZone
1191 (
1192 name,
1193 cellZones[name],
1194 cnt,
1195 mesh.cellZones()
1196 );
1197 }
1198 }
1199 if (faceZones.size())
1200 {
1201 const labelList& own = mesh.faceOwner();
1202 const labelList& nei = mesh.faceNeighbour();
1203 const pointField& centers = mesh.faceCentres();
1204 const pointField& points = mesh.points();
1205
1206 forAll(faceZones.toc(), cnt)
1207 {
1208 word name = faceZones.toc()[cnt];
1209 const labelList& oldIndizes = faceZones[name];
1210 labelList indizes(oldIndizes.size());
1211
1212 Info<< " Face Zone " << name << " " << tab
1213 << oldIndizes.size() << endl;
1214
1215 forAll(indizes, i)
1216 {
1217 const label old = oldIndizes[i];
1218 label noveau = -1;
1219 label c1 = -1, c2 = -1;
1220 if (faceToCell[0].found(old))
1221 {
1222 c1 = faceToCell[0][old];
1223 }
1224 if (faceToCell[1].found(old))
1225 {
1226 c2 = faceToCell[1][old];
1227 }
1228 if (c1 < c2)
1229 {
1230 label tmp = c1;
1231 c1 = c2;
1232 c2 = tmp;
1233 }
1234 if (c2 == -1)
1235 {
1236 // Boundary face is part of the faceZone
1237 forAll(own, j)
1238 {
1239 if (own[j] == c1)
1240 {
1241 const face& f = boundaryFaces[old];
1242 if (mag(centers[j]- f.centre(points)) < SMALL)
1243 {
1244 noveau = j;
1245 break;
1246 }
1247 }
1248 }
1249 }
1250 else
1251 {
1252 forAll(nei, j)
1253 {
1254 if
1255 (
1256 (c1 == own[j] && c2 == nei[j])
1257 || (c2 == own[j] && c1 == nei[j])
1258 )
1259 {
1260 noveau = j;
1261 break;
1262 }
1263 }
1264 }
1265 assert(noveau > -1);
1266 indizes[i] = noveau;
1267 }
1268 fZones[cnt] = new faceZone
1269 (
1270 faceZones.toc()[cnt],
1271 indizes,
1272 false, // none are flipped
1273 cnt,
1274 mesh.faceZones()
1275 );
1276 }
1277 }
1278 mesh.addZones(pZones, fZones, cZones);
1279
1280 Info << endl;
1281 }
1282
1283 // Set the precision of the points data to 10
1284 IOstream::defaultPrecision(max(10u, IOstream::defaultPrecision()));
1285
1286 mesh.write();
1287
1288 Info<< "End\n" << endl;
1289
1290 return 0;
1291}
1292
1293
1294// ************************************************************************* //
bool found
Y[inertIndex] max(0.0)
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects.
Definition: DynamicList.H:72
DynamicList< T, SizeMin > & shrink()
Shrink the allocated space to the number of elements used.
Definition: DynamicListI.H:434
void transfer(List< T > &list)
Transfer contents of the argument List into this.
Definition: DynamicListI.H:467
void setSize(const label n)
Same as resize()
Definition: DynamicList.H:221
void append(const T &val)
Copy append an element to the end of this list.
Definition: DynamicListI.H:503
bool insert(const Key &key)
Insert a new entry, not overwriting existing entries.
Definition: HashSet.H:191
A HashTable similar to std::unordered_map.
Definition: HashTable.H:123
List< Key > toc() const
The table of contents (the keys) in unsorted order.
Definition: HashTable.C:122
bool found(const Key &key) const
Return true if hashed entry is found in table.
Definition: HashTableI.H:100
bool insert(const Key &key, const T &obj)
Copy insert a new entry, not overwriting existing entries.
Definition: HashTableI.H:180
label size() const noexcept
The number of elements in table.
Definition: HashTableI.H:52
Input from file stream, using an ISstream.
Definition: IFstream.H:57
Defines the attributes of an object for which implicit objectRegistry management is supported,...
Definition: IOobject.H:170
label lineNumber() const noexcept
Const access to the current stream line number.
Definition: IOstream.H:318
bool good() const noexcept
True if next operation might succeed.
Definition: IOstream.H:233
ISstream & getLine(std::string &str, char delim='\n')
Raw, low-level getline (until delimiter) into a string.
Definition: ISstreamI.H:76
Input from string buffer, using a ISstream. Always UNCOMPRESSED.
Definition: StringStream.H:112
void transfer(List< T > &list)
Definition: List.C:447
void setSize(const label n)
Alias for resize()
Definition: List.H:218
void append(const T &val)
Append an element at the end of the list.
Definition: ListI.H:175
A HashTable to objects of type <T> with a label key.
Definition: Map.H:60
bool found(const T &val, label pos=0) const
True if the value if found in the list.
Definition: UListI.H:265
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
T & last()
Return the last element of the list.
Definition: UListI.H:216
T get(const label index) const
Get a value from the argument at index.
Definition: argListI.H:278
bool found(const word &optName) const
Return true if the named option is found.
Definition: argListI.H:178
Maps a geometry to a set of cell primitives.
Definition: cellModel.H:73
An analytical geometric cellShape.
Definition: cellShape.H:72
faceList faces() const
Faces of this cell.
Definition: cellShapeI.H:227
A subset of mesh cells.
Definition: cellZone.H:65
A topoSetCellSource to select all cells based on usage in given faceSet(s), e.g. select cells that ar...
Definition: faceToCell.H:182
A subset of mesh faces organised as a primitive patch.
Definition: faceZone.H:67
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
A class for handling file names.
Definition: fileName.H:76
A line primitive.
Definition: line.H:68
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:81
A class for managing temporary objects.
Definition: tmp.H:65
A class for handling words, derived from Foam::string.
Definition: word.H:68
const cellModel & hex
dynamicFvMesh & mesh
engineTime & runTime
IOporosityModelList pZones(mesh)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
const labelList nFaces(UPstream::listGatherValues< label >(aMesh.nFaces()))
const pointField & points
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.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))
#define IOWarningInFunction(ios)
Report an IO warning using Foam::Warning.
#define WarningInFunction
Report a warning using Foam::Warning.
void inplaceRenumber(const labelUList &oldToNew, IntListType &lists)
Inplace renumber the values (not the indices) of a list of lists.
Definition: ensightOutput.H:90
constexpr const char *const group
Group name for atomic constants.
const dimensionedScalar c2
Second radiation constant: default SI units: [m.K].
const dimensionedScalar c1
First radiation constant: default SI units: [W/m2].
Namespace for OpenFOAM.
void component(FieldField< Field, typename FieldField< Field, Type >::cmptType > &sf, const FieldField< Field, Type > &f, const direction d)
label readLabel(const char *buf)
Parse entire buffer as a label, skipping leading/trailing whitespace.
Definition: label.H:66
messageStream Info
Information stream (stdout output on master, null elsewhere)
IntListType renumber(const labelUList &oldToNew, const IntListType &input)
Renumber the values (not the indices) of a list.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
pointField vertices(const blockVertexList &bvl)
void sort(UList< T > &list)
Sort the list.
Definition: UList.C:342
errorManip< error > abort(error &err)
Definition: errorManip.H:144
labelList invert(const label len, const labelUList &map)
Create an inverse one-to-one mapping.
Definition: ListOps.C:36
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
Definition: exprTraits.C:59
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:130
Ostream & flush(Ostream &os)
Flush stream.
Definition: Ostream.H:364
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
constexpr char tab
The tab '\t' character(0x09)
Definition: Ostream.H:52
wordList patchNames(nPatches)
labelList f(nPoints)
Foam::argList args(argc, argv)
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333