added new comparison operator to D3vector. Fixed bug in direct interpoaltor :...

added new comparison operator to D3vector. Fixed bug in direct interpoaltor : didnt work, if the gridsize was non constant in all blocks. Also added Markers to mark surfaces properly.

program/source/grid/blockgrid.cc

@@ -42,6 +42,26 @@
 int Blockgrid::id_count_grid = 0;
 
 
+//remove later
+void writePoint(D3vector v, std::string str, int Counter)
+{
+    std::stringstream ss;
+    ss  << Counter;
+    std::string num = ss.str();
+    std::string s = str + num;
+    ofstream myfile;
+      myfile.open ("C:/Users/rall/Documents/UG_Blocks_thermal/testergebnisse/" +  s + ".vtk");
+      myfile << "# vtk DataFile Version 4.2\n";
+      myfile << "vtk output\n";
+      myfile << "ASCII\n";
+      myfile << "DATASET POLYDATA\n";
+      myfile << "POINTS 1 float\n";
+      myfile << v.x << " " << v.y << " " << v.z << "\n" ;
+      myfile << "VERTICES 1 2\n";
+      myfile << "1 0 \n";
+    myfile.close();
+}
+
 
 Blockgrid::Blockgrid() {
   id_of_grid = id_count_grid; ++id_count_grid;
@@ -453,17 +473,8 @@ inline std::vector<int> find_p_quad ( dir3D quad, Hexahedron_el *hex )
 D3vector Blockgrid::Give_coord_hexahedron ( int id_hex,
     int i, int j, int k ) const
 {
-    if (id_hex < 0 || id_hex > 4)
-    {
-        cout << "doubt x " << endl;
-    }
 
-
-//    if (blockgrid_hexa_coordinates_set.size() > 0)
-//    {
-        //cout << "id_hex " << id_hex << " i " << i << " j " << j << " k " << k << endl;
-        //if (blockgrid_hexa_coordinates_set.at(id_hex).at(i + Give_Nx_hexahedron ( id_hex ) * (j  + k * Give_Ny_hexahedron ( id_hex ))  ) == true)
-        if (bg_coord != NULL && bg_coord->blockgrid_edge_coordinates_calculated)
+    if (bg_coord != NULL && bg_coord->blockgrid_edge_coordinates_calculated)
         {
             int Nx = Give_Nx_hexahedron ( id_hex ) +1;
             int Ny = Give_Ny_hexahedron ( id_hex ) +1;
@@ -471,7 +482,7 @@ D3vector Blockgrid::Give_coord_hexahedron ( int id_hex,
             //cout << "index " << i + Nx* (j  + k * (Ny))<< endl;
             return bg_coord->blockgrid_hexa_coordinates.at(id_hex).at(i + Nx* (j  + k * (Ny))  );
         }
-  //}
+
 
   double eta,xi,phi;
   D3vector Psi[12];
@@ -1014,6 +1025,8 @@ void Blockgrid_coordinates::init_blockgrid_coordinates()
 {
     blockgrid_edge_coordinates_calculated = false;
     blockgrid_hexa_coordinates_calculated = false;
+
+    blockgrid_hexa_boundaries_calculated = false;
     blockgrid_hexa_coordinates.resize(bg->Give_unstructured_grid()->Give_number_hexahedra());
     for (int id_hex = 0 ; id_hex < bg->Give_unstructured_grid()->Give_number_hexahedra() ; id_hex++)
     {
@@ -1071,6 +1084,7 @@ void Blockgrid_coordinates::init_blockgrid_coordinates_boundary()
                 {
 
                     //if ((i > 0 || i < Nx-1) ^ (j > 0 || j < Ny-1) ^ (k > 0 || k < Nz-1))
+                    //remove true here?
                     if (true || (i == 0 || i == Nx-1) || (j== 0 || j == Ny-1) || (k == 0 || k == Nz-1))
                     {
                         //inner loop
@@ -1101,6 +1115,8 @@ void Blockgrid_coordinates::init_blockgrid_coordinates_boundary()
                                             blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(iInner);
                                             blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(jInner);
                                             blockgrid_hexa_boundary.at(id_hex).at(convertedLocalIndex).push_back(kInner);
+
+
                                         }
 
                                     }

program/source/grid/examples_ug_optics.cc

@@ -148,7 +148,7 @@ D3vector transform_left_lens_inner_quad ( double t1, double t2, double* pointer_
     double r_global_data              = pointer_global_data[1];
     double curvatureLeft_global_data  = pointer_global_data[2];
     double curvatureRight_global_data = pointer_global_data[3];
-    double thickness_global_data          = pointer_global_data[4];
+    double thickness_global_data      = pointer_global_data[4];
     double offsetZ_global_data        = pointer_global_data[5];
     double z_left_inner_global_data   = pointer_global_data[6];
     double z_left_outer_global_data   = pointer_global_data[7];
@@ -162,7 +162,7 @@ D3vector transform_left_lens_inner_quad ( double t1, double t2, double* pointer_
 
 
 
-    double sign = (curvatureRight_global_data > 0) ? 1 : ((curvatureRight_global_data < 0) ? -1 : 0)  ;
+    double sign = (curvatureLeft_global_data > 0) ? 1 : ((curvatureLeft_global_data < 0) ? -1 : 0)  ;
 
 
     double actualX = + r_global_data * t1 - r_global_data * t2;
@@ -215,7 +215,8 @@ D3vector transform_left_lens_diag_NW_quad ( double t1, double t2, double* pointe
 
     actualX = -actualX;
     double radiusSquared = ((x+actualX)*(x+actualX)+(y+actualY)*(y+actualY));
-    double sign = (curvatureRight_global_data > 0) ? 1 : ((curvatureRight_global_data < 0) ? -1 : 0)  ;
+
+    double sign = (curvatureLeft_global_data > 0) ? 1 : ((curvatureLeft_global_data < 0) ? -1 : 0)  ;
 
     double z =offsetZ_global_data+sign*(-1*( sqrt(pow(curvatureLeft_global_data,2)-radiusSquared) +sign *(-curvatureLeft_global_data +( (1-t)*z_left_inner_global_data + t * z_left_outer_global_data))));
      if (std::isnan(z))
@@ -288,7 +289,8 @@ D3vector transform_left_lens_diag_SW_quad ( double t1, double t2, double* pointe
     double actualY = r_global_data -  r_global_data * t1 + (R_global_data-r_global_data) * t2 - (R_global_data-r_global_data) * t1 * t2;
     actualY = -actualY;
     actualX = -actualX;
-    double sign = (curvatureRight_global_data > 0) ? 1 : ((curvatureRight_global_data < 0) ? -1 : 0)  ;
+
+    double sign = (curvatureLeft_global_data > 0) ? 1 : ((curvatureLeft_global_data < 0) ? -1 : 0)  ;
 
     double radiusSquared = ((x+actualX)*(x+actualX)+(y+actualY)*(y+actualY));
 
@@ -323,7 +325,8 @@ D3vector transform_left_lens_diag_SE_quad ( double t1, double t2, double* pointe
     double actualX = 0             +  r_global_data * t1 + 0                             * t2 + (R_global_data-r_global_data) * t1 * t2;
     double actualY = r_global_data -  r_global_data * t1 + (R_global_data-r_global_data) * t2 - (R_global_data-r_global_data) * t1 * t2;
     actualY = -actualY;
-    double sign = (curvatureRight_global_data > 0) ? 1 : ((curvatureRight_global_data < 0) ? -1 : 0)  ;
+
+    double sign = (curvatureLeft_global_data > 0) ? 1 : ((curvatureLeft_global_data < 0) ? -1 : 0)  ;
     double radiusSquared = ((x+actualX)*(x+actualX)+(y+actualY)*(y+actualY));
     double z =offsetZ_global_data+sign*(-1*( sqrt(pow(curvatureLeft_global_data,2)-radiusSquared) +sign *(-curvatureLeft_global_data +( (1-t)*z_left_inner_global_data + t * z_left_outer_global_data))));
      if (std::isnan(z))
@@ -395,7 +398,8 @@ D3vector transform_left_lens_diag_NE_quad ( double t1, double t2, double* pointe
 
     double actualX = 0             +  r_global_data * t1 + 0                             * t2 + (R_global_data-r_global_data) * t1 * t2;
     double actualY = r_global_data -  r_global_data * t1 + (R_global_data-r_global_data) * t2 - (R_global_data-r_global_data) * t1 * t2;
-    double sign = (curvatureRight_global_data > 0) ? 1 : ((curvatureRight_global_data < 0) ? -1 : 0)  ;
+
+    double sign = (curvatureLeft_global_data > 0) ? 1 : ((curvatureLeft_global_data < 0) ? -1 : 0)  ;
 
     double radiusSquared = ((x+actualX)*(x+actualX)+(y+actualY)*(y+actualY));
 
@@ -431,6 +435,7 @@ D3vector transform_diag_inner_faces_NE_quad( double t1, double t2, double* point
 
 
     double sign = (curvatureRight_global_data > 0) ? 1 : ((curvatureRight_global_data < 0) ? -1 : 0)  ;
+    double signLeft = (curvatureLeft_global_data > 0) ? 1 : ((curvatureLeft_global_data < 0) ? -1 : 0)  ;
 
     double radiusSquared = r_global_data + (R_global_data-r_global_data)* (t2);
 
@@ -443,9 +448,9 @@ D3vector transform_diag_inner_faces_NE_quad( double t1, double t2, double* point
                           ( radiusSquared))
                       + sign*(-curvatureRight_global_data +( (1-t)*(thickness_global_data-z_right_inner_global_data) + t * (thickness_global_data-z_right_outer_global_data)))));
     double zLeft =   offsetZ_global_data
-                    +sign*(-1*( sqrt(pow(curvatureLeft_global_data,2)-
+                    +signLeft*(-1*( sqrt(pow(curvatureLeft_global_data,2)-
                             radiusSquared)
-                        +sign *(-curvatureLeft_global_data +( (1-t)*z_left_inner_global_data + t * z_left_outer_global_data))));
+                        +signLeft *(-curvatureLeft_global_data +( (1-t)*z_left_inner_global_data + t * z_left_outer_global_data))));
 
 
 
@@ -485,6 +490,7 @@ D3vector transform_inner_faces_NE_quad( double t1, double t2, double* pointer_gl
 
 
     double sign = (curvatureRight_global_data > 0) ? 1 : ((curvatureRight_global_data < 0) ? -1 : 0)  ;
+    double signLeft = (curvatureLeft_global_data > 0) ? 1 : ((curvatureLeft_global_data < 0) ? -1 : 0)  ;
 
     double x = r_global_data* (t1-1);
     double y = r_global_data* (t1);
@@ -494,8 +500,8 @@ D3vector transform_inner_faces_NE_quad( double t1, double t2, double* pointer_gl
                 +sign*(( sqrt(pow(curvatureRight_global_data,2)- radiusSquared) //pow(r_global_data,2) * (t*t + (1-t)*(1-t)))
                          +sign*(- curvatureRight_global_data + (thickness_global_data-z_right_inner_global_data)) ));
     double zLeft = offsetZ_global_data
-                +sign*(-1*( sqrt(pow(curvatureLeft_global_data,2)- radiusSquared)
-                            - sign* (curvatureLeft_global_data - z_left_inner_global_data)) );
+                +signLeft*(-1*( sqrt(pow(curvatureLeft_global_data,2)- radiusSquared)
+                            - signLeft* (curvatureLeft_global_data - z_left_inner_global_data)) );
 
     //works for one side only:
     double z = zRight * (t2) + zLeft * (1-t2);
@@ -678,9 +684,13 @@ D3vector transform_lens_SE_square ( double t, double* pointer_global_data ) {
 }
 
 
-Lens_Geometry::Lens_Geometry (double Radius, double radius, double thickness, double curvatureLeft, double curvatureRight, double offsetX, double offsetY, double offsetZ, bool inner_grid_arched)
+Lens_Geometry::Lens_Geometry (double Radius, double thickness, double curvatureLeft, double curvatureRight, double offsetX, double offsetY, double offsetZ, bool inner_grid_arched, double radius)
 {
-
+    if (radius == 0)
+    {
+        radius = Radius / 2.0;
+    }
+    assert(Radius > radius);
     size_pointer_global_data = 10;
     if (fabs(curvatureLeft) >1e10)
     {

program/source/grid/examples_ug_optics.h

@@ -8,7 +8,7 @@
 
 class Lens_Geometry : public Unstructured_grid {
  public:
-  Lens_Geometry(double Radius, double radius, double thickness, double curvatureLeft, double curvatureRight, double offsetX, double offsetY, double offsetZ, bool inner_grid_arched = false);
+  Lens_Geometry(double Radius, double thickness, double curvatureLeft, double curvatureRight, double offsetX, double offsetY, double offsetZ, bool inner_grid_arched = false, double radius = 0.0);
   ~Lens_Geometry(){};
 
   double Radius;

program/source/grid/marker.cc

@@ -213,6 +213,227 @@ void Marker::selectQuads(unsigned int firstQuad, unsigned int numQuads)
     }
 }
 
+void Marker::operatorAnd(Marker *mark)
+{
+    for ( int id=0; id<ug->Give_number_hexahedra(); ++id )
+    {
+        if (Give_marker<hexahedronEl> ( id ) == mark->Give_marker<hexahedronEl> ( id ))
+        {
+            //Set_marker<hexahedronEl> ( id, yes_mark );
+        }
+        else
+        {
+            Set_marker<hexahedronEl> ( id, no_mark );
+        }
+    }
+
+    for ( int id=0; id<ug->Give_number_edges(); ++id )
+    {
+        if (Give_marker<edgeEl> ( id ) == mark->Give_marker<edgeEl> ( id ))
+        {
+            //Set_marker<edgeEl> ( id, yes_mark );
+        }
+        else
+        {
+            Set_marker<edgeEl> ( id, no_mark );
+        }
+    }
+
+    for ( int id=0; id<ug->Give_number_points(); ++id )
+    {
+        if (Give_marker<pointEl> ( id ) == mark->Give_marker<pointEl> ( id ))
+        {
+            //Set_marker<pointEl> ( id, yes_mark );
+        }
+        else
+        {
+            Set_marker<pointEl> ( id, no_mark );
+        }
+    }
+
+    for ( int id=0; id<ug->Give_number_quadrangles(); ++id )
+    {
+        if (Give_marker<quadrangleEl> ( id ) == mark->Give_marker<quadrangleEl> ( id ))
+        {
+            //Set_marker<quadrangleEl> ( id, yes_mark );
+        }
+        else
+        {
+            Set_marker<quadrangleEl> ( id, no_mark );
+        }
+    }
+}
+
+void Marker::operatorOr(Marker *mark)
+{
+
+    for ( int id=0; id<ug->Give_number_hexahedra(); ++id )
+    {
+        if (Give_marker<hexahedronEl> ( id ) != mark->Give_marker<hexahedronEl> ( id ) || Give_marker<hexahedronEl> ( id ) == yes_mark)
+        {
+            Set_marker<hexahedronEl> ( id, yes_mark );
+        }
+        else
+        {
+            Set_marker<hexahedronEl> ( id, no_mark );
+        }
+    }
+
+    for ( int id=0; id<ug->Give_number_edges(); ++id )
+    {
+        if (Give_marker<edgeEl> ( id ) != mark->Give_marker<edgeEl> ( id ) || Give_marker<edgeEl> ( id ) == yes_mark)
+        {
+            Set_marker<edgeEl> ( id, yes_mark );
+        }
+        else
+        {
+            Set_marker<edgeEl> ( id, no_mark );
+        }
+    }
+
+    for ( int id=0; id<ug->Give_number_points(); ++id )
+    {
+        if (Give_marker<pointEl> ( id ) != mark->Give_marker<pointEl> ( id ) || Give_marker<pointEl> ( id ) == yes_mark)
+        {
+            Set_marker<pointEl> ( id, yes_mark );
+        }
+        else
+        {
+            Set_marker<pointEl> ( id, no_mark );
+        }
+    }
+
+    for ( int id=0; id<ug->Give_number_quadrangles(); ++id )
+    {
+        if (Give_marker<quadrangleEl> ( id ) != mark->Give_marker<quadrangleEl> ( id ) || Give_marker<quadrangleEl> ( id ) == yes_mark)
+        {
+            Set_marker<quadrangleEl> ( id, yes_mark );
+        }
+        else
+        {
+            Set_marker<quadrangleEl> ( id, no_mark );
+        }
+    }
+}
+
+
+void Marker::operatorXOr(Marker *mark)
+{
+
+    for ( int id=0; id<ug->Give_number_hexahedra(); ++id )
+    {
+        if (Give_marker<hexahedronEl> ( id ) != mark->Give_marker<hexahedronEl> ( id ))
+        {
+            Set_marker<hexahedronEl> ( id, yes_mark );
+        }
+        else
+        {
+            Set_marker<hexahedronEl> ( id, no_mark );
+        }
+    }
+
+    for ( int id=0; id<ug->Give_number_edges(); ++id )
+    {
+        if (Give_marker<edgeEl> ( id ) != mark->Give_marker<edgeEl> ( id ))
+        {
+            Set_marker<edgeEl> ( id, yes_mark );
+        }
+        else
+        {
+            Set_marker<edgeEl> ( id, no_mark );
+        }
+    }
+
+    for ( int id=0; id<ug->Give_number_points(); ++id )
+    {
+        if (Give_marker<pointEl> ( id ) != mark->Give_marker<pointEl> ( id ))
+        {
+            Set_marker<pointEl> ( id, yes_mark );
+        }
+        else
+        {
+            Set_marker<pointEl> ( id, no_mark );
+        }
+    }
+
+    for ( int id=0; id<ug->Give_number_quadrangles(); ++id )
+    {
+        if (Give_marker<quadrangleEl> ( id ) != mark->Give_marker<quadrangleEl> ( id ))
+        {
+            Set_marker<quadrangleEl> ( id, yes_mark );
+        }
+        else
+        {
+            Set_marker<quadrangleEl> ( id, no_mark );
+        }
+    }
+}
+
+
+void Marker::complement ( Marker* marker ) {
+    for ( int id=0; id<ug->Give_number_hexahedra(); ++id ) {
+        if ( marker->Give_marker<hexahedronEl> ( id ) == no_mark )
+            Set_marker<hexahedronEl> ( id, yes_mark );
+        else Set_marker<hexahedronEl> ( id, no_mark );
+    }
+    for ( int id=0; id<ug->Give_number_quadrangles(); ++id ) {
+        if ( marker->Give_marker<quadrangleEl> ( id ) == no_mark )
+            Set_marker<quadrangleEl> ( id, yes_mark );
+        else Set_marker<quadrangleEl> ( id, no_mark );
+    }
+    for ( int id=0; id<ug->Give_number_edges(); ++id ) {
+        if ( marker->Give_marker<edgeEl> ( id ) == no_mark )
+            Set_marker<edgeEl> ( id, yes_mark );
+        else Set_marker<edgeEl> ( id, no_mark );
+    }
+    for ( int id=0; id<ug->Give_number_points(); ++id ) {
+        if ( marker->Give_marker<pointEl> ( id ) == no_mark )
+            Set_marker<pointEl> ( id, yes_mark );
+        else Set_marker<pointEl> ( id, no_mark );
+    }
+}
+
+
+void Marker::interior_boundary ( Marker* marker ) {
+  // kopieren
+  for ( int id=0; id<ug->Give_number_hexahedra(); ++id ) {
+     if ( marker->Give_marker<hexahedronEl> ( id ) == no_mark )
+    Set_marker<hexahedronEl> ( id, no_mark );
+     else Set_marker<hexahedronEl> ( id, yes_mark );
+  }
+    for ( int id=0; id<ug->Give_number_quadrangles(); ++id ) {
+        if ( marker->Give_marker<quadrangleEl> ( id ) == no_mark )
+            Set_marker<quadrangleEl> ( id, no_mark );
+        else Set_marker<quadrangleEl> ( id, yes_mark );
+    }
+    for ( int id=0; id<ug->Give_number_edges(); ++id ) {
+        if ( marker->Give_marker<edgeEl> ( id ) == no_mark )
+            Set_marker<edgeEl> ( id, no_mark );
+        else Set_marker<edgeEl> ( id, yes_mark );
+    }
+    for ( int id=0; id<ug->Give_number_points(); ++id ) {
+        if ( marker->Give_marker<pointEl> ( id ) == no_mark )
+            Set_marker<pointEl> ( id, no_mark );
+        else Set_marker<pointEl> ( id, yes_mark );
+    }
+
+
+  //aendern
+  for ( int id_quad=0; id_quad<ug->Give_number_quadrangles(); ++id_quad ) {
+     if ( marker->Give_marker<quadrangleEl> ( id_quad ) == no_mark ) {
+
+       for(int e=0;e < 4;++e) {
+          Set_marker<edgeEl> (ug->Give_quadrangle ( id_quad )->Give_id_edge ( (dir2D)e ),
+                   no_mark );
+       }
+       for(int c=0;c < 4;++c) {
+          Set_marker<pointEl> (ug->Give_quadrangle ( id_quad )->Give_id_corner ( (dir2D_sons)c ),
+                no_mark );
+       }
+     }
+  }
+}
+
 Marker::~Marker() {
 	if ( marker_hexahedra != NULL )
 		delete [] marker_hexahedra;
@@ -274,6 +495,25 @@ void Boundary_Marker::markBoundaryDir3D(dir3D dir)
     }
 }
 
+void Boundary_Marker::unmarkBoundaryDir3D(dir3D dir)
+{
+    Quadrangle_el* quad;
+    Hexahedron_el* hex ;
+    for ( int id = 0;id < ug->Give_number_hexahedra();++id ) {
+            hex = ug->Give_hexahedron ( id );
+
+            if ( ug->Give_quadrangle ( hex->Give_id_quadrangle ( dir ) )->Give_exists_exterior() == false ) {
+
+            quad = ug->Give_quadrangle ( hex->Give_id_quadrangle ( dir ) );
+            Set_marker<quadrangleEl> ( hex->Give_id_quadrangle ( dir ), no_mark );
+            for ( int i=0; i<4; ++i )
+                Set_marker<edgeEl> ( quad->Give_id_edge ( ( dir2D ) i ), no_mark );
+            for ( int i=0; i<4; ++i )
+                Set_marker<pointEl> ( quad->Give_id_corner ( ( dir2D_sons ) i ), no_mark );
+        }
+    }
+}
+
 
 Unstructured_grid_Marker::Unstructured_grid_Marker ( Unstructured_grid* ug_ ) : Marker ( ug_ ) {
 	for ( int id=0; id<ug->Give_number_hexahedra(); ++id ) {
@@ -297,70 +537,6 @@ Blockgrid_Marker::Blockgrid_Marker ( Unstructured_grid* ug_ ) : Marker ( ug_ ) {
 
 Blockgrid_Marker::~Blockgrid_Marker() {}
 
-void Unstructured_grid_Marker::complement ( Marker* marker ) {
-	for ( int id=0; id<ug->Give_number_hexahedra(); ++id ) {
-		if ( marker->Give_marker<hexahedronEl> ( id ) == no_mark )
-			Set_marker<hexahedronEl> ( id, yes_mark );
-		else Set_marker<hexahedronEl> ( id, no_mark );
-	}
-	for ( int id=0; id<ug->Give_number_quadrangles(); ++id ) {
-		if ( marker->Give_marker<quadrangleEl> ( id ) == no_mark )
-			Set_marker<quadrangleEl> ( id, yes_mark );
-		else Set_marker<quadrangleEl> ( id, no_mark );
-	}
-	for ( int id=0; id<ug->Give_number_edges(); ++id ) {
-		if ( marker->Give_marker<edgeEl> ( id ) == no_mark )
-			Set_marker<edgeEl> ( id, yes_mark );
-		else Set_marker<edgeEl> ( id, no_mark );
-	}
-	for ( int id=0; id<ug->Give_number_points(); ++id ) {
-		if ( marker->Give_marker<pointEl> ( id ) == no_mark )
-			Set_marker<pointEl> ( id, yes_mark );
-		else Set_marker<pointEl> ( id, no_mark );
-	}
-}
-
-
-void Unstructured_grid_Marker::interior_boundary ( Marker* marker ) {
-  // kopieren
-  for ( int id=0; id<ug->Give_number_hexahedra(); ++id ) {
-     if ( marker->Give_marker<hexahedronEl> ( id ) == no_mark )
-	Set_marker<hexahedronEl> ( id, no_mark );
-     else Set_marker<hexahedronEl> ( id, yes_mark );
-  }
-  	for ( int id=0; id<ug->Give_number_quadrangles(); ++id ) {
-		if ( marker->Give_marker<quadrangleEl> ( id ) == no_mark )
-			Set_marker<quadrangleEl> ( id, no_mark );
-		else Set_marker<quadrangleEl> ( id, yes_mark );
-	}
-	for ( int id=0; id<ug->Give_number_edges(); ++id ) {
-		if ( marker->Give_marker<edgeEl> ( id ) == no_mark )
-			Set_marker<edgeEl> ( id, no_mark );
-		else Set_marker<edgeEl> ( id, yes_mark );
-	}
-	for ( int id=0; id<ug->Give_number_points(); ++id ) {
-		if ( marker->Give_marker<pointEl> ( id ) == no_mark )
-			Set_marker<pointEl> ( id, no_mark );
-		else Set_marker<pointEl> ( id, yes_mark );
-	}
-	
-	
-  //aendern
-  for ( int id_quad=0; id_quad<ug->Give_number_quadrangles(); ++id_quad ) {
-     if ( marker->Give_marker<quadrangleEl> ( id_quad ) == no_mark ) {
-       
-       for(int e=0;e < 4;++e) {
-          Set_marker<edgeEl> (ug->Give_quadrangle ( id_quad )->Give_id_edge ( (dir2D)e ),
-			       no_mark );
-       }
-       for(int c=0;c < 4;++c) {
-          Set_marker<pointEl> (ug->Give_quadrangle ( id_quad )->Give_id_corner ( (dir2D_sons)c ),
-				no_mark );
-       } 
-     }
-  }
-}
-
 
 
 void Unstructured_grid_Marker::operator= ( bool ( *Marker_formula ) ( double x, double y, double z ) ) {

program/source/grid/marker.h

@@ -56,6 +56,13 @@ public:
 		else off();
 	}
 
+    void operatorAnd(Marker * mark);
+    void operatorOr(Marker * mark);
+    void operatorXOr(Marker * mark);
+
+    void complement (Marker* marker);
+    void interior_boundary (Marker* marker);
+
 	void Mark_hexahedra_only( void ); 		///< Marks hexahedra and unmarks everything else
 	void Mark_quadrangles_only( void );	///< Marks quadrangles and unmarks everything else
 	void Mark_edges_only( void );			///< Marks edges and unmarks everything else
@@ -87,15 +94,14 @@ public:
 	~Boundary_Marker();
     //* mark boundary in dirction dir3D (Ndir3D, Edir3D... )
     void markBoundaryDir3D(dir3D dir);
+    void unmarkBoundaryDir3D(dir3D dir);
 
 };
 
 class Unstructured_grid_Marker :  public Marker {
 public:
 	Unstructured_grid_Marker(Unstructured_grid* ug_);
-	~Unstructured_grid_Marker();
-	void complement (Marker* marker);
-	void interior_boundary (Marker* marker);
+    ~Unstructured_grid_Marker();
 	void operator= (bool  (*Marker_formula)(double x, double y, double z));
 	void operator=(bool setValue) {
 		if(setValue) on();

program/source/interpol/interpol.h

@@ -342,7 +342,9 @@ class Interpolate_direct {
       ,kPrevPrevPrev(-1)
       ,counterFast(0)
       ,counterFastest(0)
+      ,counterSamePoint(0)
       ,counterSecondTry(0)
+      ,counterThirdTry(0)
       ,counterSlow(0)
       ,counterHexa(0)
       ,checkCounter(0)
@@ -356,7 +358,20 @@ class Interpolate_direct {
       ,typCounter5(0){}
 
     void init();
-    void initBoundary();
+    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);
+