boundaryMesh.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-2016 OpenFOAM Foundation
9 Copyright (C) 2020-2022 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
27\*---------------------------------------------------------------------------*/
28
29#include "boundaryMesh.H"
30#include "Time.H"
31#include "polyMesh.H"
33#include "faceList.H"
34#include "indexedOctree.H"
36#include "triSurface.H"
37#include "SortableList.H"
38#include "OFstream.H"
39#include "primitivePatch.H"
41
42// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
43
44namespace Foam
45{
47}
48
49// Normal along which to divide faces into categories (used in getNearest)
50const Foam::vector Foam::boundaryMesh::splitNormal_(3, 2, 1);
51
52// Distance to face tolerance for getNearest
53const Foam::scalar Foam::boundaryMesh::distanceTol_ = 1e-2;
54
55
56// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
57
58// Returns number of feature edges connected to pointi
59Foam::label Foam::boundaryMesh::nFeatureEdges(label pointi) const
60{
61 label nFeats = 0;
62
63 const labelList& pEdges = mesh().pointEdges()[pointi];
64
65 forAll(pEdges, pEdgeI)
66 {
67 label edgeI = pEdges[pEdgeI];
68
69 if (edgeToFeature_[edgeI] != -1)
70 {
71 nFeats++;
72 }
73 }
74 return nFeats;
75}
76
77
78// Returns next feature edge connected to pointi
79Foam::label Foam::boundaryMesh::nextFeatureEdge
80(
81 const label edgeI,
82 const label vertI
83) const
84{
85 const labelList& pEdges = mesh().pointEdges()[vertI];
86
87 forAll(pEdges, pEdgeI)
88 {
89 label nbrEdgeI = pEdges[pEdgeI];
90
91 if (nbrEdgeI != edgeI)
92 {
93 label featI = edgeToFeature_[nbrEdgeI];
94
95 if (featI != -1)
96 {
97 return nbrEdgeI;
98 }
99 }
100 }
101
102 return -1;
103}
104
105
106// Finds connected feature edges, starting from startPointi and returns
107// feature labels (not edge labels). Marks feature edges handled in
108// featVisited.
109Foam::labelList Foam::boundaryMesh::collectSegment
110(
111 const boolList& isFeaturePoint,
112 const label startEdgeI,
113 boolList& featVisited
114) const
115{
116 // Find starting feature point on edge.
117
118 label edgeI = startEdgeI;
119
120 const edge& e = mesh().edges()[edgeI];
121
122 label vertI = e.start();
123
124 while (!isFeaturePoint[vertI])
125 {
126 // Step to next feature edge
127
128 edgeI = nextFeatureEdge(edgeI, vertI);
129
130 if ((edgeI == -1) || (edgeI == startEdgeI))
131 {
132 break;
133 }
134
135 // Step to next vertex on edge
136
137 const edge& e = mesh().edges()[edgeI];
138
139 vertI = e.otherVertex(vertI);
140 }
141
142 //
143 // Now we have:
144 // edgeI : first edge on this segment
145 // vertI : one of the endpoints of this segment
146 //
147 // Start walking other way and storing edges as we go along.
148 //
149
150 // Untrimmed storage for current segment
151 labelList featLabels(featureEdges_.size());
152
153 label featLabelI = 0;
154
155 label initEdgeI = edgeI;
156
157 do
158 {
159 // Mark edge as visited
160 label featI = edgeToFeature_[edgeI];
161
162 if (featI == -1)
163 {
165 << "Problem" << abort(FatalError);
166 }
167 featLabels[featLabelI++] = featI;
168
169 featVisited[featI] = true;
170
171 // Step to next vertex on edge
172
173 const edge& e = mesh().edges()[edgeI];
174
175 vertI = e.otherVertex(vertI);
176
177 // Step to next feature edge
178
179 edgeI = nextFeatureEdge(edgeI, vertI);
180
181 if ((edgeI == -1) || (edgeI == initEdgeI))
182 {
183 break;
184 }
185 }
186 while (!isFeaturePoint[vertI]);
187
188
189 // Trim to size
190 featLabels.setSize(featLabelI);
191
192 return featLabels;
193}
194
195
196void Foam::boundaryMesh::markEdges
197(
198 const label maxDistance,
199 const label edgeI,
200 const label distance,
201 labelList& minDistance,
202 DynamicList<label>& visited
203) const
204{
205 if (distance < maxDistance)
206 {
207 // Don't do anything if reached beyond maxDistance.
208
209 if (minDistance[edgeI] == -1)
210 {
211 // First visit of edge. Store edge label.
212 visited.append(edgeI);
213 }
214 else if (minDistance[edgeI] <= distance)
215 {
216 // Already done this edge
217 return;
218 }
219
220 minDistance[edgeI] = distance;
221
222 const edge& e = mesh().edges()[edgeI];
223
224 // Do edges connected to e.start
225 const labelList& startEdges = mesh().pointEdges()[e.start()];
226
227 forAll(startEdges, pEdgeI)
228 {
229 markEdges
230 (
231 maxDistance,
232 startEdges[pEdgeI],
233 distance+1,
234 minDistance,
235 visited
236 );
237 }
238
239 // Do edges connected to e.end
240 const labelList& endEdges = mesh().pointEdges()[e.end()];
241
242 forAll(endEdges, pEdgeI)
243 {
244 markEdges
245 (
246 maxDistance,
247 endEdges[pEdgeI],
248 distance+1,
249 minDistance,
250 visited
251 );
252 }
253 }
254}
255
256
257Foam::label Foam::boundaryMesh::findPatchID
258(
259 const polyPatchList& patches,
260 const word& patchName
261) const
262{
263 forAll(patches, patchi)
264 {
265 if (patches[patchi].name() == patchName)
266 {
267 return patchi;
268 }
269 }
270
271 return -1;
272}
273
274
276{
277 wordList names(patches_.size());
278
279 forAll(patches_, patchi)
280 {
281 names[patchi] = patches_[patchi].name();
282 }
283 return names;
284}
285
286
287Foam::label Foam::boundaryMesh::whichPatch
288(
289 const polyPatchList& patches,
290 const label facei
291) const
292{
293 forAll(patches, patchi)
294 {
295 const polyPatch& pp = patches[patchi];
296
297 if ((facei >= pp.start()) && (facei < (pp.start() + pp.size())))
298 {
299 return patchi;
300 }
301 }
302 return -1;
303}
304
305
306// Gets labels of changed faces and propagates them to the edges. Returns
307// labels of edges changed.
309(
310 const boolList& regionEdge,
311 const label region,
312 const labelList& changedFaces,
313 labelList& edgeRegion
314) const
315{
316 labelList changedEdges(mesh().nEdges(), -1);
317 label changedI = 0;
318
319 forAll(changedFaces, i)
320 {
321 label facei = changedFaces[i];
322
323 const labelList& fEdges = mesh().faceEdges()[facei];
324
325 forAll(fEdges, fEdgeI)
326 {
327 label edgeI = fEdges[fEdgeI];
328
329 if (!regionEdge[edgeI] && (edgeRegion[edgeI] == -1))
330 {
331 edgeRegion[edgeI] = region;
332
333 changedEdges[changedI++] = edgeI;
334 }
335 }
336 }
337
338 changedEdges.setSize(changedI);
339
340 return changedEdges;
341}
342
343
344// Reverse of faceToEdge: gets edges and returns faces
346(
347 const label region,
348 const labelList& changedEdges,
349 labelList& faceRegion
350) const
351{
352 labelList changedFaces(mesh().size(), -1);
353 label changedI = 0;
354
355 forAll(changedEdges, i)
356 {
357 label edgeI = changedEdges[i];
358
359 const labelList& eFaces = mesh().edgeFaces()[edgeI];
360
361 forAll(eFaces, eFacei)
362 {
363 label facei = eFaces[eFacei];
364
365 if (faceRegion[facei] == -1)
366 {
367 faceRegion[facei] = region;
368
369 changedFaces[changedI++] = facei;
370 }
371 }
372 }
373
374 changedFaces.setSize(changedI);
375
376 return changedFaces;
377}
378
379
380// Finds area, starting at facei, delimited by borderEdge
381void Foam::boundaryMesh::markZone
382(
383 const boolList& borderEdge,
384 label facei,
385 label currentZone,
386 labelList& faceZone
387) const
388{
389 faceZone[facei] = currentZone;
390
391 // List of faces whose faceZone has been set.
392 labelList changedFaces(1, facei);
393 // List of edges whose faceZone has been set.
394 labelList changedEdges;
395
396 // Zones on all edges.
397 labelList edgeZone(mesh().nEdges(), -1);
398
399 while (true)
400 {
401 changedEdges = faceToEdge
402 (
403 borderEdge,
404 currentZone,
405 changedFaces,
406 edgeZone
407 );
408
409 if (debug)
410 {
411 Pout<< "From changedFaces:" << changedFaces.size()
412 << " to changedEdges:" << changedEdges.size()
413 << endl;
414 }
415
416 if (changedEdges.empty())
417 {
418 break;
419 }
420
421 changedFaces = edgeToFace(currentZone, changedEdges, faceZone);
422
423 if (debug)
424 {
425 Pout<< "From changedEdges:" << changedEdges.size()
426 << " to changedFaces:" << changedFaces.size()
427 << endl;
428 }
429
430 if (changedFaces.empty())
431 {
432 break;
433 }
434 }
435}
436
437
438// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
439
441:
442 meshPtr_(nullptr),
443 patches_(),
444 meshFace_(),
445 featurePoints_(),
446 featureEdges_(),
447 featureToEdge_(),
448 edgeToFeature_(),
449 featureSegments_(),
450 extraEdges_()
451{}
452
453
454// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
455
457{
458 meshPtr_.reset(nullptr);
459}
460
461
463{
464 patches_.clear();
465
466 patches_.setSize(mesh.boundaryMesh().size());
467
468 // Number of boundary faces
469 const label nBFaces = mesh.nBoundaryFaces();
470
471 faceList bFaces(nBFaces);
472
473 meshFace_.setSize(nBFaces);
474
475 label bFacei = 0;
476
477 // Collect all boundary faces.
478 forAll(mesh.boundaryMesh(), patchi)
479 {
480 const polyPatch& pp = mesh.boundaryMesh()[patchi];
481
482 patches_.set
483 (
484 patchi,
485 new boundaryPatch
486 (
487 pp.name(),
488 patchi,
489 pp.size(),
490 bFacei,
491 pp.type()
492 )
493 );
494
495 // Collect all faces in global numbering.
496 forAll(pp, patchFacei)
497 {
498 meshFace_[bFacei] = pp.start() + patchFacei;
499
500 bFaces[bFacei] = pp[patchFacei];
501
502 bFacei++;
503 }
504 }
505
506
507 if (debug)
508 {
509 Pout<< "read : patches now:" << endl;
510
511 forAll(patches_, patchi)
512 {
513 const boundaryPatch& bp = patches_[patchi];
514
515 Pout<< " name : " << bp.name() << endl
516 << " size : " << bp.size() << endl
517 << " start : " << bp.start() << endl
518 << " type : " << bp.physicalType() << endl
519 << endl;
520 }
521 }
522
523 //
524 // Construct single patch for all of boundary
525 //
526
527 // Temporary primitivePatch to calculate compact points & faces.
528 primitiveFacePatch globalPatch(bFaces, mesh.points());
529
530 // Store in local(compact) addressing
531 clearOut();
532
533 meshPtr_.reset
534 (
535 new bMesh(globalPatch.localFaces(), globalPatch.localPoints())
536 );
537
538 if (debug & 2)
539 {
540 const bMesh& msh = *meshPtr_;
541
542 Pout<< "** Start of Faces **" << endl;
543
544 forAll(msh, facei)
545 {
546 const face& f = msh[facei];
547
548 point ctr(Zero);
549
550 forAll(f, fp)
551 {
552 ctr += msh.points()[f[fp]];
553 }
554 ctr /= f.size();
555
556 Pout<< " " << facei
557 << " ctr:" << ctr
558 << " verts:" << f
559 << endl;
560 }
561
562 Pout<< "** End of Faces **" << endl;
563
564 Pout<< "** Start of Points **" << endl;
565
566 forAll(msh.points(), pointi)
567 {
568 Pout<< " " << pointi
569 << " coord:" << msh.points()[pointi]
570 << endl;
571 }
572
573 Pout<< "** End of Points **" << endl;
574 }
575
576 // Clear edge storage
577 featurePoints_.clear();
578 featureEdges_.clear();
579
580 featureToEdge_.clear();
581 edgeToFeature_.resize(meshPtr_->nEdges());
582 edgeToFeature_ = -1;
583
584 featureSegments_.clear();
585 extraEdges_.clear();
586}
587
588
590{
591 triSurface surf(fName);
592
593 if (surf.empty())
594 {
595 return;
596 }
597
598 // Sort according to region
599 SortableList<label> regions(surf.size());
600
601 forAll(surf, triI)
602 {
603 regions[triI] = surf[triI].region();
604 }
605 regions.sort();
606
607 // Determine region mapping.
608 Map<label> regionToBoundaryPatch;
609
610 label oldRegion = -1111;
611 label boundPatch = 0;
612
613 forAll(regions, i)
614 {
615 if (regions[i] != oldRegion)
616 {
617 regionToBoundaryPatch.insert(regions[i], boundPatch);
618
619 oldRegion = regions[i];
620 boundPatch++;
621 }
622 }
623
624 const geometricSurfacePatchList& surfPatches = surf.patches();
625
626 patches_.clear();
627
628 if (surfPatches.size() == regionToBoundaryPatch.size())
629 {
630 // There are as many surface patches as region numbers in triangles
631 // so use the surface patches
632
633 patches_.setSize(surfPatches.size());
634
635 // Take over patches, setting size to 0 for now.
636 forAll(surfPatches, patchi)
637 {
638 const geometricSurfacePatch& surfPatch = surfPatches[patchi];
639
640 patches_.set
641 (
642 patchi,
643 new boundaryPatch
644 (
645 surfPatch.name(),
646 patchi,
647 0,
648 0,
649 surfPatch.geometricType()
650 )
651 );
652 }
653 }
654 else
655 {
656 // There are not enough surface patches. Make up my own.
657
658 patches_.setSize(regionToBoundaryPatch.size());
659
660 forAll(patches_, patchi)
661 {
662 patches_.set
663 (
664 patchi,
665 new boundaryPatch
666 (
668 patchi,
669 0,
670 0,
672 )
673 );
674 }
675 }
676
677 //
678 // Copy according into bFaces according to regions
679 //
680
681 const labelList& indices = regions.indices();
682
683 faceList bFaces(surf.size());
684
685 meshFace_.setSize(surf.size());
686
687 label bFacei = 0;
688
689 // Current region number
690 label surfRegion = regions[0];
691 label foamRegion = regionToBoundaryPatch[surfRegion];
692
693 Pout<< "Surface region " << surfRegion << " becomes boundary patch "
694 << foamRegion << " with name " << patches_[foamRegion].name() << endl;
695
696
697 // Index in bFaces of start of current patch
698 label startFacei = 0;
699
700 forAll(indices, indexI)
701 {
702 label triI = indices[indexI];
703
704 const labelledTri& tri = surf.localFaces()[triI];
705
706 if (tri.region() != surfRegion)
707 {
708 // Change of region. We now know the size of the previous one.
709 boundaryPatch& bp = patches_[foamRegion];
710
711 bp.size() = bFacei - startFacei;
712 bp.start() = startFacei;
713
714 surfRegion = tri.region();
715 foamRegion = regionToBoundaryPatch[surfRegion];
716
717 Pout<< "Surface region " << surfRegion << " becomes boundary patch "
718 << foamRegion << " with name " << patches_[foamRegion].name()
719 << endl;
720
721 startFacei = bFacei;
722 }
723
724 meshFace_[bFacei] = triI;
725
726 bFaces[bFacei++] = face(tri);
727 }
728
729 // Final region
730 boundaryPatch& bp = patches_[foamRegion];
731
732 bp.size() = bFacei - startFacei;
733 bp.start() = startFacei;
734
735 //
736 // Construct single primitivePatch for all of boundary
737 //
738
739 clearOut();
740
741 // Store compact.
742 meshPtr_.reset(new bMesh(bFaces, surf.localPoints()));
743
744 // Clear edge storage
745 featurePoints_.clear();
746 featureEdges_.clear();
747
748 featureToEdge_.clear();
749 edgeToFeature_.resize(meshPtr_->nEdges());
750 edgeToFeature_ = -1;
751
752 featureSegments_.clear();
753 extraEdges_.clear();
754}
755
756
758{
759 geometricSurfacePatchList surfPatches(patches_.size());
760
761 forAll(patches_, patchi)
762 {
763 const boundaryPatch& bp = patches_[patchi];
764
765 surfPatches[patchi] =
767 (
768 bp.name(),
769 patchi,
770 bp.physicalType()
771 );
772 }
773
774 //
775 // Simple triangulation.
776 //
777
778 // Get number of triangles per face
779 labelList nTris(mesh().size());
780
781 label totalNTris = getNTris(0, mesh().size(), nTris);
782
783 // Determine per face the starting triangle.
784 labelList startTri(mesh().size());
785
786 label triI = 0;
787
788 forAll(mesh(), facei)
789 {
790 startTri[facei] = triI;
791
792 triI += nTris[facei];
793 }
794
795 // Triangulate
796 labelList triVerts(3*totalNTris);
797
798 triangulate(0, mesh().size(), totalNTris, triVerts);
799
800 // Convert to labelledTri
801
802 List<labelledTri> tris(totalNTris);
803
804 triI = 0;
805
806 forAll(patches_, patchi)
807 {
808 const boundaryPatch& bp = patches_[patchi];
809
810 forAll(bp, patchFacei)
811 {
812 label facei = bp.start() + patchFacei;
813
814 label triVertI = 3*startTri[facei];
815
816 for (label faceTriI = 0; faceTriI < nTris[facei]; faceTriI++)
817 {
818 label v0 = triVerts[triVertI++];
819 label v1 = triVerts[triVertI++];
820 label v2 = triVerts[triVertI++];
821
822 tris[triI++] = labelledTri(v0, v1, v2, patchi);
823 }
824 }
825 }
826
827 triSurface surf(tris, surfPatches, mesh().points());
828
829 OFstream surfStream(fName);
830
831 surf.write(surfStream);
832}
833
834
835// Get index in this (boundaryMesh) of face nearest to each boundary face in
836// pMesh.
837// Originally all triangles/faces of boundaryMesh would be bunged into
838// one big octree. Problem was that faces on top of each other, differing
839// only in sign of normal, could not be found separately. It would always
840// find only one. We could detect that it was probably finding the wrong one
841// (based on normal) but could not 'tell' the octree to retrieve the other
842// one (since they occupy exactly the same space)
843// So now faces get put into different octrees depending on normal.
844// !It still will not be possible to differentiate between two faces on top
845// of each other having the same normal
847(
848 const primitiveMesh& pMesh,
849 const vector& searchSpan
850) const
851{
852
853 // Divide faces into two bins acc. to normal
854 // - left of splitNormal
855 // - right ,,
856 DynamicList<label> leftFaces(mesh().size()/2);
857 DynamicList<label> rightFaces(mesh().size()/2);
858
859 forAll(mesh(), bFacei)
860 {
861 scalar sign = mesh().faceNormals()[bFacei] & splitNormal_;
862
863 if (sign > -1e-5)
864 {
865 rightFaces.append(bFacei);
866 }
867 if (sign < 1e-5)
868 {
869 leftFaces.append(bFacei);
870 }
871 }
872
873 leftFaces.shrink();
874 rightFaces.shrink();
875
876 if (debug)
877 {
878 Pout<< "getNearest :"
879 << " rightBin:" << rightFaces.size()
880 << " leftBin:" << leftFaces.size()
881 << endl;
882 }
883
885 (
886 UIndirectList<face>(mesh(), leftFaces),
887 mesh().points()
888 );
889 uindirectPrimitivePatch rightPatch
890 (
891 UIndirectList<face>(mesh(), rightFaces),
892 mesh().points()
893 );
894
895
896 // Overall bb
897 treeBoundBox overallBb(mesh().localPoints());
898
899 // Extend domain slightly (also makes it 3D if was 2D)
900 // Note asymmetry to avoid having faces align with octree cubes.
901 scalar tol = 1e-6 * overallBb.avgDim();
902
903 point& bbMin = overallBb.min();
904 bbMin.x() -= tol;
905 bbMin.y() -= tol;
906 bbMin.z() -= tol;
907
908 point& bbMax = overallBb.max();
909 bbMax.x() += 2*tol;
910 bbMax.y() += 2*tol;
911 bbMax.z() += 2*tol;
912
913 const scalar planarTol =
915 perturbTol();
916
917
918 // Create the octrees
920 <
922 > leftTree
923 (
925 (
926 false, // cacheBb
927 leftPatch,
928 planarTol
929 ),
930 overallBb,
931 10, // maxLevel
932 10, // leafSize
933 3.0 // duplicity
934 );
936 <
938 > rightTree
939 (
941 (
942 false, // cacheBb
943 rightPatch,
944 planarTol
945 ),
946 overallBb,
947 10, // maxLevel
948 10, // leafSize
949 3.0 // duplicity
950 );
951
952 if (debug)
953 {
954 Pout<< "getNearest : built trees" << endl;
955 }
956
957
958 const vectorField& ns = mesh().faceNormals();
959
960
961 //
962 // Search nearest triangle centre for every polyMesh boundary face
963 //
964
965 labelList nearestBFacei(pMesh.nBoundaryFaces());
966
967 treeBoundBox tightest;
968
969 const scalar searchDimSqr = magSqr(searchSpan);
970
971 forAll(nearestBFacei, patchFacei)
972 {
973 label meshFacei = pMesh.nInternalFaces() + patchFacei;
974
975 const point& ctr = pMesh.faceCentres()[meshFacei];
976
977 if (debug && (patchFacei % 1000) == 0)
978 {
979 Pout<< "getNearest : patchFace:" << patchFacei
980 << " meshFacei:" << meshFacei << " ctr:" << ctr << endl;
981 }
982
983
984 // Get normal from area vector
985 vector n = pMesh.faceAreas()[meshFacei];
986 scalar area = mag(n);
987 n /= area;
988
989 scalar typDim = -GREAT;
990 const face& f = pMesh.faces()[meshFacei];
991
992 forAll(f, fp)
993 {
994 typDim = max(typDim, mag(pMesh.points()[f[fp]] - ctr));
995 }
996
997 // Search right tree
998 pointIndexHit rightInfo = rightTree.findNearest(ctr, searchDimSqr);
999
1000 // Search left tree. Note: could start from rightDist bounding box
1001 // instead of starting from top.
1002 pointIndexHit leftInfo = leftTree.findNearest(ctr, searchDimSqr);
1003
1004 if (rightInfo.hit())
1005 {
1006 if (leftInfo.hit())
1007 {
1008 // Found in both trees. Compare normals.
1009 label rightFacei = rightFaces[rightInfo.index()];
1010 label leftFacei = leftFaces[leftInfo.index()];
1011
1012 label rightDist = mag(rightInfo.hitPoint()-ctr);
1013 label leftDist = mag(leftInfo.hitPoint()-ctr);
1014
1015 scalar rightSign = n & ns[rightFacei];
1016 scalar leftSign = n & ns[leftFacei];
1017
1018 if
1019 (
1020 (rightSign > 0 && leftSign > 0)
1021 || (rightSign < 0 && leftSign < 0)
1022 )
1023 {
1024 // Both same sign. Choose nearest.
1025 if (rightDist < leftDist)
1026 {
1027 nearestBFacei[patchFacei] = rightFacei;
1028 }
1029 else
1030 {
1031 nearestBFacei[patchFacei] = leftFacei;
1032 }
1033 }
1034 else
1035 {
1036 // Differing sign.
1037 // - if both near enough choose one with correct sign
1038 // - otherwise choose nearest.
1039
1040 // Get local dimension as max of distance between ctr and
1041 // any face vertex.
1042
1043 typDim *= distanceTol_;
1044
1045 if (rightDist < typDim && leftDist < typDim)
1046 {
1047 // Different sign and nearby. Choosing matching normal
1048 if (rightSign > 0)
1049 {
1050 nearestBFacei[patchFacei] = rightFacei;
1051 }
1052 else
1053 {
1054 nearestBFacei[patchFacei] = leftFacei;
1055 }
1056 }
1057 else
1058 {
1059 // Different sign but faraway. Choosing nearest.
1060 if (rightDist < leftDist)
1061 {
1062 nearestBFacei[patchFacei] = rightFacei;
1063 }
1064 else
1065 {
1066 nearestBFacei[patchFacei] = leftFacei;
1067 }
1068 }
1069 }
1070 }
1071 else
1072 {
1073 // Found in right but not in left. Choose right regardless if
1074 // correct sign. Note: do we want this?
1075 label rightFacei = rightFaces[rightInfo.index()];
1076 nearestBFacei[patchFacei] = rightFacei;
1077 }
1078 }
1079 else
1080 {
1081 // No face found in right tree.
1082
1083 if (leftInfo.hit())
1084 {
1085 // Found in left but not in right. Choose left regardless if
1086 // correct sign. Note: do we want this?
1087 nearestBFacei[patchFacei] = leftFaces[leftInfo.index()];
1088 }
1089 else
1090 {
1091 // No face found in left tree.
1092 nearestBFacei[patchFacei] = -1;
1093 }
1094 }
1095 }
1096
1097 return nearestBFacei;
1098}
1099
1100
1102(
1103 const labelList& nearest,
1104 const polyBoundaryMesh& oldPatches,
1105 polyMesh& newMesh
1106) const
1107{
1108
1109 // 2 cases to be handled:
1110 // A- patches in boundaryMesh patches_
1111 // B- patches not in boundaryMesh patches_ but in polyMesh
1112
1113 // Create maps from patch name to new patch index.
1114 HashTable<label> nameToIndex(2*patches_.size());
1115
1116 Map<word> indexToName(2*patches_.size());
1117
1118
1119 label nNewPatches = patches_.size();
1120
1121 forAll(oldPatches, oldPatchi)
1122 {
1123 const polyPatch& patch = oldPatches[oldPatchi];
1124 const label newPatchi = findPatchID(patch.name());
1125
1126 if (newPatchi != -1)
1127 {
1128 nameToIndex.insert(patch.name(), newPatchi);
1129 indexToName.insert(newPatchi, patch.name());
1130 }
1131 }
1132
1133 // Include all boundaryPatches not yet in nameToIndex (i.e. not in old
1134 // patches)
1135 forAll(patches_, bPatchi)
1136 {
1137 const boundaryPatch& bp = patches_[bPatchi];
1138
1139 if (!nameToIndex.found(bp.name()))
1140 {
1141 nameToIndex.insert(bp.name(), bPatchi);
1142 indexToName.insert(bPatchi, bp.name());
1143 }
1144 }
1145
1146 // Pass1:
1147 // Copy names&type of patches (with zero size) from old mesh as far as
1148 // possible. First patch created gets all boundary faces; all others get
1149 // zero faces (repatched later on). Exception is coupled patches which
1150 // keep their size.
1151
1152 List<polyPatch*> newPatchPtrList(nNewPatches);
1153
1154 label meshFacei = newMesh.nInternalFaces();
1155
1156 // First patch gets all non-coupled faces
1157 label facesToBeDone = newMesh.nBoundaryFaces();
1158
1159 forAll(patches_, bPatchi)
1160 {
1161 const boundaryPatch& bp = patches_[bPatchi];
1162
1163 const label newPatchi = nameToIndex[bp.name()];
1164
1165 // Find corresponding patch in polyMesh
1166 const label oldPatchi = findPatchID(oldPatches, bp.name());
1167
1168 if (oldPatchi == -1)
1169 {
1170 // Newly created patch. Gets all or zero faces.
1171 if (debug)
1172 {
1173 Pout<< "patchify : Creating new polyPatch:" << bp.name()
1174 << " type:" << bp.physicalType() << endl;
1175 }
1176
1177 newPatchPtrList[newPatchi] = polyPatch::New
1178 (
1179 bp.physicalType(),
1180 bp.name(),
1181 facesToBeDone,
1182 meshFacei,
1183 newPatchi,
1184 newMesh.boundaryMesh()
1185 ).ptr();
1186
1187 meshFacei += facesToBeDone;
1188
1189 // first patch gets all boundary faces; all others get 0.
1190 facesToBeDone = 0;
1191 }
1192 else
1193 {
1194 // Existing patch. Gets all or zero faces.
1195 const polyPatch& oldPatch = oldPatches[oldPatchi];
1196
1197 if (debug)
1198 {
1199 Pout<< "patchify : Cloning existing polyPatch:"
1200 << oldPatch.name() << endl;
1201 }
1202
1203 newPatchPtrList[newPatchi] = oldPatch.clone
1204 (
1205 newMesh.boundaryMesh(),
1206 newPatchi,
1207 facesToBeDone,
1208 meshFacei
1209 ).ptr();
1210
1211 meshFacei += facesToBeDone;
1212
1213 // first patch gets all boundary faces; all others get 0.
1214 facesToBeDone = 0;
1215 }
1216 }
1217
1218
1219 if (debug)
1220 {
1221 Pout<< "Patchify : new polyPatch list:" << endl;
1222
1223 forAll(newPatchPtrList, patchi)
1224 {
1225 const polyPatch& newPatch = *newPatchPtrList[patchi];
1226
1227 if (debug)
1228 {
1229 Pout<< "polyPatch:" << newPatch.name() << endl
1230 << " type :" << newPatch.typeName << endl
1231 << " size :" << newPatch.size() << endl
1232 << " start:" << newPatch.start() << endl
1233 << " index:" << patchi << endl;
1234 }
1235 }
1236 }
1237
1238 // Actually add new list of patches
1239 repatchPolyTopoChanger polyMeshRepatcher(newMesh);
1240 polyMeshRepatcher.changePatches(newPatchPtrList);
1241
1242
1243 // Pass2:
1244 // Change patch type for face
1245
1246 if (newPatchPtrList.size())
1247 {
1248 List<DynamicList<label>> patchFaces(nNewPatches);
1249
1250 // Give reasonable estimate for size of patches
1251 label nAvgFaces = newMesh.nBoundaryFaces() / nNewPatches;
1252
1253 forAll(patchFaces, newPatchi)
1254 {
1255 patchFaces[newPatchi].setCapacity(nAvgFaces);
1256 }
1257
1258 //
1259 // Sort faces acc. to new patch index. Can loop over all old patches
1260 // since will contain all faces.
1261 //
1262
1263 forAll(oldPatches, oldPatchi)
1264 {
1265 const polyPatch& patch = oldPatches[oldPatchi];
1266
1267 forAll(patch, patchFacei)
1268 {
1269 // Put face into region given by nearest boundary face
1270
1271 label meshFacei = patch.start() + patchFacei;
1272
1273 label bFacei = meshFacei - newMesh.nInternalFaces();
1274
1275 patchFaces[whichPatch(nearest[bFacei])].append(meshFacei);
1276 }
1277 }
1278
1279 forAll(patchFaces, newPatchi)
1280 {
1281 patchFaces[newPatchi].shrink();
1282 }
1283
1284
1285 // Change patch > 0. (since above we put all faces into the zeroth
1286 // patch)
1287
1288 for (label newPatchi = 1; newPatchi < patchFaces.size(); newPatchi++)
1289 {
1290 const labelList& pFaces = patchFaces[newPatchi];
1291
1292 forAll(pFaces, pFacei)
1293 {
1294 polyMeshRepatcher.changePatchID(pFaces[pFacei], newPatchi);
1295 }
1296 }
1297
1298 polyMeshRepatcher.repatch();
1299 }
1300}
1301
1302
1304{
1305 edgeToFeature_.setSize(mesh().nEdges());
1306
1307 edgeToFeature_ = -1;
1308
1309 // 1. Mark feature edges
1310
1311 // Storage for edge labels that are features. Trim later.
1312 featureToEdge_.setSize(mesh().nEdges());
1313
1314 label featureI = 0;
1315
1316 if (minCos >= 0.9999)
1317 {
1318 // Select everything
1319 forAll(mesh().edges(), edgeI)
1320 {
1321 edgeToFeature_[edgeI] = featureI;
1322 featureToEdge_[featureI++] = edgeI;
1323 }
1324 }
1325 else
1326 {
1327 forAll(mesh().edges(), edgeI)
1328 {
1329 const labelList& eFaces = mesh().edgeFaces()[edgeI];
1330
1331 if (eFaces.size() == 2)
1332 {
1333 label face0I = eFaces[0];
1334
1335 label face1I = eFaces[1];
1336
1339 //if (whichPatch(face0I) != whichPatch(face1I))
1340 //{
1341 // edgeToFeature_[edgeI] = featureI;
1342 // featureToEdge_[featureI++] = edgeI;
1343 //}
1344 //else
1345 {
1346 const vector& n0 = mesh().faceNormals()[face0I];
1347
1348 const vector& n1 = mesh().faceNormals()[face1I];
1349
1350 float cosAng = n0 & n1;
1351
1352 if (cosAng < minCos)
1353 {
1354 edgeToFeature_[edgeI] = featureI;
1355 featureToEdge_[featureI++] = edgeI;
1356 }
1357 }
1358 }
1359 else
1360 {
1361 //Should not occur: 0 or more than two faces
1362 edgeToFeature_[edgeI] = featureI;
1363 featureToEdge_[featureI++] = edgeI;
1364 }
1365 }
1366 }
1367
1368 // Trim featureToEdge_ to actual number of edges.
1369 featureToEdge_.setSize(featureI);
1370
1371 //
1372 // Compact edges i.e. relabel vertices.
1373 //
1374
1375 featureEdges_.setSize(featureI);
1376 featurePoints_.setSize(mesh().nPoints());
1377
1378 labelList featToMeshPoint(mesh().nPoints(), -1);
1379
1380 label featPtI = 0;
1381
1382 forAll(featureToEdge_, fEdgeI)
1383 {
1384 label edgeI = featureToEdge_[fEdgeI];
1385
1386 const edge& e = mesh().edges()[edgeI];
1387
1388 label start = featToMeshPoint[e.start()];
1389
1390 if (start == -1)
1391 {
1392 featToMeshPoint[e.start()] = featPtI;
1393
1394 featurePoints_[featPtI] = mesh().points()[e.start()];
1395
1396 start = featPtI;
1397
1398 featPtI++;
1399 }
1400
1401 label end = featToMeshPoint[e.end()];
1402
1403 if (end == -1)
1404 {
1405 featToMeshPoint[e.end()] = featPtI;
1406
1407 featurePoints_[featPtI] = mesh().points()[e.end()];
1408
1409 end = featPtI;
1410
1411 featPtI++;
1412 }
1413
1414 // Store with renumbered vertices.
1415 featureEdges_[fEdgeI] = edge(start, end);
1416 }
1417
1418 // Compact points
1419 featurePoints_.setSize(featPtI);
1420
1421
1422 //
1423 // 2. Mark endpoints of feature segments. These are points with
1424 // != 2 feature edges connected.
1425 // Note: can add geometric constraint here as well that if 2 feature
1426 // edges the angle between them should be less than xxx.
1427 //
1428
1429 boolList isFeaturePoint(mesh().nPoints(), false);
1430
1431 forAll(featureToEdge_, featI)
1432 {
1433 label edgeI = featureToEdge_[featI];
1434
1435 const edge& e = mesh().edges()[edgeI];
1436
1437 if (nFeatureEdges(e.start()) != 2)
1438 {
1439 isFeaturePoint[e.start()] = true;
1440 }
1441
1442 if (nFeatureEdges(e.end()) != 2)
1443 {
1444 isFeaturePoint[e.end()] = true;
1445 }
1446 }
1447
1448
1449 //
1450 // 3: Split feature edges into segments:
1451 // find point with not 2 feature edges -> start of feature segment
1452 //
1453
1454 DynamicList<labelList> segments;
1455
1456
1457 boolList featVisited(featureToEdge_.size(), false);
1458
1459 do
1460 {
1461 label startFeatI = -1;
1462
1463 forAll(featVisited, featI)
1464 {
1465 if (!featVisited[featI])
1466 {
1467 startFeatI = featI;
1468
1469 break;
1470 }
1471 }
1472
1473 if (startFeatI == -1)
1474 {
1475 // No feature lines left.
1476 break;
1477 }
1478
1479 segments.append
1480 (
1481 collectSegment
1482 (
1483 isFeaturePoint,
1484 featureToEdge_[startFeatI],
1485 featVisited
1486 )
1487 );
1488 }
1489 while (true);
1490
1491
1492 //
1493 // Store in *this
1494 //
1495 featureSegments_.setSize(segments.size());
1496
1497 forAll(featureSegments_, segmentI)
1498 {
1499 featureSegments_[segmentI] = segments[segmentI];
1500 }
1501}
1502
1503
1505{
1506 labelList minDistance(mesh().nEdges(), -1);
1507
1508 // All edge labels encountered
1509 DynamicList<label> visitedEdges;
1510
1511 // Floodfill from edgeI starting from distance 0. Stop at distance.
1512 markEdges(8, edgeI, 0, minDistance, visitedEdges);
1513
1514 // Set edge labels to display
1515 extraEdges_.transfer(visitedEdges);
1516}
1517
1518
1519Foam::label Foam::boundaryMesh::whichPatch(const label facei) const
1520{
1521 forAll(patches_, patchi)
1522 {
1523 const boundaryPatch& bp = patches_[patchi];
1524
1525 if ((facei >= bp.start()) && (facei < (bp.start() + bp.size())))
1526 {
1527 return patchi;
1528 }
1529 }
1530
1532 << "Cannot find face " << facei << " in list of boundaryPatches "
1533 << patches_
1534 << abort(FatalError);
1535
1536 return -1;
1537}
1538
1539
1540Foam::label Foam::boundaryMesh::findPatchID(const word& patchName) const
1541{
1542 forAll(patches_, patchi)
1543 {
1544 if (patches_[patchi].name() == patchName)
1545 {
1546 return patchi;
1547 }
1548 }
1549
1550 return -1;
1551}
1552
1553
1555{
1556 patches_.setSize(patches_.size() + 1);
1557
1558 // Add empty patch at end of patch list.
1559
1560 label patchi = patches_.size()-1;
1561
1562 boundaryPatch* bpPtr = new boundaryPatch
1563 (
1564 patchName,
1565 patchi,
1566 0,
1567 mesh().size(),
1568 "empty"
1569 );
1570
1571 patches_.set(patchi, bpPtr);
1572
1573 if (debug)
1574 {
1575 Pout<< "addPatch : patches now:" << endl;
1576
1577 forAll(patches_, patchi)
1578 {
1579 const boundaryPatch& bp = patches_[patchi];
1580
1581 Pout<< " name : " << bp.name() << endl
1582 << " size : " << bp.size() << endl
1583 << " start : " << bp.start() << endl
1584 << " type : " << bp.physicalType() << endl
1585 << endl;
1586 }
1587 }
1588}
1589
1590
1592{
1593 const label delPatchi = findPatchID(patchName);
1594
1595 if (delPatchi == -1)
1596 {
1598 << "Can't find patch named " << patchName
1599 << abort(FatalError);
1600 }
1601
1602 if (patches_[delPatchi].size())
1603 {
1605 << "Trying to delete non-empty patch " << patchName
1606 << endl << "Current size:" << patches_[delPatchi].size()
1607 << abort(FatalError);
1608 }
1609
1610 PtrList<boundaryPatch> newPatches(patches_.size() - 1);
1611
1612 for (label patchi = 0; patchi < delPatchi; patchi++)
1613 {
1614 newPatches.set(patchi, patches_[patchi].clone());
1615 }
1616
1617 // Move patches down, starting from delPatchi.
1618
1619 for (label patchi = delPatchi + 1; patchi < patches_.size(); patchi++)
1620 {
1621 newPatches.set(patchi - 1, patches_[patchi].clone());
1622 }
1623
1624 patches_.clear();
1625
1626 patches_ = newPatches;
1627
1628 if (debug)
1629 {
1630 Pout<< "deletePatch : patches now:" << endl;
1631
1632 forAll(patches_, patchi)
1633 {
1634 const boundaryPatch& bp = patches_[patchi];
1635
1636 Pout<< " name : " << bp.name() << endl
1637 << " size : " << bp.size() << endl
1638 << " start : " << bp.start() << endl
1639 << " type : " << bp.physicalType() << endl
1640 << endl;
1641 }
1642 }
1643}
1644
1645
1647(
1648 const word& patchName,
1649 const word& patchType
1650)
1651{
1652 const label changeI = findPatchID(patchName);
1653
1654 if (changeI == -1)
1655 {
1657 << "Can't find patch named " << patchName
1658 << abort(FatalError);
1659 }
1660
1661
1662 // Cause we can't reassign to individual PtrList elems ;-(
1663 // work on copy
1664
1665
1666 PtrList<boundaryPatch> newPatches(patches_.size());
1667
1668 forAll(patches_, patchi)
1669 {
1670 if (patchi == changeI)
1671 {
1672 // Create copy but for type
1673 const boundaryPatch& oldBp = patches_[patchi];
1674
1675 boundaryPatch* bpPtr = new boundaryPatch
1676 (
1677 oldBp.name(),
1678 oldBp.index(),
1679 oldBp.size(),
1680 oldBp.start(),
1681 patchType
1682 );
1683
1684 newPatches.set(patchi, bpPtr);
1685 }
1686 else
1687 {
1688 // Create copy
1689 newPatches.set(patchi, patches_[patchi].clone());
1690 }
1691 }
1692
1693 patches_ = newPatches;
1694}
1695
1696
1698(
1699 const labelList& patchIDs,
1700 labelList& oldToNew
1701)
1702{
1703 if (patchIDs.size() != mesh().size())
1704 {
1706 << "List of patchIDs not equal to number of faces." << endl
1707 << "PatchIDs size:" << patchIDs.size()
1708 << " nFaces::" << mesh().size()
1709 << abort(FatalError);
1710 }
1711
1712 // Count number of faces for each patch
1713
1714 labelList nFaces(patches_.size(), Zero);
1715
1716 forAll(patchIDs, facei)
1717 {
1718 label patchID = patchIDs[facei];
1719
1720 if (patchID < 0 || patchID >= patches_.size())
1721 {
1723 << "PatchID " << patchID << " out of range"
1724 << abort(FatalError);
1725 }
1726 nFaces[patchID]++;
1727 }
1728
1729
1730 // Determine position in faces_ for each patch
1731
1732 labelList startFace(patches_.size());
1733
1734 startFace[0] = 0;
1735
1736 for (label patchi = 1; patchi < patches_.size(); patchi++)
1737 {
1738 startFace[patchi] = startFace[patchi-1] + nFaces[patchi-1];
1739 }
1740
1741 // Update patch info
1742 PtrList<boundaryPatch> newPatches(patches_.size());
1743
1744 forAll(patches_, patchi)
1745 {
1746 const boundaryPatch& bp = patches_[patchi];
1747
1748 newPatches.set
1749 (
1750 patchi,
1751 new boundaryPatch
1752 (
1753 bp.name(),
1754 patchi,
1755 nFaces[patchi],
1756 startFace[patchi],
1757 bp.physicalType()
1758 )
1759 );
1760 }
1761 patches_ = newPatches;
1762
1763 if (debug)
1764 {
1765 Pout<< "changeFaces : patches now:" << endl;
1766
1767 forAll(patches_, patchi)
1768 {
1769 const boundaryPatch& bp = patches_[patchi];
1770
1771 Pout<< " name : " << bp.name() << endl
1772 << " size : " << bp.size() << endl
1773 << " start : " << bp.start() << endl
1774 << " type : " << bp.physicalType() << endl
1775 << endl;
1776 }
1777 }
1778
1779
1780 // Construct face mapping array
1781 oldToNew.setSize(patchIDs.size());
1782
1783 forAll(patchIDs, facei)
1784 {
1785 int patchID = patchIDs[facei];
1786
1787 oldToNew[facei] = startFace[patchID]++;
1788 }
1789
1790 // Copy faces into correct position and maintain label of original face
1791 faceList newFaces(mesh().size());
1792
1793 labelList newMeshFace(mesh().size());
1794
1795 forAll(oldToNew, facei)
1796 {
1797 newFaces[oldToNew[facei]] = mesh()[facei];
1798 newMeshFace[oldToNew[facei]] = meshFace_[facei];
1799 }
1800
1801 // Reconstruct 'mesh' from new faces and (copy of) existing points.
1802
1803 unique_ptr<bMesh> newMeshPtr(new bMesh(newFaces, mesh().points()));
1804
1805 // Reset meshFace_ to new ordering.
1806 meshFace_.transfer(newMeshFace);
1807
1808
1809 // Remove old PrimitivePatch on meshPtr_.
1810 clearOut();
1811
1812 // And insert new 'mesh'.
1813 meshPtr_ = std::move(newMeshPtr);
1814}
1815
1816
1817Foam::label Foam::boundaryMesh::getNTris(const label facei) const
1818{
1819 const face& f = mesh()[facei];
1820
1821 return f.nTriangles(mesh().points());
1822}
1823
1824
1826(
1827 const label startFacei,
1828 const label nFaces,
1829 labelList& nTris
1830) const
1831{
1832 label totalNTris = 0;
1833
1834 nTris.setSize(nFaces);
1835
1836 for (label i = 0; i < nFaces; i++)
1837 {
1838 label faceNTris = getNTris(startFacei + i);
1839
1840 nTris[i] = faceNTris;
1841
1842 totalNTris += faceNTris;
1843 }
1844 return totalNTris;
1845}
1846
1847
1848// Simple triangulation of face subset. Stores vertices in tris[] as three
1849// consecutive vertices per triangle.
1851(
1852 const label startFacei,
1853 const label nFaces,
1854 const label totalNTris,
1855 labelList& triVerts
1856) const
1857{
1858 // Triangulate faces.
1859 triVerts.setSize(3*totalNTris);
1860
1861 label vertI = 0;
1862
1863 for (label i = 0; i < nFaces; i++)
1864 {
1865 label facei = startFacei + i;
1866
1867 const face& f = mesh()[facei];
1868
1869 // Have face triangulate itself (results in faceList)
1870 faceList triFaces(f.nTriangles(mesh().points()));
1871
1872 label nTri = 0;
1873
1874 f.triangles(mesh().points(), nTri, triFaces);
1875
1876 // Copy into triVerts
1877
1878 forAll(triFaces, triFacei)
1879 {
1880 const face& triF = triFaces[triFacei];
1881
1882 triVerts[vertI++] = triF[0];
1883 triVerts[vertI++] = triF[1];
1884 triVerts[vertI++] = triF[2];
1885 }
1886 }
1887}
1888
1889
1890// Number of local points in subset.
1892(
1893 const label startFacei,
1894 const label nFaces
1895) const
1896{
1897 primitivePatch patch
1898 (
1899 SubList<face>(mesh(), nFaces, startFacei),
1900 mesh().points()
1901 );
1902
1903 return patch.nPoints();
1904}
1905
1906
1907// Triangulation of face subset in local coords.
1909(
1910 const label startFacei,
1911 const label nFaces,
1912 const label totalNTris,
1913 labelList& triVerts,
1914 labelList& localToGlobal
1915) const
1916{
1917 primitivePatch patch
1918 (
1919 SubList<face>(mesh(), nFaces, startFacei),
1920 mesh().points()
1921 );
1922
1923 // Map from local to mesh().points() addressing
1924 localToGlobal = patch.meshPoints();
1925
1926 // Triangulate patch faces.
1927 triVerts.setSize(3*totalNTris);
1928
1929 label vertI = 0;
1930
1931 for (label i = 0; i < nFaces; i++)
1932 {
1933 // Local face
1934 const face& f = patch.localFaces()[i];
1935
1936 // Have face triangulate itself (results in faceList)
1937 faceList triFaces(f.nTriangles(patch.localPoints()));
1938
1939 label nTri = 0;
1940
1941 f.triangles(patch.localPoints(), nTri, triFaces);
1942
1943 // Copy into triVerts
1944
1945 forAll(triFaces, triFacei)
1946 {
1947 const face& triF = triFaces[triFacei];
1948
1949 triVerts[vertI++] = triF[0];
1950 triVerts[vertI++] = triF[1];
1951 triVerts[vertI++] = triF[2];
1952 }
1953 }
1954}
1955
1956
1958(
1959 const labelList& protectedEdges,
1960 const label seedFacei,
1961 boolList& visited
1962) const
1963{
1964 boolList protectedEdge(mesh().nEdges(), false);
1965
1966 forAll(protectedEdges, i)
1967 {
1968 protectedEdge[protectedEdges[i]] = true;
1969 }
1970
1971
1972 // Initialize zone for all faces to -1
1973 labelList currentZone(mesh().size(), -1);
1974
1975 // Mark with 0 all faces reachable from seedFacei
1976 markZone(protectedEdge, seedFacei, 0, currentZone);
1977
1978 // Set in visited all reached ones.
1979 visited.setSize(mesh().size());
1980
1981 forAll(currentZone, facei)
1982 {
1983 if (currentZone[facei] == 0)
1984 {
1985 visited[facei] = true;
1986 }
1987 else
1988 {
1989 visited[facei] = false;
1990 }
1991 }
1992}
1993
1994
1995// ************************************************************************* //
label n
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 append(const T &val)
Copy append an element to the end of this list.
Definition: DynamicListI.H:503
A HashTable similar to std::unordered_map.
Definition: HashTable.H:123
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
std::enable_if< std::is_same< bool, TypeT >::value, bool >::type set(const label i, bool val=true)
A bitSet::set() method for a list of bool.
Definition: List.H:330
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
void clear()
Clear the list, i.e. set size to zero.
Definition: ListI.H:116
A HashTable to objects of type <T> with a label key.
Definition: Map.H:60
virtual label triangulate()
Triangulate in-place, returning the number of triangles added.
Output to file stream, using an OSstream.
Definition: OFstream.H:57
virtual const fileName & name() const
Get the name of the stream.
Definition: OSstream.H:107
label edgeToFace()
Propagate from edge to face.
label faceToEdge()
Propagate from face to edge.
This class describes the interaction of (usually) a face and a point. It carries the info of a succes...
Definition: PointIndexHit.H:66
label index() const noexcept
Return the hit index.
bool hit() const noexcept
Is there a hit?
const point_type & hitPoint() const
Return hit point. Fatal if not hit.
A list of faces which address into the list of points.
const labelListList & pointEdges() const
Return point-edge addressing.
const Field< point_type > & localPoints() const
Return pointField of points in patch.
const Field< point_type > & points() const noexcept
Return reference to global points.
const List< face_type > & localFaces() const
Return patch faces addressing into local point list.
A list of pointers to objects of type <T>, with allocation/deallocation management of the pointers....
Definition: PtrList.H:73
const T * set(const label i) const
Definition: PtrList.H:138
virtual bool read()
Re-read model coefficients if they have changed.
A list that is sorted upon construction or when explicitly requested with the sort() method.
Definition: SortableList.H:63
const labelList & indices() const noexcept
Return the list of sorted indices. Updated every sort.
Definition: SortableList.H:108
void sort()
Forward (stable) sort the list (if changed after construction).
Definition: SortableList.C:124
A List obtained as a section of another List.
Definition: SubList.H:70
static autoPtr< Time > New()
Construct (dummy) Time - no functionObjects or libraries.
Definition: Time.C:717
A List with indirect addressing. Like IndirectList but does not store addressing.
Definition: IndirectList.H:79
bool empty() const noexcept
True if the UList is empty (ie, size() is zero)
Definition: UListI.H:427
iterator end() noexcept
Return an iterator to end traversing the UList.
Definition: UListI.H:350
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
label size() const noexcept
The number of elements in the list.
Definition: UPtrListI.H:106
const Cmpt & z() const
Access to the vector z component.
Definition: VectorI.H:85
const Cmpt & y() const
Access to the vector y component.
Definition: VectorI.H:79
const Cmpt & x() const
Access to the vector x component.
Definition: VectorI.H:73
const point & min() const
Minimum describing the bounding box.
Definition: boundBoxI.H:91
const point & max() const
Maximum describing the bounding box.
Definition: boundBoxI.H:97
scalar avgDim() const
Average length/height/width dimension.
Definition: boundBoxI.H:157
Addressing for all faces on surface of mesh. Can either be read from polyMesh or from triSurface....
Definition: boundaryMesh.H:63
void addPatch(const word &patchName)
Add to back of patch list.
void writeTriSurface(const fileName &) const
Write to file.
Definition: boundaryMesh.C:757
void setFeatureEdges(const scalar minCos)
Set featureEdges, edgeToFeature, featureSegments according.
label getNTris(const label facei) const
Simple triangulation of face subset. Returns number of triangles.
void setExtraEdges(const label edgeI)
Set extraEdges to edges 'near' to edgeI. Uses point-edge walk.
void readTriSurface(const fileName &)
Read from triSurface.
Definition: boundaryMesh.C:589
void deletePatch(const word &patchName)
Delete from patch list.
wordList patchNames() const
Get names of patches.
Definition: boundaryMesh.C:275
boundaryMesh()
Default construct.
Definition: boundaryMesh.C:440
void triangulateLocal(const label startFacei, const label nFaces, const label totalNTris, labelList &triVerts, labelList &localToGlobal) const
Same as triangulate but in local vertex numbering.
labelList getNearest(const primitiveMesh &pMesh, const vector &searchSpan) const
Get bMesh index of nearest face for every boundary face in.
Definition: boundaryMesh.C:847
const bMesh & mesh() const
Definition: boundaryMesh.H:206
void changeFaces(const labelList &patchIDs, labelList &oldToNew)
Recalculate face ordering and patches. Return old to new.
label getNPoints(const label startFacei, const label nFaces) const
Number of points used in face subset.
void markFaces(const labelList &protectedEdges, const label facei, boolList &visited) const
void patchify(const labelList &nearest, const polyBoundaryMesh &oldPatches, polyMesh &newMesh) const
Take over patches onto polyMesh from nearest face in *this.
void changePatchType(const word &patchName, const word &type)
Change patch.
Like polyPatch but without reference to mesh. Used in boundaryMesh to hold data on patches....
Definition: boundaryPatch.H:57
label start() const
label size() const
static word defaultName
The default cloud name: defaultCloud.
Definition: cloud.H:90
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Definition: edge.H:66
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
Identifies a surface patch/zone by name and index, with geometric type.
const word & geometricType() const noexcept
The geometric type of the patch/zone.
static constexpr const char *const emptyType
The name for an 'empty' type.
const word & name() const noexcept
The patch/zone name.
Non-pointer based hierarchical recursive searching.
Definition: indexedOctree.H:74
A triFace with additional (region) index.
Definition: labelledTri.H:60
label region() const noexcept
Return the region index.
Definition: labelledTri.H:147
const word & physicalType() const noexcept
The (optional) physical type of the patch.
label index() const noexcept
The index of this patch in the boundaryMesh.
const word & name() const noexcept
The patch name.
A polyBoundaryMesh is a polyPatch list with additional search methods and registered IO.
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:81
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:456
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:1083
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:75
virtual autoPtr< polyPatch > clone(const labelList &faceCells) const
Construct and return a clone, setting faceCells.
Definition: polyPatch.H:238
label start() const
Return start label of this patch in the polyMesh face list.
Definition: polyPatch.H:364
Cell-face mesh analysis engine.
Definition: primitiveMesh.H:79
label nBoundaryFaces() const noexcept
Number of boundary faces (== nFaces - nInternalFaces)
const labelListList & pointEdges() const
const edgeList & edges() const
Return mesh edges. Uses calcEdges.
const vectorField & faceCentres() const
virtual const faceList & faces() const =0
Return faces.
label nInternalFaces() const noexcept
Number of internal faces.
const labelListList & edgeFaces() const
const labelListList & faceEdges() const
const vectorField & faceAreas() const
virtual const pointField & points() const =0
Return mesh points.
A mesh which allows changes in the patch distribution of the boundary faces. The change in patching i...
void changePatchID(const label faceID, const label patchID)
Change patch ID for a boundary face. Note: patchID should be in new.
void changePatches(const List< polyPatch * > &patches)
Change patches.
void repatch()
Re-patch the mesh.
Standard boundBox with extra functionality for use in octree.
Definition: treeBoundBox.H:89
Encapsulation of data needed to search on PrimitivePatches.
Triangulated surface description with patch information.
Definition: triSurface.H:79
void write(Ostream &os) const
Write to Ostream in simple OpenFOAM format.
Definition: triSurfaceIO.C:336
const geometricSurfacePatchList & patches() const noexcept
Definition: triSurface.H:399
A class for handling words, derived from Foam::string.
Definition: word.H:68
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
Definition: className.H:121
const polyBoundaryMesh & patches
dynamicFvMesh & mesh
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
const labelList nFaces(UPstream::listGatherValues< label >(aMesh.nFaces()))
const labelList nEdges(UPstream::listGatherValues< label >(aMesh.nEdges()))
const pointField & points
label nPoints
Namespace for OpenFOAM.
scalar distance(const vector &p1, const vector &p2)
Definition: curveTools.C:12
PtrList< polyPatch > polyPatchList
Store lists of polyPatch as a PtrList.
Definition: polyPatch.H:63
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:47
dimensionedScalar sign(const dimensionedScalar &ds)
List< label > labelList
A List of labels.
Definition: List.H:66
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
errorManip< error > abort(error &err)
Definition: errorManip.H:144
PrimitivePatch< List< face >, const pointField > bMesh
Holder of faceList and points. (v.s. e.g. primitivePatch which references points)
Definition: bMesh.H:48
List< bool > boolList
A List of bools.
Definition: List.H:64
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
error FatalError
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
Definition: exprTraits.C:59
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
Info<< "Finished reading KIVA file"<< endl;cellShapeList cellShapes(nPoints);labelList cellZoning(nPoints, -1);const cellModel &hex=cellModel::ref(cellModel::HEX);labelList hexLabels(8);label activeCells=0;labelList pointMap(nPoints);forAll(pointMap, i){ pointMap[i]=i;}for(label i=0;i< nPoints;i++){ if(f[i] > 0.0) { hexLabels[0]=i;hexLabels[1]=i1tab[i];hexLabels[2]=i3tab[i1tab[i]];hexLabels[3]=i3tab[i];hexLabels[4]=i8tab[i];hexLabels[5]=i1tab[i8tab[i]];hexLabels[6]=i3tab[i1tab[i8tab[i]]];hexLabels[7]=i3tab[i8tab[i]];cellShapes[activeCells].reset(hex, hexLabels);edgeList edges=cellShapes[activeCells].edges();forAll(edges, ei) { if(edges[ei].mag(points)< SMALL) { label start=pointMap[edges[ei].start()];while(start !=pointMap[start]) { start=pointMap[start];} label end=pointMap[edges[ei].end()];while(end !=pointMap[end]) { end=pointMap[end];} label minLabel=min(start, end);pointMap[start]=pointMap[end]=minLabel;} } cellZoning[activeCells]=idreg[i];activeCells++;}}cellShapes.setSize(activeCells);cellZoning.setSize(activeCells);forAll(cellShapes, celli){ cellShape &cs=cellShapes[celli];forAll(cs, i) { cs[i]=pointMap[cs[i]];} cs.collapse();}label bcIDs[11]={-1, 0, 2, 4, -1, 5, -1, 6, 7, 8, 9};const label nBCs=12;const word *kivaPatchTypes[nBCs]={ &wallPolyPatch::typeName, &wallPolyPatch::typeName, &wallPolyPatch::typeName, &wallPolyPatch::typeName, &symmetryPolyPatch::typeName, &wedgePolyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &symmetryPolyPatch::typeName, &oldCyclicPolyPatch::typeName};enum patchTypeNames{ PISTON, VALVE, LINER, CYLINDERHEAD, AXIS, WEDGE, INFLOW, OUTFLOW, PRESIN, PRESOUT, SYMMETRYPLANE, CYCLIC};const char *kivaPatchNames[nBCs]={ "piston", "valve", "liner", "cylinderHead", "axis", "wedge", "inflow", "outflow", "presin", "presout", "symmetryPlane", "cyclic"};List< SLList< face > > pFaces[nBCs]
Definition: readKivaGrid.H:235
labelList f(nPoints)
volScalarField & e
Definition: createFields.H:11
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333