Optimized patch creation in mesh optimizer

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));
       }
     }
   }
@@ -130,7 +190,8 @@ static std::set<MElement*> getSurroundingPatch(MElement *el, int minLayers,
           ((*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.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
         }
       }