point out of quad fix

libsrc/meshing/meshclass.cpp

@@ -4299,7 +4299,6 @@ namespace netgen
     //SZ 
     if(SurfaceElement(element).GetType()==QUAD)
       {
-        //cout << "check quad element" << endl;
         const Element2d & el = SurfaceElement(element); 
 
         const Point3d & p1 = Point(el.PNum(1)); 
@@ -4314,6 +4313,7 @@ namespace netgen
         Vec3d c = p4 - a; 
         Vec3d d = p3 - a - b - c;
 
+     
         double dxb = d.X()*b.Y()-d.Y()*b.X();
         double dxc = d.X()*c.Y()-d.Y()*c.X(); 
         double dxa = d.X()*a.Y()-d.Y()*a.X(); 
@@ -4329,16 +4329,19 @@ namespace netgen
 	double d1 = dp13.Length2();
 	double d2 = dp24.Length2();
 
-        // if(fabs(d.X()) <= eps && fabs(d.Y())<= eps)
+	// if(fabs(d.X()) <= eps && fabs(d.Y())<= eps)
 	//if (d.Length2() < sqr(eps))
         if (d.Length2() < sqr(eps)*d1 && d.Length2() < sqr(eps)*d2)
           {
 	    //Solve Linear System
-
+	    Vec2d sol;
-            Vec2d sol;
             SolveLinearSystemLS (b, c, p-a, sol);
             lami[0] = sol.X();
             lami[1] = sol.Y();
+
+	    if(lami[1]<=1.+eps && lami[1]>=0.-eps && lami[0]<=1.+eps && lami[0]>=0.-eps)
+	      return true;
+	    
             /*
             lami[0]=(c.Y()*(p.X()-a.X())-c.X()*(p.Y()-a.Y()))/
               (b.X()*c.Y() -b.Y()*c.X()); 
@@ -4347,66 +4350,103 @@ namespace netgen
             */
           } 
         else
-          //if(fabs(dxb) <= eps)
+          if(fabs(dxb) <= eps*fabs(dxc))
-	  if(fabs(dxb) <= eps*fabs(dxc))
             {
 	      lami[1] = (dxp-dxa)/dxc;
-              //if(fabs(b.X()-d.X()*lami[1])>=eps)
+              if(fabs(b.X()+d.X()*lami[1])>=fabs(b.Y()+d.Y()*lami[1]))
-	      if(fabs(b.X()+d.X()*lami[1])>=fabs(b.Y()+d.Y()*lami[1]))
                 lami[0] = (p.X()-a.X() - c.X()*lami[1])/(b.X()+d.X()*lami[1]); 
               else
                 lami[0] = (p.Y()-a.Y() - c.Y()*lami[1])/(b.Y()+d.Y()*lami[1]);
 
+	      if(lami[1]<=1.+eps && lami[1]>=0.-eps && lami[0]<=1.+eps && lami[0]>=0.-eps)
+		return true;
             }
           else
-            //if(fabs(dxc) <= eps)
+            if(fabs(dxc) <= eps*fabs(dxb))
-	    if(fabs(dxc) <= eps*fabs(dxb))
               {
 		lami[0] = (dxp-dxa)/dxb;
-                //if(fabs(c.X()-d.X()*lami[0])>=eps)
+                if(fabs(c.X()+d.X()*lami[0])>=fabs(c.Y()+d.Y()*lami[0]))
-		if(fabs(c.X()+d.X()*lami[0])>=fabs(c.Y()+d.Y()*lami[0]))
                   lami[1] = (p.X()-a.X() - b.X()*lami[0])/(c.X()+d.X()*lami[0]); 
                 else
                   lami[1] = (p.Y()-a.Y() - b.Y()*lami[0])/(c.Y()+d.Y()*lami[0]);
 
+		if(lami[1]<=1.+eps && lami[1]>=0.-eps && lami[0]<=1.+eps && lami[0]>=0.-eps)
+		  return true;
               }
             else //Solve quadratic equation
               {
-		//if(fabs(d.X()) >= eps)
+		c2 = -d.X()*dxb;
-		if(fabs(d.X()) >= fabs(d.Y()))
+		c1 = b.X()*dxc - c.X()*dxb + d.X()*(dxp-dxa);
-                  {
+		c0 = c.X()*(dxp-dxa) + (a.X()-p.X())*dxc;
-                    c2 = d.X()*dxc;
+		double rt =  c1*c1 - 4*c2*c0;
-                    //c1 = d.X()*dxc - c.X()*dxb - d.X()*(dxp-dxa);
-		    c1 = b.X()*dxc - c.X()*dxb - d.X()*(dxp-dxa);
-                    c0 = -b.X()*(dxp -dxa) - (a.X()-p.X())*dxb;
-		  }
-                else 
-                  {
-                    c2 = d.Y()*dxc;
-                    //c1 = d.Y()*dxc - c.Y()*dxb - d.Y()*(dxp-dxa);
-		    c1 = b.Y()*dxc - c.Y()*dxb - d.Y()*(dxp-dxa);
-                    c0 = -b.Y()*(dxp -dxa) - (a.Y()-p.Y())*dxb;
-                  }
 		
-                double rt =  c1*c1 - 4*c2*c0;
+		if (rt < 0.) return false; 
+		lami[1] = (-c1 + sqrt(rt))/2/c2;
 
-                if (rt < 0.) return false; 
-                lami[1] = (-c1 + sqrt(rt))/2/c2;
 
 		if(lami[1]<=1.+eps && lami[1]>=0.-eps)
-                  {
+		  {
-                    lami[0] = (dxp - dxa -dxc*lami[1])/dxb;
+		    lami[0] = (dxp - dxa -dxb*lami[1])/dxc;
+		    
+		    if(lami[0]<=1.+eps && lami[0]>=0.-eps)
+		      return true;
+		  }
+		lami[1] = (-c1 - sqrt(rt))/2/c2;
+
+		lami[0] = (dxp - dxa -dxb*lami[1])/dxc;
+
+		if(lami[1]<=1.+eps && lami[1]>=0.-eps && lami[0]<=1.+eps && lami[0]>=0.-eps)
+		  return true;
+
+		c2 = d.Y()*dxb;
+		c1 = b.Y()*dxc - c.Y()*dxb + d.Y()*(dxp-dxa);
+		c0 = c.Y()*(dxp -dxa) + (a.Y()-p.Y())*dxc;
+		rt =  c1*c1 - 4*c2*c0;
+		
+		if (rt < 0.) return false; 
+		lami[1] = (-c1 + sqrt(rt))/2/c2;
+
+		if(lami[1]<=1.+eps && lami[1]>=0.-eps)
+		  {
+		    lami[0] = (dxp - dxa -dxb*lami[1])/dxc;
+
+		    if(lami[0]<=1.+eps && lami[0]>=0.-eps)
+		      return true;
+		  }
+		lami[1] = (-c1 - sqrt(rt))/2/c2;
+
+		lami[0] = (dxp - dxa -dxb*lami[1])/dxc;
+
+		if(lami[1]<=1.+eps && lami[1]>=0.-eps && lami[0]<=1.+eps && lami[0]>=0.-eps)
+		  return true;
+
+		c2 = -d.X()*dxc;
+		c1 = -b.X()*dxc + c.X()*dxb + d.X()*(dxp-dxa);
+		c0 = b.X()*(dxp -dxa) + (a.X()-p.X())*dxb;
+		rt =  c1*c1 - 4*c2*c0;
+		
+		if (rt < 0.) return false; 
+		lami[1] = (-c1 + sqrt(rt))/2/c2;
+
+		if(lami[1]<=1.+eps && lami[1]>=0.-eps)
+		  {
+		    lami[0] = (dxp - dxa -dxc*lami[1])/dxb;
+
+		    if(lami[0]<=1.+eps && lami[0]>=0.-eps)
+		      return true;
+		  }
+		lami[1] = (-c1 - sqrt(rt))/2/c2;
+
+		lami[0] = (dxp - dxa -dxc*lami[1])/dxb;
+
+		if(lami[1]<=1.+eps && lami[1]>=0.-eps && lami[0]<=1.+eps && lami[0]>=0.-eps)
+		  return true;
+	      }
 
-                    if(lami[0]<=1.+eps && lami[0]>=0.-eps)
-                      return true;
-                  }
       
-                lami[1] = (-c1 - sqrt(rt))/2/c2;
-                lami[0] = (dxp - dxa -dxc*lami[1])/dxb;
-              }
         //cout << "lam0,1 = " << lami[0] << ", " << lami[1] << endl;
         
-        if( lami[0] <= 1.+eps  && lami[0] >= -eps && lami[1]<=1.+eps && lami[1]>=-eps)
+        /*if( lami[0] <= 1.+eps  && lami[0] >= -eps && lami[1]<=1.+eps && lami[1]>=-eps)
           {
             if(consider3D)
               {
@@ -4419,7 +4459,7 @@ namespace netgen
               }
             else
               return true;
-          }
+	      }*/
 
         return false;