SplitImprove2 - cleanup, new point at min dist of edges

libsrc/gprim/geomtest3d.cpp

@@ -1035,6 +1035,39 @@ double MinDistLP2 (const Point3d & lp1, const Point3d & lp2, const Point3d & p)
 }
 
 
+double MinDistLP2 (const Point3d & lp1, const Point3d & lp2, const Point3d & p, double & lam)
+{
+  Vec3d v(lp1, lp2);
+  Vec3d vlp(lp1, p);
+
+  // dist(lam) = \| vlp \|^2 - 2 lam (v1p, v) + lam^2 \| v \|^2
+
+  // lam = (v * vlp) / (v * v);
+  // if (lam < 0) lam = 0;
+  // if (lam > 1) lam = 1;
+
+  double num = v*vlp;
+  double den = v*v;
+
+  if (num <= 0) 
+  {
+    lam = 0.0;
+    return Dist2 (lp1, p);
+  }
+
+  if (num >= den) 
+  {
+    lam = 1.0;
+    return Dist2 (lp2, p);
+  }
+  
+  lam = num/den;
+  if (den > 0)
+      return vlp.Length2() - num * num /den;
+  else
+    return vlp.Length2();
+}
+
 
 double MinDistTP2 (const Point3d & tp1, const Point3d & tp2, 
 		   const Point3d & tp3, const Point3d & p)
@@ -1102,7 +1135,7 @@ double MinDistTP2 (const Point3d & tp1, const Point3d & tp2,
 
 // 0 checks !!!
 double MinDistLL2 (const Point3d & l1p1, const Point3d & l1p2,
-		  const Point3d & l2p1, const Point3d & l2p2)
+		  const Point3d & l2p1, const Point3d & l2p2, double & lam1, double & lam2 )
 {
   // dist(lam1,lam2) = \| l2p1+lam2v2 - (l1p1+lam1 v1) \|
   // min !
@@ -1112,7 +1145,7 @@ double MinDistLL2 (const Point3d & l1p1, const Point3d & l1p2,
   Vec3d v2 (l2p1, l2p2);
 
   double a11, a12, a22, rs1, rs2;
-  double lam1, lam2, det;
+  double det;
 
   a11 = v1*v1;
   a12 = -(v1*v2);
@@ -1138,14 +1171,27 @@ double MinDistLL2 (const Point3d & l1p1, const Point3d & l1p2,
     }
 
   double minv, hv;
-  minv = MinDistLP2 (l1p1, l1p2, l2p1);
+  minv = MinDistLP2 (l1p1, l1p2, l2p1, lam1);
-  hv =  MinDistLP2 (l1p1, l1p2, l2p2);
+  lam2 = 0.;
-  if (hv < minv) minv = hv;
+  hv =  MinDistLP2 (l1p1, l1p2, l2p2, lam1);
+  if (hv < minv)
+  {
+    lam2 = 1.;
+    minv = hv;
+  }
 
-  hv =  MinDistLP2 (l2p1, l2p2, l1p1);
+  hv =  MinDistLP2 (l2p1, l2p2, l1p1, lam2);
-  if (hv < minv) minv = hv;
+  if (hv < minv)
-  hv =  MinDistLP2 (l2p1, l2p2, l1p2);
+  {
-  if (hv < minv) minv = hv;
+    lam1 = 0.;
+    minv = hv;
+  }
+  hv =  MinDistLP2 (l2p1, l2p2, l1p2, lam2);
+  if (hv < minv)
+  {
+    lam1 = 1.;
+    minv = hv;
+  }
 
   return minv;
 }

libsrc/gprim/geomtest3d.hpp

@@ -91,6 +91,9 @@ extern double MinDistTP2 (const Point3d & tp1, const Point3d & tp2,
 extern double MinDistLL2 (const Point3d & l1p1, const Point3d & l1p2,
 			  const Point3d & l2p1, const Point3d & l2p2);
 
+extern double MinDistLL2 (const Point3d & l1p1, const Point3d & l1p2,
+		  const Point3d & l2p1, const Point3d & l2p2, double & lam1, double & lam2 );
+
 }
 
 #endif

libsrc/meshing/improve3.cpp

@@ -3920,12 +3920,7 @@ void MeshOptimize3d :: SplitImprove2 (Mesh & mesh, OPTIMIZEGOAL goal)
     {
       if(mp.only3D_domain_nr && mp.only3D_domain_nr != mesh[ei].GetIndex())
 	continue;
-//       el_badness[ei] = CalcBad (mesh.Points(), mesh[ei], 0);
+      el_badness[ei] = CalcBad (mesh.Points(), mesh[ei], 0);
-      const auto & el = mesh[ei];
-      if(el.GetNP()==4)
-        el_badness[ei] = CalcTetBadness (mesh[el[0]], mesh[el[1]], mesh[el[2]], mesh[el[3]], 0, mp);
-      else
-        el_badness[ei] = 0.0;
     }
   });
 
@@ -3953,6 +3948,7 @@ void MeshOptimize3d :: SplitImprove2 (Mesh & mesh, OPTIMIZEGOAL goal)
 
     int minedge = -1;
     double mindist = 1e99;
+    double minlam0, minlam1;
 
     for (int i : Range(3))
     {
@@ -3967,11 +3963,14 @@ void MeshOptimize3d :: SplitImprove2 (Mesh & mesh, OPTIMIZEGOAL goal)
 //       Point3d p2 = mesh[pi2];
 //       Point3d p3 = mesh[pi3];
 
-      double dist = MinDistLL2(mesh[pi0], mesh[pi1], mesh[pi2], mesh[pi3]);
+      double lam0, lam1;
+      double dist = MinDistLL2(mesh[pi0], mesh[pi1], mesh[pi2], mesh[pi3], lam0, lam1 );
       if(dist<mindist)
       {
         mindist = dist;
         minedge = i;
+        minlam0 = lam0;
+        minlam1 = lam1;
       }
     }
 
@@ -3995,7 +3994,7 @@ void MeshOptimize3d :: SplitImprove2 (Mesh & mesh, OPTIMIZEGOAL goal)
     ArrayMem<ElementIndex, 50> has_both_points1;
 
     Point3d p[4] = { mesh[el[0]], mesh[el[1]], mesh[el[2]], mesh[el[3]] };
-    auto center = Center(p[0], p[1], p[2], p[3]);
+    auto center = Center(p[0]+minlam0*(p[1]-p[0]), p[2]+minlam1*(p[3]-p[2]));
     MeshPoint pnew;
 
     pnew(0) = center.X();
@@ -4074,13 +4073,8 @@ void MeshOptimize3d :: SplitImprove2 (Mesh & mesh, OPTIMIZEGOAL goal)
           if(newel1[i] == pi0) newel1[i] = pinew;
           if(newel2[i] == pi1) newel2[i] = pinew;
         }
-        if (!mesh.LegalTet (newel1)) cout << "newel1 illegal tet" << endl;
         mesh.AddVolumeElement (newel1);
-        if (!mesh.LegalTet (newel2)) cout << "newel2 illegal tet" << endl;
         mesh.AddVolumeElement (newel2);
-//         cout << "del element " << ei1 << ':' << oldel[0] << ',' << oldel[1] << ',' << oldel[2] << ',' << oldel[3] << endl;
-//         cout << "add element " << newel1[0] << ',' << newel1[1] << ',' << newel1[2] << ',' << newel1[3] << endl;
-//         cout << "add element " << newel2[0] << ',' << newel2[1] << ',' << newel2[2] << ',' << newel2[3] << endl;
       }
       for (auto ei1 : has_both_points1)
       {
@@ -4096,18 +4090,11 @@ void MeshOptimize3d :: SplitImprove2 (Mesh & mesh, OPTIMIZEGOAL goal)
           if(newel1[i] == pi2) newel1[i] = pinew;
           if(newel2[i] == pi3) newel2[i] = pinew;
         }
-//         cout << "del element " << ei1 << ':' << oldel[0] << ',' << oldel[1] << ',' << oldel[2] << ',' << oldel[3] << endl;
-//         cout << "add element " << newel1[0] << ',' << newel1[1] << ',' << newel1[2] << ',' << newel1[3] << endl;
-//         cout << "add element " << newel2[0] << ',' << newel2[1] << ',' << newel2[2] << ',' << newel2[3] << endl;
-        if (!mesh.LegalTet (newel1)) cout << "newel1 illegal tet" << endl;
         mesh.AddVolumeElement (newel1);
-        if (!mesh.LegalTet (newel2)) cout << "newel2 illegal tet" << endl;
         mesh.AddVolumeElement (newel2);
       }
     }
   }
-  MeshTopology mt(mesh);
-  mt.Update();
   mesh.Compress();
   bad1 = mesh.CalcTotalBad (mp);
   cout << "Total badness after splitimprove2= " << bad1 << endl;

libsrc/meshing/improve3.hpp

@@ -50,6 +50,14 @@ public:
     return 0;
   }
 
+  double 
+  CalcBadPow (const Mesh::T_POINTS & points, const Element & elem, double h)
+  {
+    if (elem.GetType() == TET)
+      return pow (max2(CalcTetBadness (points[elem[0]], points[elem[1]],  
+			     points[elem[2]], points[elem[3]], h, mp), 1e-10) , 1.0/mp.opterrpow);
+    return 0;
+  }
 
   double CalcTotalBad (const Mesh::T_POINTS & points, 
 		       const Array<Element> & elements)

libsrc/meshing/meshclass.cpp

@@ -3620,8 +3620,6 @@ namespace netgen
     //  FindOpenElements();
     timestamp = NextTimeStamp();
     lock.UnLock();
-    FindOpenElements();
-    CheckVolumeMesh ();
   }
 
   void Mesh :: OrderElements()