implementing trilinear interpolation

program/source/interpol/interpol.cc

@@ -1551,6 +1551,7 @@ void Interpolate_direct::init()
                 cEST = blockgrid->Give_coord_hexahedron(id_hex,i+1,j  ,k+1);
                 cWNT = blockgrid->Give_coord_hexahedron(id_hex,i,  j+1,k+1);
                 cENT = blockgrid->Give_coord_hexahedron(id_hex,i+1,j+1,k+1);
+                //isPlane()
 
                     // bounding box calculation
                 boxWSD.x = MIN(MIN(MIN(cWSD.x,cESD.x),MIN(cWND.x,cEND.x)),
@@ -2891,11 +2892,11 @@ int Interpolate_direct::checkBox(int id_Hex, int i, int j, int k, D3vector v)
         return -1;
     }
     checkCounter++;
-    D3vector cWSD, cESD;
-    D3vector cWND, cEND;
+//    D3vector cWSD, cESD;
+//    D3vector cWND, cEND;
 
-    D3vector cWST, cEST;
-    D3vector cWNT, cENT;
+//    D3vector cWST, cEST;
+//    D3vector cWNT, cENT;
 
     D3vector boxWSD, boxENT;
 
@@ -2934,10 +2935,43 @@ int Interpolate_direct::checkBox(int id_Hex, int i, int j, int k, D3vector v)
         cWNT = blockgrid->Give_coord_hexahedron(id_Hex,i,  j+1,k+1);
         cENT = blockgrid->Give_coord_hexahedron(id_Hex,i+1,j+1,k+1);
 
+//        std::cout << "delete line below " << std::endl;
+//        k = 1;
+//        i = 3;
+//        j = 3;
+//        cWSD = blockgrid->Give_coord_hexahedron(0,i,  j,  k  );
+//        cESD = blockgrid->Give_coord_hexahedron(0,i+1,j  ,k  );
+//        cWND = blockgrid->Give_coord_hexahedron(0,i,  j+1,k  );
+//        cEND = blockgrid->Give_coord_hexahedron(0,i+1,j+1,k  );
+
+//        cWST = blockgrid->Give_coord_hexahedron(0,i,  j,  k+1);
+//        cEST = blockgrid->Give_coord_hexahedron(0,i+1,j  ,k+1);
+//        cWNT = blockgrid->Give_coord_hexahedron(0,i,  j+1,k+1);
+//        cENT = blockgrid->Give_coord_hexahedron(0,i+1,j+1,k+1);
+
+
+        cWSD.Print();std::cout<<std::endl;
+        cESD.Print();std::cout<<std::endl;
+        cWND.Print();std::cout<<std::endl;
+        cEND.Print();std::cout<<std::endl;
+        cWST.Print();std::cout<<std::endl;
+        cEST.Print();std::cout<<std::endl;
+        cWNT.Print();std::cout<<std::endl;
+        cENT.Print();std::cout<<std::endl;
+
+
+
+
+
+
         D3vector lam, lamTemp ;
 
 
 
+        D3vector lamTrilinar = trilinarInterpolation(v,id_Hex,i,  j,  k  );
+
+
+
         lam = lambda_of_p_in_tet(v,cWND,cWNT,cWST,cEST);
         lamTemp = lam;
         if (contained_in_tet_strong(lam) && new_lam_better(lambda,lam)){typ=0;typCounter0++;}
@@ -3365,6 +3399,132 @@ void Interpolate_direct::writeInterpolationEfficiency()
 
 }
 
+D3vector Interpolate_direct::trilinarInterpolation(D3vector X, int id_Hex, int i, int j, int k)
+{
+    std::cout << "to be implemented" << std::endl;
+    D3vector A = cEND;        // 000
+    D3vector B = cESD - cEND; // 100 - 000
+    D3vector C = cWND - cEND; // 010 - 000
+    D3vector D = cENT - cEND; // 001 - 000
+    D3vector E = cWSD - cESD - cWND + cEND;        // 110 - 100 - 010 + 000
+    D3vector F = cWNT - cWND - cENT + cEND;        // 011 - 010 - 001 + 000
+    D3vector G = cEST - cESD - cENT + cEND;        // 101 - 100 - 001 + 000
+    D3vector H = cWST + cESD + cWND + cENT - cEND - cWSD - cWNT - cEST;        // 111 + 100 + 010 + 001 - 000 - 110 - 011 - 101
+//    D3vector cEST = {2,0,2};      // 1,0,1 : x6
+//    D3vector cWNT = {0,1,1};      // 0,1,1 : x4
+//    D3vector cENT = {-1,0,2};      // 0,0,1 : x7
+    D3vector normal1 = cross_product(B,C);
+    D3vector normal1opp = cross_product(cENT-cEST,cENT - cWNT);
+
+    D3vector J = cross_product(E / D3VectorNorm(E),H / D3VectorNorm(H));
+    D3vector K = cross_product(F / D3VectorNorm(F),H / D3VectorNorm(H));
+    D3vector L = cross_product(G / D3VectorNorm(G),H / D3VectorNorm(H));
+
+    if (J != 0)
+    {
+        J = J / D3VectorNorm(J);
+    }
+    if (K != 0)
+    {
+        K = K / D3VectorNorm(K);
+    }
+    if (L != 0)
+    {
+        L = L / D3VectorNorm(L);
+    }
+
+
+    if (J == 0 || K == 0 || L == 0 )
+    {
+        std::cout << "seems to be orthogonal!" << std::endl;
+    }
+
+    double eta = 0.5;
+    double xi  = 0.5;
+    double phi = 0.5;
+    D3vector coord(eta,xi,phi);
+
+    coord.Print();std::cout<<std::endl;
+    D3vector R = A    + B * coord.x + C * coord.y + D * coord.z + E * coord.x*coord.y + F * coord.y*coord.z + G * coord.x * coord.z + H * coord.x*coord.y*coord.z;
+    R.Print();std::cout<<std::endl;
+
+//    std::cout << " product(A,J) " << product(A,J) << std::endl;
+//    std::cout << " product(B,J) " << product(B,J) << std::endl;
+//    std::cout << " product(C,J) " << product(C,J) << std::endl;
+//    std::cout << " product(D,J) " << product(D,J) << std::endl;
+//    std::cout << " product(E,J) " << product(E,J) << std::endl;
+//    std::cout << " product(F,J) " << product(F,J) << std::endl;
+//    std::cout << " product(G,J) " << product(G,J) << std::endl;
+//    std::cout << " product(H,J) " << product(H,J) << std::endl;
+
+//    std::cout << " product(A,K) " << product(A,K) << std::endl;
+//    std::cout << " product(B,K) " << product(B,K) << std::endl;
+//    std::cout << " product(C,K) " << product(C,K) << std::endl;
+//    std::cout << " product(D,K) " << product(D,K) << std::endl;
+//    std::cout << " product(E,K) " << product(E,K) << std::endl;
+//    std::cout << " product(F,K) " << product(F,K) << std::endl;
+//    std::cout << " product(G,K) " << product(G,K) << std::endl;
+//    std::cout << " product(H,K) " << product(H,K) << std::endl;
+
+//    std::cout << " product(A,L) " << product(A,L) << std::endl;
+//    std::cout << " product(B,L) " << product(B,L) << std::endl;
+//    std::cout << " product(C,L) " << product(C,L) << std::endl;
+//    std::cout << " product(D,L) " << product(D,L) << std::endl;
+//    std::cout << " product(E,L) " << product(E,L) << std::endl;
+//    std::cout << " product(F,L) " << product(F,L) << std::endl;
+//    std::cout << " product(G,L) " << product(G,L) << std::endl;
+//    std::cout << " product(H,L) " << product(H,L) << std::endl;
+
+    double f = product(B,J) * coord.x + product(C,J) * coord.y + product(D,J) * coord.z                                    + product(F,J) * coord.y * coord.z + product(G,J) * coord.x * coord.z + product(A,J) - product(X,J);
+    double g = product(B,K) * coord.x + product(C,K) * coord.y + product(D,K) * coord.z + product(E,K) * coord.x * coord.y                                    + product(G,K) * coord.x * coord.z + product(A,K) - product(X,K);
+    double h = product(B,L) * coord.x + product(C,L) * coord.y + product(D,L) * coord.z + product(E,L) * coord.x * coord.y + product(F,L) * coord.y * coord.z                                    + product(A,L) - product(X,L);
+
+    for (int iter = 0 ; iter < 10 ; iter++)
+    {
+        double dxdx = product(B,J) + product(G,J) *  coord.z ;
+        double dydx = product(B,K) +  product(E,K)  * coord.y + product(G,K) * coord.z ;
+        double dzdx = product(B,L) + product(E,L) * coord.y ;
+        double dxdy = product(C,J) +  product(F,J) * coord.z ;
+        double dydy = product(C,K) + product(E,K) * coord.x;
+        double dzdy =  product(C,L) +  product(E,L) * coord.x  + product(F,L) * coord.z;
+        double dxdz =  product(D,J) + product(F,J) * coord.y + product(G,J) * coord.x ;
+        double dydz = product(D,K) +  product(G,K) * coord.x ;
+        double dzdz = product(D,L)  + product(F,L) * coord.y ;
+
+        D3vector cx(dxdx,dydx,dzdx); // first column
+        D3vector cy(dxdy,dydy,dzdy); // second column
+        D3vector cz(dxdz,dydz,dzdz); // third column
+
+
+
+        D3matrix jac(cx,cy,cz);
+        //jac.Print();
+        D3matrix Temp(cx,cy,cz);
+        jac.invert_gauss_elimination();
+        //jac.matrixMultiply(Temp).Print();
+
+        double f = product(B,J) * coord.x + product(C,J) * coord.y + product(D,J) * coord.z + product(F,J) * coord.y * coord.z + product(G,J) * coord.x * coord.z + product(A,J) - product(X,J);
+        double g = product(B,K) * coord.x + product(C,K) * coord.y + product(D,K) * coord.z + product(E,K) * coord.x * coord.y + product(G,K) * coord.x * coord.z + product(A,K) - product(X,K);
+        double h = product(B,L) * coord.x + product(C,L) * coord.y + product(D,L) * coord.z + product(E,L) * coord.x * coord.y + product(F,L) * coord.y * coord.z + product(A,L) - product(X,L);
+        D3vector f_xn = {f,g,h};
+
+        std::cout << "vector F "; f_xn.Print();std::cout<<std::endl;
+        if (D3VectorNormSquared(f_xn) < 1e-15)
+        {
+            iter = 1000;
+        }
+        coord = coord - jac.vectorMultiply(f_xn);
+        coord.Print();std::cout<<std::endl;
+
+    }
+    D3vector x = A    + B * coord.x + C * coord.y + D * coord.z + E * coord.x*coord.y + F * coord.y*coord.z + G * coord.x * coord.z + H * coord.x*coord.y*coord.z;
+    X.Print();
+    x.Print();
+    std::cout << std::endl;
+
+
+}
+
 
 int Interpolate_direct::checkForHexaNeighbours(int idHex, int i, int j, int k, D3vector v)
 {

program/source/interpol/interpol.h

@@ -389,6 +389,7 @@ class Interpolate_direct {
     D3vector vNow{1e10,1e10,1e10}, vPrev, vPrevPrev, vPrevPrevPrev;
 
 
+    D3vector trilinarInterpolation(D3vector v, int id_Hex,int i,  int j,  int k );
         bool vectorInBox(D3vector vWSD, D3vector vENT, D3vector v, double eps = 1e-10);
         int checkBox(int idHex, int i, int j, int k, D3vector v);
     private:
@@ -401,7 +402,11 @@ class Interpolate_direct {
       double lamUpperLimit{1.1};
       int typ_tet;
 
+      D3vector cWSD, cESD;
+      D3vector cWND, cEND;
 
+      D3vector cWST, cEST;
+      D3vector cWNT, cENT;
       Blockgrid    *blockgrid;
       std::vector<std::vector<std::vector<D3vector> > > arrayBoxWSDENT;
 

program/source/math_lib/math_lib.h

@@ -63,7 +63,7 @@ class D3vector {
   //bool operator==(const D3vector& v) {if(fabs(v.x - x) <1e-10 && fabs(v.y - y) <1e-10  && fabs(v.z - z) <1e-10 ){return true;} else {return false;} }
   bool operator==(const D3vector& v) {if( (v.x == x)  && (v.y == y) && (v.z == z) ){return true;} else {return false;} }
   bool operator==(const double v) {if( (v == x)  && (v == y) && (v == z) ){return true;} else {return false;} }
-  bool operator!=(const double v) {if( (v != x)  && (v != y) && (v != z) ){return true;} else {return false;} }
+  bool operator!=(const double v) {if( (v != x)  || (v != y) || (v != z) ){return true;} else {return false;} }
   bool operator<(const D3vector& v) {if(v.x > x && v.y > y && v.z > z){return true;} else {return false;} }
   bool operator>(const D3vector& v) {if(v.x < x && v.y < y && v.z < z){return true;} else {return false;} }
 //  bool operator<=(const D3vector& v) {if(v.x >= x && v.y >= y && v.z >= z){return true;} else {return false;} }
@@ -203,9 +203,9 @@ class D3matrix {
   }
 
   inline void Print() {
-    std::cout << "Matrix: " << std::endl;
-    std::cout << x1 << ", " << y1 << ", " << z1 << ";" <<std::endl;
-    std::cout << x2 << ", " << y2 << ", " << z2 << ";" <<std::endl;
+    std::cout << "Matrix: " << ";\n";
+    std::cout << x1 << ", " << y1 << ", " << z1 << ";\n" ;
+    std::cout << x2 << ", " << y2 << ", " << z2 << ";\n" ;
     std::cout << x3 << ", " << y3 << ", " << z3 << ";" <<std::endl;
   }
 
@@ -400,8 +400,8 @@ inline double product(const D3vector& v, const D3vector& w) {
 
 inline D3vector cross_product(const D3vector& v, const D3vector& w) {
   return D3vector(v.y*w.z-v.z*w.y,
-		  v.z*w.x-v.x*w.z,
-		  v.x*w.y-v.y*w.x);
+                  v.z*w.x-v.x*w.z,
+                  v.x*w.y-v.y*w.x);
 }
 
 inline D3vector normal_vector_of_triangle(const D3vector& va,
@@ -475,7 +475,7 @@ double calc_maximal_face_angle(const D3vector& va,
 // ----------
 
 inline void D3vector::Print() {
-  std::cout << "Coordinate: " << x << ", " << y << ", " << z << ";";
+  std::cout << "Coordinate: " << x << ", " << y << ", " << z << ";\n";
 }
 
 inline void D3vector::Print(std::ofstream *Datei) {