00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028 #include <meshTools/triSurfaceTools.H>
00029
00030 #include <triSurface/triSurface.H>
00031 #include <OpenFOAM/OFstream.H>
00032 #include <OpenFOAM/mergePoints.H>
00033 #include <OpenFOAM/polyMesh.H>
00034 #include <OpenFOAM/plane.H>
00035 #include <meshTools/geompack.H>
00036
00037
00038
00039
00040 const Foam::label Foam::triSurfaceTools::ANYEDGE = -1;
00041 const Foam::label Foam::triSurfaceTools::NOEDGE = -2;
00042 const Foam::label Foam::triSurfaceTools::COLLAPSED = -3;
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055 void Foam::triSurfaceTools::calcRefineStatus
00056 (
00057 const triSurface& surf,
00058 const label faceI,
00059 List<refineType>& refine
00060 )
00061 {
00062 if (refine[faceI] == RED)
00063 {
00064
00065 }
00066 else
00067 {
00068
00069 refine[faceI] = RED;
00070
00071 const labelList& myNeighbours = surf.faceFaces()[faceI];
00072
00073 forAll(myNeighbours, myNeighbourI)
00074 {
00075 label neighbourFaceI = myNeighbours[myNeighbourI];
00076
00077 if (refine[neighbourFaceI] == GREEN)
00078 {
00079
00080 calcRefineStatus(surf, neighbourFaceI, refine);
00081 }
00082 else if (refine[neighbourFaceI] == NONE)
00083 {
00084 refine[neighbourFaceI] = GREEN;
00085 }
00086 }
00087 }
00088 }
00089
00090
00091
00092 void Foam::triSurfaceTools::greenRefine
00093 (
00094 const triSurface& surf,
00095 const label faceI,
00096 const label edgeI,
00097 const label newPointI,
00098 DynamicList<labelledTri>& newFaces
00099 )
00100 {
00101 const labelledTri& f = surf.localFaces()[faceI];
00102 const edge& e = surf.edges()[edgeI];
00103
00104
00105
00106 label fp0 = findIndex(f, e[0]);
00107 label fp1 = f.fcIndex(fp0);
00108 label fp2 = f.fcIndex(fp1);
00109
00110 if (f[fp1] == e[1])
00111 {
00112
00113 newFaces.append
00114 (
00115 labelledTri
00116 (
00117 f[fp0],
00118 newPointI,
00119 f[fp2],
00120 f.region()
00121 )
00122 );
00123 newFaces.append
00124 (
00125 labelledTri
00126 (
00127 newPointI,
00128 f[fp1],
00129 f[fp2],
00130 f.region()
00131 )
00132 );
00133 }
00134 else
00135 {
00136 newFaces.append
00137 (
00138 labelledTri
00139 (
00140 f[fp2],
00141 newPointI,
00142 f[fp1],
00143 f.region()
00144 )
00145 );
00146 newFaces.append
00147 (
00148 labelledTri
00149 (
00150 newPointI,
00151 f[fp0],
00152 f[fp1],
00153 f.region()
00154 )
00155 );
00156 }
00157 }
00158
00159
00160
00161 Foam::triSurface Foam::triSurfaceTools::doRefine
00162 (
00163 const triSurface& surf,
00164 const List<refineType>& refineStatus
00165 )
00166 {
00167
00168 DynamicList<point> newPoints(surf.nPoints());
00169 forAll(surf.localPoints(), pointI)
00170 {
00171 newPoints.append(surf.localPoints()[pointI]);
00172 }
00173 label newVertI = surf.nPoints();
00174
00175
00176 DynamicList<labelledTri> newFaces(surf.size());
00177
00178
00179
00180 labelList edgeMid(surf.nEdges(), -1);
00181
00182 forAll(refineStatus, faceI)
00183 {
00184 if (refineStatus[faceI] == RED)
00185 {
00186
00187 const labelList& fEdges = surf.faceEdges()[faceI];
00188
00189 forAll(fEdges, i)
00190 {
00191 label edgeI = fEdges[i];
00192
00193 if (edgeMid[edgeI] == -1)
00194 {
00195 const edge& e = surf.edges()[edgeI];
00196
00197
00198 newPoints.append
00199 (
00200 0.5
00201 * (
00202 surf.localPoints()[e.start()]
00203 + surf.localPoints()[e.end()]
00204 )
00205 );
00206 edgeMid[edgeI] = newVertI++;
00207 }
00208 }
00209
00210
00211
00212
00213 const edgeList& edges = surf.edges();
00214
00215
00216 newFaces.append
00217 (
00218 labelledTri
00219 (
00220 edgeMid[fEdges[0]],
00221 edges[fEdges[0]].commonVertex(edges[fEdges[1]]),
00222 edgeMid[fEdges[1]],
00223 surf[faceI].region()
00224 )
00225 );
00226
00227 newFaces.append
00228 (
00229 labelledTri
00230 (
00231 edgeMid[fEdges[1]],
00232 edges[fEdges[1]].commonVertex(edges[fEdges[2]]),
00233 edgeMid[fEdges[2]],
00234 surf[faceI].region()
00235 )
00236 );
00237
00238 newFaces.append
00239 (
00240 labelledTri
00241 (
00242 edgeMid[fEdges[2]],
00243 edges[fEdges[2]].commonVertex(edges[fEdges[0]]),
00244 edgeMid[fEdges[0]],
00245 surf[faceI].region()
00246 )
00247 );
00248
00249
00250 newFaces.append
00251 (
00252 labelledTri
00253 (
00254 edgeMid[fEdges[0]],
00255 edgeMid[fEdges[1]],
00256 edgeMid[fEdges[2]],
00257 surf[faceI].region()
00258 )
00259 );
00260
00261
00262
00263 forAll(fEdges, i)
00264 {
00265 const label edgeI = fEdges[i];
00266
00267 label otherFaceI = otherFace(surf, faceI, edgeI);
00268
00269 if ((otherFaceI != -1) && (refineStatus[otherFaceI] == GREEN))
00270 {
00271 greenRefine
00272 (
00273 surf,
00274 otherFaceI,
00275 edgeI,
00276 edgeMid[edgeI],
00277 newFaces
00278 );
00279 }
00280 }
00281 }
00282 }
00283
00284
00285 forAll(refineStatus, faceI)
00286 {
00287 if (refineStatus[faceI] == NONE)
00288 {
00289 newFaces.append(surf.localFaces()[faceI]);
00290 }
00291 }
00292
00293 newFaces.shrink();
00294 newPoints.shrink();
00295
00296
00297
00298 pointField allPoints;
00299 allPoints.transfer(newPoints);
00300 newPoints.clear();
00301
00302 return triSurface(newFaces, surf.patches(), allPoints, true);
00303 }
00304
00305
00306
00307
00308 Foam::scalar Foam::triSurfaceTools::faceCosAngle
00309 (
00310 const point& pStart,
00311 const point& pEnd,
00312 const point& pLeft,
00313 const point& pRight
00314 )
00315 {
00316 const vector common(pEnd - pStart);
00317 const vector base0(pLeft - pStart);
00318 const vector base1(pRight - pStart);
00319
00320 vector n0(common ^ base0);
00321 n0 /= Foam::mag(n0);
00322
00323 vector n1(base1 ^ common);
00324 n1 /= Foam::mag(n1);
00325
00326 return n0 & n1;
00327 }
00328
00329
00330
00331
00332
00333
00334 void Foam::triSurfaceTools::protectNeighbours
00335 (
00336 const triSurface& surf,
00337 const label vertI,
00338 labelList& faceStatus
00339 )
00340 {
00341
00342
00343
00344
00345
00346
00347
00348
00349
00350
00351
00352 const labelList& myEdges = surf.pointEdges()[vertI];
00353 forAll(myEdges, i)
00354 {
00355 const labelList& myFaces = surf.edgeFaces()[myEdges[i]];
00356
00357 forAll(myFaces, myFaceI)
00358 {
00359 label faceI = myFaces[myFaceI];
00360
00361 if ((faceStatus[faceI] == ANYEDGE) || (faceStatus[faceI] >= 0))
00362 {
00363 faceStatus[faceI] = NOEDGE;
00364 }
00365 }
00366 }
00367 }
00368
00369
00370
00371
00372
00373
00374
00375 Foam::labelHashSet Foam::triSurfaceTools::getCollapsedFaces
00376 (
00377 const triSurface& surf,
00378 label edgeI
00379 )
00380 {
00381 const edge& e = surf.edges()[edgeI];
00382 label v1 = e.start();
00383 label v2 = e.end();
00384
00385
00386 const labelList& myFaces = surf.edgeFaces()[edgeI];
00387
00388 labelHashSet facesToBeCollapsed(2*myFaces.size());
00389
00390 forAll(myFaces, myFaceI)
00391 {
00392 facesToBeCollapsed.insert(myFaces[myFaceI]);
00393 }
00394
00395
00396
00397
00398
00399 const labelList& v1Faces = surf.pointFaces()[v1];
00400
00401 forAll(v1Faces, v1FaceI)
00402 {
00403 label face1I = v1Faces[v1FaceI];
00404
00405 label otherEdgeI = oppositeEdge(surf, face1I, v1);
00406
00407
00408 label face2I = otherFace(surf, face1I, otherEdgeI);
00409
00410 if (face2I != -1)
00411 {
00412
00413 if (oppositeVertex(surf, face2I, otherEdgeI) == v2)
00414 {
00415
00416
00417 facesToBeCollapsed.insert(face1I);
00418 facesToBeCollapsed.insert(face2I);
00419 }
00420 }
00421 }
00422
00423 return facesToBeCollapsed;
00424 }
00425
00426
00427
00428 Foam::label Foam::triSurfaceTools::vertexUsesFace
00429 (
00430 const triSurface& surf,
00431 const labelHashSet& faceUsed,
00432 const label vertI
00433 )
00434 {
00435 const labelList& myFaces = surf.pointFaces()[vertI];
00436
00437 forAll(myFaces, myFaceI)
00438 {
00439 label face1I = myFaces[myFaceI];
00440
00441 if (faceUsed.found(face1I))
00442 {
00443 return face1I;
00444 }
00445 }
00446 return -1;
00447 }
00448
00449
00450
00451 void Foam::triSurfaceTools::getMergedEdges
00452 (
00453 const triSurface& surf,
00454 const label edgeI,
00455 const labelHashSet& collapsedFaces,
00456 HashTable<label, label, Hash<label> >& edgeToEdge,
00457 HashTable<label, label, Hash<label> >& edgeToFace
00458 )
00459 {
00460 const edge& e = surf.edges()[edgeI];
00461 label v1 = e.start();
00462 label v2 = e.end();
00463
00464 const labelList& v1Faces = surf.pointFaces()[v1];
00465 const labelList& v2Faces = surf.pointFaces()[v2];
00466
00467
00468 labelHashSet v2FacesHash(v2Faces.size());
00469
00470 forAll(v2Faces, v2FaceI)
00471 {
00472 if (!collapsedFaces.found(v2Faces[v2FaceI]))
00473 {
00474 v2FacesHash.insert(v2Faces[v2FaceI]);
00475 }
00476 }
00477
00478
00479 forAll(v1Faces, v1FaceI)
00480 {
00481 label face1I = v1Faces[v1FaceI];
00482
00483 if (collapsedFaces.found(face1I))
00484 {
00485 continue;
00486 }
00487
00488
00489
00490
00491
00492 label vert1I = -1;
00493 label vert2I = -1;
00494 otherVertices
00495 (
00496 surf,
00497 face1I,
00498 v1,
00499 vert1I,
00500 vert2I
00501 );
00502
00503
00504
00505
00506
00507 label commonVert = vert1I;
00508 label face2I = vertexUsesFace(surf, v2FacesHash, commonVert);
00509 if (face2I == -1)
00510 {
00511 commonVert = vert2I;
00512 face2I = vertexUsesFace(surf, v2FacesHash, commonVert);
00513 }
00514
00515 if (face2I != -1)
00516 {
00517
00518 label edge1I = getEdge(surf, v1, commonVert);
00519 label edge2I = getEdge(surf, v2, commonVert);
00520
00521 edgeToEdge.insert(edge1I, edge2I);
00522 edgeToEdge.insert(edge2I, edge1I);
00523
00524 edgeToFace.insert(edge1I, face2I);
00525 edgeToFace.insert(edge2I, face1I);
00526 }
00527 }
00528 }
00529
00530
00531
00532
00533 Foam::scalar Foam::triSurfaceTools::edgeCosAngle
00534 (
00535 const triSurface& surf,
00536 const label v1,
00537 const point& pt,
00538 const labelHashSet& collapsedFaces,
00539 const HashTable<label, label, Hash<label> >& edgeToEdge,
00540 const HashTable<label, label, Hash<label> >& edgeToFace,
00541 const label faceI,
00542 const label edgeI
00543 )
00544 {
00545 const pointField& localPoints = surf.localPoints();
00546
00547 label A = surf.edges()[edgeI].start();
00548 label B = surf.edges()[edgeI].end();
00549 label C = oppositeVertex(surf, faceI, edgeI);
00550
00551 label D = -1;
00552
00553 label face2I = -1;
00554
00555 if (edgeToEdge.found(edgeI))
00556 {
00557
00558 label edge2I = edgeToEdge[edgeI];
00559 face2I = edgeToFace[edgeI];
00560
00561 D = oppositeVertex(surf, face2I, edge2I);
00562 }
00563 else
00564 {
00565
00566 face2I = otherFace(surf, faceI, edgeI);
00567
00568 if ((face2I != -1) && !collapsedFaces.found(face2I))
00569 {
00570 D = oppositeVertex(surf, face2I, edgeI);
00571 }
00572 }
00573
00574 scalar cosAngle = 1;
00575
00576 if (D != -1)
00577 {
00578 if (A == v1)
00579 {
00580 cosAngle = faceCosAngle
00581 (
00582 pt,
00583 localPoints[B],
00584 localPoints[C],
00585 localPoints[D]
00586 );
00587 }
00588 else if (B == v1)
00589 {
00590 cosAngle = faceCosAngle
00591 (
00592 localPoints[A],
00593 pt,
00594 localPoints[C],
00595 localPoints[D]
00596 );
00597 }
00598 else if (C == v1)
00599 {
00600 cosAngle = faceCosAngle
00601 (
00602 localPoints[A],
00603 localPoints[B],
00604 pt,
00605 localPoints[D]
00606 );
00607 }
00608 else if (D == v1)
00609 {
00610 cosAngle = faceCosAngle
00611 (
00612 localPoints[A],
00613 localPoints[B],
00614 localPoints[C],
00615 pt
00616 );
00617 }
00618 else
00619 {
00620 FatalErrorIn("edgeCosAngle")
00621 << "face " << faceI << " does not use vertex "
00622 << v1 << " of collapsed edge" << abort(FatalError);
00623 }
00624 }
00625 return cosAngle;
00626 }
00627
00628
00629
00630
00631 Foam::scalar Foam::triSurfaceTools::collapseMinCosAngle
00632 (
00633 const triSurface& surf,
00634 const label v1,
00635 const point& pt,
00636 const labelHashSet& collapsedFaces,
00637 const HashTable<label, label, Hash<label> >& edgeToEdge,
00638 const HashTable<label, label, Hash<label> >& edgeToFace
00639 )
00640 {
00641 const labelList& v1Faces = surf.pointFaces()[v1];
00642
00643 scalar minCos = 1;
00644
00645 forAll(v1Faces, v1FaceI)
00646 {
00647 label faceI = v1Faces[v1FaceI];
00648
00649 if (collapsedFaces.found(faceI))
00650 {
00651 continue;
00652 }
00653
00654 const labelList& myEdges = surf.faceEdges()[faceI];
00655
00656 forAll(myEdges, myEdgeI)
00657 {
00658 label edgeI = myEdges[myEdgeI];
00659
00660 minCos =
00661 min
00662 (
00663 minCos,
00664 edgeCosAngle
00665 (
00666 surf,
00667 v1,
00668 pt,
00669 collapsedFaces,
00670 edgeToEdge,
00671 edgeToFace,
00672 faceI,
00673 edgeI
00674 )
00675 );
00676 }
00677 }
00678
00679 return minCos;
00680 }
00681
00682
00683
00684
00685 bool Foam::triSurfaceTools::collapseCreatesFold
00686 (
00687 const triSurface& surf,
00688 const label v1,
00689 const point& pt,
00690 const labelHashSet& collapsedFaces,
00691 const HashTable<label, label, Hash<label> >& edgeToEdge,
00692 const HashTable<label, label, Hash<label> >& edgeToFace,
00693 const scalar minCos
00694 )
00695 {
00696 const labelList& v1Faces = surf.pointFaces()[v1];
00697
00698 forAll(v1Faces, v1FaceI)
00699 {
00700 label faceI = v1Faces[v1FaceI];
00701
00702 if (collapsedFaces.found(faceI))
00703 {
00704 continue;
00705 }
00706
00707 const labelList& myEdges = surf.faceEdges()[faceI];
00708
00709 forAll(myEdges, myEdgeI)
00710 {
00711 label edgeI = myEdges[myEdgeI];
00712
00713 if
00714 (
00715 edgeCosAngle
00716 (
00717 surf,
00718 v1,
00719 pt,
00720 collapsedFaces,
00721 edgeToEdge,
00722 edgeToFace,
00723 faceI,
00724 edgeI
00725 )
00726 < minCos
00727 )
00728 {
00729 return true;
00730 }
00731 }
00732 }
00733
00734 return false;
00735 }
00736
00737
00740
00741
00742
00743
00744
00745
00746
00747
00748
00749
00750
00751
00752
00753
00754
00755
00756
00757
00758
00759
00760
00761
00762
00763
00764
00765
00766
00767
00768
00769
00770
00771
00772
00773
00774
00775
00776
00777
00778
00779
00780
00781
00782
00783
00784
00785
00786
00787
00788
00789
00790
00791
00792
00793
00794
00795
00796
00797
00798
00799
00800
00801
00802
00803
00804
00805
00806
00807
00808
00809
00810
00811
00812
00813
00814
00815
00816
00817
00818
00819
00820
00821
00822
00823
00824
00825
00826
00827
00828
00829
00830
00831
00832
00833
00834
00835
00836
00837
00838
00839
00840
00841
00842
00843
00844 Foam::surfaceLocation Foam::triSurfaceTools::cutEdge
00845 (
00846 const triSurface& s,
00847 const label triI,
00848 const label excludeEdgeI,
00849 const label excludePointI,
00850
00851 const point& triPoint,
00852 const plane& cutPlane,
00853 const point& toPoint
00854 )
00855 {
00856 const pointField& points = s.points();
00857 const labelledTri& f = s[triI];
00858 const labelList& fEdges = s.faceEdges()[triI];
00859
00860
00861 FixedList<scalar, 3> d;
00862
00863 scalar norm = 0;
00864 forAll(d, fp)
00865 {
00866 d[fp] = (points[f[fp]]-cutPlane.refPoint()) & cutPlane.normal();
00867 norm += mag(d[fp]);
00868 }
00869
00870
00871 forAll(d, i)
00872 {
00873 d[i] /= norm;
00874
00875 if (mag(d[i]) < 1E-6)
00876 {
00877 d[i] = 0.0;
00878 }
00879 }
00880
00881
00882 surfaceLocation cut;
00883
00884 if (excludePointI != -1)
00885 {
00886
00887
00888 label fp0 = findIndex(s.localFaces()[triI], excludePointI);
00889
00890 if (fp0 == -1)
00891 {
00892 FatalErrorIn("cutEdge(..)") << "excludePointI:" << excludePointI
00893 << " localF:" << s.localFaces()[triI] << abort(FatalError);
00894 }
00895
00896 label fp1 = f.fcIndex(fp0);
00897 label fp2 = f.fcIndex(fp1);
00898
00899
00900 if (d[fp1] == 0.0)
00901 {
00902 cut.setHit();
00903 cut.setPoint(points[f[fp1]]);
00904 cut.elementType() = triPointRef::POINT;
00905 cut.setIndex(s.localFaces()[triI][fp1]);
00906 }
00907 else if (d[fp2] == 0.0)
00908 {
00909 cut.setHit();
00910 cut.setPoint(points[f[fp2]]);
00911 cut.elementType() = triPointRef::POINT;
00912 cut.setIndex(s.localFaces()[triI][fp2]);
00913 }
00914 else if
00915 (
00916 (d[fp1] < 0 && d[fp2] < 0)
00917 || (d[fp1] > 0 && d[fp2] > 0)
00918 )
00919 {
00920
00921
00922 }
00923 else
00924 {
00925 cut.setHit();
00926 cut.setPoint
00927 (
00928 (d[fp2]*points[f[fp1]] - d[fp1]*points[f[fp2]])
00929 / (d[fp2] - d[fp1])
00930 );
00931 cut.elementType() = triPointRef::EDGE;
00932 cut.setIndex(fEdges[fp1]);
00933 }
00934 }
00935 else
00936 {
00937
00938 FixedList<surfaceLocation, 2> inters;
00939 label interI = 0;
00940
00941 forAll(f, fp0)
00942 {
00943 label fp1 = f.fcIndex(fp0);
00944
00945 if (d[fp0] == 0)
00946 {
00947 if (interI >= 2)
00948 {
00949 FatalErrorIn("cutEdge(..)")
00950 << "problem : triangle has three intersections." << nl
00951 << "triangle:" << f.tri(points)
00952 << " d:" << d << abort(FatalError);
00953 }
00954 inters[interI].setHit();
00955 inters[interI].setPoint(points[f[fp0]]);
00956 inters[interI].elementType() = triPointRef::POINT;
00957 inters[interI].setIndex(s.localFaces()[triI][fp0]);
00958 interI++;
00959 }
00960 else if
00961 (
00962 (d[fp0] < 0 && d[fp1] > 0)
00963 || (d[fp0] > 0 && d[fp1] < 0)
00964 )
00965 {
00966 if (interI >= 2)
00967 {
00968 FatalErrorIn("cutEdge(..)")
00969 << "problem : triangle has three intersections." << nl
00970 << "triangle:" << f.tri(points)
00971 << " d:" << d << abort(FatalError);
00972 }
00973 inters[interI].setHit();
00974 inters[interI].setPoint
00975 (
00976 (d[fp0]*points[f[fp1]] - d[fp1]*points[f[fp0]])
00977 / (d[fp0] - d[fp1])
00978 );
00979 inters[interI].elementType() = triPointRef::EDGE;
00980 inters[interI].setIndex(fEdges[fp0]);
00981 interI++;
00982 }
00983 }
00984
00985
00986 if (interI == 0)
00987 {
00988
00989 }
00990 else if (interI == 1)
00991 {
00992
00993 cut = inters[0];
00994 }
00995 else if (interI == 2)
00996 {
00997
00998 if
00999 (
01000 inters[0].elementType() == triPointRef::EDGE
01001 && inters[0].index() == excludeEdgeI
01002 )
01003 {
01004 cut = inters[1];
01005 }
01006 else if
01007 (
01008 inters[1].elementType() == triPointRef::EDGE
01009 && inters[1].index() == excludeEdgeI
01010 )
01011 {
01012 cut = inters[0];
01013 }
01014 else
01015 {
01016
01017 if
01018 (
01019 magSqr(inters[0].rawPoint() - toPoint)
01020 < magSqr(inters[1].rawPoint() - toPoint)
01021 )
01022 {
01023 cut = inters[0];
01024 }
01025 else
01026 {
01027 cut = inters[1];
01028 }
01029 }
01030 }
01031 }
01032 return cut;
01033 }
01034
01035
01036
01037 void Foam::triSurfaceTools::snapToEnd
01038 (
01039 const triSurface& s,
01040 const surfaceLocation& end,
01041 surfaceLocation& current
01042 )
01043 {
01044 if (end.elementType() == triPointRef::NONE)
01045 {
01046 if (current.elementType() == triPointRef::NONE)
01047 {
01048
01049 if (current.index() == end.index())
01050 {
01051
01052
01053
01054
01055
01056 current = end;
01057 current.setHit();
01058 }
01059 }
01060
01061 }
01062 else if (end.elementType() == triPointRef::EDGE)
01063 {
01064 if (current.elementType() == triPointRef::NONE)
01065 {
01066
01067 const labelList& fEdges = s.faceEdges()[current.index()];
01068
01069 if (findIndex(fEdges, end.index()) != -1)
01070 {
01071
01072
01073
01074
01075
01076 current = end;
01077 current.setHit();
01078 }
01079 }
01080 else if (current.elementType() == triPointRef::EDGE)
01081 {
01082
01083 if (current.index() == end.index())
01084 {
01085
01086
01087
01088
01089
01090 current = end;
01091 current.setHit();
01092 }
01093 }
01094 else
01095 {
01096
01097 const edge& e = s.edges()[end.index()];
01098
01099 if (current.index() == e[0] || current.index() == e[1])
01100 {
01101
01102
01103
01104
01105
01106 current = end;
01107 current.setHit();
01108 }
01109 }
01110 }
01111 else
01112 {
01113 if (current.elementType() == triPointRef::NONE)
01114 {
01115
01116 const labelledTri& f = s.localFaces()[current.index()];
01117
01118 if (findIndex(f, end.index()) != -1)
01119 {
01120
01121
01122
01123
01124
01125 current = end;
01126 current.setHit();
01127 }
01128 }
01129 else if (current.elementType() == triPointRef::EDGE)
01130 {
01131
01132 const edge& e = s.edges()[current.index()];
01133
01134 if (end.index() == e[0] || end.index() == e[1])
01135 {
01136
01137
01138
01139
01140
01141 current = end;
01142 current.setHit();
01143 }
01144 }
01145 else
01146 {
01147
01148 if (current.index() == end.index())
01149 {
01150
01151
01152
01153
01154
01155 current = end;
01156 current.setHit();
01157 }
01158 }
01159 }
01160 }
01161
01162
01163
01164
01165
01166
01167 Foam::surfaceLocation Foam::triSurfaceTools::visitFaces
01168 (
01169 const triSurface& s,
01170 const labelList& eFaces,
01171 const surfaceLocation& start,
01172 const label excludeEdgeI,
01173 const label excludePointI,
01174 const surfaceLocation& end,
01175 const plane& cutPlane
01176 )
01177 {
01178 surfaceLocation nearest;
01179
01180 scalar minDistSqr = Foam::sqr(GREAT);
01181
01182 forAll(eFaces, i)
01183 {
01184 label triI = eFaces[i];
01185
01186
01187 if (triI != start.triangle())
01188 {
01189 if (end.elementType() == triPointRef::NONE && end.index() == triI)
01190 {
01191
01192 nearest = end;
01193 nearest.setHit();
01194 nearest.triangle() = triI;
01195 break;
01196 }
01197 else
01198 {
01199
01200
01201 surfaceLocation cutInfo = cutEdge
01202 (
01203 s,
01204 triI,
01205 excludeEdgeI,
01206 excludePointI,
01207 start.rawPoint(),
01208 cutPlane,
01209 end.rawPoint()
01210 );
01211
01212
01213 if (excludeEdgeI != -1 && !cutInfo.hit())
01214 {
01215 FatalErrorIn("triSurfaceTools::visitFaces(..)")
01216 << "Triangle:" << triI
01217 << " excludeEdge:" << excludeEdgeI
01218 << " point:" << start.rawPoint()
01219 << " plane:" << cutPlane
01220 << " . No intersection!" << abort(FatalError);
01221 }
01222
01223 if (cutInfo.hit())
01224 {
01225 scalar distSqr = magSqr(cutInfo.rawPoint()-end.rawPoint());
01226
01227 if (distSqr < minDistSqr)
01228 {
01229 minDistSqr = distSqr;
01230 nearest = cutInfo;
01231 nearest.triangle() = triI;
01232 nearest.setMiss();
01233 }
01234 }
01235 }
01236 }
01237 }
01238
01239 if (nearest.triangle() == -1)
01240 {
01241
01242
01243
01244 }
01245
01246 return nearest;
01247 }
01248
01249
01250
01251
01252
01253
01254 void Foam::triSurfaceTools::writeOBJ
01255 (
01256 const fileName& fName,
01257 const pointField& pts
01258 )
01259 {
01260 OFstream outFile(fName);
01261
01262 forAll(pts, pointI)
01263 {
01264 const point& pt = pts[pointI];
01265
01266 outFile<< "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << endl;
01267 }
01268 Pout<< "Written " << pts.size() << " vertices to file " << fName << endl;
01269 }
01270
01271
01272
01273 void Foam::triSurfaceTools::writeOBJ
01274 (
01275 const triSurface& surf,
01276 const fileName& fName,
01277 const boolList& markedVerts
01278 )
01279 {
01280 OFstream outFile(fName);
01281
01282 label nVerts = 0;
01283 forAll(markedVerts, vertI)
01284 {
01285 if (markedVerts[vertI])
01286 {
01287 const point& pt = surf.localPoints()[vertI];
01288
01289 outFile<< "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << endl;
01290
01291 nVerts++;
01292 }
01293 }
01294 Pout<< "Written " << nVerts << " vertices to file " << fName << endl;
01295 }
01296
01297
01298
01299
01300
01301
01302
01303 void Foam::triSurfaceTools::getVertexTriangles
01304 (
01305 const triSurface& surf,
01306 const label edgeI,
01307 labelList& edgeTris
01308 )
01309 {
01310 const edge& e = surf.edges()[edgeI];
01311 const labelList& myFaces = surf.edgeFaces()[edgeI];
01312
01313 label face1I = myFaces[0];
01314 label face2I = -1;
01315 if (myFaces.size() == 2)
01316 {
01317 face2I = myFaces[1];
01318 }
01319
01320 const labelList& startFaces = surf.pointFaces()[e.start()];
01321 const labelList& endFaces = surf.pointFaces()[e.end()];
01322
01323
01324
01325 edgeTris.setSize(startFaces.size() + endFaces.size() - myFaces.size());
01326
01327 label nTris = 0;
01328 forAll(startFaces, startFaceI)
01329 {
01330 edgeTris[nTris++] = startFaces[startFaceI];
01331 }
01332
01333 forAll(endFaces, endFaceI)
01334 {
01335 label faceI = endFaces[endFaceI];
01336
01337 if ((faceI != face1I) && (faceI != face2I))
01338 {
01339 edgeTris[nTris++] = faceI;
01340 }
01341 }
01342 }
01343
01344
01345
01346 Foam::labelList Foam::triSurfaceTools::getVertexVertices
01347 (
01348 const triSurface& surf,
01349 const edge& e
01350 )
01351 {
01352 const edgeList& edges = surf.edges();
01353
01354 label v1 = e.start();
01355 label v2 = e.end();
01356
01357
01358 labelHashSet vertexNeighbours;
01359
01360 const labelList& v1Edges = surf.pointEdges()[v1];
01361
01362 forAll(v1Edges, v1EdgeI)
01363 {
01364 const edge& e = edges[v1Edges[v1EdgeI]];
01365 vertexNeighbours.insert(e.otherVertex(v1));
01366 }
01367
01368 const labelList& v2Edges = surf.pointEdges()[v2];
01369
01370 forAll(v2Edges, v2EdgeI)
01371 {
01372 const edge& e = edges[v2Edges[v2EdgeI]];
01373
01374 label vertI = e.otherVertex(v2);
01375
01376 vertexNeighbours.insert(vertI);
01377 }
01378 return vertexNeighbours.toc();
01379 }
01380
01381
01383
01384
01385
01386
01387
01388
01389
01390
01391
01392
01393
01394
01395
01396
01397
01398
01399
01400
01401
01402
01403
01404
01405
01406
01407
01408
01409
01410
01411
01412
01413
01414
01415
01416
01417
01418
01419
01420
01421
01422
01423
01424
01425
01426
01427
01428
01429
01430
01431 Foam::label Foam::triSurfaceTools::otherFace
01432 (
01433 const triSurface& surf,
01434 const label faceI,
01435 const label edgeI
01436 )
01437 {
01438 const labelList& myFaces = surf.edgeFaces()[edgeI];
01439
01440 if (myFaces.size() != 2)
01441 {
01442 return -1;
01443 }
01444 else
01445 {
01446 if (faceI == myFaces[0])
01447 {
01448 return myFaces[1];
01449 }
01450 else
01451 {
01452 return myFaces[0];
01453 }
01454 }
01455 }
01456
01457
01458
01459 void Foam::triSurfaceTools::otherEdges
01460 (
01461 const triSurface& surf,
01462 const label faceI,
01463 const label edgeI,
01464 label& e1,
01465 label& e2
01466 )
01467 {
01468 const labelList& eFaces = surf.faceEdges()[faceI];
01469
01470 label i0 = findIndex(eFaces, edgeI);
01471
01472 if (i0 == -1)
01473 {
01474 FatalErrorIn
01475 (
01476 "otherEdges"
01477 "(const triSurface&, const label, const label,"
01478 " label&, label&)"
01479 ) << "Edge " << surf.edges()[edgeI] << " not in face "
01480 << surf.localFaces()[faceI] << abort(FatalError);
01481 }
01482
01483 label i1 = eFaces.fcIndex(i0);
01484 label i2 = eFaces.fcIndex(i1);
01485
01486 e1 = eFaces[i1];
01487 e2 = eFaces[i2];
01488 }
01489
01490
01491
01492 void Foam::triSurfaceTools::otherVertices
01493 (
01494 const triSurface& surf,
01495 const label faceI,
01496 const label vertI,
01497 label& vert1I,
01498 label& vert2I
01499 )
01500 {
01501 const labelledTri& f = surf.localFaces()[faceI];
01502
01503 if (vertI == f[0])
01504 {
01505 vert1I = f[1];
01506 vert2I = f[2];
01507 }
01508 else if (vertI == f[1])
01509 {
01510 vert1I = f[2];
01511 vert2I = f[0];
01512 }
01513 else if (vertI == f[2])
01514 {
01515 vert1I = f[0];
01516 vert2I = f[1];
01517 }
01518 else
01519 {
01520 FatalErrorIn
01521 (
01522 "otherVertices"
01523 "(const triSurface&, const label, const label,"
01524 " label&, label&)"
01525 ) << "Vertex " << vertI << " not in face " << f << abort(FatalError);
01526 }
01527 }
01528
01529
01530
01531 Foam::label Foam::triSurfaceTools::oppositeEdge
01532 (
01533 const triSurface& surf,
01534 const label faceI,
01535 const label vertI
01536 )
01537 {
01538 const labelList& myEdges = surf.faceEdges()[faceI];
01539
01540 forAll(myEdges, myEdgeI)
01541 {
01542 label edgeI = myEdges[myEdgeI];
01543
01544 const edge& e = surf.edges()[edgeI];
01545
01546 if ((e.start() != vertI) && (e.end() != vertI))
01547 {
01548 return edgeI;
01549 }
01550 }
01551
01552 FatalErrorIn
01553 (
01554 "oppositeEdge"
01555 "(const triSurface&, const label, const label)"
01556 ) << "Cannot find vertex " << vertI << " in edges of face " << faceI
01557 << abort(FatalError);
01558
01559 return -1;
01560 }
01561
01562
01563
01564 Foam::label Foam::triSurfaceTools::oppositeVertex
01565 (
01566 const triSurface& surf,
01567 const label faceI,
01568 const label edgeI
01569 )
01570 {
01571 const labelledTri& f = surf.localFaces()[faceI];
01572
01573 const edge& e = surf.edges()[edgeI];
01574
01575 forAll(f, fp)
01576 {
01577 label vertI = f[fp];
01578
01579 if (vertI != e.start() && vertI != e.end())
01580 {
01581 return vertI;
01582 }
01583 }
01584
01585 FatalErrorIn("triSurfaceTools::oppositeVertex")
01586 << "Cannot find vertex opposite edge " << edgeI << " vertices " << e
01587 << " in face " << faceI << " vertices " << f << abort(FatalError);
01588
01589 return -1;
01590 }
01591
01592
01593
01594 Foam::label Foam::triSurfaceTools::getEdge
01595 (
01596 const triSurface& surf,
01597 const label v1,
01598 const label v2
01599 )
01600 {
01601 const labelList& v1Edges = surf.pointEdges()[v1];
01602
01603 forAll(v1Edges, v1EdgeI)
01604 {
01605 label edgeI = v1Edges[v1EdgeI];
01606
01607 const edge& e = surf.edges()[edgeI];
01608
01609 if ((e.start() == v2) || (e.end() == v2))
01610 {
01611 return edgeI;
01612 }
01613 }
01614 return -1;
01615 }
01616
01617
01618
01619 Foam::label Foam::triSurfaceTools::getTriangle
01620 (
01621 const triSurface& surf,
01622 const label e0I,
01623 const label e1I,
01624 const label e2I
01625 )
01626 {
01627 if ((e0I == e1I) || (e0I == e2I) || (e1I == e2I))
01628 {
01629 FatalErrorIn
01630 (
01631 "getTriangle"
01632 "(const triSurface&, const label, const label,"
01633 " const label)"
01634 ) << "Duplicate edge labels : e0:" << e0I << " e1:" << e1I
01635 << " e2:" << e2I
01636 << abort(FatalError);
01637 }
01638
01639 const labelList& eFaces = surf.edgeFaces()[e0I];
01640
01641 forAll(eFaces, eFaceI)
01642 {
01643 label faceI = eFaces[eFaceI];
01644
01645 const labelList& myEdges = surf.faceEdges()[faceI];
01646
01647 if
01648 (
01649 (myEdges[0] == e1I)
01650 || (myEdges[1] == e1I)
01651 || (myEdges[2] == e1I)
01652 )
01653 {
01654 if
01655 (
01656 (myEdges[0] == e2I)
01657 || (myEdges[1] == e2I)
01658 || (myEdges[2] == e2I)
01659 )
01660 {
01661 return faceI;
01662 }
01663 }
01664 }
01665 return -1;
01666 }
01667
01668
01669
01670 Foam::triSurface Foam::triSurfaceTools::collapseEdges
01671 (
01672 const triSurface& surf,
01673 const labelList& collapsableEdges
01674 )
01675 {
01676 pointField edgeMids(surf.nEdges());
01677
01678 forAll(edgeMids, edgeI)
01679 {
01680 const edge& e = surf.edges()[edgeI];
01681
01682 edgeMids[edgeI] =
01683 0.5
01684 * (
01685 surf.localPoints()[e.start()]
01686 + surf.localPoints()[e.end()]
01687 );
01688 }
01689
01690
01691 labelList faceStatus(surf.size(), ANYEDGE);
01692
01694
01695
01696
01697
01698
01699
01700
01701
01702
01703
01704
01705
01706
01707
01708
01709
01710
01711
01712
01713
01714
01715
01716
01717
01718 return collapseEdges(surf, collapsableEdges, edgeMids, faceStatus);
01719 }
01720
01721
01722
01723 Foam::triSurface Foam::triSurfaceTools::collapseEdges
01724 (
01725 const triSurface& surf,
01726 const labelList& collapseEdgeLabels,
01727 const pointField& edgeMids,
01728 labelList& faceStatus
01729 )
01730 {
01731 const labelListList& edgeFaces = surf.edgeFaces();
01732 const pointField& localPoints = surf.localPoints();
01733 const edgeList& edges = surf.edges();
01734
01735
01736 pointField newPoints(localPoints);
01737
01738
01739 labelList pointMap(localPoints.size());
01740 forAll(localPoints, pointI)
01741 {
01742 pointMap[pointI] = pointI;
01743 }
01744
01745
01746
01747
01748 forAll(collapseEdgeLabels, collapseEdgeI)
01749 {
01750 const label edgeI = collapseEdgeLabels[collapseEdgeI];
01751
01752 if ((edgeI < 0) || (edgeI >= surf.nEdges()))
01753 {
01754 FatalErrorIn("collapseEdges")
01755 << "Edge label outside valid range." << endl
01756 << "edge label:" << edgeI << endl
01757 << "total number of edges:" << surf.nEdges() << endl
01758 << abort(FatalError);
01759 }
01760
01761 const labelList& neighbours = edgeFaces[edgeI];
01762
01763 if (neighbours.size() == 2)
01764 {
01765 const label stat0 = faceStatus[neighbours[0]];
01766 const label stat1 = faceStatus[neighbours[1]];
01767
01768
01769 if
01770 (
01771 ((stat0 == ANYEDGE) || (stat0 == edgeI))
01772 && ((stat1 == ANYEDGE) || (stat1 == edgeI))
01773 )
01774 {
01775 const edge& e = edges[edgeI];
01776
01777
01778 if
01779 (
01780 (pointMap[e.start()] != e.start())
01781 || (pointMap[e.end()] != e.end())
01782 )
01783 {
01784 FatalErrorIn("collapseEdges")
01785 << "points already mapped. Double collapse." << endl
01786 << "edgeI:" << edgeI
01787 << " start:" << e.start()
01788 << " end:" << e.end()
01789 << " pointMap[start]:" << pointMap[e.start()]
01790 << " pointMap[end]:" << pointMap[e.end()]
01791 << abort(FatalError);
01792 }
01793
01794 const label minVert = min(e.start(), e.end());
01795 pointMap[e.start()] = minVert;
01796 pointMap[e.end()] = minVert;
01797
01798
01799 newPoints[minVert] = edgeMids[edgeI];
01800
01801
01802 protectNeighbours(surf, e.start(), faceStatus);
01803 protectNeighbours(surf, e.end(), faceStatus);
01804 protectNeighbours
01805 (
01806 surf,
01807 oppositeVertex(surf, neighbours[0], edgeI),
01808 faceStatus
01809 );
01810 protectNeighbours
01811 (
01812 surf,
01813 oppositeVertex(surf, neighbours[1], edgeI),
01814 faceStatus
01815 );
01816
01817
01818 labelList collapseFaces =
01819 getCollapsedFaces
01820 (
01821 surf,
01822 edgeI
01823 ).toc();
01824
01825 forAll(collapseFaces, collapseI)
01826 {
01827 faceStatus[collapseFaces[collapseI]] = COLLAPSED;
01828 }
01829 }
01830 }
01831 }
01832
01833
01834
01835 List<labelledTri> newTris(surf.size());
01836 label newTriI = 0;
01837
01838 const List<labelledTri>& localFaces = surf.localFaces();
01839
01840
01841
01842 forAll(localFaces, faceI)
01843 {
01844 const labelledTri& f = localFaces[faceI];
01845
01846 const label a = pointMap[f[0]];
01847 const label b = pointMap[f[1]];
01848 const label c = pointMap[f[2]];
01849
01850 if
01851 (
01852 (a != b) && (a != c) && (b != c)
01853 && (faceStatus[faceI] != COLLAPSED)
01854 )
01855 {
01856
01857 newTris[newTriI++] = labelledTri(a, b, c, f.region());
01858 }
01859 else
01860 {
01861
01862
01863 }
01864 }
01865 newTris.setSize(newTriI);
01866
01867
01868
01869
01870
01871 triSurface tempSurf(newTris, surf.patches(), newPoints);
01872
01873 return
01874 triSurface
01875 (
01876 tempSurf.localFaces(),
01877 tempSurf.patches(),
01878 tempSurf.localPoints()
01879 );
01880 }
01881
01882
01883
01884
01885 Foam::triSurface Foam::triSurfaceTools::redGreenRefine
01886 (
01887 const triSurface& surf,
01888 const labelList& refineFaces
01889 )
01890 {
01891 List<refineType> refineStatus(surf.size(), NONE);
01892
01893
01894 forAll(refineFaces, refineFaceI)
01895 {
01896 calcRefineStatus(surf, refineFaces[refineFaceI], refineStatus);
01897 }
01898
01899
01900 return doRefine(surf, refineStatus);
01901 }
01902
01903
01904 Foam::triSurface Foam::triSurfaceTools::greenRefine
01905 (
01906 const triSurface& surf,
01907 const labelList& refineEdges
01908 )
01909 {
01910
01911 List<refineType> refineStatus(surf.size(), NONE);
01912
01913
01914 DynamicList<labelledTri> newFaces(0);
01915
01916 pointField newPoints(surf.localPoints());
01917 newPoints.setSize(surf.nPoints() + surf.nEdges());
01918 label newPointI = surf.nPoints();
01919
01920
01921
01922 forAll(refineEdges, refineEdgeI)
01923 {
01924 label edgeI = refineEdges[refineEdgeI];
01925
01926 const labelList& myFaces = surf.edgeFaces()[edgeI];
01927
01928 bool neighbourIsRefined= false;
01929
01930 forAll(myFaces, myFaceI)
01931 {
01932 if (refineStatus[myFaces[myFaceI]] != NONE)
01933 {
01934 neighbourIsRefined = true;
01935 }
01936 }
01937
01938
01939 if (!neighbourIsRefined)
01940 {
01941
01942 const edge& e = surf.edges()[edgeI];
01943
01944 point mid =
01945 0.5
01946 * (
01947 surf.localPoints()[e.start()]
01948 + surf.localPoints()[e.end()]
01949 );
01950
01951 newPoints[newPointI] = mid;
01952
01953
01954 forAll(myFaces, myFaceI)
01955 {
01956
01957 greenRefine
01958 (
01959 surf,
01960 myFaces[myFaceI],
01961 edgeI,
01962 newPointI,
01963 newFaces
01964 );
01965
01966
01967 refineStatus[myFaces[myFaceI]] = GREEN;
01968 }
01969
01970 newPointI++;
01971 }
01972 }
01973
01974
01975 forAll(surf.localFaces(), faceI)
01976 {
01977 if (refineStatus[faceI] == NONE)
01978 {
01979 newFaces.append(surf.localFaces()[faceI]);
01980 }
01981 }
01982
01983 newFaces.shrink();
01984 newPoints.setSize(newPointI);
01985
01986 return triSurface(newFaces, surf.patches(), newPoints, true);
01987 }
01988
01989
01990
01991
01992
01993
01994
01995
01996 Foam::label Foam::triSurfaceTools::minEdge
01997 (
01998 const triSurface& surf,
01999 const labelList& edgeIndices
02000 )
02001 {
02002 scalar minLength = GREAT;
02003 label minIndex = -1;
02004 forAll(edgeIndices, i)
02005 {
02006 const edge& e = surf.edges()[edgeIndices[i]];
02007
02008 scalar length =
02009 mag
02010 (
02011 surf.localPoints()[e.end()]
02012 - surf.localPoints()[e.start()]
02013 );
02014
02015 if (length < minLength)
02016 {
02017 minLength = length;
02018 minIndex = i;
02019 }
02020 }
02021 return edgeIndices[minIndex];
02022 }
02023
02024
02025
02026 Foam::label Foam::triSurfaceTools::maxEdge
02027 (
02028 const triSurface& surf,
02029 const labelList& edgeIndices
02030 )
02031 {
02032 scalar maxLength = -GREAT;
02033 label maxIndex = -1;
02034 forAll(edgeIndices, i)
02035 {
02036 const edge& e = surf.edges()[edgeIndices[i]];
02037
02038 scalar length =
02039 mag
02040 (
02041 surf.localPoints()[e.end()]
02042 - surf.localPoints()[e.start()]
02043 );
02044
02045 if (length > maxLength)
02046 {
02047 maxLength = length;
02048 maxIndex = i;
02049 }
02050 }
02051 return edgeIndices[maxIndex];
02052 }
02053
02054
02055
02056 Foam::triSurface Foam::triSurfaceTools::mergePoints
02057 (
02058 const triSurface& surf,
02059 const scalar mergeTol
02060 )
02061 {
02062 pointField newPoints(surf.nPoints());
02063
02064 labelList pointMap(surf.nPoints());
02065
02066 bool hasMerged = Foam::mergePoints
02067 (
02068 surf.localPoints(),
02069 mergeTol,
02070 false,
02071 pointMap,
02072 newPoints
02073 );
02074
02075 if (hasMerged)
02076 {
02077
02078
02079
02080 List<labelledTri> newTriangles(surf.size());
02081 label newTriangleI = 0;
02082
02083 forAll(surf, faceI)
02084 {
02085 const labelledTri& f = surf.localFaces()[faceI];
02086
02087 label newA = pointMap[f[0]];
02088 label newB = pointMap[f[1]];
02089 label newC = pointMap[f[2]];
02090
02091 if ((newA != newB) && (newA != newC) && (newB != newC))
02092 {
02093 newTriangles[newTriangleI++] =
02094 labelledTri(newA, newB, newC, f.region());
02095 }
02096 }
02097 newTriangles.setSize(newTriangleI);
02098
02099 return triSurface
02100 (
02101 newTriangles,
02102 surf.patches(),
02103 newPoints
02104 );
02105 }
02106 else
02107 {
02108 return surf;
02109 }
02110 }
02111
02112
02113
02114 Foam::vector Foam::triSurfaceTools::surfaceNormal
02115 (
02116 const triSurface& surf,
02117 const label nearestFaceI,
02118 const point& nearestPt
02119 )
02120 {
02121 const labelledTri& f = surf[nearestFaceI];
02122 const pointField& points = surf.points();
02123
02124 label nearType;
02125 label nearLabel;
02126 triPointRef
02127 (
02128 points[f[0]],
02129 points[f[1]],
02130 points[f[2]]
02131 ).classify(nearestPt, 1E-6, nearType, nearLabel);
02132
02133 if (nearType == triPointRef::NONE)
02134 {
02135
02136 return surf.faceNormals()[nearestFaceI];
02137 }
02138 else if (nearType == triPointRef::EDGE)
02139 {
02140
02141 label edgeI = surf.faceEdges()[nearestFaceI][nearLabel];
02142
02143
02144 const labelList& eFaces = surf.edgeFaces()[edgeI];
02145
02146 vector edgeNormal(vector::zero);
02147
02148 forAll(eFaces, i)
02149 {
02150 edgeNormal += surf.faceNormals()[eFaces[i]];
02151 }
02152 return edgeNormal/(mag(edgeNormal) + VSMALL);
02153 }
02154 else
02155 {
02156
02157 const labelledTri& localF = surf.localFaces()[nearestFaceI];
02158 return surf.pointNormals()[localF[nearLabel]];
02159 }
02160 }
02161
02162
02163 Foam::triSurfaceTools::sideType Foam::triSurfaceTools::edgeSide
02164 (
02165 const triSurface& surf,
02166 const point& sample,
02167 const point& nearestPoint,
02168 const label edgeI
02169 )
02170 {
02171 const labelList& eFaces = surf.edgeFaces()[edgeI];
02172
02173 if (eFaces.size() != 2)
02174 {
02175
02176 return UNKNOWN;
02177 }
02178 else
02179 {
02180 const vectorField& faceNormals = surf.faceNormals();
02181
02182
02183
02184
02185 vector n = 0.5*(faceNormals[eFaces[0]] + faceNormals[eFaces[1]]);
02186
02187 if (((sample - nearestPoint) & n) > 0)
02188 {
02189 return OUTSIDE;
02190 }
02191 else
02192 {
02193 return INSIDE;
02194 }
02195 }
02196 }
02197
02198
02199
02200 Foam::triSurfaceTools::sideType Foam::triSurfaceTools::surfaceSide
02201 (
02202 const triSurface& surf,
02203 const point& sample,
02204 const label nearestFaceI,
02205 const point& nearestPoint,
02206 const scalar tol
02207 )
02208 {
02209 const labelledTri& f = surf[nearestFaceI];
02210 const pointField& points = surf.points();
02211
02212
02213
02214 label nearType, nearLabel;
02215 triPointRef
02216 (
02217 points[f[0]],
02218 points[f[1]],
02219 points[f[2]]
02220 ).classify(nearestPoint, tol, nearType, nearLabel);
02221
02222 if (nearType == triPointRef::NONE)
02223 {
02224
02225 scalar c = (sample - nearestPoint) & surf.faceNormals()[nearestFaceI];
02226
02227 if (c > 0)
02228 {
02229 return OUTSIDE;
02230 }
02231 else
02232 {
02233 return INSIDE;
02234 }
02235 }
02236 else if (nearType == triPointRef::EDGE)
02237 {
02238
02239
02240
02241
02242 label edgeI = surf.faceEdges()[nearestFaceI][nearLabel];
02243
02244
02245
02246
02247
02248
02249
02250
02251
02252
02253
02254
02255
02256
02257
02258
02259
02260
02261
02262
02263
02264
02265
02266 return edgeSide(surf, sample, nearestPoint, edgeI);
02267 }
02268 else
02269 {
02270
02271
02272
02273
02274
02275 const labelledTri& localF = surf.localFaces()[nearestFaceI];
02276 label nearPointI = localF[nearLabel];
02277
02278 const edgeList& edges = surf.edges();
02279 const pointField& localPoints = surf.localPoints();
02280 const point& base = localPoints[nearPointI];
02281
02282 const labelList& pEdges = surf.pointEdges()[nearPointI];
02283
02284 scalar minDistSqr = Foam::sqr(GREAT);
02285 label minEdgeI = -1;
02286
02287 forAll(pEdges, i)
02288 {
02289 label edgeI = pEdges[i];
02290
02291 const edge& e = edges[edgeI];
02292
02293 label otherPointI = e.otherVertex(nearPointI);
02294
02295
02296 vector eVec(localPoints[otherPointI] - base);
02297 scalar magEVec = mag(eVec);
02298
02299 if (magEVec > VSMALL)
02300 {
02301 eVec /= magEVec;
02302
02303
02304 const point perturbPoint = base + eVec;
02305
02306 scalar distSqr = Foam::magSqr(sample - perturbPoint);
02307
02308 if (distSqr < minDistSqr)
02309 {
02310 minDistSqr = distSqr;
02311 minEdgeI = edgeI;
02312 }
02313 }
02314 }
02315
02316 if (minEdgeI == -1)
02317 {
02318 FatalErrorIn("treeDataTriSurface::getSide")
02319 << "Problem: did not find edge closer than " << minDistSqr
02320 << abort(FatalError);
02321 }
02322
02323 return edgeSide(surf, sample, nearestPoint, minEdgeI);
02324 }
02325 }
02326
02327
02328
02329 Foam::triSurface Foam::triSurfaceTools::triangulate
02330 (
02331 const polyBoundaryMesh& bMesh,
02332 const labelHashSet& includePatches,
02333 const bool verbose
02334 )
02335 {
02336 const polyMesh& mesh = bMesh.mesh();
02337
02338
02339 DynamicList<labelledTri> triangles
02340 (
02341 mesh.nFaces() - mesh.nInternalFaces()
02342 );
02343
02344 label newPatchI = 0;
02345
02346 for
02347 (
02348 labelHashSet::const_iterator iter = includePatches.begin();
02349 iter != includePatches.end();
02350 ++iter
02351 )
02352 {
02353 label patchI = iter.key();
02354
02355 const polyPatch& patch = bMesh[patchI];
02356
02357 const pointField& points = patch.points();
02358
02359 label nTriTotal = 0;
02360
02361 forAll(patch, patchFaceI)
02362 {
02363 const face& f = patch[patchFaceI];
02364
02365 faceList triFaces(f.nTriangles(points));
02366
02367 label nTri = 0;
02368
02369 f.triangles(points, nTri, triFaces);
02370
02371 forAll(triFaces, triFaceI)
02372 {
02373 const face& f = triFaces[triFaceI];
02374
02375 triangles.append(labelledTri(f[0], f[1], f[2], newPatchI));
02376
02377 nTriTotal++;
02378 }
02379 }
02380
02381 if (verbose)
02382 {
02383 Pout<< patch.name() << " : generated " << nTriTotal
02384 << " triangles from " << patch.size() << " faces with"
02385 << " new patchid " << newPatchI << endl;
02386 }
02387
02388 newPatchI++;
02389 }
02390 triangles.shrink();
02391
02392
02393 triSurface rawSurface(triangles, mesh.points());
02394
02395
02396 triSurface surface
02397 (
02398 rawSurface.localFaces(),
02399 rawSurface.localPoints()
02400 );
02401
02402
02403 surface.patches().setSize(newPatchI);
02404
02405 newPatchI = 0;
02406
02407 for
02408 (
02409 labelHashSet::const_iterator iter = includePatches.begin();
02410 iter != includePatches.end();
02411 ++iter
02412 )
02413 {
02414 label patchI = iter.key();
02415
02416 const polyPatch& patch = bMesh[patchI];
02417
02418 surface.patches()[newPatchI].name() = patch.name();
02419
02420 surface.patches()[newPatchI].geometricType() = patch.type();
02421
02422 newPatchI++;
02423 }
02424
02425 return surface;
02426 }
02427
02428
02429
02430 Foam::triSurface Foam::triSurfaceTools::triangulateFaceCentre
02431 (
02432 const polyBoundaryMesh& bMesh,
02433 const labelHashSet& includePatches,
02434 const bool verbose
02435 )
02436 {
02437 const polyMesh& mesh = bMesh.mesh();
02438
02439
02440 const pointField& points = mesh.points();
02441 const pointField& faceCentres = mesh.faceCentres();
02442
02443 pointField newPoints(points.size() + faceCentres.size());
02444
02445 label newPointI = 0;
02446
02447 forAll(points, pointI)
02448 {
02449 newPoints[newPointI++] = points[pointI];
02450 }
02451 forAll(faceCentres, faceI)
02452 {
02453 newPoints[newPointI++] = faceCentres[faceI];
02454 }
02455
02456
02457
02458 DynamicList<labelledTri> triangles
02459 (
02460 mesh.nFaces() - mesh.nInternalFaces()
02461 );
02462
02463 label newPatchI = 0;
02464
02465 for
02466 (
02467 labelHashSet::const_iterator iter = includePatches.begin();
02468 iter != includePatches.end();
02469 ++iter
02470 )
02471 {
02472 label patchI = iter.key();
02473
02474 const polyPatch& patch = bMesh[patchI];
02475
02476 label nTriTotal = 0;
02477
02478 forAll(patch, patchFaceI)
02479 {
02480
02481 const face& f = patch[patchFaceI];
02482
02483
02484 label fc = points.size() + patchFaceI + patch.start();
02485
02486 forAll(f, fp)
02487 {
02488 label fp1 = (fp + 1) % f.size();
02489
02490 triangles.append(labelledTri(f[fp], f[fp1], fc, newPatchI));
02491
02492 nTriTotal++;
02493 }
02494 }
02495
02496 if (verbose)
02497 {
02498 Pout<< patch.name() << " : generated " << nTriTotal
02499 << " triangles from " << patch.size() << " faces with"
02500 << " new patchid " << newPatchI << endl;
02501 }
02502
02503 newPatchI++;
02504 }
02505 triangles.shrink();
02506
02507
02508
02509 triSurface rawSurface(triangles, newPoints);
02510
02511
02512 triSurface surface
02513 (
02514 rawSurface.localFaces(),
02515 rawSurface.localPoints()
02516 );
02517
02518
02519 surface.patches().setSize(newPatchI);
02520
02521 newPatchI = 0;
02522
02523 for
02524 (
02525 labelHashSet::const_iterator iter = includePatches.begin();
02526 iter != includePatches.end();
02527 ++iter
02528 )
02529 {
02530 label patchI = iter.key();
02531
02532 const polyPatch& patch = bMesh[patchI];
02533
02534 surface.patches()[newPatchI].name() = patch.name();
02535
02536 surface.patches()[newPatchI].geometricType() = patch.type();
02537
02538 newPatchI++;
02539 }
02540
02541 return surface;
02542 }
02543
02544
02545 Foam::triSurface Foam::triSurfaceTools::delaunay2D(const List<vector2D>& pts)
02546 {
02547
02548
02549 List<doubleScalar> geompackVertices(2*pts.size());
02550 label doubleI = 0;
02551 forAll(pts, i)
02552 {
02553 geompackVertices[doubleI++] = pts[i][0];
02554 geompackVertices[doubleI++] = pts[i][1];
02555 }
02556
02557
02558 int m2 = 3;
02559 List<int> triangle_node(m2*3*pts.size());
02560 List<int> triangle_neighbor(m2*3*pts.size());
02561
02562
02563 int nTris = 0;
02564 dtris2
02565 (
02566 pts.size(),
02567 geompackVertices.begin(),
02568 &nTris,
02569 triangle_node.begin(),
02570 triangle_neighbor.begin()
02571 );
02572
02573
02574 triangle_node.setSize(3*nTris);
02575 triangle_neighbor.setSize(3*nTris);
02576
02577
02578 List<labelledTri> faces(nTris);
02579
02580 forAll(faces, i)
02581 {
02582 faces[i] = labelledTri
02583 (
02584 triangle_node[3*i]-1,
02585 triangle_node[3*i+1]-1,
02586 triangle_node[3*i+2]-1,
02587 0
02588 );
02589 }
02590
02591 pointField points(pts.size());
02592 forAll(pts, i)
02593 {
02594 points[i][0] = pts[i][0];
02595 points[i][1] = pts[i][1];
02596 points[i][2] = 0.0;
02597 }
02598
02599 return triSurface(faces, points);
02600 }
02601
02602
02604
02605
02606
02607
02608
02609
02610
02611
02612
02613
02614
02615
02616
02617
02618
02619
02620
02621
02622
02623
02624
02625
02626
02627
02628
02629
02630
02631
02632
02633
02634
02635
02636
02637
02638
02639
02640
02641
02642
02643
02644
02645
02646
02647
02648
02649
02650
02651
02652
02653
02654
02655
02656
02657
02658
02659
02660
02661
02662
02663
02664
02665
02666
02667
02668
02669
02670
02671
02672
02673
02674
02675
02676
02677
02678
02679
02680
02681 void Foam::triSurfaceTools::calcInterpolationWeights
02682 (
02683 const triPointRef& tri,
02684 const point& p,
02685 FixedList<scalar, 3>& weights
02686 )
02687 {
02688
02689
02690 FixedList<vector, 3> edge;
02691 edge[0] = tri.c()-tri.b();
02692 edge[1] = tri.a()-tri.c();
02693 edge[2] = tri.b()-tri.a();
02694
02695 vector triangleFaceNormal = edge[1] ^ edge[2];
02696
02697
02698 FixedList<vector, 3> normal;
02699 for(label i=0; i<3; i++)
02700 {
02701 normal[i] = triangleFaceNormal ^ edge[i];
02702 normal[i] /= mag(normal[i]) + VSMALL;
02703 }
02704
02705 weights[0] = ((p-tri.b()) & normal[0]) / max(VSMALL, normal[0] & edge[1]);
02706 weights[1] = ((p-tri.c()) & normal[1]) / max(VSMALL, normal[1] & edge[2]);
02707 weights[2] = ((p-tri.a()) & normal[2]) / max(VSMALL, normal[2] & edge[0]);
02708 }
02709
02710
02711
02712
02713 void Foam::triSurfaceTools::calcInterpolationWeights
02714 (
02715 const triSurface& s,
02716 const pointField& samplePts,
02717 List<FixedList<label, 3> >& allVerts,
02718 List<FixedList<scalar, 3> >& allWeights
02719 )
02720 {
02721 allVerts.setSize(samplePts.size());
02722 allWeights.setSize(samplePts.size());
02723
02724 const pointField& points = s.points();
02725
02726 forAll(samplePts, i)
02727 {
02728 const point& samplePt = samplePts[i];
02729
02730
02731 FixedList<label, 3>& verts = allVerts[i];
02732 FixedList<scalar, 3>& weights = allWeights[i];
02733
02734 scalar minDistance = GREAT;
02735
02736 forAll(s, faceI)
02737 {
02738 const labelledTri& f = s[faceI];
02739
02740 triPointRef tri(f.tri(points));
02741
02742 pointHit nearest = tri.nearestPoint(samplePt);
02743
02744 if (nearest.hit())
02745 {
02746
02747 verts[0] = f[0];
02748 verts[1] = f[1];
02749 verts[2] = f[2];
02750
02751 calcInterpolationWeights(tri, nearest.rawPoint(), weights);
02752
02753
02754
02755
02756
02757
02758
02759 break;
02760 }
02761 else if (nearest.distance() < minDistance)
02762 {
02763 minDistance = nearest.distance();
02764
02765
02766 label nearType, nearLabel;
02767 tri.classify
02768 (
02769 nearest.rawPoint(),
02770 1E-6,
02771 nearType,
02772 nearLabel
02773 );
02774
02775 if (nearType == triPointRef::POINT)
02776 {
02777 verts[0] = f[nearLabel];
02778 weights[0] = 1;
02779 verts[1] = -1;
02780 weights[1] = -GREAT;
02781 verts[2] = -1;
02782 weights[2] = -GREAT;
02783
02784
02785
02786
02787
02788 }
02789 else if (nearType == triPointRef::EDGE)
02790 {
02791 verts[0] = f[nearLabel];
02792 verts[1] = f[f.fcIndex(nearLabel)];
02793 verts[2] = -1;
02794
02795 const point& p0 = points[verts[0]];
02796 const point& p1 = points[verts[1]];
02797
02798 scalar s = min
02799 (
02800 1,
02801 max
02802 (
02803 0,
02804 mag(nearest.rawPoint()-p0)/mag(p1-p0)
02805 )
02806 );
02807
02808
02809 weights[0] = 1-s;
02810 weights[1] = s;
02811 weights[2] = -GREAT;
02812
02813
02814
02815
02816
02817 }
02818 else
02819 {
02820
02821 verts[0] = f[0];
02822 verts[1] = f[1];
02823 verts[2] = f[2];
02824
02825 calcInterpolationWeights(tri, nearest.rawPoint(), weights);
02826
02827 break;
02828 }
02829 }
02830 }
02831 }
02832 }
02833
02834
02835
02836
02837
02838
02839
02840 Foam::surfaceLocation Foam::triSurfaceTools::classify
02841 (
02842 const triSurface& s,
02843 const label triI,
02844 const point& trianglePoint
02845 )
02846 {
02847 surfaceLocation nearest;
02848
02849
02850 label index, elemType;
02851
02852 triPointRef(s[triI].tri(s.points())).classify
02853 (
02854 trianglePoint,
02855 1E-6,
02856 elemType,
02857 index
02858 );
02859
02860 nearest.setPoint(trianglePoint);
02861
02862 if (elemType == triPointRef::NONE)
02863 {
02864 nearest.setHit();
02865 nearest.setIndex(triI);
02866 nearest.elementType() = triPointRef::NONE;
02867 }
02868 else if (elemType == triPointRef::EDGE)
02869 {
02870 nearest.setMiss();
02871 nearest.setIndex(s.faceEdges()[triI][index]);
02872 nearest.elementType() = triPointRef::EDGE;
02873 }
02874 else
02875 {
02876 nearest.setMiss();
02877 nearest.setIndex(s.localFaces()[triI][index]);
02878 nearest.elementType() = triPointRef::POINT;
02879 }
02880
02881 return nearest;
02882 }
02883
02884
02885 Foam::surfaceLocation Foam::triSurfaceTools::trackToEdge
02886 (
02887 const triSurface& s,
02888 const surfaceLocation& start,
02889 const surfaceLocation& end,
02890 const plane& cutPlane
02891 )
02892 {
02893
02894 surfaceLocation nearest = start;
02895 nearest.setMiss();
02896
02897
02898 snapToEnd(s, end, nearest);
02899
02900 if (!nearest.hit())
02901 {
02902
02903
02904 if (start.elementType() == triPointRef::NONE)
02905 {
02906
02907
02908
02909
02910
02911
02912 nearest = cutEdge
02913 (
02914 s,
02915 start.index(),
02916 -1,
02917 -1,
02918 start.rawPoint(),
02919 cutPlane,
02920 end.rawPoint()
02921 );
02922 nearest.elementType() = triPointRef::EDGE;
02923 nearest.triangle() = start.index();
02924 nearest.setMiss();
02925 }
02926 else if (start.elementType() == triPointRef::EDGE)
02927 {
02928
02929 const labelList& eFaces = s.edgeFaces()[start.index()];
02930
02931 nearest = visitFaces
02932 (
02933 s,
02934 eFaces,
02935 start,
02936 start.index(),
02937 -1,
02938 end,
02939 cutPlane
02940 );
02941 }
02942 else
02943 {
02944 const labelList& pFaces = s.pointFaces()[start.index()];
02945
02946 nearest = visitFaces
02947 (
02948 s,
02949 pFaces,
02950 start,
02951 -1,
02952 start.index(),
02953 end,
02954 cutPlane
02955 );
02956 }
02957 snapToEnd(s, end, nearest);
02958 }
02959 return nearest;
02960 }
02961
02962
02963 void Foam::triSurfaceTools::track
02964 (
02965 const triSurface& s,
02966 const surfaceLocation& endInfo,
02967 const plane& cutPlane,
02968 surfaceLocation& hitInfo
02969 )
02970 {
02971
02972
02973
02974
02975
02976
02977 while (true)
02978 {
02979
02980
02981
02982
02983 hitInfo = trackToEdge
02984 (
02985 s,
02986 hitInfo,
02987 endInfo,
02988 cutPlane
02989 );
02990
02991
02992
02993
02994
02995
02996
02997
02998 if (hitInfo.hit() || hitInfo.triangle() == -1)
02999 {
03000 break;
03001 }
03002 }
03003 }
03004
03005
03006