Merge branch 'master' of i10git.cs.fau.de:er96apow/UGBlocks_V3

program/source/extemp/print_var_vtk_cc.h

@@ -221,6 +221,13 @@ void Variable<DTyp>::Print_VTK(std::ostream& Datei, double (*convert)(DTyp x),do
 template <class DTyp>
 void Variable<DTyp>::QPrint_VTK(QString DateiName, double (*convert)(DTyp x),double stretch_z, QString title,  
 				Unstructured_grid_Marker * marker, bool cutSmallData) {
+
+    if (this==NULL)
+    {
+        std::cout << "could not write file, var == NULL" << std::endl;
+        return;
+    }
+
   D3vector v;
   int Nx, Ny, Nz;
   int num_total = 0;

program/source/extemp/variable.h

@@ -296,6 +296,9 @@ class Variable : public Expr<Variable< DTyp > >, public Object_based_on_ug {
 		
 		template <elementTyp TYP_EL>
 		inline DTyp Give_data ( params_in ) const;
+
+        template <elementTyp TYP_EL>
+        inline void Set_data ( params_in, DTyp value );
 		
 		
 		// fuer debug-zwecke mit positionsrueckgabe
@@ -303,7 +306,7 @@ class Variable : public Expr<Variable< DTyp > >, public Object_based_on_ug {
 		inline DTyp Give_data ( params_in, D3vector &posVector ) const ;
 		
 		void UpdateHexahedra(); ///> Updated alle Randdaten in allen Hexahedra
-		
+        void UpdateHexahedraBack();
 		/// todo update nur falls notwendig
 		template <elementTyp TYP_EL>
 		void Update ( int id ) const;
@@ -420,6 +423,13 @@ void Variable<DTyp>::UpdateHexahedra() {
      }
 }
 
+template <class DTyp>
+void Variable<DTyp>::UpdateHexahedraBack() {
+     for(int id_hex=0;id_hex<ug->Give_number_hexahedra();++id_hex) {
+         Update_back(id_hex);
+     }
+}
+
 template <class DTyp>
 int Variable<DTyp>::checkThreshold ( double threshold, int kz ) {
 			int id_hex = 1;
@@ -501,6 +511,12 @@ inline double Variable<double>::Give_data<hexahedronEl> ( params_in ) const {
 	return data_hexahedra[id][i+ ( Nx+1 ) * ( j+ ( Ny+1 ) *k ) ];
 }
 
+template <>
+template <>
+inline void Variable<double>::Set_data<hexahedronEl> ( params_in, double value )  {
+    data_hexahedra[id][i+ ( Nx+1 ) * ( j+ ( Ny+1 ) *k ) ] = value;
+}
+
 // fuer debug zwecke koordinatenrueckgabe
 template <>
 template <>

program/source/extemp/variable_cc.h

@@ -989,7 +989,7 @@ assert(	rank_rec == 0);
 			  //  data_hexahedra[id_hex][i*Ni+j*Nj+Nconst] =
 			  //			   data_quadrangles[id_quad][i+ ( Nxquad+1 ) *j];
 			  data_quadrangles[id_quad][i+ ( Nxquad+1 ) *j] =
-			       data_hexahedra[id_hex][i*Ni+j*Nj+Nconst];
+                   data_hexahedra[id_hex][i*Ni+j*Nj+Nconst];
 		      }
 	       }
 	       else {
@@ -1034,7 +1034,7 @@ assert(	rank_rec == 0);
 			   for(i=1; i<Nxedge; ++i ) {
 			       // vorwaerts
 			       //data_hexahedra[id_hex][i*Ni+Nconst] = data_edges[id_edge][i];
-			       data_edges[id_edge][i] = data_hexahedra[id_hex][i*Ni+Nconst];
+                   data_edges[id_edge][i] = data_hexahedra[id_hex][i*Ni+Nconst];
 			   }
 			} 
 			else {
@@ -1064,7 +1064,9 @@ assert(	rank_rec == 0);
 			if(rank_send == rank_rec )      // for parallel
                            // vorwaerts
 			   //data_hexahedra[id_hex][Nconst] = data_points[id_point][0];
-			   data_points[id_point][0] = data_hexahedra[id_hex][Nconst];
+              {
+                data_points[id_point][0] = data_hexahedra[id_hex][Nconst];
+            }
 			else {
 			  assert(false);
 			}

program/source/grid/blockgrid.cc

@@ -101,6 +101,7 @@ Blockgrid::Blockgrid ( Unstructured_grid *ug_, int Nx, int Ny, int Nz ) {
 
 
 Blockgrid::Blockgrid ( Unstructured_grid *ug_, int Nx, int Ny, int Nz, int Nr ) {
+    cout << "Nx " << Nx <<  "Ny " << Ny <<  "Nz " << Nz <<  "Nr " << Nr << endl;
   id_of_grid = id_count_grid; ++id_count_grid;
     phase_shift = NULL;
   bg_coord = NULL;

program/source/interpol/interpol.cc

@@ -1607,10 +1607,6 @@ void Interpolate_direct::init()
             j = 0;
             for (i = 0 ; i < Nx+1; i++) for (k = 0 ; k < Nz+1; k++)
             {
-                if (i == 0 && j == 0 && k == 4)
-                {
-                    cout << "debuigf " << endl;
-                }
                 if (!indexAllFacesOutside.at(3).empty())
                 {calculateNeighbourIndex(correctRelation.at(3), indexAllFacesOutside.at(3).at(0).at(0), id_hex, i,j,k,Nx,Ny,Nz);}
                 correctRelation.at(3);
@@ -3182,16 +3178,17 @@ int Interpolate_direct::checkEdge(int idHex, int i, int j, int k, D3vector v)
 
 bool Interpolate_direct::checkOverlapOfBoxes(D3vector vWSD, D3vector vENT, D3vector wWSD, D3vector wENT)
 {
-    if ( (( (vWSD.x <= wENT.x) && (vWSD.x >= wWSD.x) ) || ( (vENT.x <= wENT.x) && (vENT.x >= wWSD.x) ) )
-      && (( (vWSD.y <= wENT.y) && (vWSD.y >= wWSD.y) ) || ( (vENT.y <= wENT.y) && (vENT.y >= wWSD.y) ) )
-      && (( (vWSD.z <= wENT.z) && (vWSD.z >= wWSD.z) ) || ( (vENT.z <= wENT.z) && (vENT.z >= wWSD.z) ) ) )
+    double epsilon = 1e-10;
+    if ( (( (vWSD.x <= wENT.x + epsilon) && (vWSD.x + epsilon >= wWSD.x) ) || ( (vENT.x <= wENT.x + epsilon) && (vENT.x + epsilon >= wWSD.x) ) )
+      && (( (vWSD.y <= wENT.y + epsilon) && (vWSD.y + epsilon >= wWSD.y) ) || ( (vENT.y <= wENT.y + epsilon) && (vENT.y + epsilon >= wWSD.y) ) )
+      && (( (vWSD.z <= wENT.z + epsilon) && (vWSD.z + epsilon >= wWSD.z) ) || ( (vENT.z <= wENT.z + epsilon) && (vENT.z + epsilon >= wWSD.z) ) ) )
     {
         return true;
     }
 
-    if ( (( (wWSD.x <= vENT.x) && (wWSD.x >= vWSD.x) ) || ( (wENT.x <= vENT.x) && (wENT.x >= vWSD.x) ) )
-      && (( (wWSD.y <= vENT.y) && (wWSD.y >= vWSD.y) ) || ( (wENT.y <= vENT.y) && (wENT.y >= vWSD.y) ) )
-      && (( (wWSD.z <= vENT.z) && (wWSD.z >= vWSD.z) ) || ( (wENT.z <= vENT.z) && (wENT.z >= vWSD.z) ) ) )
+    if ( (( (wWSD.x <= vENT.x + epsilon) && (wWSD.x + epsilon >= vWSD.x) ) || ( (wENT.x <= vENT.x + epsilon) && (wENT.x + epsilon >= vWSD.x) ) )
+      && (( (wWSD.y <= vENT.y + epsilon) && (wWSD.y + epsilon >= vWSD.y) ) || ( (wENT.y <= vENT.y + epsilon) && (wENT.y + epsilon >= vWSD.y) ) )
+      && (( (wWSD.z <= vENT.z + epsilon) && (wWSD.z + epsilon >= vWSD.z) ) || ( (wENT.z <= vENT.z + epsilon) && (wENT.z + epsilon >= vWSD.z) ) ) )
     {
         return true;
     }

program/source/interpol/interpol.h

@@ -28,9 +28,10 @@
 /////////////////////////////////////////////////////////////
 // 1. Interpolate from  blockgrid to rectangular blockgrid
 // 2. Interpolate from  blockgrid  to  blockgrid
+// 3. Interpolate from Variable on a blockgrid to any point using structured intermediate grid
+// 4. Interpolate from blockgrid direct to any point
 /////////////////////////////////////////////////////////////
 
-
 /////////////////////////////////////////////////////////////
 // 1. Interpolate from  blockgrid to rectangular blockgrid
 /////////////////////////////////////////////////////////////
@@ -328,10 +329,19 @@ private:
 };
 /* @} */
 
-//template <class DTyp>
+
+/////////////////////////////////////////////////////////////
+// 4. Interpolate from blockgrid direct to any point
+/////////////////////////////////////////////////////////////
+
+
+/** \addtogroup InterpolationOperators **/
+/* @{ */
+/**
+ * Interpolate from blockgrid direct to any point
+***/
 class Interpolate_direct {
  public:
-
     Interpolate_direct (Blockgrid* bg):blockgrid(bg)
       ,idHexPrev(-1)
       ,iPrev(-1)
@@ -362,37 +372,31 @@ class Interpolate_direct {
       ,typCounter4(0)
       ,typCounter5(0){}
 
+    /**
+     * preparation for interpolation : calculates the neighbour index for each cell.
+     * also  calculates the neighbour cell, if in other blockgrid.
+    **/
     void init();
-    std::vector<std::vector<int> > calculateNeighbourIndexRelation(std::vector<std::vector<int> > inner, std::vector<std::vector<int> > outer);
-    std::vector<int> calculateNeighbourIndex(std::vector<std::vector<int> > relation, int id_hex_outside,int id_hex_inside, int i, int j, int k, int Nx, int Ny, int Nz);
-    std::vector<std::vector<int> > filterCorrectNeighbours(std::vector<std::vector<int> > outer);
-    bool compareIndicies(std::vector<std::vector<int> > inner, std::vector<std::vector<int> > outer, int notCheck);
-
-
-    std::vector<std::vector<int> > switchIJ(std::vector<std::vector<int> > v);
-    std::vector<std::vector<int> > switchIK(std::vector<std::vector<int> > v);
-    std::vector<std::vector<int> > switchJK(std::vector<std::vector<int> > v);
-    std::vector<std::vector<int> > invertI (std::vector<std::vector<int> > v);
-    std::vector<std::vector<int> > invertJ (std::vector<std::vector<int> > v);
-    std::vector<std::vector<int> > invertK (std::vector<std::vector<int> > v);
-
-
 
+    /**
+     * Iterates throuh all cells : if D3vector v is inside of a cell, the weighting of the 8 cell points are calculated (D3vector lambda).
+     * Uses previously cell indexes and its neighbours, since it's assumed that two consecutive find_cell(v) calls are close to each other.
+     * Saves : idNow, idPrev, ifPrevPrev, idPrevPrevPrev
+     *
+    **/
     void find_cell(D3vector v);
-    void setPrevIndex();
-    void setPrevPrevIndex();
-    void setPrevPrevPrevIndex();
-    int checkBox(int idHex, int i, int j, int k, D3vector v);
-    int checkBoxSurrounding(int idHex, int i, int j, int k, D3vector v);
-    int checkBoxSurroundingOptimized(int idHex, int i, int j, int k, D3vector v);
-    int checkCorner(int idHex, int i, int j, int k, D3vector v);
-    int checkEdge(int idHex, int i, int j, int k, D3vector v);
-    bool checkOverlapOfBoxes(D3vector vWSD, D3vector vENT, D3vector wWSD, D3vector wENT);
-    double calculateOverlapOfBoxes(D3vector vWSD, D3vector vENT, D3vector wWSD, D3vector wENT);
+    /**
+     * Print a vtk file with a box, surrounding the cell.
+    **/
     void writeBox(D3vector vWSD, D3vector vENT, std::string str);
+    /**
+     * Print a vtk file with a point.
+    **/
     static void writePoint(D3vector v, std::string str, int Counter);
+    /**
+     * Writes interpolation efficiency: how many direct hits, second try, iterate through all, ect...
+    **/
    void writeInterpolationEfficiency();
-    int checkForHexaNeighbours(int idHex, int i, int j, int k, D3vector v);
    //   void interpolate(Variable<DTyp>& u, DTyp* data, DTyp defaultInterpolation);
     template <class DTyp>
     DTyp evaluate(Variable<DTyp>& u);
@@ -417,11 +421,42 @@ class Interpolate_direct {
       std::vector<std::vector<std::vector<D3vector> > > arrayBoxWSDENT;
 
       std::vector<std::vector<std::vector<int> > > array_box_boundary;
-    //  template <class DTyp>
-    //  Variable<DTyp>* u;
-};
+private:
+      /**
+       * Necessary to calculate the neighbour relation.
+      **/
+      std::vector<std::vector<int> > calculateNeighbourIndexRelation(std::vector<std::vector<int> > inner, std::vector<std::vector<int> > outer);
+      std::vector<int> calculateNeighbourIndex(std::vector<std::vector<int> > relation, int id_hex_outside,int id_hex_inside, int i, int j, int k, int Nx, int Ny, int Nz);
+      std::vector<std::vector<int> > filterCorrectNeighbours(std::vector<std::vector<int> > outer);
+      bool compareIndicies(std::vector<std::vector<int> > inner, std::vector<std::vector<int> > outer, int notCheck);
+      std::vector<std::vector<int> > switchIJ(std::vector<std::vector<int> > v);
+      std::vector<std::vector<int> > switchIK(std::vector<std::vector<int> > v);
+      std::vector<std::vector<int> > switchJK(std::vector<std::vector<int> > v);
+      std::vector<std::vector<int> > invertI (std::vector<std::vector<int> > v);
+      std::vector<std::vector<int> > invertJ (std::vector<std::vector<int> > v);
+      std::vector<std::vector<int> > invertK (std::vector<std::vector<int> > v);
+
+      /**
+       * Necessary to evaluate the neighbour cell relation.
+      **/
+      void setPrevIndex();
+      void setPrevPrevIndex();
+      void setPrevPrevPrevIndex();
+      int checkBox(int idHex, int i, int j, int k, D3vector v);
+      int checkForHexaNeighbours(int idHex, int i, int j, int k, D3vector v);
+      int checkBoxSurrounding(int idHex, int i, int j, int k, D3vector v);
+      int checkBoxSurroundingOptimized(int idHex, int i, int j, int k, D3vector v);
+      int checkCorner(int idHex, int i, int j, int k, D3vector v);
+      int checkEdge(int idHex, int i, int j, int k, D3vector v);
+      bool checkOverlapOfBoxes(D3vector vWSD, D3vector vENT, D3vector wWSD, D3vector wENT);
+      double calculateOverlapOfBoxes(D3vector vWSD, D3vector vENT, D3vector wWSD, D3vector wENT);
 
 
+
+
+};
+
+/* @} */
 template<class DTyp>
 DTyp Interpolate_direct::evaluate(Variable<DTyp> &u)
 {

program/source/math_lib/math_lib.h

@@ -56,6 +56,10 @@ class D3vector {
   void Print();
   void Print(std::ofstream *Datei);
   void operator=(const D3vector& v) { x=v.x; y=v.y; z=v.z; }
+  void operator+=(const D3vector& v) { x+=v.x; y+=v.y; z+=v.z; }
+  void operator-=(const D3vector& v) { x-=v.x; y-=v.y; z-=v.z; }
+  void operator*=(const double v) { x*=v; y*=v; z*=v; }
+  void operator/=(const double v) { x/=v; y/=v; z/=v; }
   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 D3vector& v) {if(v.x < x && v.y < y && v.z < z){return true;} else {return false;} }
@@ -149,6 +153,7 @@ class D3matrix {
   double Determinante() {
     return   (x1 * (y2*z3 - y3*z2) - y1 * (x2*z3 - x3*z2) + z1 * (x2*y3 - x3*y2));
   }
+
     /*
   D3matrix(D3vector cx, D3vector cy, D3vector cz) {
     x1 = cx.x; y1 = cy.x;  z1 = cz.x;
@@ -178,7 +183,7 @@ class D3matrix {
   D3vector invert_apply(const D3vector& v) {
     double x,y,z;
     double det =  Determinante(); //x1 * (y2*z3 - y3*z2) - y1 * (x2*z3 - x3*z2) + z1 * (x2*y3 - x3*y2);
-
+   // std::cout << "det "<< det << std::endl;
     /*
     x =  (  (y2*z3 - y3*z2) * v.x - (x2*z3 - x3*z2) * v.y + (x2*y3 - x3*y2) * v.z ) / det;
     y =  (- (y1*z3 - y3*z1) * v.x + (x1*z3 - x3*z1) * v.y - (x1*y3 - x3*y1) * v.z ) / det;
@@ -192,6 +197,176 @@ class D3matrix {
     return D3vector(x,y,z);
   }
 
+  inline void Print() {
+    std::cout << "Matrix: " << std::endl;
+    std::cout << x1 << ", " << y1 << ", " << z1 << ";" <<std::endl;
+    std::cout << x2 << ", " << y2 << ", " << z2 << ";" <<std::endl;
+    std::cout << x3 << ", " << y3 << ", " << z3 << ";" <<std::endl;
+  }
+
+  D3matrix matrixMultiply(D3matrix m)
+  {
+      D3vector cx = {x1 * m.x1 + y1 * m.x2 +z1 * m.x3,
+                      x2 * m.x1 + y2 * m.x2 +z2 * m.x3,
+                      x3 * m.x1 + y3 * m.x2 +z3 * m.x3};
+      D3vector cy = {x1 * m.y1 + y1 * m.y2 +z1 * m.y3,
+                      x2 * m.y1 + y2 * m.y2 +z2 * m.y3,
+                      x3 * m.y1 + y3 * m.y2 +z3 * m.y3};
+      D3vector cz = {x1 * m.z1 + y1 * m.z2 +z1 * m.z3,
+                      x2 * m.z1 + y2 * m.z2 +z2 * m.z3,
+                      x3 * m.z1 + y3 * m.z2 +z3 * m.z3};
+
+
+      return D3matrix(cx,cy,cz);
+
+  }
+  D3vector vectorMultiply(D3vector m)
+  {
+//      D3vector res = {x1 * m.x + x2 * m.x + x3 * m.x,
+//                      y1 * m.y + y2 * m.y + y3 * m.y,
+//                      z1 * m.z + z2 * m.z + z3 * m.z};
+//      D3vector res = {x1 * m.x + x2 * m.y + x3 * m.z,
+//                      y1 * m.x + y2 * m.y + y3 * m.z,
+//                      z1 * m.x + z2 * m.y + z3 * m.z};
+      D3vector res = {x1 * m.x + y1 * m.y + z1 * m.z,
+                      x2 * m.x + y2 * m.y + z2 * m.z,
+                      x3 * m.x + y3 * m.y + z3 * m.z};
+
+
+      return res;
+
+  }
+
+  void transpose()
+  {
+      double temp;
+      temp = x2; x2 = y1 ; y1 = temp;
+      temp = x3; x3 = z1 ; z1 = temp;
+      temp = y3; y3 = z2 ; z2 = temp;
+
+  }
+  void invert_gauss_elimination()
+  {
+      int n;
+         double a[3][6];
+         for (int i = 0 ; i<3;i++)
+         {
+             for (int j = 0 ; j<6;j++)
+             {
+                 //std::cout << "i j " << i << " " << j << std::endl;
+                 a[i][j] = 0;
+             }
+         }
+         int order = 3;
+         n = 3;
+         a[0][0] = x1;
+         a[1][0] = y1;
+         a[2][0] = z1;
+         a[0][1] = x2;
+         a[1][1] = y2;
+         a[2][1] = z2;
+         a[0][2] = x3;
+         a[1][2] = y3;
+         a[2][2] = z3;
+
+//         a[0][0] = x1;
+//         a[1][0] = x2;
+//         a[2][0] = x3;
+//         a[0][1] = y1;
+//         a[1][1] = y2;
+//         a[2][1] = y3;
+//         a[0][2] = z1;
+//         a[1][2] = z2;
+//         a[2][2] = z3;
+
+
+             double temp;
+
+
+             for (int i = 0; i < order; i++) {
+
+                 for (int j = 0; j < 2 * order; j++) {
+
+                     // Add '1' at the diagonal places of
+                     // the matrix to create a identity matirx
+                     if (j == (i + order))
+                         a[i][j] = 1;
+                 }
+             }
+
+             for (int i = order - 1; i > 0; i--) {
+
+                  //Swapping each and every element of the two rows
+                  if (a[i - 1][0] < a[i][0])
+                  for (int j = 0; j < 2 * order; j++) {
+
+                         // Swapping of the row, if above
+                         // condition satisfied.
+                  temp = a[i][j];
+                  a[i][j] = a[i - 1][j];
+                  a[i - 1][j] = temp;
+                     }
+             }
+
+
+
+             // Replace a row by sum of itself and a
+             // constant multiple of another row of the matrix
+             for (int i = 0; i < order; i++) {
+
+                 for (int j = 0; j < order; j++) {
+
+                     if (j != i) {
+
+                         temp = a[j][i] / a[i][i];
+                         for (int k = 0; k < 2 * order; k++) {
+
+                             a[j][k] -= a[i][k] * temp;
+                         }
+                     }
+                 }
+             }
+
+             // Multiply each row by a nonzero integer.
+             // Divide row element by the diagonal element
+             for (int i = 0; i < order; i++) {
+
+                 temp = a[i][i];
+                 for (int j = 0; j < 2 * order; j++) {
+
+                     a[i][j] = a[i][j] / temp;
+                 }
+             }
+//             std::cout << "inverted matrix : " << std::endl;
+//             for (int i = 0; i < n; i++) {
+//                 for (int j = order; j < 6; j++) {
+//                     std::cout << a[i][j] << "  ";
+//                 }
+//                 std::cout<<std::endl;
+//             }
+
+             x1 = a[0][0+order];
+             y1 = a[1][0+order];
+             z1 = a[2][0+order];
+             x2 = a[0][1+order];
+             y2 = a[1][1+order];
+             z2 = a[2][1+order];
+             x3 = a[0][2+order];
+             y3 = a[1][2+order];
+             z3 = a[2][2+order];
+
+//             x1 = a[0][0+order];
+//             x2 = a[1][0+order];
+//             x3 = a[2][0+order];
+//             y1 = a[0][1+order];
+//             y2 = a[1][1+order];
+//             y3 = a[2][1+order];
+//             z1 = a[0][2+order];
+//             z2 = a[1][2+order];
+//             z3 = a[2][2+order];
+             return;
+  }
+
 
 };
 

program/source/pde_op/stenop_cc.h

@@ -89,11 +89,12 @@ void Sten_matrix<TYPE2, STENCIL_TYP>::Allocate_data()
 			if (sten_typ.Not_const_y(id)==false) Ny = 2;
 			if (sten_typ.Not_const_z(id)==false) Nz = 2;
 
-			//cout << "Nx: " << Nx << " Ny: " << Ny << " Nz: " << Nz << endl;
+            //cout << "Hex: "<<id<< "Nx: " << Nx << " Ny: " << Ny << " Nz: " << Nz << endl;
 
-			N_total = (Nx-1) * (Ny-1) * (Nz-1) * 27;
-			sten_m_hexahedra[id] = new TYPE2[N_total];
-			for (i=0;i<N_total;++i) sten_m_hexahedra[id][i] = (TYPE2)0;
+            N_total = (Nx-1) * (Ny-1) * (Nz-1) * 27;
+            sten_m_hexahedra[id] = new TYPE2[N_total];
+            for (i=0;i<N_total;++i)
+            {sten_m_hexahedra[id][i] = (TYPE2)0;}
 		} else sten_m_hexahedra[id] = NULL;
 	}