53 void Foam::edgeIntersections::checkEdges(
const triSurface& surf)
55 const pointField& localPoints = surf.localPoints();
56 const edgeList& edges = surf.edges();
59 treeBoundBox bb(localPoints);
61 scalar minSize = SMALL * bb.minDim();
65 const edge&
e = edges[edgeI];
67 scalar eMag = e.mag(localPoints);
73 "Foam::edgeIntersections::checkEdges(const triSurface& surf)"
74 ) <<
"Edge " << edgeI <<
" vertices " << e
75 <<
" coords:" << localPoints[e[0]] <<
' '
76 << localPoints[e[1]] <<
" is very small compared to bounding"
77 <<
" box dimensions " << bb <<
endl
78 <<
"This might lead to problems in intersection"
82 if (edgeFaces[edgeI].
size() == 1)
86 "Foam::edgeIntersections::checkEdges(const triSurface& surf)"
87 ) <<
"Edge " << edgeI <<
" vertices " << e
88 <<
" coords:" << localPoints[e[0]] <<
' '
89 << localPoints[e[1]] <<
" has only one face connected to it:"
90 << edgeFaces[edgeI] <<
endl
91 <<
"This might lead to problems in intersection"
99 void Foam::edgeIntersections::intersectEdges
101 const triSurface& surf1,
103 const triSurfaceSearch& querySurf2,
108 const triSurface& surf2 = querySurf2.surface();
111 const labelList& meshPoints = surf1.meshPoints();
115 Pout<<
"Calculating intersection of " << edgeLabels.size() <<
" edges"
116 <<
" out of " << surf1.nEdges() <<
" with " << surf2.size()
117 <<
" triangles ..." <<
endl;
121 const indexedOctree<treeDataTriSurface>& tree = querySurf2.tree();
130 label edgeI = edgeLabels[i];
132 if (debug && (i % 1000 == 0))
134 Pout<<
"Intersecting edge " << edgeI <<
" with surface" <<
endl;
137 const edge& e = surf1.edges()[edgeI];
139 const point& pStart = points1[meshPoints[e.start()]];
140 const point& pEnd = points1[meshPoints[e.end()]];
142 const vector eVec(pEnd - pStart);
143 const scalar eMag =
mag(eVec);
144 const vector n(eVec/(eMag + VSMALL));
147 const point tolVec = 1e-6*eVec;
154 point p0 = pStart - 0.5*surf1PointTol[e[0]]*n;
155 const point p1 = pEnd + 0.5*surf1PointTol[e[1]]*n;
156 const scalar maxS =
mag(p1 - pStart);
160 DynamicList<pointIndexHit> currentIntersections(100);
161 DynamicList<label> currentIntersectionTypes(100);
171 currentIntersections.append(pHit);
176 if (
mag(pHit.hitPoint() - pStart) < surf1PointTol[e[0]])
180 else if (
mag(pHit.hitPoint() - pEnd) < surf1PointTol[e[1]])
184 else if (
mag(n & normals2[pHit.index()]) < alignedCos_)
186 Pout<<
"Flat angle edge:" << edgeI
187 <<
" face:" << pHit.index()
188 <<
" cos:" <<
mag(n & normals2[pHit.index()])
193 currentIntersectionTypes.append(edgeEnd);
203 p0 = pHit.hitPoint() + tolVec;
205 if (((p0-pStart) & n) >= maxS)
220 operator[](edgeI).transfer(currentIntersections);
221 classification_[edgeI].transfer(currentIntersectionTypes);
226 Pout<<
"Found " << nHits <<
" intersections of edges with surface ..."
239 bool Foam::edgeIntersections::inlinePerturb
241 const triSurface& surf1,
249 bool hasPerturbed =
false;
255 const labelList& edgeEnds = classification_[edgeI];
259 bool perturbStart =
false;
260 bool perturbEnd =
false;
263 if (edgeEnds[0] == 0)
268 if (edgeEnds[edgeEnds.size()-1] == 1)
274 if (perturbStart || perturbEnd)
276 const edge& e = surf1.edges()[edgeI];
278 label v0 = surf1.meshPoints()[e[0]];
279 label
v1 = surf1.meshPoints()[e[1]];
281 vector eVec(points1[v1] - points1[v0]);
287 scalar t = 4.0*(rndGen.scalar01() - 0.5);
288 points1[v0] += t*surf1PointTol[e[0]]*n;
290 const labelList& pEdges = surf1.pointEdges()[e[0]];
294 affectedEdges[pEdges[i]] =
true;
300 scalar t = 4.0*(rndGen.scalar01() - 0.5);
301 points1[
v1] += t*surf1PointTol[e[1]]*n;
303 const labelList& pEdges = surf1.pointEdges()[e[1]];
307 affectedEdges[pEdges[i]] =
true;
319 bool Foam::edgeIntersections::rotatePerturb
321 const triSurface& surf1,
330 const labelList& meshPoints = surf1.meshPoints();
332 const labelList& edgeEnds = classification_[edgeI];
334 bool hasPerturbed =
false;
338 if (edgeEnds[i] == 2)
340 const edge& e = surf1.edges()[edgeI];
343 label pointI = e[rndGen.bit()];
347 vector rndVec = rndGen.vector01() -
vector(0.5, 0.5, 0.5);
350 vector n(points1[meshPoints[e[1]]] - points1[meshPoints[e[0]]]);
351 scalar magN =
mag(n) + VSMALL;
354 rndVec -= n*(n & rndVec);
357 rndVec /=
mag(rndVec) + VSMALL;
362 Pout<<
"rotating: shifting endpoint " << meshPoints[pointI]
363 <<
" of edge:" << edgeI <<
" verts:"
364 << points1[meshPoints[e[0]]] <<
' '
365 << points1[meshPoints[e[1]]]
367 <<
" tol:" << surf1PointTol[pointI] <<
endl;
369 points1[meshPoints[pointI]] += rndVec;
372 const labelList& pEdges = surf1.pointEdges()[pointI];
376 affectedEdges[pEdges[i]] =
true;
392 bool Foam::edgeIntersections::offsetPerturb
394 const triSurface& surf1,
395 const triSurface& surf2,
403 const labelList& meshPoints = surf1.meshPoints();
405 const edge& e = surf1.edges()[edgeI];
407 const List<pointIndexHit>& hits = operator[](edgeI);
410 bool hasPerturbed =
false;
418 label surf2FaceI = pHit.index();
420 const labelledTri& f2 = surf2.localFaces()[surf2FaceI];
422 const pointField& surf2Pts = surf2.localPoints();
437 scalar tolDim = 0.001*
mag(tri.a() - ctr);
439 tri.classify(pHit.hitPoint(), tolDim, nearType, nearLabel);
444 vector offset = 0.01*rndGen.scalar01()*(ctr - pHit.hitPoint());
447 points1[meshPoints[e[0]]] += offset;
450 const labelList& pEdges0 = surf1.pointEdges()[e[0]];
454 affectedEdges[pEdges0[i]] =
true;
458 points1[meshPoints[e[1]]] += offset;
461 const labelList& pEdges1 = surf1.pointEdges()[e[1]];
465 affectedEdges[pEdges1[i]] =
true;
498 classification_(surf1.
nEdges())
532 classification_(classification)
548 const labelList& pEdges = pointEdges[pointI];
550 scalar minDist = GREAT;
554 minDist =
min(minDist, edges[pEdges[i]].
mag(localPoints));
557 minLen[pointI] = minDist;
590 for (; iter < nIters; iter++)
602 label edgeI = edgesToTest[i];
605 if (!affectedEdges[edgeI])
609 bool shiftedEdgeEndPoints =
620 nShifted += (shiftedEdgeEndPoints ? 1 : 0);
622 if (!shiftedEdgeEndPoints)
635 nRotated += (rotatedEdge ? 1 : 0);
643 bool offsetEdgePoints =
654 nOffset += (offsetEdgePoints ? 1 : 0);
662 Pout<<
"Edges to test : " <<
nl
663 <<
" total:" << edgesToTest.size() <<
nl
664 <<
" resolved by:" <<
nl
665 <<
" shifting : " << nShifted <<
nl
666 <<
" rotating : " << nRotated <<
nl
667 <<
" offsetting : " << nOffset <<
nl
672 if (nShifted == 0 && nRotated == 0 && nOffset == 0)
683 forAll(affectedEdges, edgeI)
685 if (affectedEdges[edgeI])
687 newEdgesToTest[newEdgeI++] = edgeI;
690 newEdgesToTest.setSize(newEdgeI);
694 Pout<<
"Edges to test:" <<
nl
695 <<
" was : " << edgesToTest.size() <<
nl
696 <<
" is : " << newEdgesToTest.size() <<
nl
701 edgesToTest.transfer(newEdgesToTest);
703 if (edgesToTest.empty())