diff --git a/example/ALL_test.cpp b/example/ALL_test.cpp
index dfd9876cc54e17d5efd9a9adb3ee5cb0a7c934e4..919e3a54768a8aefb54c71bf35ca7639bbe29233 100644
--- a/example/ALL_test.cpp
+++ b/example/ALL_test.cpp
@@ -62,8 +62,6 @@
 #define MAX_NEIG 1024
 #define ALL_HISTOGRAM_DEFAULT_WIDTH 1.0
 
-#undef ALL_VTK_OUTPUT
-
 #ifdef ALL_VORONOI
 #define ALL_VORONOI_PREP_STEPS 50
 #endif
@@ -294,20 +292,24 @@ void print_points(
     auto writer = vtkSmartPointer<vtkXMLPolyDataWriter>::New();
     writer->SetInputData(polydata);
     writer->SetFileName(ss_local.str().c_str());
+    writer->SetDataModeToAscii();
     writer->Write();
 
     std::ostringstream ss_para;
     ss_para << "vtk_points/ALL_vtk_points_" 
         << std::setw(7) << std::setfill('0') << step << ".pvtp";
 
-    auto parallel_writer = vtkSmartPointer<vtkXMLPPolyDataWriter>::New();
-    parallel_writer->SetFileName(ss_para.str().c_str());
-    parallel_writer->SetNumberOfPieces(n_ranks);
-    parallel_writer->SetStartPiece(rank);
-    parallel_writer->SetEndPiece(rank);
-    parallel_writer->SetInputData(polydata);
-    parallel_writer->Write();
 
+        auto parallel_writer = vtkSmartPointer<vtkXMLPPolyDataWriter>::New();
+        parallel_writer->SetFileName(ss_para.str().c_str());
+        parallel_writer->SetNumberOfPieces(n_ranks);
+        parallel_writer->SetStartPiece(rank);
+        parallel_writer->SetEndPiece(rank);
+        parallel_writer->SetInputData(polydata);
+        parallel_writer->SetDataModeToAscii();
+        parallel_writer->Update();
+    MPI_Barrier(comm);
+        parallel_writer->Write();
 }
 #endif
 
@@ -1020,6 +1022,9 @@ int main(int argc, char** argv)
         double limit_efficiency = 0.5;
         // create ALL object
         ALL<double,double> lb_obj(chosen_method,sys_dim,vertices,gamma);
+        lb_obj.set_proc_grid_params(local_coords,global_dim);
+        lb_obj.set_communicator(MPI_COMM_WORLD);
+        lb_obj.setup();
         for (int i_loop = 0; i_loop < max_loop; ++i_loop)
         {
             MPI_Barrier(cart_comm);
@@ -1209,17 +1214,14 @@ int main(int argc, char** argv)
             MPI_Barrier(cart_comm);
             if( local_rank == 0 ) std::cout << "Setting process grid information..." << std::endl;
 #endif
-            lb_obj.set_proc_grid_params(local_coords,global_dim);
 #ifdef ALL_DEBUG_ENABLED
             MPI_Barrier(cart_comm);
             if( local_rank == 0 ) std::cout << "Setting communicator..." << std::endl;
 #endif
-            lb_obj.set_communicator(MPI_COMM_WORLD);
 #ifdef ALL_DEBUG_ENABLED
             MPI_Barrier(cart_comm);
             if( local_rank == 0 ) std::cout << "Setting up method..." << std::endl;
 #endif
-            lb_obj.setup();
 #ifdef ALL_DEBUG_ENABLED
             MPI_Barrier(cart_comm);
             if( local_rank == 0 ) std::cout << "Setting method data..." << std::endl;
@@ -1257,8 +1259,8 @@ int main(int argc, char** argv)
             MPI_Barrier(cart_comm);
             if( local_rank == 0 ) std::cout << "Updating vertices..." << std::endl;
 #endif
-            new_vertices = lb_obj.get_result_vertices();
-            old_vertices = vertices;
+            new_vertices = lb_obj.get_vertices();
+            old_vertices = lb_obj.get_prev_vertices();
             vertices = new_vertices;
 
             if (chosen_method == ALL_LB_t::VORONOI)
diff --git a/include/ALL.hpp b/include/ALL.hpp
index 8532992b11c4de120b0737332b0eddfc591ea3f5..c60da2ed776ce1343f1f187a1bd5c83c5a739e7c 100644
--- a/include/ALL.hpp
+++ b/include/ALL.hpp
@@ -183,7 +183,7 @@ template <class T, class W> class ALL
 
         // vertices 
         std::vector<ALL_Point<T>>& get_vertices() {return balancer->get_vertices();};
-        std::vector<ALL_Point<T>>& get_result_vertices() {return balancer->get_shifted_vertices();};
+        std::vector<ALL_Point<T>>& get_prev_vertices() {return balancer->get_prev_vertices();};
 
         // dimension
         int get_dimension() {return balancer->get_dimension();}
@@ -570,7 +570,7 @@ template <class T, class W> std::vector<ALL_Point<T>>& ALL<T,W>::balance(bool in
         std::cout << e.what() << std::endl;
         MPI_Abort(MPI_COMM_WORLD,-1);
     }
-    return balancer->get_shifted_vertices();
+    return balancer->get_vertices();
 
 }
 
diff --git a/include/ALL_Histogram.hpp b/include/ALL_Histogram.hpp
index 6e9edd77aad9a1c22571f89a1654708825e8767b..42d0cb2f36cb519aaf0ac9dbb2a0ab4946e8d027 100644
--- a/include/ALL_Histogram.hpp
+++ b/include/ALL_Histogram.hpp
@@ -359,10 +359,10 @@ template<class T, class W> void ALL_Histogram_LB<T,W>::balance(int step)
     // size of one slice
     T size = this->sys_size[2*i+1] / (T)work_collection.size();
 
-    this->shifted_vertices = this->vertices;
+    this->prev_vertices = this->vertices;
 
-    this->shifted_vertices.at(0).set_coordinate(i, (T)down * size);
-    this->shifted_vertices.at(1).set_coordinate(i, (T)up * size);
+    this->vertices.at(0).set_coordinate(i, (T)down * size);
+    this->vertices.at(1).set_coordinate(i, (T)up * size);
 
     // compute min / max work in current dimension
     MPI_Allreduce(work_new_dimension.data()+this->local_coords.at(i), 
@@ -432,8 +432,8 @@ template<class T, class W> void ALL_Histogram_LB<T,W>::find_neighbors()
     MPI_Cart_shift(this->global_comm,1,1,&rank_left,&rank_right);
 
     // collect border information from local column
-    vertices_loc[0] = this->shifted_vertices.at(0).x(0);
-    vertices_loc[1] = this->shifted_vertices.at(1).x(0);
+    vertices_loc[0] = this->vertices.at(0).x(0);
+    vertices_loc[1] = this->vertices.at(1).x(0);
     MPI_Allgather(vertices_loc,
                   2,
                   mpi_data_type_T,
@@ -553,10 +553,10 @@ template<class T, class W> void ALL_Histogram_LB<T,W>::find_neighbors()
     MPI_Cart_shift(this->global_comm,2,1,&rank_left,&rank_right);
 
     // collect border information from local column
-    vertices_loc[0] = this->shifted_vertices.at(0).x(0);
-    vertices_loc[1] = this->shifted_vertices.at(1).x(0);
-    vertices_loc[2] = this->shifted_vertices.at(0).x(1);
-    vertices_loc[3] = this->shifted_vertices.at(1).x(1);
+    vertices_loc[0] = this->vertices.at(0).x(0);
+    vertices_loc[1] = this->vertices.at(1).x(0);
+    vertices_loc[2] = this->vertices.at(0).x(1);
+    vertices_loc[3] = this->vertices.at(1).x(1);
 
     MPI_Barrier(this->global_comm);
 
diff --git a/include/ALL_LB.hpp b/include/ALL_LB.hpp
index 1db74299dfe011c83bfca308962518571062eaaf..59dcea277abaca7189a7ed2fa4f1088faeaf2612 100644
--- a/include/ALL_LB.hpp
+++ b/include/ALL_LB.hpp
@@ -85,7 +85,7 @@ template <class T, class W> class ALL_LB {
     
     //virtual void set_vertices(std::vector<T>&);
     virtual std::vector<ALL_Point<T>>& get_vertices();
-    virtual std::vector<ALL_Point<T>>& get_shifted_vertices();
+    virtual std::vector<ALL_Point<T>>& get_prev_vertices();
     
     virtual void set_communicator(MPI_Comm comm) 
     {
@@ -95,8 +95,6 @@ template <class T, class W> class ALL_LB {
         MPI_Comm_rank(comm,&local_rank);
     }
     
-    void resize_shifted_vertices(int);
-    
     int get_n_vertices() {return vertices.size();}
     
     virtual void set_data(void*) = 0;
@@ -118,25 +116,28 @@ protected:
     std::vector<W> work;
     // vertices
     std::vector<ALL_Point<T>> vertices;
-    // shifted vertices
-    std::vector<ALL_Point<T>> shifted_vertices;
+    // previous (original) vertices
+    std::vector<ALL_Point<T>> prev_vertices;
     // global dimension of the system in each coordinate
     std::vector<int> global_dims;
     // coordinates of the local domain in the system
     std::vector<int> local_coords;
     // periodicity of the MPI communicator / system
     std::vector<int> periodicity;
+    // resize the vector to store the vertices
+    void resize_vertices(int);
+    
 };
 
 // set the actual vertices (unsafe due to array copy)
 template <class T, class W> void ALL_LB<T,W>::set_vertices(std::vector<ALL_Point<T>>& vertices_in)
 {
-    vertices = vertices_in;
-    // as a basic assumption the number of shifted vertices
+    // as a basic assumption the number of resulting vertices
     // is equal to the number of input vertices:
     // exceptions:
     //      VORONOI
-    shifted_vertices = vertices;
+    vertices = vertices_in;
+    prev_vertices = vertices_in;
 }
 
 template <class T, class W> std::vector<ALL_Point<T>>& ALL_LB<T,W>::get_vertices()
@@ -144,13 +145,13 @@ template <class T, class W> std::vector<ALL_Point<T>>& ALL_LB<T,W>::get_vertices
     return vertices;
 }
 
-template <class T, class W> std::vector<ALL_Point<T>>& ALL_LB<T,W>::get_shifted_vertices()
+template <class T, class W> std::vector<ALL_Point<T>>& ALL_LB<T,W>::get_prev_vertices()
 {
-    return shifted_vertices;
+    return prev_vertices;
 }
 
-template <class T, class W> void ALL_LB<T,W>::resize_shifted_vertices(int new_size)
+template <class T, class W> void ALL_LB<T,W>::resize_vertices(int new_size)
 {
-    shifted_vertices.resize(new_size);
+    vertices.resize(new_size);
 }
 #endif  // ALL_LB_HEADER_INCLUDED
diff --git a/include/ALL_Staggered.hpp b/include/ALL_Staggered.hpp
index 5b159842bf2f5eadcc9853adb9a399656fd8ec0a..1619463c5b7b5d33dbb07ec2bc1a817b71e1ae3e 100644
--- a/include/ALL_Staggered.hpp
+++ b/include/ALL_Staggered.hpp
@@ -187,6 +187,7 @@ template <class T, class W> void ALL_Staggered_LB<T,W>::setup()
 
 template<class T, class W> void ALL_Staggered_LB<T,W>::balance(int)
 {
+    this->prev_vertices = this->vertices;
     int dimension = this->get_dimension();
     
     // loop over all available dimensions
@@ -280,25 +281,25 @@ template<class T, class W> void ALL_Staggered_LB<T,W>::balance(int)
         // if a left neighbor exists: shift left border
                 // if a left neighbor exists: shift left border 
         if (rank_left != MPI_PROC_NULL && this->local_coords[i] != 0)
-            this->shifted_vertices.at(0).set_coordinate(i, this->vertices.at(0).x(i) + remote_shift);
+            this->vertices.at(0).set_coordinate(i, this->prev_vertices.at(0).x(i) + remote_shift);
         else
-            this->shifted_vertices.at(0).set_coordinate(i, this->vertices.at(0).x(i));
+            this->vertices.at(0).set_coordinate(i, this->prev_vertices.at(0).x(i));
 
         // if a right neighbor exists: shift right border
         if (rank_right != MPI_PROC_NULL && this->local_coords[i] != this->global_dims[i]-1)
-            this->shifted_vertices.at(1).set_coordinate(i, this->vertices.at(1).x(i) + shift);
+            this->vertices.at(1).set_coordinate(i, this->prev_vertices.at(1).x(i) + shift);
         else
-            this->shifted_vertices.at(1).set_coordinate(i, this->vertices.at(1).x(i));
+            this->vertices.at(1).set_coordinate(i, this->prev_vertices.at(1).x(i));
 
 
         // check if vertices are crossed and throw exception if something went wrong
-        if (this->shifted_vertices.at(1).x(i) < this->shifted_vertices.at(1).x(i))
+        if (this->vertices.at(1).x(i) < this->vertices.at(1).x(i))
         {
             std::cout << "ERROR on process: " << this->local_rank << " "
+                                              << this->prev_vertices.at(0).x(i) << " "
+                                              << this->prev_vertices.at(1).x(i) << " "
                                               << this->vertices.at(0).x(i) << " "
                                               << this->vertices.at(1).x(i) << " "
-                                              << this->shifted_vertices.at(0).x(i) << " "
-                                              << this->shifted_vertices.at(1).x(i) << " "
                                               << remote_shift << " "
                                               << shift << " "
                                               << remote_size << " "
@@ -326,7 +327,6 @@ template<class T, class W> void ALL_Staggered_LB<T,W>::find_neighbors()
     int dimension = this->get_dimension();
     auto work = this->get_work();
     auto vertices = this->get_vertices();
-    auto shifted_vertices = this->get_shifted_vertices();
     
     neighbors.clear();
     // collect work from right neighbor plane
@@ -360,8 +360,8 @@ template<class T, class W> void ALL_Staggered_LB<T,W>::find_neighbors()
     MPI_Cart_shift(this->global_comm,1,1,&rank_left,&rank_right);
 
     // collect border information from local column
-    vertices_loc[0] = shifted_vertices.at(0).x(0);
-    vertices_loc[1] = shifted_vertices.at(1).x(0);
+    vertices_loc[0] = vertices.at(0).x(0);
+    vertices_loc[1] = vertices.at(1).x(0);
     MPI_Allgather(vertices_loc,2,mpi_data_type_T,vertices_rem+2*this->global_dims[0],2,mpi_data_type_T,communicators[1]);
 
     // exchange local column information with upper neighbor in Y direction (cart grid)
@@ -404,7 +404,7 @@ template<class T, class W> void ALL_Staggered_LB<T,W>::find_neighbors()
     MPI_Barrier(this->global_comm);
 
     // exchange local column information with lower neighbor in Y direction (cart grid)
-    MPI_Irecv(vertices_rem                 ,2*this->global_dims[0],mpi_data_type_T,rank_right, 0,this->global_comm,&rreq);
+    MPI_Irecv(vertices_rem                       ,2*this->global_dims[0],mpi_data_type_T,rank_right, 0,this->global_comm,&rreq);
     MPI_Isend(vertices_rem+2*this->global_dims[0],2*this->global_dims[0],mpi_data_type_T,rank_left,  0,this->global_comm,&sreq);
 
     // determine the offset in ranks
@@ -446,17 +446,17 @@ template<class T, class W> void ALL_Staggered_LB<T,W>::find_neighbors()
     MPI_Cart_shift(this->global_comm,2,1,&rank_left,&rank_right);
 
     // collect border information from local column
-    vertices_loc[0] = shifted_vertices.at(0).x(0);
-    vertices_loc[1] = shifted_vertices.at(1).x(0);
-    vertices_loc[2] = shifted_vertices.at(0).x(1);
-    vertices_loc[3] = shifted_vertices.at(1).x(1);
+    vertices_loc[0] = vertices.at(0).x(0);
+    vertices_loc[1] = vertices.at(1).x(0);
+    vertices_loc[2] = vertices.at(0).x(1);
+    vertices_loc[3] = vertices.at(1).x(1);
 
     MPI_Barrier(this->global_comm);
 
     MPI_Allgather(vertices_loc,4,mpi_data_type_T,vertices_rem+4*this->global_dims[0]*this->global_dims[1],4,mpi_data_type_T,communicators[2]);
 
     // exchange local column information with upper neighbor in Z direction (cart grid)
-    MPI_Irecv(vertices_rem                                ,4*this->global_dims[0]*this->global_dims[1],mpi_data_type_T,rank_left, 0,this->global_comm,&rreq);
+    MPI_Irecv(vertices_rem                                            ,4*this->global_dims[0]*this->global_dims[1],mpi_data_type_T,rank_left, 0,this->global_comm,&rreq);
     MPI_Isend(vertices_rem+4*this->global_dims[0]*this->global_dims[1],4*this->global_dims[0]*this->global_dims[1],mpi_data_type_T,rank_right,0,this->global_comm,&sreq);
 
     // determine the offset in ranks
@@ -509,7 +509,7 @@ template<class T, class W> void ALL_Staggered_LB<T,W>::find_neighbors()
     MPI_Barrier(this->global_comm);
 
     // exchange local column information with upper neighbor in Y direction (cart grid)
-    MPI_Irecv(vertices_rem                                ,4*this->global_dims[0]*this->global_dims[1],mpi_data_type_T,rank_right, 0,this->global_comm,&rreq);
+    MPI_Irecv(vertices_rem                                            ,4*this->global_dims[0]*this->global_dims[1],mpi_data_type_T,rank_right, 0,this->global_comm,&rreq);
     MPI_Isend(vertices_rem+4*this->global_dims[0]*this->global_dims[1],4*this->global_dims[0]*this->global_dims[1],mpi_data_type_T,rank_left,0,this->global_comm,&sreq);
 
     // determine the offset in ranks
diff --git a/include/ALL_Tensor.hpp b/include/ALL_Tensor.hpp
index 192de8ebc643b5d15c1a5b7aaabc871917760ae0..039158ac2c98c2b20f0f08925430a9e6d3d9e14f 100644
--- a/include/ALL_Tensor.hpp
+++ b/include/ALL_Tensor.hpp
@@ -155,7 +155,7 @@ template<class T, class W> void ALL_Tensor_LB<T,W>::balance(int)
 {
     int dim = this->get_dimension();
     std::vector<ALL_Point<T>>& vertices = this->get_vertices();
-    std::vector<ALL_Point<T>>& shifted_vertices = this->get_shifted_vertices();
+    std::vector<ALL_Point<T>>& prev_vertices = this->get_vertices();
     
     // loop over all available dimensions
     for (int i = 0; i < dim; ++i)
@@ -185,20 +185,23 @@ template<class T, class W> void ALL_Tensor_LB<T,W>::balance(int)
 
         // collect size in dimension from right neighbor plane
 
-        local_size = vertices.at(1).x(i) - vertices.at(0).x(i);
+        local_size = prev_vertices.at(1).x(i) - prev_vertices.at(0).x(i);
         
         MPI_Irecv(&remote_size,1,mpi_data_type_T,rank_right,0,this->global_comm,&rreq);
         MPI_Isend(&local_size,1,mpi_data_type_T,rank_left,0,this->global_comm,&sreq);
         MPI_Wait(&sreq,&sstat);
         MPI_Wait(&rreq,&rstat);
 
+        // automatic selection of gamma to guarantee stability of method (choice of user gamma ignored)
+        this->gamma = std::max(4.1, 1.1 * (1.0 + std::max(local_size,remote_size) / std::min(local_size,remote_size)));
+
         // calculate shift of borders:
         // s = 0.5 * gamma * (W_r - W_l) / (W_r + W_l) * (d_r + d_l)
         // *_r = * of neighbor (remote)
         // *_l = * of local domain
         T shift = (T)0;
         if (rank_right != MPI_PROC_NULL && this->local_coords.at(i) != this->global_dims.at(i)-1 && !( remote_work == (T)0 && work_local_plane == (T)0) )
-            shift = 1.0 / this->get_gamma() * 0.5 * (remote_work - work_local_plane) / (remote_work + work_local_plane) * (local_size + remote_size);
+            shift = 1.0 / this->gamma * 0.5 * (remote_work - work_local_plane) / (remote_work + work_local_plane) * (local_size + remote_size);
 
         // determine a maximum move in order to avoid overlapping borders and to guarantee
         // a stable shift of the domain borders (avoid the loss of a domain due to too large
@@ -224,15 +227,15 @@ template<class T, class W> void ALL_Tensor_LB<T,W>::balance(int)
 
         // if a left neighbor exists: shift left border 
         if (rank_left != MPI_PROC_NULL && this->local_coords[i] != 0)
-            this->shifted_vertices.at(0).set_coordinate(i, vertices.at(0).x(i) + remote_shift);
+            this->vertices.at(0).set_coordinate(i, prev_vertices.at(0).x(i) + remote_shift);
         else
-            this->shifted_vertices.at(0).set_coordinate(i, vertices.at(0).x(i));
+            this->vertices.at(0).set_coordinate(i, prev_vertices.at(0).x(i));
 
         // if a right neighbor exists: shift right border
         if (rank_right != MPI_PROC_NULL && this->local_coords[i] != this->global_dims[i]-1)
-            this->shifted_vertices.at(1).set_coordinate(i, vertices.at(1).x(i) + shift);
+            this->vertices.at(1).set_coordinate(i, prev_vertices.at(1).x(i) + shift);
         else
-            this->shifted_vertices.at(1).set_coordinate(i, vertices.at(1).x(i));
+            this->vertices.at(1).set_coordinate(i, prev_vertices.at(1).x(i));
     }
 }
 
diff --git a/include/ALL_Unstructured.hpp b/include/ALL_Unstructured.hpp
index 76cb14a327bdcababa26b1d08451d8ccb5a81fe7..6cde2fb71e875ca34d66b03acf3f55bd50835682 100644
--- a/include/ALL_Unstructured.hpp
+++ b/include/ALL_Unstructured.hpp
@@ -482,6 +482,9 @@ template <class T, class W> void ALL_Unstructured_LB<T,W>::setup()
 
 template <class T, class W> void ALL_Unstructured_LB<T,W>::balance(int /*step*/)
 {
+
+    this->prev_vertices = this->vertices;
+
     // geometrical center and work of each neighbor for each vertex
     // work is cast to vertex data type, since in the end a shift
     // of the vertex is computed
@@ -505,7 +508,7 @@ template <class T, class W> void ALL_Unstructured_LB<T,W>::balance(int /*step*/)
     {
         for (int d = 0; d < this->dimension; ++d)
         {
-            center[d] += ( (this->vertices.at(v)).x(d) / (T)n_vertices );
+            center[d] += ( (this->prev_vertices.at(v)).x(d) / (T)n_vertices );
         }
         local_info[v * (this->dimension+2) + this->dimension] = (T)this->work.at(0);
     }
@@ -513,7 +516,7 @@ template <class T, class W> void ALL_Unstructured_LB<T,W>::balance(int /*step*/)
     {
         for (int d = 0; d < this->dimension; ++d)
         {
-            local_info[v * (this->dimension+2) + d] = center[d] - (this->vertices.at(v)).x(d);
+            local_info[v * (this->dimension+2) + d] = center[d] - (this->prev_vertices.at(v)).x(d);
         }
     }
 
@@ -588,28 +591,28 @@ template <class T, class W> void ALL_Unstructured_LB<T,W>::balance(int /*step*/)
     // compute for each vertex the maximum movement towards the center
 
     local_info[ 0 * (this->dimension+2) + this->dimension+1 ] =
-        distance( &(this->vertices), 0, 1, 2, 4, this->local_rank);
+        distance( &(this->prev_vertices), 0, 1, 2, 4, this->local_rank);
 
     local_info[ 1 * (this->dimension+2) + this->dimension+1 ] =
-        distance( &(this->vertices), 1, 0, 3, 5, this->local_rank);
+        distance( &(this->prev_vertices), 1, 0, 3, 5, this->local_rank);
 
     local_info[ 2 * (this->dimension+2) + this->dimension+1 ] =
-        distance( &(this->vertices), 2, 3, 0, 6, this->local_rank);
+        distance( &(this->prev_vertices), 2, 3, 0, 6, this->local_rank);
 
     local_info[ 3 * (this->dimension+2) + this->dimension+1 ] =
-        distance( &(this->vertices), 3, 2, 1, 7, this->local_rank);
+        distance( &(this->prev_vertices), 3, 2, 1, 7, this->local_rank);
 
     local_info[ 4 * (this->dimension+2) + this->dimension+1 ] =
-        distance( &(this->vertices), 4, 5, 6, 0, this->local_rank);
+        distance( &(this->prev_vertices), 4, 5, 6, 0, this->local_rank);
 
     local_info[ 5 * (this->dimension+2) + this->dimension+1 ] =
-        distance( &(this->vertices), 5, 4, 7, 1, this->local_rank);
+        distance( &(this->prev_vertices), 5, 4, 7, 1, this->local_rank);
 
     local_info[ 6 * (this->dimension+2) + this->dimension+1 ] =
-        distance( &(this->vertices), 6, 7, 4, 2, this->local_rank);
+        distance( &(this->prev_vertices), 6, 7, 4, 2, this->local_rank);
 
     local_info[ 7 * (this->dimension+2) + this->dimension+1 ] =
-        distance( &(this->vertices), 7, 6, 5, 3, this->local_rank);
+        distance( &(this->prev_vertices), 7, 6, 5, 3, this->local_rank);
 
     /*
     if (this->local_rank == 4)
@@ -719,8 +722,8 @@ template <class T, class W> void ALL_Unstructured_LB<T,W>::balance(int /*step*/)
     }
     MPI_Waitall(n_vertices,request,status);
 
-    this->shifted_vertices = this->vertices;
     /*
+    this->vertices = this->prev_vertices;
     if (this->local_rank == 7)
     {
         std::cout << "center: ";
@@ -736,8 +739,8 @@ template <class T, class W> void ALL_Unstructured_LB<T,W>::balance(int /*step*/)
         //if (this->local_rank == 4) std::cout << "vertex " << v << ": ";
         for (int d = 0; d < this->dimension; ++d)
         {
-            this->shifted_vertices.at(v).set_coordinate(d, this->vertices.at(v).x(d) + shift_vectors[v * this->dimension + d]);
-            if (this->shifted_vertices.at(v).x(d) != this->shifted_vertices.at(v).x(d))
+            this->vertices.at(v).set_coordinate(d, this->prev_vertices.at(v).x(d) + shift_vectors[v * this->dimension + d]);
+            if (this->vertices.at(v).x(d) != this->vertices.at(v).x(d))
             {
                 std::cout << this->local_rank << " " << v << " " << d << std::endl;
                 MPI_Abort(MPI_COMM_WORLD,-200);
diff --git a/include/ALL_Voronoi.hpp b/include/ALL_Voronoi.hpp
index 6c4b742e6870577c4c501ad94cf168856dcf4acf..7ca1df40ddb72f85ec66f6caca9abf0acda993f8 100644
--- a/include/ALL_Voronoi.hpp
+++ b/include/ALL_Voronoi.hpp
@@ -157,14 +157,14 @@ template<class T, class W> void ALL_Voronoi_LB<T,W>::balance(int loadbalancing_s
     // collect local information in array
     int info_length = 2*dimension+1;
     std::vector<T> local_info(info_length);
+    std::vector<ALL_Point<T>>& prev_vertices = this->get_vertices();
     std::vector<ALL_Point<T>>& vertices = this->get_vertices();
-    std::vector<ALL_Point<T>>& shifted_vertices = this->get_shifted_vertices();
     std::vector<W>& work = this->get_work();
     
     for (auto i = 0; i < 2; ++i)
         for (auto d = 0; d < dimension; ++d)
         {
-            local_info.at(d) = vertices.at(i).x(d);
+            local_info.at(d) = prev_vertices.at(i).x(d);
         }
     local_info.at(2*dimension) = (T)work.at(0);
 
@@ -547,7 +547,7 @@ template<class T, class W> void ALL_Voronoi_LB<T,W>::balance(int loadbalancing_s
 #endif
 
     for (int d = 0; d < dimension; ++d)
-       shifted_vertices[d] = local_info.at(d);
+       vertices[d] = local_info.at(d);
 
     // compute new neighboring cells and generator points
 
diff --git a/src/ALL_fortran.cpp b/src/ALL_fortran.cpp
index 6f20f0c959af7088b0aef6272253b0af47bc369f..1199bc8edc352d005dabbc51aff65391a0a05d72 100644
--- a/src/ALL_fortran.cpp
+++ b/src/ALL_fortran.cpp
@@ -92,7 +92,7 @@ extern "C"
     // wrapper to get number of new vertices
     void ALL_get_new_number_of_vertices_f(ALL<double,double>* all_obj, int* n_vertices)
     {
-        *n_vertices = (all_obj->get_result_vertices()).size(); 
+        *n_vertices = (all_obj->get_vertices()).size(); 
     }
 
     // wrapper to return new vertices
@@ -111,16 +111,16 @@ extern "C"
     }
 
     // wrapper to return new vertices
-    void ALL_get_new_vertices_f(ALL<double,double>* all_obj, int* n_vertices, double* new_vertices)
+    void ALL_get_prev_vertices_f(ALL<double,double>* all_obj, int* n_vertices, double* prev_vertices)
     {
-        std::vector<ALL_Point<double>> tmp_vertices = all_obj->get_result_vertices();
+        std::vector<ALL_Point<double>> tmp_vertices = all_obj->get_prev_vertices();
         int dimension = all_obj->get_dimension();
         *n_vertices = tmp_vertices.size();
         for (int i = 0; i < *n_vertices; ++i)
         {
             for (int j = 0; j < dimension; ++j)
             {
-                new_vertices[i * dimension + j] = tmp_vertices.at(i).x(j);
+                prev_vertices[i * dimension + j] = tmp_vertices.at(i).x(j);
             }
         }
     }
diff --git a/src/ALL_module.f90 b/src/ALL_module.f90
index e6f5d49d866e2f2539251850be87eef04666d9c2..d7ce908818bcbbcfeee13b6b6cd700eacce2762a 100644
--- a/src/ALL_module.f90
+++ b/src/ALL_module.f90
@@ -102,8 +102,8 @@ MODULE ALL_module
             INTEGER(c_int)                  ::  n
             REAL(c_double)                  ::  vertices(*) 
         END SUBROUTINE
-        SUBROUTINE ALL_get_new_vertices_int(obj,n,vertices) &
-            BIND(C,NAME="ALL_get_new_vertices_f")
+        SUBROUTINE ALL_get_prev_vertices_int(obj,n,vertices) &
+            BIND(C,NAME="ALL_get_prev_vertices_f")
             USE ISO_C_BINDING
             TYPE(c_ptr), VALUE              ::  obj
             INTEGER(c_int)                  ::  n
@@ -121,7 +121,7 @@ MODULE ALL_module
             TYPE(ALL_t)     :: obj
             INTEGER         :: method
             INTEGER         :: dim
-            REAL(8)         :: gamma
+            REAL(c_double)  :: gamma
             obj%object = ALL_init_int(INT(method,c_int),INT(dim,c_int),REAL(gamma,c_double))
         END SUBROUTINE
         SUBROUTINE ALL_set_proc_grid_params(obj, loc, ranks)
@@ -131,21 +131,21 @@ MODULE ALL_module
             call ALL_set_proc_grid_params_int(obj%object, size(loc,1), loc, size(ranks,1), ranks)
         END SUBROUTINE
         SUBROUTINE ALL_set_min_domain_size(obj, dim, domain_size)
-            INTEGER                  ::  dim
-            REAL(8), DIMENSION(dim)  ::  domain_size
-            TYPE(ALL_t)              ::  obj
+            INTEGER                         ::  dim
+            REAL(c_double), DIMENSION(dim)  ::  domain_size
+            TYPE(ALL_t)                     ::  obj
             CALL ALL_set_min_domain_size_int(obj%object, dim, domain_size)
         END SUBROUTINE
         SUBROUTINE ALL_set_work(obj, work)
             TYPE(ALL_t)     ::  obj
-            REAL(8)         ::  work
+            REAL(c_double)  ::  work
             CALL ALL_set_work_int(obj%object, work)
         END SUBROUTINE
         SUBROUTINE ALL_set_vertices(obj, n, dim, vertices)
             TYPE(ALL_t)     ::  obj
             INTEGER         ::  n
             INTEGER         ::  dim
-            REAL(8)         ::  vertices(n*dim)
+            REAL(c_double)  ::  vertices(n*dim)
             CALL ALL_set_vertices_int(obj%object, n, dim, vertices)
         END SUBROUTINE
         SUBROUTINE ALL_set_communicator(obj, comm)
@@ -171,14 +171,14 @@ MODULE ALL_module
         SUBROUTINE ALL_get_vertices(obj,n,vertices)
             TYPE(ALL_t)         ::  obj
             INTEGER             ::  n
-            REAL(8)             ::  vertices(*)
+            REAL(c_double)      ::  vertices(*)
             CALL ALL_get_vertices_int(obj%object,n,vertices)
         END SUBROUTINE
-        SUBROUTINE ALL_get_new_vertices(obj,n,vertices)
+        SUBROUTINE ALL_get_prev_vertices(obj,n,vertices)
             TYPE(ALL_t)         ::  obj
             INTEGER             ::  n
-            REAL(8)             ::  vertices(*)
-            CALL ALL_get_new_vertices_int(obj%object,n,vertices)
+            REAL(c_double)      ::  vertices(*)
+            CALL ALL_get_prev_vertices_int(obj%object,n,vertices)
         END SUBROUTINE
         SUBROUTINE ALL_print_vtk_outlines(obj,step)
             TYPE(ALL_t)    ::  obj