bug fix

.gitignore

@@ -5,4 +5,5 @@ cgrid3d.cpp
 cmesh3d.cpp
 build/
 debug_tt_mex/
+test_lelievre11.py
 

ttcr/Grid3Dui.h

@@ -621,7 +621,11 @@ namespace ttcr {
                 std::swap(iA, iC);
                 std::swap(vertexA, vertexC);
             }
-
+            if ( vertexB->getTT(threadNo) > vertexC->getTT(threadNo) ) {
+                std::swap(iB, iC);
+                std::swap(vertexB, vertexC);
+            }
+
             if ( vertexA->getTT(threadNo) == std::numeric_limits<T1>::max() ) {
                 continue;
             }
@@ -641,6 +645,7 @@ namespace ttcr {
                     vertexB->getY() - vertexA->getY(),
                     vertexB->getZ() - vertexA->getZ() };
 
+                // vector normal to plane
                 sxyz<T1> v_n = cross(v_b, v_c);
 
                 T1 b = norm( v_b );
@@ -651,44 +656,30 @@ namespace ttcr {
 
                 T1 phi = c*b*sin(alpha);
 
+                // TODO check for negative value
                 T1 w_tilde = sqrt( vertexD->getNodeSlowness()*vertexD->getNodeSlowness()*phi*phi -
                                   u*u*b*b - v*v*c*c + 2.*u*v*d2 );
 
+                // Point (ξ_0 , ζ_0 ) is the normalized projection of node D onto face ABC
                 // project D on plane
 
                 T1 d_tmp = -vertexA->getX()*v_n.x - vertexA->getY()*v_n.y - vertexA->getZ()*v_n.z;
 
                 T1 k = -(d_tmp + v_n.x*vertexD->getX() + v_n.y*vertexD->getY() + v_n.z*vertexD->getZ())/
                 norm2(v_n);
 
-                sxyz<T1> pt;
+                sxyz<T1> pt;   // -> Point (ξ_0 , ζ_0 ) 
                 pt.x = vertexD->getX() + k*v_n.x;
                 pt.y = vertexD->getY() + k*v_n.y;
                 pt.z = vertexD->getZ() + k*v_n.z;
 
                 T1 rho0 = vertexD->getDistance( pt );
 
                 sxyz<T1> v_pt = {pt.x-vertexA->getX(), pt.y-vertexA->getY(), pt.z-vertexA->getZ()};
-
-                // project point on AB
-
-//                k = dot(v_pt,v_c)/dot(v_c,v_c);
-//                pt.x = vertexA->getX() + k*v_c.x;
-//                pt.y = vertexA->getY() + k*v_c.y;
-//                pt.z = vertexA->getZ() + k*v_c.z;
-//                
-//                T1 xi0 = vertexA->getDistance( pt )/c;
-                T1 xi0 = dot(v_pt,v_c)/dot(v_c,v_c);
-
-                // project point on AC
-
-//                k = dot(v_pt,v_b)/dot(v_b,v_b);
-//                pt.x = vertexA->getX() + k*v_b.x;
-//                pt.y = vertexA->getY() + k*v_b.y;
-//                pt.z = vertexA->getZ() + k*v_b.z;
-//                
-//                T1 zeta0 = vertexA->getDistance( pt )/b;
-                T1 zeta0 = dot(v_pt,v_b)/dot(v_b,v_b);
+
+                T1 xi0;
+                T1 zeta0;
+                projNorm(v_b/b, v_c/c, v_pt, xi0, zeta0);
 
                 T1 beta = u*b*b - v*d2;
                 T1 gamma = v*c*c - u*d2;

ttcr/ttcr_t.h

@@ -562,7 +562,6 @@ namespace ttcr {
         return v3;
     }
 
-
     template<typename T>
     T norm(const sxz<T>& v) {
         return sqrt( v.x*v.x + v.z*v.z );
@@ -583,6 +582,47 @@ namespace ttcr {
         return v.x*v.x + v.y*v.y + v.z*v.z;
     }
 
+    template<typename T>
+    void projNorm(const sxyz<T>& b, const sxyz<T>& c, const sxyz<T>& p, T& xi0, T& zeta0) {
+        //                B
+        //                o
+        //              /   \
+        //           /          \
+        //        /       o P       \
+        //    A o----------------------o  C
+        //
+        // b is _unit_ vector AC
+        // c is _unit_ vector AB
+        // p is vector AP
+        // xi0 & zeta0 are normalized coordinates of P in plane ABC (Lelièvre et al, 2011, GJI 184, 885-896)
+
+        // solved using xi*c + zeta*b = p
+
+        if ( c.x != 0.0 && c.x*b.y != c.y*b.x) {
+            zeta0 = (p.y - c.y/c.x*p.x) / (b.y - c.y/c.x*b.x);
+            xi0 = (p.x - zeta0*b.x)/c.x;
+        } else if ( c.x != 0.0 && c.x*b.z != c.z*b.x) {
+            zeta0 = (p.z - c.z/c.x*p.x) / (b.z - c.z/c.x*b.x);
+            xi0 = (p.x - zeta0*b.x)/c.x;
+        } else if ( c.y != 0.0 && c.y*b.x != c.x*b.y ) {
+            zeta0 = (p.x - c.x/c.y*p.y) / (p.x - c.x/c.y*p.y);
+            xi0 = (p.y - zeta0*b.y)/c.y;
+        } else if ( c.y != 0.0 && c.y*b.z != c.z*b.y ) {
+            zeta0 = (p.z - c.z/c.y*p.y) / (p.z - c.z/c.y*p.y);
+            xi0 = (p.y - zeta0*b.y)/c.y;
+        } else if ( c.z != 0.0 && c.z*b.x != c.x*b.z ) {
+            zeta0 = (p.x - c.x/c.z*p.z) / (b.x - c.x/c.z*b.z);
+            xi0 = (p.z - zeta0*b.z)/c.z;
+        } else if ( c.z != 0.0 && c.z*b.y != c.y*b.z ) {
+            zeta0 = (p.y - c.y/c.z*p.z) / (b.y - c.y/c.z*b.z);
+            xi0 = (p.z - zeta0*b.z)/c.z;
+        } else {
+            // error -> return negative values for both xi & zeta
+            xi0 = -1.0;
+            zeta0 = -1.0;
+        }
+    }
+
 #ifndef _MSC_VER
     // following 3 fct from
     // http://stackoverflow.com/questions/1903954/is-there-a-standard-sign-function-signum-sgn-in-c-c