stable pseudo-inverse, improve edge-analysis for revolution

libsrc/csg/revolution.cpp

@@ -948,6 +948,67 @@ namespace netgen
 
     return VecInSolid(p,v1+0.01*v2,eps);
   }
+
+  void Revolution ::
+  GetTangentialVecSurfaceIndices2 (const Point<3> & p, const Vec<3> & v1, const Vec<3> & v2,
+                                   NgArray<int> & surfind, double eps) const
+  {
+    *testout << "tangentialvecsurfind2, p = " << p << endl;
+    for (int i = 0; i < faces.Size(); i++)
+      if (faces[i]->PointInFace (p, eps))
+        {
+          *testout << "check face " << i << endl;
+          Point<2> p2d;
+          Vec<2> v12d;
+          faces[i]->CalcProj(p,p2d,v1,v12d);
+          *testout << "v12d = " << v12d << endl;
+          auto & spline = faces[i]->GetSpline();
+          if (Dist2 (spline.StartPI(), p2d) < sqr(eps))
+            {
+              *testout << "start pi" << endl;
+              Vec<2> tang = spline.GetTangent(0);
+              double ip = tang*v12d;
+              *testout << "ip = " << ip << endl;
+              if (ip > eps)
+                surfind.Append(GetSurfaceId(i));
+              else if (ip > -eps)
+                {
+                  Vec<2> v22d;
+                  faces[i]->CalcProj(p,p2d,v2,v22d);
+                  double ip2 = tang*v22d;
+                  *testout << "ip2 = " << ip2 << endl;
+                  if (ip2 > -eps)
+                    surfind.Append(GetSurfaceId(i));                    
+                }
+            }
+          else if (Dist2 (faces[i]->GetSpline().EndPI(), p2d) < sqr(eps))
+            {
+              *testout << "end pi" << endl;
+              
+              Vec<2> tang = spline.GetTangent(1);
+              double ip = tang*v12d;
+              *testout << "ip = " << ip << endl;
+              if (ip < -eps)
+                surfind.Append(GetSurfaceId(i));
+              else if (ip < eps)
+                {
+                  Vec<2> v22d;
+                  faces[i]->CalcProj(p,p2d,v2,v22d);
+                  double ip2 = tang*v22d;
+                  *testout << "ip2 = " << ip2 << endl;                  
+                  if (ip2 < eps)
+                    surfind.Append(GetSurfaceId(i));                    
+                }
+            }
+          else
+            {
+              *testout << "inner point" << endl;
+              surfind.Append(GetSurfaceId(i));
+            }
+        }
+  }
+
+
   
   int Revolution :: GetNSurfaces() const
   {

libsrc/csg/revolution.hpp

@@ -158,7 +158,10 @@ namespace netgen
 				      const Vec<3> & v2,
 				      double eps) const;
 
-  
+    virtual void GetTangentialVecSurfaceIndices2 (const Point<3> & p, const Vec<3> & v1, const Vec<3> & v2,
+						  NgArray<int> & surfind, double eps) const;
+
+    
     virtual int GetNSurfaces() const;
     virtual Surface & GetSurface (int i = 0);
     virtual const Surface & GetSurface (int i = 0) const;

libsrc/csg/specpoin.cpp

@@ -1908,8 +1908,11 @@ namespace netgen
 		    }
 
 
-		  if (Abs2 (t) < 1e-8)
+		  if (Abs2 (t) < 1e-16)
-		    continue;
+                    {
+                      cerr << "normal vectors degenerated" << endl;
+                      continue;
+                    }
 
 #ifdef DEVELOP
 		  *testout << "           tangential vector " << t << endl;
@@ -1950,10 +1953,15 @@ namespace netgen
 		      rhs(0) = -t * (hessej * t);
 		      rhs(1) = -t * (hessek * t);
 
-		      CalcInverse (mat, inv);
+                      CalcInverse (mat, inv);
-		      t2 = inv * rhs;
+                      t2 = inv * rhs;
+                      // t2 = StableSolve (mat, rhs);
+#ifdef DEVELOP
+                      *testout << "t = " << t << ", t2 = " << t2 << endl;
+#endif
 
-		      
+
+                      
 		      /*
 		      ageometry.GetIndependentSurfaceIndices 
 			(locsol, p, t, surfind2);
@@ -2010,7 +2018,9 @@ namespace netgen
 
 			  // locsol2 -> TangentialSolid2 (p, m1, locsol3, surfind3, 1e-9*geomsize);
 			  locsol -> TangentialEdgeSolid (p, t, t2, m1, locsol3, surfind3, ideps*geomsize);
-
+#ifdef DEVELOP
+                          (*testout) << "m1 = " << m1 << ", surfind3 = " << surfind3 << endl;
+#endif
 			  //ageometry.GetIndependentSurfaceIndices (surfind3);
 
 			  if (surfind3.Contains(surfind2[l]))
@@ -2019,6 +2029,9 @@ namespace netgen
 			  
 			  // locsol2 -> TangentialSolid2 (p, m2, locsol3, surfind3, 1e-9*geomsize);
 			  locsol -> TangentialEdgeSolid (p, t, t2, m2, locsol3, surfind3, ideps*geomsize); 
+#ifdef DEVELOP
+                          (*testout) << "m2 = " << m2 << ", surfind3 = " << surfind3 << endl;
+#endif
 
 			  // ageometry.GetIndependentSurfaceIndices (surfind3);
 

libsrc/csg/surface.hpp

@@ -293,9 +293,13 @@ namespace netgen
 				      const Vec<3> & m,
 				      double eps) const;
 
+    // for a point p in the surface, into which (closed) surfaces does v point into ? 
     virtual void GetTangentialVecSurfaceIndices (const Point<3> & p, const Vec<3> & v,
 						 NgArray<int> & surfind, double eps) const;
 
+    // a point p in the surface, and v a tangential vector
+    // for arbitrary small, but positive t consider q := Project(p+t*v)
+    // into which (closed) surfaces does v2 point into, when starting from q ?     
     virtual void GetTangentialVecSurfaceIndices2 (const Point<3> & p, const Vec<3> & v1, const Vec<3> & v2,
 						  NgArray<int> & surfind, double eps) const;
 

libsrc/gprim/geomfuncs.hpp

@@ -142,10 +142,24 @@ namespace netgen
 
   inline void CalcInverse (const Mat<2,3> & m, Mat<3,2> & inv)
   {
+    Vec<3> a0 = m.Row(0);
+    Vec<3> a1 = m.Row(1);
+    Vec<3> n = Cross(a0, a1);
+    Vec<3> d0 = Cross(a1, n);
+    Vec<3> d1 = Cross(a0, n);
+    double s0 = 1.0/(a0*d0);
+    double s1 = 1.0/(a1*d1);
+    for (int i = 0; i < 3; i++)
+      {
+        inv(i,0) = s0*d0(i);
+        inv(i,1) = s1*d1(i);
+      }
+    /*
     Mat<2,2> a = m * Trans (m);
     Mat<2,2> ainv;
     CalcInverse (a, ainv);
     inv = Trans (m) * ainv;
+    */
   }
 
   void CalcInverse (const Mat<3,2> & m, Mat<2,3> & inv);
@@ -166,6 +180,31 @@ namespace netgen
   void EigenValues (const Mat<3,3> & m, Vec<3> & ev);
   void EigenValues (const Mat<2,2> & m, Vec<3> & ev);
 
+
+  template <typename T> 
+  Vec<3,T> StableSolve (Mat<2,3,T> mat, Vec<2,T> rhs)
+  {
+    Vec<3> a0 = mat.Row(0);
+    Vec<3> a1 = mat.Row(1);
+    /*
+    Vec<3> d = Cross ( Cross (a0, a1), a0);
+
+    double alpha = rhs(0) / a0.Length2();
+    double beta = (rhs(1)-alpha* (a0*a1)) / (d*a1);
+    return alpha * a0 + beta * d;
+    */
+    Vec<3> n = Cross(a0, a1);
+    Vec<3> d0 = Cross(a1, n);
+    Vec<3> d1 = Cross(a0, n);
+    double alpha = rhs(0) / (a0*d0);
+    double beta = rhs(1) / (a1*d1);
+    return alpha * d0 + beta * d1;
+  }
+
+
+
+
+  
 }
 
 #endif