rayspace interpolatir improved: added new efficient interpolation method, should work better

program/source/extemp/variable.h

@@ -514,6 +514,18 @@ inline double Variable<double>::Give_data<hexahedronEl> ( params_in ) const {
 template <>
 template <>
 inline void Variable<double>::Set_data<hexahedronEl> ( params_in, double value )  {
+    if (id == 4 && i == 8 && j == 8 && k == 8)
+    {
+        std::cout << "broh " << std::endl;
+    }
+    else if (id == 5 && i == 0 && j == 0 && k == 8)
+    {
+        std::cout << "broh " << std::endl;
+    }
+    else if (id == 8 && i == 0 && j == 8 && k == 8)
+    {
+        std::cout << "broh " << std::endl;
+    }
     data_hexahedra[id][i+ ( Nx+1 ) * ( j+ ( Ny+1 ) *k ) ] = value;
 }
 

program/source/grid/blockgrid.cc

@@ -1082,7 +1082,8 @@ void Blockgrid_coordinates::init_blockgrid_coordinates_boundary()
         int Ny = bg->Give_Ny_hexahedron(id_hex)+1;
         int Nz = bg->Give_Nz_hexahedron(id_hex)+1;
         blockgrid_hexa_boundary.at(id_hex).resize(Nx*Ny*Nz);
-                for(int i=0;i<Nx;++i) for(int j=0;j<Ny;++j) for(int k=0;k<Nz;++k)
+        //for(int i=0;i<Nx;++i) for(int j=0;j<Ny;++j) for(int k=0;k<Nz;++k)
+             for(int k=0;k<Nz;++k) for(int j=0;j<Ny;++j) for(int i=0;i<Nx;++i)
                 {
 
                     //if ((i > 0 || i < Nx-1) ^ (j > 0 || j < Ny-1) ^ (k > 0 || k < Nz-1))
@@ -1102,13 +1103,14 @@ void Blockgrid_coordinates::init_blockgrid_coordinates_boundary()
                             int NxInner = bg->Give_Nx_hexahedron(id_hexInner)+1;
                             int NyInner = bg->Give_Ny_hexahedron(id_hexInner)+1;
                             int NzInner = bg->Give_Nz_hexahedron(id_hexInner)+1;
-                                for(int iInner=0;iInner<NxInner;++iInner) for(int jInner=0;jInner<NyInner;++jInner) for(int kInner=0;kInner<NzInner;++kInner)
+                            //for(int iInner=0;iInner<NxInner;++iInner) for(int jInner=0;jInner<NyInner;++jInner) for(int kInner=0;kInner<NzInner;++kInner)
+                             for(int kInner=0;kInner<NzInner;++kInner) for(int jInner=0;jInner<NyInner;++jInner) for(int iInner=0;iInner<NxInner;++iInner)
                                 {
                                     // not sure which if loop is correct:
                                     //if ((iInner > 0 || iInner < NxInner-1) ^ (jInner > 0 || jInner < NyInner-1) ^ (kInner > 0 || kInner < NzInner-1))
                                     if (true || (iInner == 0 || iInner == NxInner-1) || (jInner == 0 || jInner == NyInner-1) || (kInner == 0 || kInner == NzInner-1))
                                     {
-                                        if(bg->Give_coord_hexahedron(id_hex,i,j,k) == bg->Give_coord_hexahedron(id_hexInner,iInner,jInner,kInner) )
+                                        if(D3VectorNormSquared(bg->Give_coord_hexahedron(id_hex,i,j,k) - bg->Give_coord_hexahedron(id_hexInner,iInner,jInner,kInner)) < 1e-10 )
                                         {
 
                                             int convertedLocalIndex = i    +(Nx)*(j    +(Ny)* k);

program/source/grid/blockgrid.h

@@ -49,6 +49,7 @@ class Blockgrid_coordinates {
      * */
     void init_blockgrid_coordinates_boundary();
     void delete_blockgrid_coordinates();
+    std::vector<std::vector<std::vector<int> > > getHexaCoordiantesBoundary(){return blockgrid_hexa_boundary;}
     std::vector<std::vector<std::vector<int> > > blockgrid_hexa_boundary;
 
     std::vector< std::vector< D3vector > > blockgrid_hexa_coordinates;

program/source/grid/examples_ug_optics.cc

@@ -1886,8 +1886,18 @@ Lens_Geometry_cutted_edges::Lens_Geometry_cutted_edges(double RadiusLeft, double
     }
 
 
-    constructionParameters[10] = cutLeft;
-    constructionParameters[11] = cutRight;
+    if (cutLeft)
+    {constructionParameters[10] = Radius;}
+    else
+    {
+       constructionParameters[10] = 0;
+    }
+    if (cutRight)
+    {constructionParameters[11] = Radius;}
+    else
+    {
+       constructionParameters[11] = 0;
+    }
 
 
 //    if (curvatureRight < 0 && curvatureLeft < 0)
@@ -2042,13 +2052,6 @@ Lens_Geometry_cutted_edges::Lens_Geometry_cutted_edges(double RadiusLeft, double
   Set_transformation_face(16,19,20,23,transform_outer_boundary_SE_cut);
   Set_transformation_face(9,20,15,23,transform_right_lens_diag_SE_quad_cut);
 
-  //unteil here
-
-
-
-
-
-
     //inner block: corner ids: 0 2 4 6 8 10 12 14
 
   Set_transformation_face(0,2,4,6,transform_left_lens_inner_quad); // z=0 plane transform_left_lens_inner_quad
@@ -2095,6 +2098,8 @@ Lens_Geometry_cutted_edges::Lens_Geometry_cutted_edges(double RadiusLeft, double
 
 
 
+  // WSDdir3D, ESDdir3D, WNDdir3D, ENDdir3D, WSTdir3D, ESTdir3D, WNTdir3D, ENTdir3D
 
   construction_done();
+
 }

program/source/grid/marker.cc

@@ -476,11 +476,22 @@ Boundary_Marker::Boundary_Marker ( Unstructured_grid* ug_ ) : Marker ( ug_ )
 
 Boundary_Marker::~Boundary_Marker() {}
 
-void Boundary_Marker::markBoundaryDir3D(dir3D dir)
-{
+void Boundary_Marker::markBoundaryDir3D(dir3D dir, int blockRangeStart, int blockRangeEnd)
+    {
+    if (blockRangeStart == -1)
+    {return;}
+    if (blockRangeEnd < 0 || blockRangeEnd > ug->Give_number_hexahedra())
+    {
+      blockRangeEnd =   ug->Give_number_hexahedra();
+    }
+    if (blockRangeStart < 0 || blockRangeStart > ug->Give_number_hexahedra())
+    {
+      blockRangeStart =   0;
+    }
+
     Quadrangle_el* quad;
     Hexahedron_el* hex ;
-    for ( int id = 0;id < ug->Give_number_hexahedra();++id ) {
+    for ( int id = blockRangeStart;id < blockRangeEnd;++id ) {
             hex = ug->Give_hexahedron ( id );
 
             if ( ug->Give_quadrangle ( hex->Give_id_quadrangle ( dir ) )->Give_exists_exterior() == false ) {
@@ -492,14 +503,25 @@ void Boundary_Marker::markBoundaryDir3D(dir3D dir)
             for ( int i=0; i<4; ++i )
                 Set_marker<pointEl> ( quad->Give_id_corner ( ( dir2D_sons ) i ), yes_mark );
         }
+        }
     }
-}
 
-void Boundary_Marker::unmarkBoundaryDir3D(dir3D dir)
+void Boundary_Marker::unmarkBoundaryDir3D(dir3D dir, int blockRangeStart, int blockRangeEnd)
 {
+    if (blockRangeStart == -1)
+    {return;}
+    if (blockRangeEnd < 0 || blockRangeEnd > ug->Give_number_hexahedra())
+    {
+      blockRangeEnd =   ug->Give_number_hexahedra();
+    }
+    if (blockRangeStart < 0 || blockRangeStart > ug->Give_number_hexahedra())
+    {
+      blockRangeStart =   0;
+    }
+
     Quadrangle_el* quad;
     Hexahedron_el* hex ;
-    for ( int id = 0;id < ug->Give_number_hexahedra();++id ) {
+    for ( int id = blockRangeStart;id < blockRangeStart;++id ) {
             hex = ug->Give_hexahedron ( id );
 
             if ( ug->Give_quadrangle ( hex->Give_id_quadrangle ( dir ) )->Give_exists_exterior() == false ) {

program/source/grid/marker.h

@@ -93,8 +93,8 @@ public:
 	Boundary_Marker(Unstructured_grid* ug_);
 	~Boundary_Marker();
     //* mark boundary in dirction dir3D (Ndir3D, Edir3D... )
-    void markBoundaryDir3D(dir3D dir);
-    void unmarkBoundaryDir3D(dir3D dir);
+    void markBoundaryDir3D(dir3D dir, int blockRangeStart = 0, int blockRangeEnd = -1 );
+    void unmarkBoundaryDir3D(dir3D dir, int blockRangeStart = 0, int blockRangeEnd = -1 );
 
 };
 

program/source/interpol/interpol.cc

@@ -2472,43 +2472,6 @@ std::vector<int> Interpolate_direct::calculateNeighbourIndex(std::vector<std::ve
     int indexConverted = convertedI    +Nx*(convertedJ    +Ny* convertedK);
     //indexConverted = 0;
 
-    //int indexTemp = iTemp    +NxOuter*(jTemp    +NyOuter* convertedK); // to be saved!!!
-//    D3vector boxWSDInner =  arrayBoxWSDENT.at(id_hex_inside).at(indexConverted).at(0);
-//    D3vector boxENTInner =  arrayBoxWSDENT.at(id_hex_inside).at(indexConverted).at(1);
-
-//    D3vector cWSD, cESD;
-//    D3vector cWND, cEND;
-
-//    D3vector cWST, cEST;
-//    D3vector cWNT, cENT;
-
-//    D3vector boxWSD, boxENT;
-//    cWSD = blockgrid->Give_coord_hexahedron(id_hex_outside,iTemp,  jTemp,  kTemp  );
-//    cESD = blockgrid->Give_coord_hexahedron(id_hex_outside,iTemp+1,jTemp  ,kTemp  );
-//    cWND = blockgrid->Give_coord_hexahedron(id_hex_outside,iTemp,  jTemp+1,kTemp  );
-//    cEND = blockgrid->Give_coord_hexahedron(id_hex_outside,iTemp+1,jTemp+1,kTemp  );
-
-//    cWST = blockgrid->Give_coord_hexahedron(id_hex_outside,iTemp,  jTemp,  kTemp+1);
-//    cEST = blockgrid->Give_coord_hexahedron(id_hex_outside,iTemp+1,jTemp  ,kTemp+1);
-//    cWNT = blockgrid->Give_coord_hexahedron(id_hex_outside,iTemp,  jTemp+1,kTemp+1);
-//    cENT = blockgrid->Give_coord_hexahedron(id_hex_outside,iTemp+1,jTemp+1,kTemp+1);
-
-//        // bounding box calculation
-//    boxWSD.x = MIN(MIN(MIN(cWSD.x,cESD.x),MIN(cWND.x,cEND.x)),
-//           MIN(MIN(cWST.x,cEST.x),MIN(cWNT.x,cENT.x)));// - factor *hx;
-//    boxWSD.y = MIN(MIN(MIN(cWSD.y,cESD.y),MIN(cWND.y,cEND.y)),
-//           MIN(MIN(cWST.y,cEST.y),MIN(cWNT.y,cENT.y)));// - factor *hy;
-//    boxWSD.z = MIN(MIN(MIN(cWSD.z,cESD.z),MIN(cWND.z,cEND.z)),
-//           MIN(MIN(cWST.z,cEST.z),MIN(cWNT.z,cENT.z)));// - factor *hz;
-
-//    boxENT.x = MAX(MAX(MAX(cWSD.x,cESD.x),MAX(cWND.x,cEND.x)),
-//           MAX(MAX(cWST.x,cEST.x),MAX(cWNT.x,cENT.x)));// + factor *hx;
-//    boxENT.y = MAX(MAX(MAX(cWSD.y,cESD.y),MAX(cWND.y,cEND.y)),
-//           MAX(MAX(cWST.y,cEST.y),MAX(cWNT.y,cENT.y)));// + factor *hy;
-//    boxENT.z = MAX(MAX(MAX(cWSD.z,cESD.z),MAX(cWND.z,cEND.z)),
-//           MAX(MAX(cWST.z,cEST.z),MAX(cWNT.z,cENT.z)));// + factor *hz;
-    //D3vector boxWSDOuter = arrayBoxWSDENT.at(idHexTempTest).at(indexTemp).at(0);
-    //D3vector boxENTOuter = arrayBoxWSDENT.at(idHexTempTest).at(indexTemp).at(1);
 
     //check, if boundary already written here
     bool write = true;
@@ -2521,6 +2484,10 @@ std::vector<int> Interpolate_direct::calculateNeighbourIndex(std::vector<std::ve
     }
     if (write)
     {
+        if (array_box_boundary.at(id_hex_inside).at(indexConverted).size() > 9)
+        {
+            std::cout << "4 edges?? " << std::endl;
+        }
         array_box_boundary.at(id_hex_inside).at(indexConverted).push_back(id_hex_outside);
         array_box_boundary.at(id_hex_inside).at(indexConverted).push_back(iTemp);
         array_box_boundary.at(id_hex_inside).at(indexConverted).push_back(jTemp);
@@ -3316,7 +3283,7 @@ void Interpolate_direct::writePoint(D3vector v, std::string str, int Counter)
     ss << std::setw(5) << std::setfill('0') << Counter;
     std::string s = ss.str();
     ofstream myfile;
-      myfile.open ("C:/Users/rall/Documents/UG_Blocks_thermal/testergebnisse/" + s + str + ".vtk");
+      myfile.open ("testergebnisse/" + s + str + ".vtk");
       myfile << "# vtk DataFile Version 4.2\n";
       myfile << "vtk output\n";
       myfile << "ASCII\n";

program/source/interpol/interpol.h

@@ -373,13 +373,15 @@ class Interpolate_direct {
     template <class DTyp>
     DTyp evaluate(Variable<DTyp>& u);
 
+    std::vector<std::vector<std::vector<int> > > getBoundaryBox(){return array_box_boundary;}
+
     int counterFast{}, counterFastest{}, counterSamePoint{}, counterSecondTry{}, counterThirdTry{}, counterSlow{}, counterHexa{}, counterCorner{}, counterEdge{};
     int checkCounter{};
     int boxCounter{};
     int typCounter0{},typCounter1{},typCounter2{},typCounter3{},typCounter4{},typCounter5{};
     bool debugTest;
     D3vector lambda;
-    D3vector vNow{9124,123,52143}, vPrev, vPrevPrev, vPrevPrevPrev;
+    D3vector vNow{1e10,1e10,1e10}, vPrev, vPrevPrev, vPrevPrevPrev;
 
 
         bool vectorInBox(D3vector vWSD, D3vector vENT, D3vector v, double eps = 1e-10);
@@ -399,6 +401,7 @@ class Interpolate_direct {
       std::vector<std::vector<std::vector<D3vector> > > arrayBoxWSDENT;
 
       std::vector<std::vector<std::vector<int> > > array_box_boundary;
+      std::vector<std::vector<std::vector<int> > > array_point_boundary;
 private:
       /**
        * Necessary to calculate the neighbour relation.