diff --git a/contrib/MeshOptimizer/MeshOptimizer.cpp b/contrib/MeshOptimizer/MeshOptimizer.cpp
index 61053b12c0ffd1a31a582de470ed7b6bf1e9edd2..4cf44f2d1acb7a1db9e317c82b7f8558d5e51a77 100644
--- a/contrib/MeshOptimizer/MeshOptimizer.cpp
+++ b/contrib/MeshOptimizer/MeshOptimizer.cpp
@@ -71,36 +71,96 @@ static std::set<MVertex *> getAllBndVertices
}
+// Test intersection between sphere and segment
+static bool testSegSphereIntersect(SPoint3 A, SPoint3 B, const SPoint3& P, const double rr)
+{
+ // Test if separating plane between sphere and segment vertices
+ // For each vertex, separation if vertex is outside sphere and P on opposite side
+ // to other seg. vertex w.r.t plane of normal (vertex-P) through vertex
+ const SVector3 PA(P, A), PB(P, B);
+ const double aa = dot(PA, PA), ab = dot(PA, PB);
+ if ((aa > rr) & (ab > aa)) return false;
+ const double bb = dot(PB, PB);
+ if ((bb > rr) & (ab > bb)) return false;
+
+ // Test if separating plane between sphere and line
+ // For A, separation if projection Q of P on (AB) lies outside the sphere
+ const SVector3 AB(A, B);
+ const double d = ab - aa, e = dot(AB, AB);
+ const SVector3 PQ = PA * e - d * AB;
+ if (dot(PQ, PQ) > rr * e * e) return false;
+
+ // Return true (intersection) if no separation at all
+ return true;
+}
+
+
+// Test intersection between sphere and triangle
+// Inspired by Christer Ericson, http://realtimecollisiondetection.net/blog/?p=103
+static bool testTriSphereIntersect(SPoint3 A, SPoint3 B, SPoint3 C, const SPoint3& P, const double rr)
+{
+ // Test if separating plane between sphere and triangle plane
+ const SVector3 PA(P, A), AB(A, B), AC(A, C);
+ const SVector3 V = crossprod(AB, AC); // Normal to triangle plane
+ const double d = dot(PA, V); // Dist. from P to triangle plane times norm of V
+ const double e = dot(V, V); // Norm of V
+ if (d * d > rr * e) return false; // Test if separating plane between sphere and triangle plane
+
+ // Test if separating plane between sphere and triangle vertices
+ const SVector3 PB(P, B), PC(P, B);
+ const double aa = dot(PA, PA), ab = dot(PA, PB), ac = dot(PA, PC);
+ const double bb = dot(PB, PB), bc = dot(PB, PC), cc = dot(PC, PC);
+ if ((aa > rr) & (ab > aa) & (ac > aa)) return false; // For each triangle vertex, separation if vertex is outside sphere
+ if ((bb > rr) & (ab > bb) & (bc > bb)) return false; // and P on opposite side to other two triangle vertices w.r.t
+ if ((cc > rr) & (ac > cc) & (bc > cc)) return false; // plane of normal (vertex-P) through vertex
+
+ // Test if separating plane between sphere and triangle edges
+ const SVector3 BC(B, C);
+ const double d1 = ab - aa, d2 = bc - bb, d3 = ac - cc;
+ const double e1 = dot(AB, AB), e2 = dot(BC, BC), e3 = dot(AC, AC);
+ const SVector3 PQ1 = PA * e1 - d1 * AB; // Q1 projection of P on line (AB)
+ const SVector3 PQ2 = PB * e2 - d2 * BC; // Q2 projection of P on line (BC)
+ const SVector3 PQ3 = PC * e3 + d3 * AC; // Q3 projection of P on line (AC)
+ const SVector3 PQC = PC * e1 - PQ1;
+ const SVector3 PQA = PA * e2 - PQ2;
+ const SVector3 PQB = PB * e3 - PQ3;
+ if ((dot(PQ1, PQ1) > rr * e1 * e1) & (dot(PQ1, PQC) > 0)) return false; // For A, separation if Q lies outside the sphere and if P and C
+ if ((dot(PQ2, PQ2) > rr * e2 * e2) & (dot(PQ2, PQA) > 0)) return false; // are on opposite sides of plane through AB with normal PQ
+ if ((dot(PQ3, PQ3) > rr * e3 * e3) & (dot(PQ3, PQB) > 0)) return false; // Same for other two vertices
+
+ // Return true (intersection) if no separation at all
+ return true;
+}
+
+
// Approximate test of intersection element with circle/sphere by sampling
-static bool testElInDist(const SPoint3 p, double limDist, MElement *el)
+static bool testElInDist(const SPoint3 P, double limDist, MElement *el)
{
const double sampleLen = 0.5*limDist; // Distance between sample points
+ const double limDistSq = limDist*limDist;
if (el->getDim() == 2) { // 2D?
for (int iEd = 0; iEd < el->getNumEdges(); iEd++) { // Loop over edges of element
MEdge ed = el->getEdge(iEd);
- const int nPts = int(ed.length()/sampleLen)+2; // Nb of sample points based on edge length
- for (int iPt = 0; iPt < nPts; iPt++) { // Loop over sample points
- const SPoint3 pt = ed.interpolate(iPt/float(nPts-1));
- if (p.distance(pt) < limDist) return true;
+ const SPoint3 A = ed.getVertex(0)->point();
+ const SPoint3 B = ed.getVertex(1)->point();
+ if (testSegSphereIntersect(A, B, P, limDistSq)) {
+ return true;
}
}
}
else { // 3D
for (int iFace = 0; iFace < el->getNumFaces(); iFace++) { // Loop over faces of element
MFace face = el->getFace(iFace);
- double lMax = 0.; // Max. edge length in face
- const int nVert = face.getNumVertices();
- for (int iEd = 0; iEd < nVert; iEd++)
- lMax = std::max(lMax, face.getEdge(iEd).length());
- const int nPts = int(lMax/sampleLen)+2; // Nb of sample points based on max. edge length in face
- for (int iPt0 = 0; iPt0 < nPts; iPt0++) {
- const double u = iPt0/float(nPts-1);
- for (int iPt1 = 0; iPt1 < nPts; iPt1++) { // Loop over sample points
- const double vMax = (nVert == 3) ? 1.-u : 1.;
- const SPoint3 pt = face.interpolate(u, vMax*iPt1/float(nPts-1));
- if (p.distance(pt) < limDist) return true;
- }
+ const SPoint3 A = face.getVertex(0)->point();
+ const SPoint3 B = face.getVertex(1)->point();
+ const SPoint3 C = face.getVertex(2)->point();
+ if (face.getNumVertices() == 3)
+ return testTriSphereIntersect(A, B, C, P, limDistSq);
+ else {
+ const SPoint3 D = face.getVertex(3)->point();
+ return (testTriSphereIntersect(A, B, C, P, limDistSq) ||
+ testTriSphereIntersect(A, C, D, P, limDistSq));
}
}
}
@@ -129,9 +189,10 @@ static std::set<MElement*> getSurroundingPatch(MElement *el, int minLayers,
const std::vector<MElement*> &neighbours = vertex2elements.find
((*it)->getVertex(i))->second;
for (std::vector<MElement*>::const_iterator itN = neighbours.begin();
- itN != neighbours.end(); ++itN){
- if ((d < minLayers) || testElInDist(pnt, limDist, *itN))
- if (patch.insert(*itN).second) currentLayer.push_back(*itN); // Assume that if an el is too far, its neighbours are too far as well
+ itN != neighbours.end(); ++itN) {
+ if ((patch.find(*itN) == patch.end()) &&
+ ((d < minLayers) || testElInDist(pnt, limDist, *itN)))
+ if (patch.insert(*itN).second) currentLayer.push_back(*itN); // Assume that if an el is too far, its neighbours are too far as well
}
}
}