diff --git a/libsrc/csg/polyhedra.cpp b/libsrc/csg/polyhedra.cpp
index ce35e108..7db4f94c 100644
--- a/libsrc/csg/polyhedra.cpp
+++ b/libsrc/csg/polyhedra.cpp
@@ -540,7 +540,7 @@ namespace netgen
             double lam3 = 1-lam1-lam2;
 
             if (lam1 >= -eps_base1 && lam2 >= -eps_base1 && lam3 >= -eps_base1)
-              {  // point is close to trianlge, perturbe by alpha*v 
+              {  // point is close to trianlge, perturbe by alpha*v
                 double dlamn = face.nn*v;
 
                 if (fabs(dlamn) < 1e-8) // vec also in plane
diff --git a/libsrc/csg/solid.cpp b/libsrc/csg/solid.cpp
index 58ed069e..70156c35 100644
--- a/libsrc/csg/solid.cpp
+++ b/libsrc/csg/solid.cpp
@@ -285,6 +285,7 @@ namespace netgen
   }
 
 
+  /*
   bool Solid::VectorIn2 (const Point<3> & p, const Vec<3> & v1, 
 			const Vec<3> & v2, double eps) const
   {
@@ -317,10 +318,51 @@ namespace netgen
       }
     return 0;  
   }
+  */
 
 
+  bool Solid::VectorIn2 (const Point<3> & p, const Vec<3> & v1, 
+                         const Vec<3> & v2, double eps) const
+  {
+    switch (op)
+      {
+      case TERM: case TERM_REF:
+        {
+          auto res = prim->VecInSolid2 (p, v1, v2, eps);
+          return res != IS_OUTSIDE;
+        }
+      case SECTION:
+	return s1->VectorIn2 (p, v1, v2, eps) && s2->VectorIn2 (p, v1, v2, eps);
+      case UNION:
+	return s1->VectorIn2 (p, v1, v2, eps) || s2->VectorIn2 (p, v1, v2, eps);
+      case SUB:
+	return !s1->VectorStrictIn2 (p, v1, v2, eps);
+      case ROOT:
+	return s1->VectorIn2 (p, v1, v2, eps);
+      }
+    // return 0;  
+  }
 
-
+  bool Solid :: VectorStrictIn2 (const Point<3> & p, const Vec<3> & v1, const Vec<3> & v2,
+                                 double eps) const
+  {
+    switch (op)
+      {
+      case TERM: case TERM_REF:
+        {
+          auto res = prim->VecInSolid2 (p, v1, v2, eps);
+          return (res == IS_INSIDE);
+        }
+      case SECTION:
+	return s1->VectorStrictIn2 (p, v1, v2, eps) && s2->VectorStrictIn2 (p, v1, v2, eps);
+      case UNION:
+	return s1->VectorStrictIn2 (p, v1, v2, eps) || s2->VectorStrictIn2 (p, v1, v2, eps);
+      case SUB:
+	return !s1->VectorIn2 (p, v1, v2, eps);
+      case ROOT:
+	return s1->VectorStrictIn2 (p, v1, v2, eps);
+      }
+  }
 
 
   void Solid :: Print (ostream & str) const
diff --git a/libsrc/csg/solid.hpp b/libsrc/csg/solid.hpp
index 5f3d6688..d3364203 100644
--- a/libsrc/csg/solid.hpp
+++ b/libsrc/csg/solid.hpp
@@ -109,9 +109,12 @@ namespace netgen
   
     bool VectorIn2 (const Point<3> & p, const Vec<3> & v1, const Vec<3> & v2,
 		    double eps) const;
+    /*
     bool VectorIn2Rec (const Point<3> & p, const Vec<3> & v1, const Vec<3> & v2,
 		       double eps) const;
-
+    */
+    bool VectorStrictIn2 (const Point<3> & p, const Vec<3> & v1, const Vec<3> & v2,
+                          double eps) const;
 
     /// compute localization in point p
     unique_ptr<Solid> TangentialSolid (const Point<3> & p, NgArray<int> & surfids, double eps) const;
diff --git a/libsrc/csg/specpoin.cpp b/libsrc/csg/specpoin.cpp
index 99c33377..c69a7239 100644
--- a/libsrc/csg/specpoin.cpp
+++ b/libsrc/csg/specpoin.cpp
@@ -2076,16 +2076,43 @@ namespace netgen
 			    continue;
 			  
 			  Vec<3> s = Cross (normalvecs[m], t);
+
 			  Vec<3> t2a = t + 0.01 *s;
 			  Vec<3> t2b = t - 0.01 *s;
-			  
-			  bool isface =
+
+			  bool isfaceold =
 			    (locsol->VectorIn (p, t2a, 1e-6*geomsize) &&
 			     !locsol->VectorStrictIn (p, t2a, 1e-6*geomsize))
 			    ||
 			    (locsol->VectorIn (p, t2b, 1e-6*geomsize) &&
 			     !locsol->VectorStrictIn (p, t2b, 1e-6*geomsize));
-			  
+
+                          bool isfacenew =
+                            (locsol->VectorIn2 (p, t, s, 1e-6*geomsize) && !locsol->VectorStrictIn2 (p, t, s, 1e-6*geomsize)) || 
+                            (locsol->VectorIn2 (p, t, -s, 1e-6*geomsize) && !locsol->VectorStrictIn2 (p, t, -s, 1e-6*geomsize));
+
+                          bool isface = isfacenew;
+                          
+                          if (isfaceold != isfacenew)
+                            {
+                              *testout << "different, p = " << p << ", t = " << t << ", s = " << s << endl;
+                              *testout << "tlo = " << si << endl;
+                              *testout << "isface, old = " << isface << ", isfacenew = " << isfacenew << endl;
+                              
+                              *testout << "t2a = " << t2a << endl;
+                              *testout << "vecin(p,t2a) = " << locsol->VectorIn (p, t2a, 1e-6*geomsize) << endl;
+                              *testout << "vecstrictin(p,t2a) = " << locsol->VectorStrictIn (p, t2a, 1e-6*geomsize) << endl;
+                              *testout << "vectorin2 = " << locsol->VectorIn2 (p, t, s, 1e-6*geomsize)  << endl;
+                              *testout << "vectorstrictin2 = " << locsol->VectorStrictIn2 (p, t, s, 1e-6*geomsize)  << endl;
+                              
+                              *testout << "t2b = " << t2b << endl;
+                              *testout << "vecin(p,t2b) = " << locsol->VectorIn (p, t2b, 1e-6*geomsize) << endl;
+                              *testout << "vecstrictin(p,t2b) = " << locsol->VectorStrictIn (p, t2b, 1e-6*geomsize) << endl;
+                              *testout << "vectorin2- = " << locsol->VectorIn2 (p, t, -s, 1e-6*geomsize)  << endl;
+                              *testout << "vectorstrictin2- = " << locsol->VectorStrictIn2 (p, t, -s, 1e-6*geomsize)  << endl;
+                            }
+                          
+                          
 			  /*
 			  bool isface =
 			    (locsol->VectorIn (p, t2a) &&
@@ -2095,10 +2122,8 @@ namespace netgen
 			     !locsol->VectorStrictIn (p, t2b));
 			  */
 
-			  if (isface)
-			    {
-			      cnts++;
-			    }
+                          if (isface)
+                            cnts++;
 			}
 		      if (cnts < 2) isedge = 0;
 		    }
diff --git a/libsrc/csg/surface.cpp b/libsrc/csg/surface.cpp
index 8efe1f7a..cb6ab48b 100644
--- a/libsrc/csg/surface.cpp
+++ b/libsrc/csg/surface.cpp
@@ -427,11 +427,20 @@ VecInSolid2 (const Point<3> & p,
   if (hv1 >= eps)
     return IS_OUTSIDE;
 
+  double hv2 = v2 * hv;
+  if (hv2 <= -eps)
+    return IS_INSIDE;
+  if (hv2 >= eps)
+    return IS_OUTSIDE;
+  return DOES_INTERSECT;
+  
+  /*
   double hv2 = v2 * hv;
   if (hv2 <= 0)
     return IS_INSIDE;
   else
     return IS_OUTSIDE;
+  */
 }