triSurfaceTools.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) 2016-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 "triSurfaceTools.H"
30#include "triSurface.H"
31#include "MeshedSurface.H"
32#include "OFstream.H"
33#include "mergePoints.H"
34#include "polyMesh.H"
35#include "plane.H"
36#include "geompack.H"
37
38// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
39
40const Foam::label Foam::triSurfaceTools::ANYEDGE = -1;
41const Foam::label Foam::triSurfaceTools::NOEDGE = -2;
42const Foam::label Foam::triSurfaceTools::COLLAPSED = -3;
43
44// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
45
46/*
47 Refine by splitting all three edges of triangle ('red' refinement).
48 Neighbouring triangles (which are not marked for refinement get split
49 in half ('green') refinement. (R. Verfuerth, "A review of a posteriori
50 error estimation and adaptive mesh refinement techniques",
51 Wiley-Teubner, 1996)
52*/
53
54// Facei gets refined ('red'). Propagate edge refinement.
55void Foam::triSurfaceTools::calcRefineStatus
56(
57 const triSurface& surf,
58 const label facei,
59 List<refineType>& refine
60)
61{
62 if (refine[facei] == RED)
63 {
64 // Already marked for refinement. Do nothing.
65 }
66 else
67 {
68 // Not marked or marked for 'green' refinement. Refine.
69 refine[facei] = RED;
70
71 const labelList& myNeighbours = surf.faceFaces()[facei];
72
73 for (const label neighbourFacei : myNeighbours)
74 {
75 if (refine[neighbourFacei] == GREEN)
76 {
77 // Change to red refinement and propagate
78 calcRefineStatus(surf, neighbourFacei, refine);
79 }
80 else if (refine[neighbourFacei] == NONE)
81 {
82 refine[neighbourFacei] = GREEN;
83 }
84 }
85 }
86}
87
88
89// Split facei along edgeI at position newPointi
90void Foam::triSurfaceTools::greenRefine
91(
92 const triSurface& surf,
93 const label facei,
94 const label edgeI,
95 const label newPointi,
96 DynamicList<labelledTri>& newFaces
97)
98{
99 const labelledTri& f = surf.localFaces()[facei];
100 const edge& e = surf.edges()[edgeI];
101
102 // Find index of edge in face.
103
104 label fp0 = f.find(e[0]);
105 label fp1 = f.fcIndex(fp0);
106 label fp2 = f.fcIndex(fp1);
107
108 if (f[fp1] == e[1])
109 {
110 // Edge oriented like face
111 newFaces.append
112 (
113 labelledTri
114 (
115 f[fp0],
116 newPointi,
117 f[fp2],
118 f.region()
119 )
120 );
121 newFaces.append
122 (
123 labelledTri
124 (
125 newPointi,
126 f[fp1],
127 f[fp2],
128 f.region()
129 )
130 );
131 }
132 else
133 {
134 newFaces.append
135 (
136 labelledTri
137 (
138 f[fp2],
139 newPointi,
140 f[fp1],
141 f.region()
142 )
143 );
144 newFaces.append
145 (
146 labelledTri
147 (
148 newPointi,
149 f[fp0],
150 f[fp1],
151 f.region()
152 )
153 );
154 }
155}
156
157
158// Refine all triangles marked for refinement.
159Foam::triSurface Foam::triSurfaceTools::doRefine
160(
161 const triSurface& surf,
162 const List<refineType>& refineStatus
163)
164{
165 // Storage for new points. (start after old points)
166 label newVertI = surf.nPoints();
167
168 DynamicList<point> newPoints(newVertI);
169 newPoints.append(surf.localPoints());
170
171 // Storage for new faces
172 DynamicList<labelledTri> newFaces(surf.size());
173
174
175 // Point index for midpoint on edge
176 labelList edgeMid(surf.nEdges(), -1);
177
178 forAll(refineStatus, facei)
179 {
180 if (refineStatus[facei] == RED)
181 {
182 // Create new vertices on all edges to be refined.
183 const labelList& fEdges = surf.faceEdges()[facei];
184
185 for (const label edgei : fEdges)
186 {
187 if (edgeMid[edgei] == -1)
188 {
189 const edge& e = surf.edges()[edgei];
190
191 // Create new point on mid of edge
192 newPoints.append(e.centre(surf.localPoints()));
193 edgeMid[edgei] = newVertI++;
194 }
195 }
196
197 // Now we have new mid edge vertices for all edges on face
198 // so create triangles for RED refinement.
199
200 const edgeList& edges = surf.edges();
201
202 // Corner triangles
203 newFaces.append
204 (
205 labelledTri
206 (
207 edgeMid[fEdges[0]],
208 edges[fEdges[0]].commonVertex(edges[fEdges[1]]),
209 edgeMid[fEdges[1]],
210 surf[facei].region()
211 )
212 );
213
214 newFaces.append
215 (
216 labelledTri
217 (
218 edgeMid[fEdges[1]],
219 edges[fEdges[1]].commonVertex(edges[fEdges[2]]),
220 edgeMid[fEdges[2]],
221 surf[facei].region()
222 )
223 );
224
225 newFaces.append
226 (
227 labelledTri
228 (
229 edgeMid[fEdges[2]],
230 edges[fEdges[2]].commonVertex(edges[fEdges[0]]),
231 edgeMid[fEdges[0]],
232 surf[facei].region()
233 )
234 );
235
236 // Inner triangle
237 newFaces.append
238 (
239 labelledTri
240 (
241 edgeMid[fEdges[0]],
242 edgeMid[fEdges[1]],
243 edgeMid[fEdges[2]],
244 surf[facei].region()
245 )
246 );
247
248
249 // Create triangles for GREEN refinement.
250 for (const label edgei : fEdges)
251 {
252 label otherFacei = otherFace(surf, facei, edgei);
253
254 if ((otherFacei != -1) && (refineStatus[otherFacei] == GREEN))
255 {
256 greenRefine
257 (
258 surf,
259 otherFacei,
260 edgei,
261 edgeMid[edgei],
262 newFaces
263 );
264 }
265 }
266 }
267 }
268
269 // Copy unmarked triangles since keep original vertex numbering.
270 forAll(refineStatus, facei)
271 {
272 if (refineStatus[facei] == NONE)
273 {
274 newFaces.append(surf.localFaces()[facei]);
275 }
276 }
277
278 newFaces.shrink();
279 newPoints.shrink();
280
281
282 // Transfer DynamicLists to straight ones.
284 allPoints.transfer(newPoints);
285 newPoints.clear();
286
287 return triSurface(newFaces, surf.patches(), allPoints, true);
288}
289
290
291// Angle between two neighbouring triangles,
292// angle between shared-edge end points and left and right face end points
293Foam::scalar Foam::triSurfaceTools::faceCosAngle
294(
295 const point& pStart,
296 const point& pEnd,
297 const point& pLeft,
298 const point& pRight
299)
300{
301 const vector common(pEnd - pStart);
302 const vector base0(pLeft - pStart);
303 const vector base1(pRight - pStart);
304
305 const vector n0 = normalised(common ^ base0);
306 const vector n1 = normalised(base1 ^ common);
307
308 return n0 & n1;
309}
310
311
312// Protect edges around vertex from collapsing.
313// Moving vertI to a different position
314// - affects obviously the cost of the faces using it
315// - but also their neighbours since the edge between the faces changes
316void Foam::triSurfaceTools::protectNeighbours
317(
318 const triSurface& surf,
319 const label vertI,
320 labelList& faceStatus
321)
322{
323// for (const label facei : surf.pointFaces()[vertI])
324// {
325// if ((faceStatus[facei] == ANYEDGE) || (faceStatus[facei] >= 0))
326// {
327// faceStatus[facei] = NOEDGE;
328// }
329// }
330
331 for (const label edgei : surf.pointEdges()[vertI])
332 {
333 for (const label facei : surf.edgeFaces()[edgei])
334 {
335 if ((faceStatus[facei] == ANYEDGE) || (faceStatus[facei] >= 0))
336 {
337 faceStatus[facei] = NOEDGE;
338 }
339 }
340 }
341}
342
343
344//
345// Edge collapse helper functions
346//
347
348// Get all faces that will get collapsed if edgeI collapses.
349Foam::labelHashSet Foam::triSurfaceTools::getCollapsedFaces
350(
351 const triSurface& surf,
352 label edgeI
353)
354{
355 const edge& e = surf.edges()[edgeI];
356 const label v1 = e.start();
357 const label v2 = e.end();
358
359 // Faces using edge will certainly get collapsed.
360 const labelList& myFaces = surf.edgeFaces()[edgeI];
361
362 labelHashSet facesToBeCollapsed(2*myFaces.size());
363 facesToBeCollapsed.insert(myFaces);
364
365 // From faces using v1 check if they share an edge with faces
366 // using v2.
367 // - share edge: are part of 'splay' tree and will collapse if edge
368 // collapses
369 const labelList& v1Faces = surf.pointFaces()[v1];
370
371 for (const label face1I : v1Faces)
372 {
373 label otherEdgeI = oppositeEdge(surf, face1I, v1);
374
375 // Step across edge to other face
376 label face2I = otherFace(surf, face1I, otherEdgeI);
377
378 if (face2I != -1)
379 {
380 // found face on other side of edge. Now check if uses v2.
381 if (oppositeVertex(surf, face2I, otherEdgeI) == v2)
382 {
383 // triangles face1I and face2I form splay tree and will
384 // collapse.
385 facesToBeCollapsed.insert(face1I);
386 facesToBeCollapsed.insert(face2I);
387 }
388 }
389 }
390
391 return facesToBeCollapsed;
392}
393
394
395// Return value of faceUsed for faces using vertI
396Foam::label Foam::triSurfaceTools::vertexUsesFace
397(
398 const triSurface& surf,
399 const labelHashSet& faceUsed,
400 const label vertI
401)
402{
403 for (const label face1I : surf.pointFaces()[vertI])
404 {
405 if (faceUsed.found(face1I))
406 {
407 return face1I;
408 }
409 }
410 return -1;
411}
412
413
414// Calculate new edge-face addressing resulting from edge collapse
415void Foam::triSurfaceTools::getMergedEdges
416(
417 const triSurface& surf,
418 const label edgeI,
419 const labelHashSet& collapsedFaces,
420 Map<label>& edgeToEdge,
421 Map<label>& edgeToFace
422)
423{
424 const edge& e = surf.edges()[edgeI];
425 const label v1 = e.start();
426 const label v2 = e.end();
427
428 const labelList& v1Faces = surf.pointFaces()[v1];
429 const labelList& v2Faces = surf.pointFaces()[v2];
430
431 // Mark all (non collapsed) faces using v2
432 labelHashSet v2FacesHash(v2Faces.size());
433
434 for (const label facei : v2Faces)
435 {
436 if (!collapsedFaces.found(facei))
437 {
438 v2FacesHash.insert(facei);
439 }
440 }
441
442
443 for (const label face1I: v1Faces)
444 {
445 if (collapsedFaces.found(face1I))
446 {
447 continue;
448 }
449
450 //
451 // Check if face1I uses a vertex connected to a face using v2
452 //
453
454 label vert1I = -1;
455 label vert2I = -1;
456 otherVertices
457 (
458 surf,
459 face1I,
460 v1,
461 vert1I,
462 vert2I
463 );
464 //Pout<< "Face:" << surf.localFaces()[face1I] << " other vertices:"
465 // << vert1I << ' ' << vert2I << endl;
466
467 // Check vert1, vert2 for usage by v2Face.
468
469 label commonVert = vert1I;
470 label face2I = vertexUsesFace(surf, v2FacesHash, commonVert);
471 if (face2I == -1)
472 {
473 commonVert = vert2I;
474 face2I = vertexUsesFace(surf, v2FacesHash, commonVert);
475 }
476
477 if (face2I != -1)
478 {
479 // Found one: commonVert is used by both face1I and face2I
480 label edge1I = getEdge(surf, v1, commonVert);
481 label edge2I = getEdge(surf, v2, commonVert);
482
483 edgeToEdge.insert(edge1I, edge2I);
484 edgeToEdge.insert(edge2I, edge1I);
485
486 edgeToFace.insert(edge1I, face2I);
487 edgeToFace.insert(edge2I, face1I);
488 }
489 }
490}
491
492
493// Calculates (cos of) angle across edgeI of facei,
494// taking into account updated addressing (resulting from edge collapse)
495Foam::scalar Foam::triSurfaceTools::edgeCosAngle
496(
497 const triSurface& surf,
498 const label v1,
499 const point& pt,
500 const labelHashSet& collapsedFaces,
501 const Map<label>& edgeToEdge,
502 const Map<label>& edgeToFace,
503 const label facei,
504 const label edgeI
505)
506{
507 const pointField& localPoints = surf.localPoints();
508
509 label A = surf.edges()[edgeI].start();
510 label B = surf.edges()[edgeI].end();
511 label C = oppositeVertex(surf, facei, edgeI);
512
513 label D = -1;
514
515 label face2I = -1;
516
517 if (edgeToEdge.found(edgeI))
518 {
519 // Use merged addressing
520 label edge2I = edgeToEdge[edgeI];
521 face2I = edgeToFace[edgeI];
522
523 D = oppositeVertex(surf, face2I, edge2I);
524 }
525 else
526 {
527 // Use normal edge-face addressing
528 face2I = otherFace(surf, facei, edgeI);
529
530 if ((face2I != -1) && !collapsedFaces.found(face2I))
531 {
532 D = oppositeVertex(surf, face2I, edgeI);
533 }
534 }
535
536 scalar cosAngle = 1;
537
538 if (D != -1)
539 {
540 if (A == v1)
541 {
542 cosAngle = faceCosAngle
543 (
544 pt,
545 localPoints[B],
546 localPoints[C],
547 localPoints[D]
548 );
549 }
550 else if (B == v1)
551 {
552 cosAngle = faceCosAngle
553 (
554 localPoints[A],
555 pt,
556 localPoints[C],
557 localPoints[D]
558 );
559 }
560 else if (C == v1)
561 {
562 cosAngle = faceCosAngle
563 (
564 localPoints[A],
565 localPoints[B],
566 pt,
567 localPoints[D]
568 );
569 }
570 else if (D == v1)
571 {
572 cosAngle = faceCosAngle
573 (
574 localPoints[A],
575 localPoints[B],
576 localPoints[C],
577 pt
578 );
579 }
580 else
581 {
583 << "face " << facei << " does not use vertex "
584 << v1 << " of collapsed edge" << abort(FatalError);
585 }
586 }
587 return cosAngle;
588}
589
590
591Foam::scalar Foam::triSurfaceTools::collapseMinCosAngle
592(
593 const triSurface& surf,
594 const label v1,
595 const point& pt,
596 const labelHashSet& collapsedFaces,
597 const Map<label>& edgeToEdge,
598 const Map<label>& edgeToFace
599)
600{
601 const labelList& v1Faces = surf.pointFaces()[v1];
602
603 scalar minCos = 1;
604
605 for (const label facei : v1Faces)
606 {
607 if (collapsedFaces.found(facei))
608 {
609 continue;
610 }
611
612 for (const label edgeI : surf.faceEdges()[facei])
613 {
614 minCos =
615 min
616 (
617 minCos,
618 edgeCosAngle
619 (
620 surf,
621 v1,
622 pt,
623 collapsedFaces,
624 edgeToEdge,
625 edgeToFace,
626 facei,
627 edgeI
628 )
629 );
630 }
631 }
632
633 return minCos;
634}
635
636
637// Calculate max dihedral angle after collapsing edge to v1 (at pt).
638// Note that all edges of all faces using v1 are affected.
639bool Foam::triSurfaceTools::collapseCreatesFold
640(
641 const triSurface& surf,
642 const label v1,
643 const point& pt,
644 const labelHashSet& collapsedFaces,
645 const Map<label>& edgeToEdge,
646 const Map<label>& edgeToFace,
647 const scalar minCos
648)
649{
650 const labelList& v1Faces = surf.pointFaces()[v1];
651
652 for (const label facei : v1Faces)
653 {
654 if (collapsedFaces.found(facei))
655 {
656 continue;
657 }
658
659 const labelList& myEdges = surf.faceEdges()[facei];
660
661 for (const label edgeI : myEdges)
662 {
663 if
664 (
665 edgeCosAngle
666 (
667 surf,
668 v1,
669 pt,
670 collapsedFaces,
671 edgeToEdge,
672 edgeToFace,
673 facei,
674 edgeI
675 )
676 < minCos
677 )
678 {
679 return true;
680 }
681 }
682 }
683
684 return false;
685}
686
687
690// a vertex.
691//bool Foam::triSurfaceTools::collapseCreatesDuplicates
692//(
693// const triSurface& surf,
694// const label edgeI,
695// const labelHashSet& collapsedFaces
696//)
697//{
698//Pout<< "duplicate : edgeI:" << edgeI << surf.edges()[edgeI]
699// << " collapsedFaces:" << collapsedFaces.toc() << endl;
700//
701// // Mark neighbours of faces to be collapsed, i.e. get the first layer
702// // of triangles outside collapsedFaces.
703// // neighbours actually contains the
704// // edge with which triangle connects to collapsedFaces.
705//
706// Map<label> neighbours;
707//
708// labelList collapsed = collapsedFaces.toc();
709//
710// for (const label facei : collapsed)
711// {
712// const labelList& myEdges = surf.faceEdges()[facei];
713//
714// Pout<< "collapsing facei:" << facei << " uses edges:" << myEdges
715// << endl;
716//
717// forAll(myEdges, myEdgeI)
718// {
719// const labelList& myFaces = surf.edgeFaces()[myEdges[myEdgeI]];
720//
721// Pout<< "Edge " << myEdges[myEdgeI] << " is used by faces "
722// << myFaces << endl;
723//
724// if ((myEdges[myEdgeI] != edgeI) && (myFaces.size() == 2))
725// {
726// // Get the other face
727// label neighbourFacei = myFaces[0];
728//
729// if (neighbourFacei == facei)
730// {
731// neighbourFacei = myFaces[1];
732// }
733//
734// // Is 'outside' face if not in collapsedFaces itself
735// if (!collapsedFaces.found(neighbourFacei))
736// {
737// neighbours.insert(neighbourFacei, myEdges[myEdgeI]);
738// }
739// }
740// }
741// }
742//
743// // Now neighbours contains first layer of triangles outside of
744// // collapseFaces
745// // There should be
746// // -two if edgeI is a boundary edge
747// // since the outside 'edge' of collapseFaces should
748// // form a triangle and the face connected to edgeI is not inserted.
749// // -four if edgeI is not a boundary edge since then the outside edge forms
750// // a diamond.
751//
752// // Check if any of the faces in neighbours share an edge. (n^2)
753//
754// labelList neighbourList = neighbours.toc();
755//
756//Pout<< "edgeI:" << edgeI << " neighbourList:" << neighbourList << endl;
757//
758//
759// forAll(neighbourList, i)
760// {
761// const labelList& faceIEdges = surf.faceEdges()[neighbourList[i]];
762//
763// for (label j = i+1; j < neighbourList.size(); i++)
764// {
765// const labelList& faceJEdges = surf.faceEdges()[neighbourList[j]];
766//
767// // Check if facei and faceJ share an edge
768// forAll(faceIEdges, fI)
769// {
770// forAll(faceJEdges, fJ)
771// {
772// Pout<< " comparing " << faceIEdges[fI] << " to "
773// << faceJEdges[fJ] << endl;
774//
775// if (faceIEdges[fI] == faceJEdges[fJ])
776// {
777// return true;
778// }
779// }
780// }
781// }
782// }
783// Pout<< "Found no match. Returning false" << endl;
784// return false;
785//}
786
787
788// Finds the triangle edge cut by the plane between a point inside / on edge
789// of a triangle and a point outside. Returns:
790// - cut through edge/point
791// - miss()
792Foam::surfaceLocation Foam::triSurfaceTools::cutEdge
793(
794 const triSurface& s,
795 const label triI,
796 const label excludeEdgeI,
797 const label excludePointi,
798
799 const point& triPoint,
800 const plane& cutPlane,
801 const point& toPoint
802)
803{
804 const pointField& points = s.points();
805 const labelledTri& f = s[triI];
806 const labelList& fEdges = s.faceEdges()[triI];
807
808 // Get normal distance to planeN
809 FixedList<scalar, 3> d;
810
811 scalar norm = 0;
812 forAll(d, fp)
813 {
814 d[fp] = cutPlane.signedDistance(points[f[fp]]);
815 norm += mag(d[fp]);
816 }
817
818 // Normalise and truncate
819 forAll(d, i)
820 {
821 d[i] /= norm;
822
823 if (mag(d[i]) < 1e-6)
824 {
825 d[i] = 0.0;
826 }
827 }
828
829 // Return information
830 surfaceLocation cut;
831
832 if (excludePointi != -1)
833 {
834 // Excluded point. Test only opposite edge.
835
836 label fp0 = s.localFaces()[triI].find(excludePointi);
837
838 if (fp0 == -1)
839 {
841 << " localF:" << s.localFaces()[triI] << abort(FatalError);
842 }
843
844 label fp1 = f.fcIndex(fp0);
845 label fp2 = f.fcIndex(fp1);
846
847
848 if (d[fp1] == 0.0)
849 {
850 cut.setHit();
851 cut.setPoint(points[f[fp1]]);
852 cut.elementType() = triPointRef::POINT;
853 cut.setIndex(s.localFaces()[triI][fp1]);
854 }
855 else if (d[fp2] == 0.0)
856 {
857 cut.setHit();
858 cut.setPoint(points[f[fp2]]);
859 cut.elementType() = triPointRef::POINT;
860 cut.setIndex(s.localFaces()[triI][fp2]);
861 }
862 else if
863 (
864 (d[fp1] < 0 && d[fp2] < 0)
865 || (d[fp1] > 0 && d[fp2] > 0)
866 )
867 {
868 // Both same sign. Not crossing edge at all.
869 // cut already set to miss().
870 }
871 else
872 {
873 cut.setHit();
874 cut.setPoint
875 (
876 (d[fp2]*points[f[fp1]] - d[fp1]*points[f[fp2]])
877 / (d[fp2] - d[fp1])
878 );
879 cut.elementType() = triPointRef::EDGE;
880 cut.setIndex(fEdges[fp1]);
881 }
882 }
883 else
884 {
885 // Find the two intersections
886 FixedList<surfaceLocation, 2> inters;
887 label interI = 0;
888
889 forAll(f, fp0)
890 {
891 label fp1 = f.fcIndex(fp0);
892
893 if (d[fp0] == 0)
894 {
895 if (interI >= 2)
896 {
898 << "problem : triangle has three intersections." << nl
899 << "triangle:" << f.tri(points)
900 << " d:" << d << abort(FatalError);
901 }
902 inters[interI].setHit();
903 inters[interI].setPoint(points[f[fp0]]);
904 inters[interI].elementType() = triPointRef::POINT;
905 inters[interI].setIndex(s.localFaces()[triI][fp0]);
906 interI++;
907 }
908 else if
909 (
910 (d[fp0] < 0 && d[fp1] > 0)
911 || (d[fp0] > 0 && d[fp1] < 0)
912 )
913 {
914 if (interI >= 2)
915 {
917 << "problem : triangle has three intersections." << nl
918 << "triangle:" << f.tri(points)
919 << " d:" << d << abort(FatalError);
920 }
921 inters[interI].setHit();
922 inters[interI].setPoint
923 (
924 (d[fp0]*points[f[fp1]] - d[fp1]*points[f[fp0]])
925 / (d[fp0] - d[fp1])
926 );
927 inters[interI].elementType() = triPointRef::EDGE;
928 inters[interI].setIndex(fEdges[fp0]);
929 interI++;
930 }
931 }
932
933
934 if (interI == 0)
935 {
936 // Return miss
937 }
938 else if (interI == 1)
939 {
940 // Only one intersection. Should not happen!
941 cut = inters[0];
942 }
943 else if (interI == 2)
944 {
945 // Handle excludeEdgeI
946 if
947 (
948 inters[0].elementType() == triPointRef::EDGE
949 && inters[0].index() == excludeEdgeI
950 )
951 {
952 cut = inters[1];
953 }
954 else if
955 (
956 inters[1].elementType() == triPointRef::EDGE
957 && inters[1].index() == excludeEdgeI
958 )
959 {
960 cut = inters[0];
961 }
962 else
963 {
964 // Two cuts. Find nearest.
965 if
966 (
967 magSqr(inters[0].rawPoint() - toPoint)
968 < magSqr(inters[1].rawPoint() - toPoint)
969 )
970 {
971 cut = inters[0];
972 }
973 else
974 {
975 cut = inters[1];
976 }
977 }
978 }
979 }
980 return cut;
981}
982
983
984// 'Snap' point on to endPoint.
985void Foam::triSurfaceTools::snapToEnd
986(
987 const triSurface& s,
988 const surfaceLocation& end,
989 surfaceLocation& current
990)
991{
992 if (end.elementType() == triPointRef::NONE)
993 {
994 if (current.elementType() == triPointRef::NONE)
995 {
996 // endpoint on triangle; current on triangle
997 if (current.index() == end.index())
998 {
999 //if (debug)
1000 //{
1001 // Pout<< "snapToEnd : snapping:" << current << " onto:"
1002 // << end << endl;
1003 //}
1004 current = end;
1005 current.setHit();
1006 }
1007 }
1008 // No need to handle current on edge/point since tracking handles this.
1009 }
1010 else if (end.elementType() == triPointRef::EDGE)
1011 {
1012 if (current.elementType() == triPointRef::NONE)
1013 {
1014 // endpoint on edge; current on triangle
1015 const labelList& fEdges = s.faceEdges()[current.index()];
1016
1017 if (fEdges.found(end.index()))
1018 {
1019 //if (debug)
1020 //{
1021 // Pout<< "snapToEnd : snapping:" << current << " onto:"
1022 // << end << endl;
1023 //}
1024 current = end;
1025 current.setHit();
1026 }
1027 }
1028 else if (current.elementType() == triPointRef::EDGE)
1029 {
1030 // endpoint on edge; current on edge
1031 if (current.index() == end.index())
1032 {
1033 //if (debug)
1034 //{
1035 // Pout<< "snapToEnd : snapping:" << current << " onto:"
1036 // << end << endl;
1037 //}
1038 current = end;
1039 current.setHit();
1040 }
1041 }
1042 else
1043 {
1044 // endpoint on edge; current on point
1045 const edge& e = s.edges()[end.index()];
1046
1047 if (current.index() == e[0] || current.index() == e[1])
1048 {
1049 //if (debug)
1050 //{
1051 // Pout<< "snapToEnd : snapping:" << current << " onto:"
1052 // << end << endl;
1053 //}
1054 current = end;
1055 current.setHit();
1056 }
1057 }
1058 }
1059 else // end.elementType() == POINT
1060 {
1061 if (current.elementType() == triPointRef::NONE)
1062 {
1063 // endpoint on point; current on triangle
1064 const triSurface::face_type& f = s.localFaces()[current.index()];
1065
1066 if (f.found(end.index()))
1067 {
1068 //if (debug)
1069 //{
1070 // Pout<< "snapToEnd : snapping:" << current << " onto:"
1071 // << end << endl;
1072 //}
1073 current = end;
1074 current.setHit();
1075 }
1076 }
1077 else if (current.elementType() == triPointRef::EDGE)
1078 {
1079 // endpoint on point; current on edge
1080 const edge& e = s.edges()[current.index()];
1081
1082 if (end.index() == e[0] || end.index() == e[1])
1083 {
1084 //if (debug)
1085 //{
1086 // Pout<< "snapToEnd : snapping:" << current << " onto:"
1087 // << end << endl;
1088 //}
1089 current = end;
1090 current.setHit();
1091 }
1092 }
1093 else
1094 {
1095 // endpoint on point; current on point
1096 if (current.index() == end.index())
1097 {
1098 //if (debug)
1099 //{
1100 // Pout<< "snapToEnd : snapping:" << current << " onto:"
1101 // << end << endl;
1102 //}
1103 current = end;
1104 current.setHit();
1105 }
1106 }
1107 }
1108}
1109
1110
1111// Start:
1112// - location
1113// - element type (triangle/edge/point) and index
1114// - triangle to exclude
1115Foam::surfaceLocation Foam::triSurfaceTools::visitFaces
1116(
1117 const triSurface& s,
1118 const labelList& eFaces,
1119 const surfaceLocation& start,
1120 const label excludeEdgeI,
1121 const label excludePointi,
1122 const surfaceLocation& end,
1123 const plane& cutPlane
1124)
1125{
1126 surfaceLocation nearest;
1127
1128 scalar minDistSqr = Foam::sqr(GREAT);
1129
1130 for (const label triI : eFaces)
1131 {
1132 // Make sure we don't revisit previous face
1133 if (triI != start.triangle())
1134 {
1135 if (end.elementType() == triPointRef::NONE && end.index() == triI)
1136 {
1137 // Endpoint is in this triangle. Jump there.
1138 nearest = end;
1139 nearest.setHit();
1140 nearest.triangle() = triI;
1141 break;
1142 }
1143 else
1144 {
1145 // Which edge is cut.
1146
1147 surfaceLocation cutInfo = cutEdge
1148 (
1149 s,
1150 triI,
1151 excludeEdgeI, // excludeEdgeI
1152 excludePointi, // excludePointi
1153 start.rawPoint(),
1154 cutPlane,
1155 end.rawPoint()
1156 );
1157
1158 // If crossing an edge we expect next edge to be cut.
1159 if (excludeEdgeI != -1 && !cutInfo.hit())
1160 {
1162 << "Triangle:" << triI
1163 << " excludeEdge:" << excludeEdgeI
1164 << " point:" << start.rawPoint()
1165 << " plane:" << cutPlane
1166 << " . No intersection!" << abort(FatalError);
1167 }
1168
1169 if (cutInfo.hit())
1170 {
1171 scalar distSqr = magSqr(cutInfo.rawPoint()-end.rawPoint());
1172
1173 if (distSqr < minDistSqr)
1174 {
1175 minDistSqr = distSqr;
1176 nearest = cutInfo;
1177 nearest.triangle() = triI;
1178 nearest.setMiss();
1179 }
1180 }
1181 }
1182 }
1183 }
1184
1185 if (nearest.triangle() == -1)
1186 {
1187 // Did not move from edge. Give warning? Return something special?
1188 // For now responsibility of caller to make sure that nothing has
1189 // moved.
1190 }
1191
1192 return nearest;
1193}
1194
1195
1196// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
1197
1198
1199// Write pointField to file
1201(
1202 const fileName& fName,
1203 const pointField& pts
1204)
1205{
1206 OFstream outFile(fName);
1207
1208 for (const point& pt : pts)
1209 {
1210 outFile<< "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << endl;
1211 }
1212 Pout<< "Written " << pts.size() << " vertices to file " << fName << endl;
1213}
1214
1215
1216// Write vertex subset to OBJ format file
1218(
1219 const triSurface& surf,
1220 const fileName& fName,
1221 const boolList& markedVerts
1222)
1223{
1224 OFstream outFile(fName);
1225
1226 label nVerts = 0;
1227 forAll(markedVerts, vertI)
1228 {
1229 if (markedVerts[vertI])
1230 {
1231 const point& pt = surf.localPoints()[vertI];
1232
1233 outFile<< "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << endl;
1234
1235 nVerts++;
1236 }
1237 }
1238 Pout<< "Written " << nVerts << " vertices to file " << fName << endl;
1239}
1240
1241
1242// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
1243// Addressing helper functions:
1244
1245
1246// Get all triangles using vertices of edge
1248(
1249 const triSurface& surf,
1250 const label edgeI,
1251 labelList& edgeTris
1252)
1253{
1254 const edge& e = surf.edges()[edgeI];
1255 const labelList& myFaces = surf.edgeFaces()[edgeI];
1256
1257 label face1I = myFaces[0];
1258 label face2I = -1;
1259 if (myFaces.size() == 2)
1260 {
1261 face2I = myFaces[1];
1262 }
1263
1264 const labelList& startFaces = surf.pointFaces()[e.start()];
1265 const labelList& endFaces = surf.pointFaces()[e.end()];
1266
1267 // Number of triangles is sum of pointfaces - common faces
1268 // (= faces using edge)
1269 edgeTris.setSize(startFaces.size() + endFaces.size() - myFaces.size());
1270
1271 label nTris = 0;
1272 for (const label facei : startFaces)
1273 {
1274 edgeTris[nTris++] = facei;
1275 }
1276
1277 for (const label facei : endFaces)
1278 {
1279 if ((facei != face1I) && (facei != face2I))
1280 {
1281 edgeTris[nTris++] = facei;
1282 }
1283 }
1284}
1285
1286
1287// Get all vertices connected to vertices of edge
1289(
1290 const triSurface& surf,
1291 const edge& e
1292)
1293{
1294 const edgeList& edges = surf.edges();
1295 const label v1 = e.start();
1296 const label v2 = e.end();
1297
1298 // Get all vertices connected to v1 or v2 through an edge
1299 labelHashSet vertexNeighbours;
1300
1301 const labelList& v1Edges = surf.pointEdges()[v1];
1302
1303 for (const label edgei : v1Edges)
1304 {
1305 vertexNeighbours.insert(edges[edgei].otherVertex(v1));
1306 }
1307
1308 const labelList& v2Edges = surf.pointEdges()[v2];
1309
1310 for (const label edgei : v2Edges)
1311 {
1312 vertexNeighbours.insert(edges[edgei].otherVertex(v2));
1313 }
1314 return vertexNeighbours.toc();
1315}
1316
1317
1318// Get the other face using edgeI
1320(
1321 const triSurface& surf,
1322 const label facei,
1323 const label edgeI
1324)
1325{
1326 const labelList& myFaces = surf.edgeFaces()[edgeI];
1327
1328 if (myFaces.size() != 2)
1329 {
1330 return -1;
1331 }
1332 else if (facei == myFaces[0])
1333 {
1334 return myFaces[1];
1335 }
1336 else
1337 {
1338 return myFaces[0];
1339 }
1340}
1341
1342
1343// Get the two edges on facei counterclockwise after edgeI
1345(
1346 const triSurface& surf,
1347 const label facei,
1348 const label edgeI,
1349 label& e1,
1350 label& e2
1351)
1352{
1353 const labelList& eFaces = surf.faceEdges()[facei];
1354
1355 label i0 = eFaces.find(edgeI);
1356
1357 if (i0 == -1)
1358 {
1360 << "Edge " << surf.edges()[edgeI] << " not in face "
1361 << surf.localFaces()[facei] << abort(FatalError);
1362 }
1363
1364 label i1 = eFaces.fcIndex(i0);
1365 label i2 = eFaces.fcIndex(i1);
1366
1367 e1 = eFaces[i1];
1368 e2 = eFaces[i2];
1369}
1370
1371
1372// Get the two vertices on facei counterclockwise vertI
1374(
1375 const triSurface& surf,
1376 const label facei,
1377 const label vertI,
1378 label& vert1I,
1379 label& vert2I
1380)
1381{
1382 const labelledTri& f = surf.localFaces()[facei];
1383
1384 if (vertI == f[0])
1385 {
1386 vert1I = f[1];
1387 vert2I = f[2];
1388 }
1389 else if (vertI == f[1])
1390 {
1391 vert1I = f[2];
1392 vert2I = f[0];
1393 }
1394 else if (vertI == f[2])
1395 {
1396 vert1I = f[0];
1397 vert2I = f[1];
1398 }
1399 else
1400 {
1402 << "Vertex " << vertI << " not in face " << f << nl
1403 << abort(FatalError);
1404 }
1405}
1406
1407
1408// Get edge opposite vertex
1410(
1411 const triSurface& surf,
1412 const label facei,
1413 const label vertI
1414)
1415{
1416 const labelList& myEdges = surf.faceEdges()[facei];
1417
1418 for (const label edgei : myEdges)
1419 {
1420 const edge& e = surf.edges()[edgei];
1421
1422 if (!e.found(vertI))
1423 {
1424 return edgei;
1425 }
1426 }
1427
1429 << "Cannot find vertex " << vertI << " in edges of face " << facei
1430 << nl
1431 << abort(FatalError);
1432
1433 return -1;
1434}
1435
1436
1437// Get vertex opposite edge
1439(
1440 const triSurface& surf,
1441 const label facei,
1442 const label edgeI
1443)
1444{
1445 const triSurface::face_type& f = surf.localFaces()[facei];
1446 const edge& e = surf.edges()[edgeI];
1447
1448 for (const label pointi : f)
1449 {
1450 if (!e.found(pointi))
1451 {
1452 return pointi;
1453 }
1454 }
1455
1457 << "Cannot find vertex opposite edge " << edgeI << " vertices " << e
1458 << " in face " << facei << " vertices " << f << abort(FatalError);
1459
1460 return -1;
1461}
1462
1463
1464// Returns edge label connecting v1, v2
1466(
1467 const triSurface& surf,
1468 const label v1,
1469 const label v2
1470)
1471{
1472 const labelList& v1Edges = surf.pointEdges()[v1];
1473
1474 for (const label edgei : v1Edges)
1475 {
1476 const edge& e = surf.edges()[edgei];
1477
1478 if (e.found(v2))
1479 {
1480 return edgei;
1481 }
1482 }
1483 return -1;
1484}
1485
1486
1487// Return index of triangle (or -1) using all three edges
1489(
1490 const triSurface& surf,
1491 const label e0I,
1492 const label e1I,
1493 const label e2I
1494)
1495{
1496 if ((e0I == e1I) || (e0I == e2I) || (e1I == e2I))
1497 {
1499 << "Duplicate edge labels : e0:" << e0I << " e1:" << e1I
1500 << " e2:" << e2I
1501 << abort(FatalError);
1502 }
1503
1504 const labelList& eFaces = surf.edgeFaces()[e0I];
1505
1506 for (const label facei : eFaces)
1507 {
1508 const labelList& myEdges = surf.faceEdges()[facei];
1509
1510 if
1511 (
1512 (myEdges[0] == e1I)
1513 || (myEdges[1] == e1I)
1514 || (myEdges[2] == e1I)
1515 )
1516 {
1517 if
1518 (
1519 (myEdges[0] == e2I)
1520 || (myEdges[1] == e2I)
1521 || (myEdges[2] == e2I)
1522 )
1523 {
1524 return facei;
1525 }
1526 }
1527 }
1528 return -1;
1529}
1530
1531
1532// Collapse indicated edges. Return new tri surface.
1534(
1535 const triSurface& surf,
1536 const labelList& collapsableEdges
1537)
1538{
1539 pointField edgeMids(surf.nEdges());
1540
1541 forAll(edgeMids, edgeI)
1542 {
1543 const edge& e = surf.edges()[edgeI];
1544 edgeMids[edgeI] = e.centre(surf.localPoints());
1545 }
1546
1547
1548 labelList faceStatus(surf.size(), ANYEDGE);
1549
1551 //forAll(edges, edgeI)
1552 //{
1553 // const labelList& neighbours = edgeFaces[edgeI];
1554 //
1555 // if ((neighbours.size() != 2) && (neighbours.size() != 1))
1556 // {
1557 // FatalErrorInFunction
1558 // << abort(FatalError);
1559 // }
1560 //
1561 // if (neighbours.size() == 2)
1562 // {
1563 // if (surf[neighbours[0]].region() != surf[neighbours[1]].region())
1564 // {
1565 // // Neighbours on different regions. For now, do not allow
1566 // // any collapse.
1567 // //Pout<< "protecting face " << neighbours[0]
1568 // // << ' ' << neighbours[1] << endl;
1569 // faceStatus[neighbours[0]] = NOEDGE;
1570 // faceStatus[neighbours[1]] = NOEDGE;
1571 // }
1572 // }
1573 //}
1574
1575 return collapseEdges(surf, collapsableEdges, edgeMids, faceStatus);
1576}
1577
1578
1579// Collapse indicated edges. Return new tri surface.
1581(
1582 const triSurface& surf,
1583 const labelList& collapseEdgeLabels,
1584 const pointField& edgeMids,
1585 labelList& faceStatus
1586)
1587{
1588 const labelListList& edgeFaces = surf.edgeFaces();
1589 const pointField& localPoints = surf.localPoints();
1590 const edgeList& edges = surf.edges();
1591
1592 // Storage for new points
1593 pointField newPoints(localPoints);
1594
1595 // Map for old to new points
1596 labelList pointMap(identity(localPoints.size()));
1597
1598 // Do actual 'collapsing' of edges
1599
1600 for (const label edgei : collapseEdgeLabels)
1601 {
1602 if (edgei < 0 || edgei >= surf.nEdges())
1603 {
1605 << "Edge label outside valid range." << endl
1606 << "edge label:" << edgei << endl
1607 << "total number of edges:" << surf.nEdges() << endl
1608 << abort(FatalError);
1609 }
1610
1611 const labelList& neighbours = edgeFaces[edgei];
1612
1613 if (neighbours.size() == 2)
1614 {
1615 const label stat0 = faceStatus[neighbours[0]];
1616 const label stat1 = faceStatus[neighbours[1]];
1617
1618 // Check faceStatus to make sure this one can be collapsed
1619 if
1620 (
1621 ((stat0 == ANYEDGE) || (stat0 == edgei))
1622 && ((stat1 == ANYEDGE) || (stat1 == edgei))
1623 )
1624 {
1625 const edge& e = edges[edgei];
1626
1627 // Set up mapping to 'collapse' points of edge
1628 if
1629 (
1630 (pointMap[e.start()] != e.start())
1631 || (pointMap[e.end()] != e.end())
1632 )
1633 {
1635 << "points already mapped. Double collapse." << endl
1636 << "edgei:" << edgei
1637 << " start:" << e.start()
1638 << " end:" << e.end()
1639 << " pointMap[start]:" << pointMap[e.start()]
1640 << " pointMap[end]:" << pointMap[e.end()]
1641 << abort(FatalError);
1642 }
1643
1644 const label minVert = min(e.start(), e.end());
1645 pointMap[e.start()] = minVert;
1646 pointMap[e.end()] = minVert;
1647
1648 // Move shared vertex to mid of edge
1649 newPoints[minVert] = edgeMids[edgei];
1650
1651 // Protect neighbouring faces
1652 protectNeighbours(surf, e.start(), faceStatus);
1653 protectNeighbours(surf, e.end(), faceStatus);
1654 protectNeighbours
1655 (
1656 surf,
1657 oppositeVertex(surf, neighbours[0], edgei),
1658 faceStatus
1659 );
1660 protectNeighbours
1661 (
1662 surf,
1663 oppositeVertex(surf, neighbours[1], edgei),
1664 faceStatus
1665 );
1666
1667 // Mark all collapsing faces
1668 labelList collapseFaces =
1669 getCollapsedFaces
1670 (
1671 surf,
1672 edgei
1673 ).toc();
1674
1675 forAll(collapseFaces, collapseI)
1676 {
1677 faceStatus[collapseFaces[collapseI]] = COLLAPSED;
1678 }
1679 }
1680 }
1681 }
1682
1683
1684 // Storage for new triangles
1685 List<labelledTri> newTriangles(surf.size());
1686 label nNewTris = 0;
1687
1688 const List<labelledTri>& localFaces = surf.localFaces();
1689
1690 // Get only non-collapsed triangles and renumber vertex labels.
1691 forAll(localFaces, facei)
1692 {
1693 if (faceStatus[facei] != COLLAPSED)
1694 {
1695 // Uncollapsed triangle
1696 labelledTri f(localFaces[facei]);
1697
1698 // inplace renumber
1699 f[0] = pointMap[f[0]];
1700 f[1] = pointMap[f[1]];
1701 f[2] = pointMap[f[2]];
1702
1703 if (f.valid())
1704 {
1705 newTriangles[nNewTris++] = f;
1706 }
1707 }
1708 }
1709 newTriangles.resize(nNewTris);
1710
1711
1712 // Pack faces
1713
1714 triSurface tempSurf(newTriangles, surf.patches(), newPoints);
1715
1716 return
1718 (
1719 tempSurf.localFaces(),
1720 tempSurf.patches(),
1721 tempSurf.localPoints()
1722 );
1723}
1724
1725
1726// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
1727
1729(
1730 const triSurface& surf,
1731 const labelList& refineFaces
1732)
1733{
1734 List<refineType> refineStatus(surf.size(), NONE);
1735
1736 // Mark & propagate refinement
1737 forAll(refineFaces, refineFacei)
1738 {
1739 calcRefineStatus(surf, refineFaces[refineFacei], refineStatus);
1740 }
1741
1742 // Do actual refinement
1743 return doRefine(surf, refineStatus);
1744}
1745
1746
1747Foam::triSurface Foam::triSurfaceTools::greenRefine
1748(
1749 const triSurface& surf,
1750 const labelList& refineEdges
1751)
1752{
1753 // Storage for marking faces
1754 List<refineType> refineStatus(surf.size(), NONE);
1755
1756 // Storage for new faces
1757 DynamicList<labelledTri> newFaces(0);
1758
1759 pointField newPoints(surf.localPoints());
1760 newPoints.setSize(surf.nPoints() + surf.nEdges());
1761 label newPointi = surf.nPoints();
1762
1763
1764 // Refine edges
1765 for (const label edgei : refineEdges)
1766 {
1767 const labelList& myFaces = surf.edgeFaces()[edgei];
1768
1769 bool neighbourIsRefined= false;
1770
1771 for (const label facei : myFaces)
1772 {
1773 if (refineStatus[facei] != NONE)
1774 {
1775 neighbourIsRefined = true;
1776 }
1777 }
1778
1779 // Only refine if none of the faces is refined
1780 if (!neighbourIsRefined)
1781 {
1782 // Refine edge
1783 const edge& e = surf.edges()[edgei];
1784
1785 newPoints[newPointi] = e.centre(surf.localPoints());
1786
1787 // Refine faces using edge
1788 for (const label facei : myFaces)
1789 {
1790 // Add faces to newFaces
1791 greenRefine
1792 (
1793 surf,
1794 facei,
1795 edgei,
1796 newPointi,
1797 newFaces
1798 );
1799
1800 // Mark as refined
1801 refineStatus[facei] = GREEN;
1802 }
1803
1804 ++newPointi;
1805 }
1806 }
1807
1808 // Add unrefined faces
1809 forAll(surf.localFaces(), facei)
1810 {
1811 if (refineStatus[facei] == NONE)
1812 {
1813 newFaces.append(surf.localFaces()[facei]);
1814 }
1815 }
1816
1817 newFaces.shrink();
1818 newPoints.setSize(newPointi);
1819
1820 return triSurface(newFaces, surf.patches(), newPoints, true);
1821}
1822
1823
1824
1825// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
1826// Geometric helper functions:
1827
1828
1829// Returns element in edgeIndices with minimum length
1831(
1832 const triSurface& surf,
1833 const labelList& edgeIndices
1834)
1835{
1836 scalar minLen = GREAT;
1837 label minEdge = -1;
1838
1839 for (const label edgei : edgeIndices)
1840 {
1841 const edge& e = surf.edges()[edgei];
1842
1843 const scalar length = e.mag(surf.localPoints());
1844
1845 if (length < minLen)
1846 {
1847 minLen = length;
1848 minEdge = edgei;
1849 }
1850 }
1851
1852 return minEdge;
1853}
1854
1855
1856// Returns element in edgeIndices with maximum length
1858(
1859 const triSurface& surf,
1860 const labelList& edgeIndices
1861)
1862{
1863 scalar maxLen = -GREAT;
1864 label maxEdge = -1;
1865
1866 for (const label edgei : edgeIndices)
1867 {
1868 const edge& e = surf.edges()[edgei];
1869
1870 const scalar length = e.mag(surf.localPoints());
1871
1872 if (length > maxLen)
1873 {
1874 maxLen = length;
1875 maxEdge = edgei;
1876 }
1877 }
1878
1879 return maxEdge;
1880}
1881
1882
1883// Merge points and reconstruct surface
1885(
1886 const triSurface& surf,
1887 const scalar mergeTol
1888)
1889{
1890 labelList pointMap;
1891 labelList uniquePoints;
1892
1893 label nChanged = Foam::mergePoints
1894 (
1895 surf.localPoints(),
1896 pointMap,
1897 uniquePoints,
1898 mergeTol,
1899 false
1900 );
1901
1902 if (nChanged)
1903 {
1904 // Pack the triangles
1905
1906 // Storage for new triangles
1907 List<labelledTri> newTriangles(surf.size());
1908 label nNewTris = 0;
1909
1910 // Iterate and work on a copy
1911 for (labelledTri f : surf.localFaces())
1912 {
1913 // inplace renumber
1914 f[0] = pointMap[f[0]];
1915 f[1] = pointMap[f[1]];
1916 f[2] = pointMap[f[2]];
1917
1918 if (f.valid())
1919 {
1920 newTriangles[nNewTris++] = f;
1921 }
1922 }
1923 newTriangles.resize(nNewTris);
1924
1925 pointField newPoints(surf.localPoints(), uniquePoints);
1926
1927 return triSurface
1928 (
1929 newTriangles,
1930 surf.patches(),
1931 newPoints,
1932 true //reuse storage
1933 );
1934 }
1935
1936 return surf;
1937}
1938
1939
1940// Calculate normal on triangle
1942(
1943 const triSurface& surf,
1944 const label nearestFacei,
1945 const point& nearestPt
1946)
1947{
1948 const triSurface::face_type& f = surf[nearestFacei];
1949 const pointField& points = surf.points();
1950
1951 label nearType, nearLabel;
1952
1953 f.nearestPointClassify(nearestPt, points, nearType, nearLabel);
1954
1955 if (nearType == triPointRef::NONE)
1956 {
1957 // Nearest to face
1958 return surf.faceNormals()[nearestFacei];
1959 }
1960 else if (nearType == triPointRef::EDGE)
1961 {
1962 // Nearest to edge. Assume order of faceEdges same as face vertices.
1963 label edgeI = surf.faceEdges()[nearestFacei][nearLabel];
1964
1965 // Calculate edge normal by averaging face normals
1966 const labelList& eFaces = surf.edgeFaces()[edgeI];
1967
1968 vector edgeNormal(Zero);
1969
1970 for (const label facei : eFaces)
1971 {
1972 edgeNormal += surf.faceNormals()[facei];
1973 }
1974 return normalised(edgeNormal);
1975 }
1976 else
1977 {
1978 // Nearest to point
1979 const triSurface::face_type& localF = surf.localFaces()[nearestFacei];
1980 return surf.pointNormals()[localF[nearLabel]];
1981 }
1982}
1983
1984
1986(
1987 const triSurface& surf,
1988 const point& sample,
1989 const point& nearestPoint,
1990 const label edgeI
1991)
1992{
1993 const labelList& eFaces = surf.edgeFaces()[edgeI];
1994
1995 if (eFaces.size() != 2)
1996 {
1997 // Surface not closed.
1998 return UNKNOWN;
1999 }
2000 else
2001 {
2002 const vectorField& faceNormals = surf.faceNormals();
2003
2004 // Compare to bisector. This is actually correct since edge is
2005 // nearest so there is a knife-edge.
2006
2007 vector n = 0.5*(faceNormals[eFaces[0]] + faceNormals[eFaces[1]]);
2008
2009 if (((sample - nearestPoint) & n) > 0)
2010 {
2011 return OUTSIDE;
2012 }
2013 else
2014 {
2015 return INSIDE;
2016 }
2017 }
2018}
2019
2020
2021// Calculate normal on triangle
2023(
2024 const triSurface& surf,
2025 const point& sample,
2026 const label nearestFacei
2027)
2028{
2029 const triSurface::face_type& f = surf[nearestFacei];
2030 const pointField& points = surf.points();
2031
2032 // Find where point is on face
2033 label nearType, nearLabel;
2034
2035 pointHit pHit = f.nearestPointClassify(sample, points, nearType, nearLabel);
2036
2037 const point& nearestPoint(pHit.rawPoint());
2038
2039 if (nearType == triPointRef::NONE)
2040 {
2041 vector sampleNearestVec = (sample - nearestPoint);
2042
2043 // Nearest to face interior. Use faceNormal to determine side
2044 scalar c = sampleNearestVec & surf.faceNormals()[nearestFacei];
2045
2046 // // If the sample is essentially on the face, do not check for
2047 // // it being perpendicular.
2048
2049 // scalar magSampleNearestVec = mag(sampleNearestVec);
2050
2051 // if (magSampleNearestVec > SMALL)
2052 // {
2053 // c /= magSampleNearestVec*mag(surf.faceNormals()[nearestFacei]);
2054
2055 // if (mag(c) < 0.99)
2056 // {
2057 // FatalErrorInFunction
2058 // << "nearestPoint identified as being on triangle face "
2059 // << "but vector from nearestPoint to sample is not "
2060 // << "perpendicular to the normal." << nl
2061 // << "sample: " << sample << nl
2062 // << "nearestPoint: " << nearestPoint << nl
2063 // << "sample - nearestPoint: "
2064 // << sample - nearestPoint << nl
2065 // << "normal: " << surf.faceNormals()[nearestFacei] << nl
2066 // << "mag(sample - nearestPoint): "
2067 // << mag(sample - nearestPoint) << nl
2068 // << "normalised dot product: " << c << nl
2069 // << "triangle vertices: " << nl
2070 // << " " << points[f[0]] << nl
2071 // << " " << points[f[1]] << nl
2072 // << " " << points[f[2]] << nl
2073 // << abort(FatalError);
2074 // }
2075 // }
2076
2077 if (c > 0)
2078 {
2079 return OUTSIDE;
2080 }
2081 else
2082 {
2083 return INSIDE;
2084 }
2085 }
2086 else if (nearType == triPointRef::EDGE)
2087 {
2088 // Nearest to edge nearLabel. Note that this can only be a knife-edge
2089 // situation since otherwise the nearest point could never be the edge.
2090
2091 // Get the edge. Assume order of faceEdges same as face vertices.
2092 label edgeI = surf.faceEdges()[nearestFacei][nearLabel];
2093
2094 // if (debug)
2095 // {
2096 // // Check order of faceEdges same as face vertices.
2097 // const edge& e = surf.edges()[edgeI];
2098 // const labelList& meshPoints = surf.meshPoints();
2099 // const edge meshEdge(meshPoints[e[0]], meshPoints[e[1]]);
2100
2101 // if
2102 // (
2103 // meshEdge
2104 // != edge(f[nearLabel], f[f.fcIndex(nearLabel)])
2105 // )
2106 // {
2107 // FatalErrorInFunction
2108 // << "Edge:" << edgeI << " local vertices:" << e
2109 // << " mesh vertices:" << meshEdge
2110 // << " not at position " << nearLabel
2111 // << " in face " << f
2112 // << abort(FatalError);
2113 // }
2114 // }
2115
2116 return edgeSide(surf, sample, nearestPoint, edgeI);
2117 }
2118 else
2119 {
2120 // Nearest to point. Could use pointNormal here but is not correct.
2121 // Instead determine which edge using point is nearest and use test
2122 // above (nearType == triPointRef::EDGE).
2123
2124
2125 const triSurface::face_type& localF = surf.localFaces()[nearestFacei];
2126 label nearPointi = localF[nearLabel];
2127
2128 const edgeList& edges = surf.edges();
2129 const pointField& localPoints = surf.localPoints();
2130 const point& base = localPoints[nearPointi];
2131
2132 const labelList& pEdges = surf.pointEdges()[nearPointi];
2133
2134 scalar minDistSqr = Foam::sqr(GREAT);
2135 label minEdgeI = -1;
2136
2137 forAll(pEdges, i)
2138 {
2139 label edgeI = pEdges[i];
2140
2141 const edge& e = edges[edgeI];
2142
2143 label otherPointi = e.otherVertex(nearPointi);
2144
2145 // Get edge normal.
2146 vector eVec(localPoints[otherPointi] - base);
2147 scalar magEVec = mag(eVec);
2148
2149 if (magEVec > VSMALL)
2150 {
2151 eVec /= magEVec;
2152
2153 // Get point along vector and determine closest.
2154 const point perturbPoint = base + eVec;
2155
2156 scalar distSqr = Foam::magSqr(sample - perturbPoint);
2157
2158 if (distSqr < minDistSqr)
2159 {
2160 minDistSqr = distSqr;
2161 minEdgeI = edgeI;
2162 }
2163 }
2164 }
2165
2166 if (minEdgeI == -1)
2167 {
2169 << "Problem: did not find edge closer than " << minDistSqr
2170 << abort(FatalError);
2171 }
2172
2173 return edgeSide(surf, sample, nearestPoint, minEdgeI);
2174 }
2175}
2176
2177
2179(
2180 const polyBoundaryMesh& bMesh,
2181 const labelHashSet& includePatches,
2183 const bool verbose
2184)
2185{
2186 const polyMesh& mesh = bMesh.mesh();
2187
2188 // Storage for surfaceMesh. Size estimate.
2189 List<labelledTri> triangles;
2190
2191 // Calculate number of triangles
2192 label nTris = 0;
2193
2194 for (const label patchi : includePatches)
2195 {
2196 const polyPatch& patch = bMesh[patchi];
2197 const pointField& points = patch.points();
2198 for (const face& f : patch)
2199 {
2200 nTris += f.nTriangles(points);
2201 }
2202 }
2203
2204 triangles.setSize(nTris);
2205 faceMap.setSize(nTris);
2206 label newPatchi = 0;
2207
2208 nTris = 0;
2209 for (const label patchi : includePatches)
2210 {
2211 const polyPatch& patch = bMesh[patchi];
2212 const pointField& points = patch.points();
2213
2214 label nTriTotal = 0;
2215
2216 label faceI = 0;
2217 for (const face& f : patch)
2218 {
2219 faceList triFaces(f.nTriangles(points));
2220
2221 label nTri = 0;
2222
2223 f.triangles(points, nTri, triFaces);
2224
2225 for (const face& f : triFaces)
2226 {
2227 faceMap[nTris] = patch.start() + faceI;
2228 triangles[nTris++] = labelledTri(f[0], f[1], f[2], newPatchi);
2229
2230 ++nTriTotal;
2231 }
2232
2233 faceI++;
2234 }
2235
2236 if (verbose)
2237 {
2238 Pout<< patch.name() << " : generated " << nTriTotal
2239 << " triangles from " << patch.size() << " faces with"
2240 << " new patchid " << newPatchi << endl;
2241 }
2242
2243 newPatchi++;
2244 }
2245 //triangles.shrink();
2246
2247 // Create globally numbered tri surface
2248 triSurface rawSurface(triangles, mesh.points());
2249
2250 // Create locally numbered tri surface
2251 triSurface surface
2252 (
2253 rawSurface.localFaces(),
2254 rawSurface.localPoints()
2255 );
2256
2257 // Add patch names to surface
2258 surface.patches().setSize(newPatchi);
2259
2260 newPatchi = 0;
2261
2262 for (const label patchi : includePatches)
2263 {
2264 const polyPatch& patch = bMesh[patchi];
2265
2266 surface.patches()[newPatchi].name() = patch.name();
2267 surface.patches()[newPatchi].geometricType() = patch.type();
2268
2269 newPatchi++;
2270 }
2271
2272 return surface;
2273}
2274
2275
2277(
2278 const polyBoundaryMesh& bMesh,
2279 const labelHashSet& includePatches,
2280 const boundBox& bBox,
2281 const bool verbose
2282)
2283{
2284 const polyMesh& mesh = bMesh.mesh();
2285
2286 // Storage for surfaceMesh. Size estimate.
2288
2289 label newPatchi = 0;
2290
2291 for (const label patchi : includePatches)
2292 {
2293 const polyPatch& patch = bMesh[patchi];
2294 const pointField& points = patch.points();
2295
2296 label nTriTotal = 0;
2297
2298 forAll(patch, patchFacei)
2299 {
2300 const face& f = patch[patchFacei];
2301
2302 if (bBox.containsAny(points, f))
2303 {
2304 faceList triFaces(f.nTriangles(points));
2305
2306 label nTri = 0;
2307
2308 f.triangles(points, nTri, triFaces);
2309
2310 forAll(triFaces, triFacei)
2311 {
2312 const face& f = triFaces[triFacei];
2313
2314 triangles.append(labelledTri(f[0], f[1], f[2], newPatchi));
2315
2316 nTriTotal++;
2317 }
2318 }
2319 }
2320
2321 if (verbose)
2322 {
2323 Pout<< patch.name() << " : generated " << nTriTotal
2324 << " triangles from " << patch.size() << " faces with"
2325 << " new patchid " << newPatchi << endl;
2326 }
2327
2328 newPatchi++;
2329 }
2330 triangles.shrink();
2331
2332 // Create globally numbered tri surface
2333 triSurface rawSurface(triangles, mesh.points());
2334
2335 // Create locally numbered tri surface
2336 triSurface surface
2337 (
2338 rawSurface.localFaces(),
2339 rawSurface.localPoints()
2340 );
2341
2342 // Add patch names to surface
2343 surface.patches().setSize(newPatchi);
2344
2345 newPatchi = 0;
2346
2347 for (const label patchi : includePatches)
2348 {
2349 const polyPatch& patch = bMesh[patchi];
2350
2351 surface.patches()[newPatchi].name() = patch.name();
2352 surface.patches()[newPatchi].geometricType() = patch.type();
2353
2354 ++newPatchi;
2355 }
2356
2357 return surface;
2358}
2359
2360
2361// triangulation of boundaryMesh
2363(
2364 const polyBoundaryMesh& bMesh,
2365 const labelHashSet& includePatches,
2366 const bool verbose
2367)
2368{
2369 const polyMesh& mesh = bMesh.mesh();
2370
2371 // Storage for new points = meshpoints + face centres.
2372 const pointField& points = mesh.points();
2373 const pointField& faceCentres = mesh.faceCentres();
2374
2375 pointField newPoints(points.size() + faceCentres.size());
2376
2377 label newPointi = 0;
2378
2379 forAll(points, pointi)
2380 {
2381 newPoints[newPointi++] = points[pointi];
2382 }
2383 forAll(faceCentres, facei)
2384 {
2385 newPoints[newPointi++] = faceCentres[facei];
2386 }
2387
2388
2389 // Count number of faces.
2391
2392 label newPatchi = 0;
2393
2394 for (const label patchi : includePatches)
2395 {
2396 const polyPatch& patch = bMesh[patchi];
2397
2398 label nTriTotal = 0;
2399
2400 forAll(patch, patchFacei)
2401 {
2402 // Face in global coords.
2403 const face& f = patch[patchFacei];
2404
2405 // Index in newPointi of face centre.
2406 label fc = points.size() + patchFacei + patch.start();
2407
2408 forAll(f, fp)
2409 {
2410 label fp1 = f.fcIndex(fp);
2411
2412 triangles.append(labelledTri(f[fp], f[fp1], fc, newPatchi));
2413
2414 nTriTotal++;
2415 }
2416 }
2417
2418 if (verbose)
2419 {
2420 Pout<< patch.name() << " : generated " << nTriTotal
2421 << " triangles from " << patch.size() << " faces with"
2422 << " new patchid " << newPatchi << endl;
2423 }
2424
2425 newPatchi++;
2426 }
2427 triangles.shrink();
2428
2429
2430 // Create globally numbered tri surface
2431 triSurface rawSurface(triangles, newPoints);
2432
2433 // Create locally numbered tri surface
2434 triSurface surface
2435 (
2436 rawSurface.localFaces(),
2437 rawSurface.localPoints()
2438 );
2439
2440 // Add patch names to surface
2441 surface.patches().setSize(newPatchi);
2442
2443 newPatchi = 0;
2444
2445 for (const label patchi : includePatches)
2446 {
2447 const polyPatch& patch = bMesh[patchi];
2448
2449 surface.patches()[newPatchi].name() = patch.name();
2450 surface.patches()[newPatchi].geometricType() = patch.type();
2451
2452 newPatchi++;
2453 }
2454
2455 return surface;
2456}
2457
2458
2460{
2461 // Vertices in geompack notation. Note that could probably just use
2462 // pts.begin() if double precision.
2463 List<doubleScalar> geompackVertices(2*pts.size());
2464 label doubleI = 0;
2465 for (const vector2D& pt : pts)
2466 {
2467 geompackVertices[doubleI++] = pt[0];
2468 geompackVertices[doubleI++] = pt[1];
2469 }
2470
2471 // Storage for triangles
2472 int m2 = 3;
2473 List<int> triangle_node(m2*3*pts.size());
2474 List<int> triangle_neighbor(m2*3*pts.size());
2475
2476 // Triangulate
2477 int nTris = 0;
2478 int err = dtris2
2479 (
2480 pts.size(),
2481 geompackVertices.begin(),
2482 &nTris,
2483 triangle_node.begin(),
2484 triangle_neighbor.begin()
2485 );
2486
2487 if (err != 0)
2488 {
2490 << "Failed dtris2 with vertices:" << pts.size()
2491 << abort(FatalError);
2492 }
2493
2494 // Trim
2495 triangle_node.setSize(3*nTris);
2496 triangle_neighbor.setSize(3*nTris);
2497
2498 // Convert to triSurface.
2499 List<labelledTri> faces(nTris);
2500
2501 forAll(faces, i)
2502 {
2503 faces[i] = labelledTri
2504 (
2505 triangle_node[3*i]-1,
2506 triangle_node[3*i+1]-1,
2507 triangle_node[3*i+2]-1,
2508 0
2509 );
2510 }
2511
2512 pointField points(pts.size());
2513 forAll(pts, i)
2514 {
2515 points[i][0] = pts[i][0];
2516 points[i][1] = pts[i][1];
2517 points[i][2] = 0.0;
2518 }
2519
2520 return triSurface(faces, points);
2521}
2522
2523
2525(
2526 const triPointRef& tri,
2527 const point& p,
2528 FixedList<scalar, 3>& weights
2529)
2530{
2531 // calculate triangle edge vectors and triangle face normal
2532 // the 'i':th edge is opposite node i
2534 edge[0] = tri.c()-tri.b();
2535 edge[1] = tri.a()-tri.c();
2536 edge[2] = tri.b()-tri.a();
2537
2538 const vector triangleFaceNormal = edge[1] ^ edge[2];
2539
2540 // calculate edge normal (pointing inwards)
2541 FixedList<vector, 3> normal;
2542 for (label i=0; i<3; i++)
2543 {
2544 normal[i] = normalised(triangleFaceNormal ^ edge[i]);
2545 }
2546
2547 weights[0] = ((p-tri.b()) & normal[0]) / max(VSMALL, normal[0] & edge[1]);
2548 weights[1] = ((p-tri.c()) & normal[1]) / max(VSMALL, normal[1] & edge[2]);
2549 weights[2] = ((p-tri.a()) & normal[2]) / max(VSMALL, normal[2] & edge[0]);
2550}
2551
2552
2553// Calculate weighting factors from samplePts to triangle it is in.
2554// Uses linear search.
2556(
2557 const triSurface& s,
2558 const pointField& samplePts,
2559 List<FixedList<label, 3>>& allVerts,
2560 List<FixedList<scalar, 3>>& allWeights
2561)
2562{
2563 allVerts.setSize(samplePts.size());
2564 allWeights.setSize(samplePts.size());
2565
2566 const pointField& points = s.points();
2567
2568 forAll(samplePts, i)
2569 {
2570 const point& samplePt = samplePts[i];
2571
2572 FixedList<label, 3>& verts = allVerts[i];
2573 FixedList<scalar, 3>& weights = allWeights[i];
2574
2575 scalar minDistance = GREAT;
2576
2577 for (const labelledTri& f : s)
2578 {
2579 triPointRef tri(f.tri(points));
2580
2581 label nearType, nearLabel;
2582
2583 pointHit nearest = tri.nearestPointClassify
2584 (
2585 samplePt,
2586 nearType,
2587 nearLabel
2588 );
2589
2590 if (nearest.hit())
2591 {
2592 // samplePt inside triangle
2593 verts[0] = f[0];
2594 verts[1] = f[1];
2595 verts[2] = f[2];
2596
2597 calcInterpolationWeights(tri, nearest.rawPoint(), weights);
2598
2599 //Pout<< "calcScalingFactors : samplePt:" << samplePt
2600 // << " inside triangle:" << facei
2601 // << " verts:" << verts
2602 // << " weights:" << weights
2603 // << endl;
2604
2605 break;
2606 }
2607 else if (nearest.distance() < minDistance)
2608 {
2609 minDistance = nearest.distance();
2610
2611 // Outside triangle. Store nearest.
2612
2613 if (nearType == triPointRef::POINT)
2614 {
2615 verts[0] = f[nearLabel];
2616 weights[0] = 1;
2617 verts[1] = -1;
2618 weights[1] = -GREAT;
2619 verts[2] = -1;
2620 weights[2] = -GREAT;
2621
2622 //Pout<< "calcScalingFactors : samplePt:" << samplePt
2623 // << " distance:" << nearest.distance()
2624 // << " from point:" << points[f[nearLabel]]
2625 // << endl;
2626 }
2627 else if (nearType == triPointRef::EDGE)
2628 {
2629 verts[0] = f[nearLabel];
2630 verts[1] = f[f.fcIndex(nearLabel)];
2631 verts[2] = -1;
2632
2633 const point& p0 = points[verts[0]];
2634 const point& p1 = points[verts[1]];
2635
2636 scalar s = min
2637 (
2638 1,
2639 max
2640 (
2641 0,
2642 mag(nearest.rawPoint() - p0)/mag(p1 - p0)
2643 )
2644 );
2645
2646 // Interpolate
2647 weights[0] = 1 - s;
2648 weights[1] = s;
2649 weights[2] = -GREAT;
2650
2651 //Pout<< "calcScalingFactors : samplePt:" << samplePt
2652 // << " distance:" << nearest.distance()
2653 // << " from edge:" << p0 << p1 << " s:" << s
2654 // << endl;
2655 }
2656 else
2657 {
2658 // triangle. Can only happen because of truncation errors.
2659 verts[0] = f[0];
2660 verts[1] = f[1];
2661 verts[2] = f[2];
2662
2663 calcInterpolationWeights(tri, nearest.rawPoint(), weights);
2664
2665 //Pout<< "calcScalingFactors : samplePt:" << samplePt
2666 // << " distance:" << nearest.distance()
2667 // << " to verts:" << verts
2668 // << " weights:" << weights
2669 // << endl;
2670
2671 break;
2672 }
2673 }
2674 }
2675 }
2676}
2677
2678
2679// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
2680// Checking:
2681
2683(
2684 const triSurface& surf,
2685 const label facei,
2686 const bool verbose
2687)
2688{
2689 typedef labelledTri FaceType;
2690 const FaceType& f = surf[facei];
2691
2692 // Simple check on indices ok.
2693 for (const label pointi : f)
2694 {
2695 if (pointi < 0 || pointi >= surf.points().size())
2696 {
2697 if (verbose)
2698 {
2700 << "triangle " << facei << " vertices " << f
2701 << " uses point indices outside point range 0.."
2702 << surf.points().size()-1 << endl;
2703 }
2704 return false;
2705 }
2706 }
2707
2708 if (f[0] == f[1] || f[0] == f[2] || f[1] == f[2])
2709 {
2710 if (verbose)
2711 {
2713 << "triangle " << facei
2714 << " uses non-unique vertices " << f
2715 << " coords:" << f.points(surf.points()) << endl;
2716 }
2717 return false;
2718 }
2719
2720 // duplicate triangle check
2721
2722 const labelList& fFaces = surf.faceFaces()[facei];
2723
2724 // Check if faceNeighbours use same points as this face.
2725 // Note: discards normal information - sides of baffle are merged.
2726 for (const label nbrFacei : fFaces)
2727 {
2728 if (nbrFacei <= facei)
2729 {
2730 // lower numbered faces already checked
2731 continue;
2732 }
2733
2734 const FaceType& nbrF = surf[nbrFacei];
2735
2736 // Same as calling triFace::compare(f, nbrF) == 1 only
2737 if
2738 (
2739 (f[0] == nbrF[0] || f[0] == nbrF[1] || f[0] == nbrF[2])
2740 && (f[1] == nbrF[0] || f[1] == nbrF[1] || f[1] == nbrF[2])
2741 && (f[2] == nbrF[0] || f[2] == nbrF[1] || f[2] == nbrF[2])
2742 )
2743 {
2744 if (verbose)
2745 {
2747 << "triangle " << facei << " vertices " << f
2748 << " has the same vertices as triangle " << nbrFacei
2749 << " vertices " << nbrF
2750 << " coords:" << f.points(surf.points()) << endl;
2751 }
2752
2753 return false;
2754 }
2755 }
2756
2757 return true;
2758}
2759
2760
2762(
2763 const MeshedSurface<face>& surf,
2764 const label facei,
2765 const bool verbose
2766)
2767{
2768 typedef face FaceType;
2769 const FaceType& f = surf[facei];
2770
2771 if (f.size() != 3)
2772 {
2773 if (verbose)
2774 {
2776 << "face " << facei
2777 << " is not a triangle, it has " << f.size()
2778 << " indices" << endl;
2779 }
2780 return false;
2781 }
2782
2783 // Simple check on indices ok.
2784 for (const label pointi : f)
2785 {
2786 if (pointi < 0 || pointi >= surf.points().size())
2787 {
2788 if (verbose)
2789 {
2791 << "triangle " << facei << " vertices " << f
2792 << " uses point indices outside point range 0.."
2793 << surf.points().size()-1 << endl;
2794 }
2795 return false;
2796 }
2797 }
2798
2799 if (f[0] == f[1] || f[0] == f[2] || f[1] == f[2])
2800 {
2801 if (verbose)
2802 {
2804 << "triangle " << facei
2805 << " uses non-unique vertices " << f
2806 << " coords:" << f.points(surf.points()) << endl;
2807 }
2808 return false;
2809 }
2810
2811 // duplicate triangle check
2812
2813 const labelList& fFaces = surf.faceFaces()[facei];
2814
2815 // Check if faceNeighbours use same points as this face.
2816 // Note: discards normal information - sides of baffle are merged.
2817 for (const label nbrFacei : fFaces)
2818 {
2819 if (nbrFacei <= facei)
2820 {
2821 // lower numbered faces already checked
2822 continue;
2823 }
2824
2825 const FaceType& nbrF = surf[nbrFacei];
2826
2827 // Same as calling triFace::compare(f, nbrF) == 1 only
2828 if
2829 (
2830 (f[0] == nbrF[0] || f[0] == nbrF[1] || f[0] == nbrF[2])
2831 && (f[1] == nbrF[0] || f[1] == nbrF[1] || f[1] == nbrF[2])
2832 && (f[2] == nbrF[0] || f[2] == nbrF[1] || f[2] == nbrF[2])
2833 )
2834 {
2835 if (verbose)
2836 {
2838 << "triangle " << facei << " vertices " << f
2839 << " has the same vertices as triangle " << nbrFacei
2840 << " vertices " << nbrF
2841 << " coords:" << f.points(surf.points()) << endl;
2842 }
2843 return false;
2844 }
2845 }
2846
2847 return true;
2848}
2849
2850
2851// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
2852// Tracking:
2853
2854// Test point on surface to see if is on face,edge or point.
2856(
2857 const triSurface& s,
2858 const label triI,
2859 const point& trianglePoint
2860)
2861{
2862 surfaceLocation nearest;
2863
2864 // Nearest point could be on point or edge. Retest.
2865 label index, elemType;
2866 //bool inside =
2867 triPointRef(s[triI].tri(s.points())).classify
2868 (
2869 trianglePoint,
2870 elemType,
2871 index
2872 );
2873
2874 nearest.setPoint(trianglePoint);
2875
2876 if (elemType == triPointRef::NONE)
2877 {
2878 nearest.setHit();
2879 nearest.setIndex(triI);
2880 nearest.elementType() = triPointRef::NONE;
2881 }
2882 else if (elemType == triPointRef::EDGE)
2883 {
2884 nearest.setMiss();
2885 nearest.setIndex(s.faceEdges()[triI][index]);
2886 nearest.elementType() = triPointRef::EDGE;
2887 }
2888 else //if (elemType == triPointRef::POINT)
2889 {
2890 nearest.setMiss();
2891 nearest.setIndex(s.localFaces()[triI][index]);
2892 nearest.elementType() = triPointRef::POINT;
2893 }
2894
2895 return nearest;
2896}
2897
2898
2900(
2901 const triSurface& s,
2902 const surfaceLocation& start,
2903 const surfaceLocation& end,
2904 const plane& cutPlane
2905)
2906{
2907 // Start off from starting point
2908 surfaceLocation nearest = start;
2909 nearest.setMiss();
2910
2911 // See if in same triangle as endpoint. If so snap.
2912 snapToEnd(s, end, nearest);
2913
2914 if (!nearest.hit())
2915 {
2916 // Not yet at end point
2917
2918 if (start.elementType() == triPointRef::NONE)
2919 {
2920 // Start point is inside triangle. Trivial cases already handled
2921 // above.
2922
2923 // end point is on edge or point so cross current triangle to
2924 // see which edge is cut.
2925
2926 nearest = cutEdge
2927 (
2928 s,
2929 start.index(), // triangle
2930 -1, // excludeEdge
2931 -1, // excludePoint
2932 start.rawPoint(),
2933 cutPlane,
2934 end.rawPoint()
2935 );
2936 nearest.elementType() = triPointRef::EDGE;
2937 nearest.triangle() = start.index();
2938 nearest.setMiss();
2939 }
2940 else if (start.elementType() == triPointRef::EDGE)
2941 {
2942 // Pick connected triangle that is most in direction.
2943 const labelList& eFaces = s.edgeFaces()[start.index()];
2944
2945 nearest = visitFaces
2946 (
2947 s,
2948 eFaces,
2949 start,
2950 start.index(), // excludeEdgeI
2951 -1, // excludePointi
2952 end,
2953 cutPlane
2954 );
2955 }
2956 else // start.elementType() == triPointRef::POINT
2957 {
2958 const labelList& pFaces = s.pointFaces()[start.index()];
2959
2960 nearest = visitFaces
2961 (
2962 s,
2963 pFaces,
2964 start,
2965 -1, // excludeEdgeI
2966 start.index(), // excludePointi
2967 end,
2968 cutPlane
2969 );
2970 }
2971 snapToEnd(s, end, nearest);
2972 }
2973 return nearest;
2974}
2975
2976
2978(
2979 const triSurface& s,
2980 const surfaceLocation& endInfo,
2981 const plane& cutPlane,
2982 surfaceLocation& hitInfo
2983)
2984{
2985 //OFstream str("track.obj");
2986 //label vertI = 0;
2987 //meshTools::writeOBJ(str, hitInfo.rawPoint());
2988 //vertI++;
2989
2990 // Track across surface.
2991 while (true)
2992 {
2993 //Pout<< "Tracking from:" << nl
2994 // << " " << hitInfo.info()
2995 // << endl;
2996
2997 hitInfo = trackToEdge
2998 (
2999 s,
3000 hitInfo,
3001 endInfo,
3002 cutPlane
3003 );
3004
3005 //meshTools::writeOBJ(str, hitInfo.rawPoint());
3006 //vertI++;
3007 //str<< "l " << vertI-1 << ' ' << vertI << nl;
3008
3009 //Pout<< "Tracked to:" << nl
3010 // << " " << hitInfo.info() << endl;
3011
3012 if (hitInfo.hit() || hitInfo.triangle() == -1)
3013 {
3014 break;
3015 }
3016 }
3017}
3018
3019
3020// ************************************************************************* //
static const Foam::dimensionedScalar A("", Foam::dimPressure, 611.21)
static const Foam::dimensionedScalar B("", Foam::dimless, 18.678)
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 1D vector of objects of type <T> with a fixed length <N>.
Definition: FixedList.H:81
Minimal example by using system/controlDict.functions:
bool insert(const Key &key)
Insert a new entry, not overwriting existing entries.
Definition: HashSet.H:191
List< Key > toc() const
The table of contents (the keys) in unsorted order.
Definition: HashTable.C:122
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 resize(const label len)
Adjust allocated size of list.
Definition: ListI.H:139
A surface geometry mesh with zone information, not to be confused with the similarly named surfaceMes...
Definition: MeshedSurface.H:99
virtual label triangulate()
Triangulate in-place, returning the number of triangles added.
Output to file stream, using an OSstream.
Definition: OFstream.H:57
Describes the interaction of a face and a point. It carries the info of a successful hit and (if succ...
Definition: PointHit.H:54
const point_type & rawPoint() const noexcept
The point, no checks.
Definition: PointHit.H:172
scalar distance() const noexcept
Return distance to hit.
Definition: PointHit.H:139
bool hit() const noexcept
Is there a hit.
Definition: PointHit.H:121
void setHit() noexcept
Set the hit status on.
const point_type & rawPoint() const noexcept
The point, no checks. Same as point()
void setIndex(const label index) noexcept
Set the index.
void setPoint(const point_type &p)
Set the point.
label index() const noexcept
Return the hit index.
void setMiss() noexcept
Set the hit status off.
bool hit() const noexcept
Is there a hit?
A list of faces which address into the list of points.
label nEdges() const
Number of edges in patch.
label nPoints() const
Number of points supporting patch faces.
const labelListList & pointEdges() const
Return point-edge addressing.
const edgeList & edges() const
Return list of edges, address into LOCAL point list.
const Field< point_type > & localPoints() const
Return pointField of points in patch.
const Field< point_type > & faceNormals() const
Return face unit normals for patch.
const Field< point_type > & pointNormals() const
Return point normals for patch.
const Field< point_type > & points() const noexcept
Return reference to global points.
const labelListList & faceFaces() const
Return face-face addressing.
const labelListList & pointFaces() const
Return point-face addressing.
const labelListList & edgeFaces() const
Return edge-face addressing.
const labelListList & faceEdges() const
Return face-edge addressing.
const List< face_type > & localFaces() const
Return patch faces addressing into local point list.
bool found(const T &val, label pos=0) const
True if the value if found in the list.
Definition: UListI.H:265
iterator begin() noexcept
Return an iterator to begin traversing the UList.
Definition: UListI.H:329
iterator end() noexcept
Return an iterator to end traversing the UList.
Definition: UListI.H:350
label find(const T &val, label pos=0) const
Find index of the first occurrence of the value.
Definition: UList.C:212
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
label fcIndex(const label i) const noexcept
Definition: UListI.H:60
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
A bounding box defined in terms of min/max extrema points.
Definition: boundBox.H:64
bool containsAny(const UList< point > &points) const
Contains any of the points? (inside or on edge)
Definition: boundBox.C:249
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
virtual void track()
Do the actual tracking to fill the track data.
Definition: streamLine.C:48
A triFace with additional (region) index.
Definition: labelledTri.H:60
Geometric class that creates a 3D plane and can return the intersection point between a line and the ...
Definition: plane.H:95
A polyBoundaryMesh is a polyPatch list with additional search methods and registered IO.
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:81
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
label nBoundaryFaces() const noexcept
Number of boundary faces (== nFaces - nInternalFaces)
const vectorField & faceCentres() const
Patchify triangles based on orientation w.r.t other (triangulated or triangulatable) surfaces.
Contains information about location on a triSurface.
triPointRef::proxType & elementType() noexcept
label & triangle() noexcept
sideType
On which side of surface.
static const label COLLAPSED
static label oppositeVertex(const triSurface &surf, const label facei, const label edgeI)
Get vertex (local numbering) opposite edge.
static bool validTri(const triSurface &, const label facei, const bool verbose=true)
Check single triangle for (topological) validity.
static triSurface redGreenRefine(const triSurface &surf, const labelList &refineFaces)
Refine face by splitting all edges. Neighbouring face is.
static triSurface collapseEdges(const triSurface &surf, const labelList &collapsableEdges)
Create new triSurface by collapsing edges to edge mids.
static label getTriangle(const triSurface &surf, const label e0I, const label e1I, const label e2I)
Return index of triangle (or -1) using all three edges.
static void otherVertices(const triSurface &surf, const label facei, const label vertI, label &vert1I, label &vert2I)
Get the two vertices (local numbering) on facei counterclockwise.
static void calcInterpolationWeights(const triPointRef &tri, const point &p, FixedList< scalar, 3 > &weights)
Calculate linear interpolation weights for point (guaranteed to be.
static triSurface mergePoints(const triSurface &surf, const scalar mergeTol)
Merge points within distance.
static const label ANYEDGE
Face collapse status.
static label getEdge(const triSurface &surf, const label vert1I, const label vert2I)
Returns edge label connecting v1, v2 (local numbering)
static void writeOBJ(const fileName &fName, const pointField &pts)
Write pointField to OBJ format file.
static surfaceLocation trackToEdge(const triSurface &, const surfaceLocation &start, const surfaceLocation &end, const plane &cutPlane)
Track on surface to get closer to point.
static surfaceLocation classify(const triSurface &, const label triI, const point &trianglePoint)
Test point on plane of triangle to see if on edge or point or inside.
static vector surfaceNormal(const triSurface &surf, const label nearestFacei, const point &nearestPt)
Triangle (unit) normal. If nearest point to triangle on edge use.
static label otherFace(const triSurface &surf, const label facei, const label edgeI)
Get face connected to edge not facei.
static void getVertexTriangles(const triSurface &surf, const label edgeI, labelList &edgeTris)
Get all triangles using edge endpoint.
triSurface triangulateFaceCentre(const polyBoundaryMesh &mBesh, const labelHashSet &includePatches, const bool verbose=false)
Face-centre triangulation of (selected patches of) boundaryMesh.
static const label NOEDGE
static label minEdge(const triSurface &surf, const labelList &edgeIndices)
Returns element in edgeIndices with minimum length.
static sideType edgeSide(const triSurface &surf, const point &sample, const point &nearestPoint, const label edgeI)
If nearest point is on edgeI, determine on which side of surface.
static label oppositeEdge(const triSurface &surf, const label facei, const label vertI)
Get edge opposite vertex (local numbering)
static sideType surfaceSide(const triSurface &surf, const point &sample, const label nearestFacei)
Given nearest point (to sample) on surface determines which side.
static label maxEdge(const triSurface &surf, const labelList &edgeIndices)
Returns element in edgeIndices with minimum length.
static void otherEdges(const triSurface &surf, const label facei, const label edgeI, label &e1, label &e2)
Get the two edges on facei counterclockwise after edgeI.
static labelList getVertexVertices(const triSurface &surf, const edge &e)
Get all vertices (local numbering) connected to vertices of edge.
static triSurface delaunay2D(const List< vector2D > &)
Do unconstrained Delaunay of points. Returns triSurface with 3D.
Triangulated surface description with patch information.
Definition: triSurface.H:79
labelledTri face_type
The face type (same as the underlying PrimitivePatch)
Definition: triSurface.H:209
const geometricSurfacePatchList & patches() const noexcept
Definition: triSurface.H:399
A triangle primitive used to calculate face normals and swept volumes.
Definition: triangle.H:80
const Point & a() const
Return first vertex.
Definition: triangleI.H:86
pointHit nearestPointClassify(const point &p, label &nearType, label &nearLabel) const
Find the nearest point to p on the triangle and classify it:
Definition: triangleI.H:522
@ POINT
Close to point.
Definition: triangle.H:96
@ EDGE
Close to edge.
Definition: triangle.H:97
@ NONE
Unknown proximity.
Definition: triangle.H:95
const Point & c() const
Return third vertex.
Definition: triangleI.H:98
const Point & b() const
Return second vertex.
Definition: triangleI.H:92
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
volScalarField & p
const volScalarField & p0
Definition: EEqn.H:36
dynamicFvMesh & mesh
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
int dtris2(int point_num, double point_xy[], int *tri_num, int tri_vert[], int tri_nabe[])
Definition: geompack.C:949
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))
surfaceVectorField faceNormals(mesh.Sf()/mesh.magSf())
Geometric merging of points. See below.
#define WarningInFunction
Report a warning using Foam::Warning.
tmp< pointField > allPoints(const Triangulation &t)
Extract all points in vertex-index order.
label otherFace(const primitiveMesh &mesh, const label celli, const label facei, const label edgeI)
Return face on cell using edgeI but not facei. Throws error.
Definition: meshTools.C:555
Pair< int > faceMap(const label facePi, const face &faceP, const label faceNi, const face &faceN)
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:47
labelList identity(const label len, label start=0)
Return an identity map of the given length with (map[i] == i)
Definition: labelList.C:38
List< label > labelList
A List of labels.
Definition: List.H:66
dimensionedSymmTensor sqr(const dimensionedVector &dv)
quaternion normalised(const quaternion &q)
Return the normalised (unit) quaternion of the given quaternion.
Definition: quaternionI.H:680
PointFrompoint toPoint(const Foam::point &p)
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:44
vector point
Point is a vector.
Definition: point.H:43
label mergePoints(const PointList &points, labelList &pointToUnique, labelList &uniquePoints, const scalar mergeTol=SMALL, const bool verbose=false)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
triangle< point, const point & > triPointRef
A triangle using referred points.
Definition: triangle.H:71
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33
errorManip< error > abort(error &err)
Definition: errorManip.H:144
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
HashSet< label, Hash< label > > labelHashSet
A HashSet of labels, uses label hasher.
Definition: HashSet.H:85
error FatalError
List< edge > edgeList
A List of edges.
Definition: edgeList.H:63
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
constexpr auto end(C &c) -> decltype(c.end())
Return iterator to the end of the container c.
Definition: stdFoam.H:158
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
label newPointi
Definition: readKivaGrid.H:496
labelList f(nPoints)
volScalarField & C
const dimensionedScalar & D
volScalarField & e
Definition: createFields.H:11
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333