/**********************************************************************************
 * Copyright 2010 Christoph Pflaum 
 * 		Department Informatik Lehrstuhl 10 - Systemsimulation
 *		Friedrich-Alexander Universität Erlangen-Nürnberg
 * 
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 * http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 **********************************************************************************/
 
// ------------------------------------------------------------
//
// variable_cc.h
//
// ------------------------------------------------------------

#ifndef VARIABLE_CC_H
#define VARIABLE_CC_H

#ifdef _OPENMP
#include <omp.h>
#endif

#include "../grid/compose_grid.h"



template <class DTyp>
Variable<DTyp>::Variable(Blockgrid& blockgrid_ ) :
	Object_based_on_ug ( blockgrid_.Give_unstructured_grid()) {
	int num_total, num_neighbor;
	blockgrid = &blockgrid_;
	ug = blockgrid->Give_unstructured_grid();

	specialVariableForOneBlock = false;

	// Kopie der UG Daten zum Löschen
	ug_number_hexahedra   = ug->Give_number_hexahedra();
	ug_number_quadrangles = ug->Give_number_quadrangles();
        ug_number_edges       = ug->Give_number_edges();
	ug_number_points      = ug->Give_number_points();
        ug_edge_number_neighbors = new int[ug_number_edges];
	idGrid = blockgrid->getId();
	
	for(int i=0;i<ug_number_edges;++i)
	    ug_edge_number_neighbors[i] = ug->Give_edge(i)->Give_number_neighbors();
	
	// for parallel
	my_rank = ug->Give_my_rank();
	comm    = ug->Give_MPI_comm();

	// hexahedra
	data_hexahedra = new DTyp*[ug->Give_number_hexahedra() ];
	for ( int id=0; id<ug->Give_number_hexahedra(); ++id ) {
		num_total = blockgrid->Give_N_total_hexahedron ( id );

		if ( ug->Give_hexahedron ( id )->my_object ( my_rank ) ) {
			data_hexahedra[id] = new DTyp[num_total];

			{

#pragma omp parallel for num_threads(UGBlocks::numThreadsToTake)  if(UGBlocks::useOpenMP)
				for ( int i=0; i<num_total; ++i )
					data_hexahedra[id][i] = ( DTyp ) 0;
			}
		} else data_hexahedra[id] = NULL;
	}

	// quadrangles
	data_quadrangles = new DTyp*[ug->Give_number_quadrangles() ];
	for ( int id=0; id<ug->Give_number_quadrangles(); ++id ) {
		num_total = blockgrid->Give_N_total_quadrangle ( id );

		if ( ug->Give_quadrangle ( id )->my_object ( my_rank ) ) {
			data_quadrangles[id] = new DTyp[num_total];

			for ( int i=0; i<num_total; ++i )
				data_quadrangles[id][i] = ( DTyp ) 0;
		} else {
			data_quadrangles[id] = NULL;
		}
	}

	data_interior_quadrangles = new DTyp*[ug->Give_number_quadrangles() ];
	for ( int id=0; id<ug->Give_number_quadrangles(); ++id ) {
		num_total = blockgrid->Give_N_total_quadrangle ( id );

		if ( ug->Give_quadrangle ( id )->my_object ( my_rank ) ) {
			data_interior_quadrangles[id] = new DTyp[num_total];

			for ( int i=0; i<num_total; ++i )
				data_interior_quadrangles[id][i] = ( DTyp ) 0;
		} else data_interior_quadrangles[id] = NULL;
	}

	data_exterior_quadrangles = new DTyp*[ug->Give_number_quadrangles() ];
	for ( int id=0; id<ug->Give_number_quadrangles(); ++id ) {
		if ( ug->Give_quadrangle ( id )->Give_exists_exterior() &&
		     ug->Give_quadrangle ( id )->my_object ( my_rank ) )        {
			num_total = blockgrid->Give_N_total_quadrangle ( id );
			data_exterior_quadrangles[id] = new DTyp[num_total];

			for ( int i=0; i<num_total; ++i )
				data_exterior_quadrangles[id][i] = ( DTyp ) 0;
		} else data_exterior_quadrangles[id] = NULL;
	}

	// edges
	data_edges = new DTyp*[ug->Give_number_edges() ];
	for ( int id=0; id<ug->Give_number_edges(); ++id ) {
		num_total = blockgrid->Give_Nx_edge ( id ) +1;
		if ( ug->Give_edge ( id )->my_object ( my_rank ) ) {
			data_edges[id] = new DTyp[num_total];

			for ( int i=0; i<num_total; ++i )
				data_edges[id][i] = ( DTyp ) 0;
		} else data_edges[id] = NULL;
	}


	data_neighbor_edges = new DTyp**[ug->Give_number_edges() ];
	for ( int id=0; id<ug->Give_number_edges(); ++id ) {
		num_neighbor = ug->Give_edge ( id )->Give_number_neighbors();

		if ( ug->Give_edge ( id )->my_object ( my_rank ) ) {
			data_neighbor_edges[id] = new DTyp*[num_neighbor];
			for ( int i=0; i<num_neighbor; ++i ) {
				num_total = blockgrid->Give_Nx_edge ( id ) +1;
				data_neighbor_edges[id][i] = new DTyp[num_total];

				for ( int j=0; j<num_total; ++j )
					data_neighbor_edges[id][i][j] = ( DTyp ) 0;
			}
		} else data_neighbor_edges[id] = NULL;
	}

	// points
	data_points = new DTyp*[ug->Give_number_points() ];
	for ( int id=0; id<ug->Give_number_points(); ++id ) {
		num_total = 1;

		if ( ug->Give_point ( id )->my_object ( my_rank ) ) {
			data_points[id] = new DTyp[num_total];

			for ( int i=0; i<num_total; ++i )
				data_points[id][i] = ( DTyp ) 0;
		} else data_points[id] = NULL;
	}
	data_neighbor_points = new DTyp*[ug->Give_number_points() ];
	for ( int id=0; id<ug->Give_number_points(); ++id ) {
		num_neighbor = ug->Give_point ( id )->Give_number_neighbors();

		if ( ug->Give_point ( id )->my_object ( my_rank ) ) {
			data_neighbor_points[id] = new DTyp[num_neighbor];

			for ( int i=0; i<num_neighbor; ++i )
				data_neighbor_points[id][i] = ( DTyp ) 0;
		} else data_neighbor_points[id] = NULL;
	}
}

template <class DTyp>
Variable<DTyp>::~Variable()
{
	if ( ug!=NULL ) {
		Unregister ( ug );	
		Delete_data();
	}
}


template <class DTyp>
void Variable<DTyp>::Delete_data() {
	int num_neighbor;
	if ( ug==NULL ) return;
	ug=NULL;

	if ( data_hexahedra != NULL && (specialVariableForOneBlock==false)) {
		for ( int i=0; i< ug_number_hexahedra; ++i )
			if ( data_hexahedra[i] != NULL )
				delete [] data_hexahedra[i];
		delete [] data_hexahedra;
	}
	data_hexahedra = NULL;

	if ( data_quadrangles != NULL ) {
		for ( int i=0; i< ug_number_quadrangles; ++i )
			if ( data_quadrangles[i] != NULL )
				delete [] data_quadrangles[i];
		delete [] data_quadrangles;
	}
	data_quadrangles = NULL;

	if ( data_interior_quadrangles != NULL ) {
		for ( int i=0; i<ug_number_quadrangles; ++i )
			if ( data_interior_quadrangles[i] != NULL )
				delete [] data_interior_quadrangles[i];
		delete [] data_interior_quadrangles;
	}
	data_interior_quadrangles = NULL;

	if ( data_exterior_quadrangles != NULL ) {
		for ( int i=0; i<ug_number_quadrangles; ++i )
			if ( data_exterior_quadrangles[i] != NULL )
				delete [] data_exterior_quadrangles[i];
		delete [] data_exterior_quadrangles;
	}
	data_exterior_quadrangles = NULL;

	if ( data_edges != NULL ) {
		for ( int i=0; i<ug_number_edges; ++i )
			if ( data_edges[i] != NULL )
				delete [] data_edges[i];
		delete [] data_edges;
	}
	data_edges = NULL;

	if ( data_neighbor_edges != NULL ) {
		for ( int i=0; i<ug_number_edges; ++i ) {
			if ( data_neighbor_edges[i] != NULL ) {
				num_neighbor = ug_edge_number_neighbors[i];
				for ( int j=0; j<num_neighbor; ++j )
					if ( data_neighbor_edges[i][j] != NULL )
						delete [] data_neighbor_edges[i][j];
			}
			delete [] data_neighbor_edges[i];
		}
		delete [] data_neighbor_edges;
	}
	data_neighbor_edges = NULL;

	if ( data_points != NULL ) {
		for ( int i=0; i<ug_number_points; ++i )
			if ( data_points[i] != NULL )
				delete [] data_points[i];
		delete [] data_points;
	}
	data_points = NULL;

	if ( data_neighbor_points != NULL ) {
		for ( int i=0; i<ug_number_points; ++i )
			if ( data_neighbor_points[i] != NULL )
				delete [] data_neighbor_points[i];
		delete [] data_neighbor_points;
	}
	data_neighbor_points = NULL;

	delete [] ug_edge_number_neighbors;
}

template <class DTyp>
void Variable<DTyp>::operator= ( DTyp x )
{
	int num_total, i;

	/*
	cout << " operator=x: my_rank: "
	 << my_rank << endl;
	*/

	// hexahedra
	for ( int id=0; id<ug->Give_number_hexahedra(); ++id ) {
		if ( ug->Give_hexahedron ( id )->my_object ( my_rank ) ) {
			num_total = blockgrid->Give_N_total_hexahedron ( id );

#pragma omp parallel for num_threads(UGBlocks::numThreadsToTake) if(UGBlocks::useOpenMP)
			for ( i=0; i<num_total; ++i )
				data_hexahedra[id][i] = x;
		}
	}

	// quadrangles
	for ( int id=0; id<ug->Give_number_quadrangles(); ++id ) {
		if ( ug->Give_quadrangle ( id )->my_object ( my_rank ) ) {
			num_total = blockgrid->Give_N_total_quadrangle ( id );

			for ( i=0; i<num_total; ++i )
				data_quadrangles[id][i] = x;
		}
	}

	// edges
	for ( int id=0; id<ug->Give_number_edges(); ++id ) {
		if ( ug->Give_edge ( id )->my_object ( my_rank ) ) {
			num_total = blockgrid->Give_Nx_edge ( id ) +1;

			for ( i=0; i<num_total; ++i )
				data_edges[id][i] = x;
		}
	}

	// points
	for ( int id=0; id<ug->Give_number_points(); ++id ) {
		if ( ug->Give_point ( id )->my_object ( my_rank ) ) {
			num_total = 1;

			for ( i=0; i<num_total; ++i )
				data_points[id][i] = x;
		}
	}
}

template <class DTyp>
void Variable<DTyp>::operator= ( const DTyp_Restriction<DTyp>& a )
{
	int num_total, i;
	DTyp   x      = a.Give_x();

	// hexahedra
	for ( int id=0; id<ug->Give_number_hexahedra(); ++id ) {
		if ( ug->Give_hexahedron ( id )->my_object ( my_rank ) &&
		     a.template Give_marker<hexahedronEl> ( id ) ==yes_mark ) {
			num_total = blockgrid->Give_N_total_hexahedron ( id );

#pragma omp parallel for num_threads(UGBlocks::numThreadsToTake) if(UGBlocks::useOpenMP)
			for ( i=0; i<num_total; ++i )
				data_hexahedra[id][i] = x;
		}
	}


	// quadrangles
	for ( int id=0; id<ug->Give_number_quadrangles(); ++id ) {
		if ( ug->Give_quadrangle ( id )->my_object ( my_rank ) &&
		     a.template Give_marker<quadrangleEl> ( id ) ==yes_mark ) {
			num_total = blockgrid->Give_N_total_quadrangle ( id );

			for ( i=0; i<num_total; ++i )
				data_quadrangles[id][i] = x;
		}
	}

	// edges
	for ( int id=0; id<ug->Give_number_edges(); ++id ) {
		if ( ug->Give_edge ( id )->my_object ( my_rank ) &&
		     a.template Give_marker<edgeEl> ( id ) ==yes_mark ) {
			num_total = blockgrid->Give_Nx_edge ( id ) +1;

			for ( i=0; i<num_total; ++i )
				data_edges[id][i] = x;
		}
	}

	// points
	for ( int id=0; id<ug->Give_number_points(); ++id ) {
		if ( ug->Give_point ( id )->my_object ( my_rank ) &&
		     a.template Give_marker<pointEl> ( id ) ==yes_mark ) {
			num_total = 1;

			for ( i=0; i<num_total; ++i )
				data_points[id][i] = x;
		}
	}
}

//
//
//template <class DTyp>
//void Variable<DTyp>::operator= ( const Expr_interpolant_cell_to_point<DTyp>& ao ) {
//  int Nx, Ny, Nz, i, j, k;
//
//  stencil_typ sten_typ = ao.Give_stencil_typ();
////  control_typ cont_typ = ao.Give_control_typ();
//
//  // hexahedra
//
//  for ( int id = 0;id < ug->Give_number_hexahedra();++id ) {
//          if ( sten_typ == yes_stencil ) {
//                  ao.Update<hexahedronEl> ( id );
//          }
//
//          if ( ug->Give_hexahedron ( id )->my_object ( my_rank ) ) {
//
//                  Nx = blockgrid->Give_Nx_hexahedron ( id );
//                  Ny = blockgrid->Give_Ny_hexahedron ( id );
//                  Nz = blockgrid->Give_Nz_hexahedron ( id );
//
//                  if ( cont_typ == thread_save ) {
//#pragma omp parallel for private(j,i) num_threads(UGBlocks::numThreadsToTake)  if(UGBlocks::useOpenMP)
//                          for ( k = 1;k < Nz;++k ) {
//                                  for ( j = 1;j < Ny;++j )
//                                          for ( i = 1;i < Nx;++i )
//                                                  data_hexahedra[id][i+ ( Nx+1 ) * ( j+ ( Ny+1 ) *k ) ] =
//                                                      ao.Give_data<hexahedronEl> ( id, i, j, k, Nx, Ny );
//                          }
//                  }
//                  else {
//                          for ( k = 1;k < Nz;++k ) {
//                                  for ( j = 1;j < Ny;++j )
//                                          for ( i = 1;i < Nx;++i )
//                                                  data_hexahedra[id][i+ ( Nx+1 ) * ( j+ ( Ny+1 ) *k ) ] =
//                                                      ao.Give_data<hexahedronEl> ( id, i, j, k, Nx, Ny );
//                          }
//                  }
//          }
//  }
//
//  // quadrangles
//  for ( int id = 0;id < ug->Give_number_quadrangles();++id ) {
//          if ( sten_typ == yes_stencil )
//                  ao.Update<quadrangleEl> ( id );
//
//          if ( ug->Give_quadrangle ( id )->my_object ( my_rank ) ) {
//
//                  Nx = blockgrid->Give_Nx_quadrangle ( id );
//                  Ny = blockgrid->Give_Ny_quadrangle ( id );
//
//                  for ( j = 1;j < Ny;++j )
//                          for ( i = 1;i < Nx;++i ) {
//                                  data_quadrangles[id][i+ ( Nx+1 ) *j] =
//                                      ao.Give_data<quadrangleEl> ( id, i, j, 0, Nx, 0 );
//                          }
//          }
//  }
//
//  // edges
//  for ( int id = 0;id < ug->Give_number_edges();++id ) {
//          if ( sten_typ == yes_stencil )
//                  ao.Update<edgeEl> ( id );
//
//          if ( ug->Give_edge ( id )->my_object ( my_rank ) ) {
//                  Nx = blockgrid->Give_Nx_edge ( id );
//
//                  for ( i = 1;i < Nx;++i )
//                          data_edges[id][i] =
//                              ao.Give_data<edgeEl> ( id, i, 0, 0, 0, 0 );
//          }
//  }
//
//  // points
//  for ( int id = 0;id < ug->Give_number_points();++id ) {
//          if ( sten_typ == yes_stencil )
//                  ao.Update<pointEl> ( id );
//
//          if ( ug->Give_point ( id )->my_object ( my_rank ) ) {
//                  data_points[id][0] =
//                      ao.Give_data<pointEl> ( id, 0, 0, 0, 0, 0 );
//          }
//  }
//}


template <class DTyp>
void Variable<DTyp>::operator= ( const Variable<DTyp>& v ) {
  if(blockgrid->getId() != v.blockgrid->getId()) {
        assert(ug->isComposeGrid());
        ComposeUg* compUg = static_cast<ComposeUg*>(ug);

        int thisIdGrid = v.blockgrid->getId();

	// hexahedra
        int startHex = compUg->getStartHex(thisIdGrid);
        int endHex   = compUg->getEndHex(thisIdGrid);	
        for(int id=startHex;id<endHex;++id) {
            int idRight = id - startHex;	
	    if(ug->Give_hexahedron ( id )->my_object ( my_rank ) ) {
	       int num_total = blockgrid->Give_N_total_hexahedron ( id );

#pragma omp parallel for num_threads(UGBlocks::numThreadsToTake) if(UGBlocks::useOpenMP)
	        for(int i=0; i<num_total; ++i )
		    data_hexahedra[id][i] = v.data_hexahedra[idRight][i];
	    }
	}

	// quadrangles
        int startQuad = compUg->getStartQuad(thisIdGrid);
        int endQuad   = compUg->getEndQuad(thisIdGrid);	
        for(int id=startQuad;id<endQuad;++id) {
            int idRight = id - startQuad;	
	    if(ug->Give_quadrangle ( id )->my_object ( my_rank ) ) {
	       int num_total = blockgrid->Give_N_total_quadrangle ( id );

	       for(int i=0; i<num_total; ++i )
		   data_quadrangles[id][i] = v.data_quadrangles[idRight][i];
	    }
	}

	// edges
        int startEdge = compUg->getStartEdge(thisIdGrid);
        int endEdge   = compUg->getEndEdge(thisIdGrid);	
        for(int id=startEdge;id<endEdge;++id) {
            int idRight = id - startEdge;	
	    if(ug->Give_edge ( id )->my_object ( my_rank ) ) {
	       int num_total = blockgrid->Give_Nx_edge ( id ) +1;

	       for(int i=0; i<num_total; ++i )
		   data_edges[id][i] = v.data_edges[idRight][i];
	    }
	}

	// points
        int startPoint = compUg->getStartPoint(thisIdGrid);
        int endPoint   = compUg->getEndPoint(thisIdGrid);	
        for(int id=startPoint;id<endPoint;++id) {
            int idRight = id - startPoint;
	    if(ug->Give_point ( id )->my_object ( my_rank ) ) {
	       int num_total = 1;

	       for(int i=0; i<num_total; ++i )
		   data_points[id][i] = v.data_points[idRight][i];
	    }
	}
        return;
  }
     
  // hexahedra
  for(int id=0; id<ug->Give_number_hexahedra(); ++id ) {
      if(ug->Give_hexahedron ( id )->my_object ( my_rank ) ) {
	 int num_total = blockgrid->Give_N_total_hexahedron ( id );

#pragma omp parallel for num_threads(UGBlocks::numThreadsToTake) if(UGBlocks::useOpenMP)
			for (int i=0; i<num_total; ++i )
				data_hexahedra[id][i] = v.data_hexahedra[id][i];

      }
  }

  // quadrangles
  for(int id=0; id<ug->Give_number_quadrangles(); ++id ) {
      if(ug->Give_quadrangle ( id )->my_object ( my_rank ) ) {
	 int num_total = blockgrid->Give_N_total_quadrangle ( id );
	 for(int i=0; i<num_total; ++i )
	     data_quadrangles[id][i] = v.data_quadrangles[id][i];
      }
  }

  // edges
  for(int id=0; id<ug->Give_number_edges(); ++id ) {
      if(ug->Give_edge ( id )->my_object ( my_rank ) ) {
	 int num_total = blockgrid->Give_Nx_edge ( id ) +1;

      for(int i=0; i<num_total; ++i )
	  data_edges[id][i] = v.data_edges[id][i];
      }
  }

  // points
  for(int id=0; id<ug->Give_number_points(); ++id ) {
      if(ug->Give_point ( id )->my_object ( my_rank ) ) {
         int num_total = 1;

	 for(int i=0; i<num_total; ++i )
	     data_points[id][i] = v.data_points[id][i];
      }
  }

//        for ( int id = 0;id < ug->Give_number_hexahedra();++id )
//          Update<hexahedronEl> ( id );
//        for ( int id = 0;id < ug->Give_number_quadrangles();++id )
//          Update<quadrangleEl> ( id );
//        for ( int id = 0;id < ug->Give_number_edges();++id )
//          Update<edgeEl> ( id );
//        for ( int id = 0;id < ug->Give_number_points();++id )
//          Update<pointEl> ( id );
}

/***
 * Zur Auswertung einer komplizierten Formel in einem System
 *         Assign_System::Compute(outputVector,functor,inputVector)
 *  Zwei versionen seriell und parallel      
 **/
class Assign_System {
public:  
 /***
 *  Klasse Func folgende Methoden haben:
 *        void Func::compute(DTyp* inputData);
 *        DTyp Func::getResult(int idata);
 **/
  template <class DTyp, class Func>
  static void Compute(const std::vector<Variable<DTyp>*>& outputVector,
	       Func& functor,
	       const std::vector<Variable<DTyp>*>& inputVector) {
    	int num_total, i;

	// input output Daten
	int sizeInputData  = inputVector.size();
	double* inputData  = new double[sizeInputData];
	
	//    inputDataP[s]  = new double[sizeInputData];

	int sizeOutputData = outputVector.size();
		
	// Kopie der IterationsDaten
	Unstructured_grid* ug = outputVector.at(0)->Give_Ug();
	Blockgrid*  blockgrid = outputVector.at(0)->Give_blockgrid();
	int my_rank           = ug->Give_my_rank();;
	
	// hexahedra
	for ( int id=0; id<ug->Give_number_hexahedra(); ++id ) {
		if ( ug->Give_hexahedron ( id )->my_object ( my_rank ) ) {
			num_total = blockgrid->Give_N_total_hexahedron ( id );
			for(i=0; i<num_total; ++i ) {
			    for(int idata=0;idata<sizeInputData;++idata) {
			        inputData[idata] = inputVector.at(idata)->Give_pointer_data_hexahedra()[id][i];
			    }
			    functor.compute(inputData);
			    for(int idata=0;idata<sizeOutputData;++idata) {
				outputVector.at(idata)->Give_pointer_data_hexahedra()[id][i] = functor.getResult(idata);
			    }
			}
		}
	}

	
	
	// quadrangles
	for(int id=0; id<ug->Give_number_quadrangles(); ++id ) {
		if ( ug->Give_quadrangle ( id )->my_object ( my_rank ) ) {
			num_total = blockgrid->Give_N_total_quadrangle ( id );
			for(i=0; i<num_total; ++i ) {
			    for(int idata=0;idata<sizeInputData;++idata) {
			        inputData[idata] = inputVector.at(idata)->Give_pointer_data_quadrangle()[id][i];
			    }
			    functor.compute(inputData);
			    for(int idata=0;idata<sizeOutputData;++idata) {
				outputVector.at(idata)->Give_pointer_data_quadrangle()[id][i] 
				                                    = functor.getResult(idata);
			    }
			}
		}
	}

	// edges
	for ( int id=0; id<ug->Give_number_edges(); ++id ) {
		if ( ug->Give_edge ( id )->my_object ( my_rank ) ) {
			num_total = blockgrid->Give_Nx_edge ( id ) +1;
			for(i=0; i<num_total; ++i ) {
			    for(int idata=0;idata<sizeInputData;++idata) {
			        inputData[idata] = inputVector.at(idata)->Give_pointer_data_edges()[id][i];
			    }
			    functor.compute(inputData);
			    for(int idata=0;idata<sizeOutputData;++idata) {
				outputVector.at(idata)->Give_pointer_data_edges()[id][i] 
				                                    = functor.getResult(idata);
			    }
			}
		}
	}

	// points
	for ( int id=0; id<ug->Give_number_points(); ++id ) {
		if ( ug->Give_point ( id )->my_object ( my_rank ) ) {
			num_total = 1;
			for ( i=0; i<num_total; ++i ) {
			    for(int idata=0;idata<sizeInputData;++idata) {
			        inputData[idata] = inputVector.at(idata)->Give_pointer_data_points()[id][i];
			    }
			    functor.compute(inputData);
			    for(int idata=0;idata<sizeOutputData;++idata) {
				outputVector.at(idata)->Give_pointer_data_points()[id][i] 
				                                    = functor.getResult(idata);
			    }				
			}
		}
	}
    
    delete[] inputData;
  }
  
 /***
 *  Klasse Func folgende Methoden haben:
 *        void Func::compute(DTyp* inputData, int tid);
 *        DTyp Func::getResult(int idata, int tid);
 *  Hierbei ist tid die Nummer der OpenMp Prozesses
 **/  
  template <class DTyp, class Func> 
  static void ComputeParallel(const std::vector<Variable<DTyp>*>& outputVector,
	       Func& functor,
	       const std::vector<Variable<DTyp>*>& inputVector) {
    	int num_total, i;

	// input output Daten
	int sizeInputData  = inputVector.size();
	double* inputData  = new double[sizeInputData];
	
	double** inputDataP  = new double*[UGBlocks::numThreadsToTake];
	for(int s=0;s<UGBlocks::numThreadsToTake;++s)
	    inputDataP[s]  = new double[sizeInputData];

	int sizeOutputData = outputVector.size();
		
	// Kopie der IterationsDaten
	Unstructured_grid* ug = outputVector.at(0)->Give_Ug();
	Blockgrid*  blockgrid = outputVector.at(0)->Give_blockgrid();
	int my_rank           = ug->Give_my_rank();;
	
	// hexahedra
	for ( int id=0; id<ug->Give_number_hexahedra(); ++id ) {
		if ( ug->Give_hexahedron ( id )->my_object ( my_rank ) ) {
			num_total = blockgrid->Give_N_total_hexahedron ( id );
#pragma omp parallel for num_threads(UGBlocks::numThreadsToTake) if(UGBlocks::useOpenMP)
			for(i=0; i<num_total; ++i ) {
			    int tid = omp_get_thread_num();
			    for(int idata=0;idata<sizeInputData;++idata) {
			        inputDataP[tid][idata] = inputVector.at(idata)->Give_pointer_data_hexahedra()[id][i];
			    }
			    functor.compute(inputDataP[tid], tid);
			    for(int idata=0;idata<sizeOutputData;++idata) {
				outputVector.at(idata)->Give_pointer_data_hexahedra()[id][i] = functor.getResult(idata, tid);
			    }
			}
		}
	}

	
	
	// quadrangles
	for(int id=0; id<ug->Give_number_quadrangles(); ++id ) {
		if ( ug->Give_quadrangle ( id )->my_object ( my_rank ) ) {
			num_total = blockgrid->Give_N_total_quadrangle ( id );
			for(i=0; i<num_total; ++i ) {
			    for(int idata=0;idata<sizeInputData;++idata) {
			        inputData[idata] = inputVector.at(idata)->Give_pointer_data_quadrangle()[id][i];
			    }
			    functor.compute(inputData,0);
			    for(int idata=0;idata<sizeOutputData;++idata) {
				outputVector.at(idata)->Give_pointer_data_quadrangle()[id][i] 
				                                    = functor.getResult(idata,0);
			    }
			}
		}
	}

	// edges
	for ( int id=0; id<ug->Give_number_edges(); ++id ) {
		if ( ug->Give_edge ( id )->my_object ( my_rank ) ) {
			num_total = blockgrid->Give_Nx_edge ( id ) +1;
			for(i=0; i<num_total; ++i ) {
			    for(int idata=0;idata<sizeInputData;++idata) {
			        inputData[idata] = inputVector.at(idata)->Give_pointer_data_edges()[id][i];
			    }
			    functor.compute(inputData,0);
			    for(int idata=0;idata<sizeOutputData;++idata) {
				outputVector.at(idata)->Give_pointer_data_edges()[id][i] 
				                                    = functor.getResult(idata,0);
			    }
			}
		}
	}

	// points
	for ( int id=0; id<ug->Give_number_points(); ++id ) {
		if ( ug->Give_point ( id )->my_object ( my_rank ) ) {
			num_total = 1;
			for ( i=0; i<num_total; ++i ) {
			    for(int idata=0;idata<sizeInputData;++idata) {
			        inputData[idata] = inputVector.at(idata)->Give_pointer_data_points()[id][i];
			    }
			    functor.compute(inputData,0);
			    for(int idata=0;idata<sizeOutputData;++idata) {
				outputVector.at(idata)->Give_pointer_data_points()[id][i] 
				                                    = functor.getResult(idata,0);
			    }				
			}
		}
	}
    
    delete[] inputData;
    for(int s=0;s<UGBlocks::numThreadsToTake;++s)
	delete[] inputDataP[s];
    delete inputDataP;
  }  
};



//////////////////////////////////////////////////////////////////////////////


template <class DTyp>
void Update_back_hex(int id_hex, Blockgrid* blockgrid, DTyp** data_hexahedra,
                            DTyp** data_quadrangles, DTyp** data_edges, DTyp** data_points,
                            int my_rank, MPI_Comm comm ) {
	int i,j;
	int id_quad, id_edge, id_point;
	int Ni,Nj,Nconst;
	int Nxedge;
	int Nxquad,Nyquad;
	int Nxhex,Nyhex,Nzhex;
	int id_SW, id_SE, id_NW, idsp;
	int id_W, id_E;
	dir3D_sons c_SW ( ( dir3D_sons ) 0 ), c_SE ( ( dir3D_sons ) 0 ), c_NW ( ( dir3D_sons ) 0 );
	dir3D_sons c_W ( ( dir3D_sons ) 0 ), c_E ( ( dir3D_sons ) 0 );
	dir3D   c_x, c_y;
	Unstructured_grid* ug = blockgrid->Give_unstructured_grid();

	int rank_send, rank_rec;  // for parallel

	Nxhex = blockgrid->Give_Nx_hexahedron ( id_hex );
	Nyhex = blockgrid->Give_Ny_hexahedron ( id_hex );
	Nzhex = blockgrid->Give_Nz_hexahedron ( id_hex );

	rank_rec = ug->Give_hexahedron ( id_hex )->Give_rank(); // for parallel

// keine MPI Paralleliserung!!!!!!!!!!!!!!
assert(	rank_rec == 0);
	
	// quadrangles
	for(int fc=0; fc<6; ++fc ) {
	    id_quad = ug->Give_hexahedron ( id_hex )->Give_id_quadrangle ( ( dir3D ) fc );
	    rank_send = ug->Give_quadrangle ( id_quad )->Give_rank(); // for parallel
	    if(rank_send==my_rank || rank_rec==my_rank ) { // for parallel
	       Nxquad = blockgrid->Give_Nx_quadrangle ( id_quad );
	       Nyquad = blockgrid->Give_Ny_quadrangle ( id_quad );

	       // a) coordinate points of quadrangle
	       id_SW = ug->Give_quadrangle ( id_quad )->Give_id_corner ( SWdir2D );
	       id_SE = ug->Give_quadrangle ( id_quad )->Give_id_corner ( SEdir2D );
	       id_NW = ug->Give_quadrangle ( id_quad )->Give_id_corner ( NWdir2D );

	       // b) Find corresponding points in hexahedron
	       for(int c=0; c<8; ++c ) {
		   idsp = ug->Give_hexahedron ( id_hex )->Give_id_corner ( ( dir3D_sons ) c );
		   if ( idsp==id_SW ) c_SW = ( dir3D_sons ) c;
		   else if ( idsp==id_SE ) c_SE = ( dir3D_sons ) c;
		   else if ( idsp==id_NW ) c_NW = ( dir3D_sons ) c;
	       }

	       c_x = direction_from_to ( c_SW,c_SE );
	       c_y = direction_from_to ( c_SW,c_NW );

	       // x - direction
	       Nconst = Give_Nconst ( ( dir3D ) c_x, Nxhex, Nyhex, Nzhex );
	       Ni     = Give_Ni (    ( dir3D ) c_x, Nxhex, Nyhex, Nzhex );

	       // y - direction
	       Nconst = Give_Nconst ( ( dir3D ) c_y, Nxhex, Nyhex, Nzhex ) + Nconst;
	       Nj     = Give_Ni (    ( dir3D ) c_y, Nxhex, Nyhex, Nzhex );

	       // empty - direction
	       Nconst = Give_Nconst ( opposite3D ( ( dir3D ) fc ),
			                       Nxhex, Nyhex, Nzhex ) + Nconst;
	       if(rank_send == rank_rec ) {    // for parallel
		  for(j=1; j<Nyquad; ++j )
		      for(i=1;i<Nxquad; ++i ) {
			  // vorwaerts:
			  //  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];
		      }
	       }
	       else {
		  assert(false); 
	       }
	    }
	}

	// edges
	for ( int ed=0; ed<12; ++ed ) {
		id_edge = ug->Give_hexahedron ( id_hex )->Give_id_edge ( ( Edges_cell ) ed );
		rank_send = ug->Give_edge ( id_edge )->Give_rank(); // for parallel
		if ( rank_send==my_rank || rank_rec==my_rank ) { // for parallel
			Nxedge = blockgrid->Give_Nx_edge ( id_edge );

			// a) coordinate points of edge
			id_W = ug->Give_edge ( id_edge )->Give_id_corner_W();
			id_E = ug->Give_edge ( id_edge )->Give_id_corner_E();

			// b) Find corresponding points in hexahedron
			for ( int c=0; c<8; ++c ) {
				idsp = ug->Give_hexahedron ( id_hex )->Give_id_corner ( ( dir3D_sons ) c );
				if ( idsp==id_W ) c_W = ( dir3D_sons ) c;
				else if ( idsp==id_E ) c_E = ( dir3D_sons ) c;
			}

			c_x = direction_from_to ( c_W,c_E );

			// x - direction
			Nconst = Give_Nconst ( ( dir3D ) c_x, Nxhex, Nyhex, Nzhex );
			Ni     = Give_Ni    ( ( dir3D ) c_x, Nxhex, Nyhex, Nzhex );

			// A. empty - direction
			Nconst = Give_Nconst ( opposite3D ( OrthoA ( ( Edges_cell ) ed ) ),
			                       Nxhex, Nyhex, Nzhex ) + Nconst;

			// B. empty - direction
			Nconst = Give_Nconst ( opposite3D ( OrthoB ( ( Edges_cell ) ed ) ),
			                       Nxhex, Nyhex, Nzhex ) + Nconst;

			if(rank_send == rank_rec ) {    // for parallel
			   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];
			   }
			} 
			else {
			   assert(false); 
			}
		}
	}

	// points
	for ( int c=0; c<8; ++c ) {
		id_point = ug->Give_hexahedron ( id_hex )->Give_id_corner ( ( dir3D_sons ) c );
		rank_send = ug->Give_point ( id_point )->Give_rank(); // for parallel
		if ( rank_send==my_rank || rank_rec==my_rank ) { // for parallel

			// x- direction
			Nconst = Give_Nconst ( opposite3D ( xCoord ( ( dir3D_sons ) c ) ),
			                       Nxhex, Nyhex, Nzhex );

			// y- direction
			Nconst = Give_Nconst ( opposite3D ( yCoord ( ( dir3D_sons ) c ) ),
			                       Nxhex, Nyhex, Nzhex ) + Nconst;

			// z- direction
			Nconst = Give_Nconst ( opposite3D ( zCoord ( ( dir3D_sons ) c ) ),
			                       Nxhex, Nyhex, Nzhex ) + Nconst;

			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];
            }
			else {
			  assert(false);
			}
		}
	}
}





//////////////////////////////////////////////////////////////////////////////











template <class DTyp>
void Update_hexahedron ( int id_hex, Blockgrid* blockgrid, DTyp** data_hexahedra,
                         DTyp** data_quadrangles, DTyp** data_edges, DTyp** data_points,
                         int my_rank, MPI_Comm comm ) {
	int id_quad, id_edge, id_point;
	int Ni,Nj,Nconst;
	int Nxedge;
	int Nxquad,Nyquad;
	int Nxhex,Nyhex,Nzhex;
	int id_SW, id_SE, id_NW, idsp;
	int id_W, id_E;
	dir3D_sons c_SW ( ( dir3D_sons ) 0 ), c_SE ( ( dir3D_sons ) 0 ), c_NW ( ( dir3D_sons ) 0 );
	dir3D_sons c_W ( ( dir3D_sons ) 0 ), c_E ( ( dir3D_sons ) 0 );
	dir3D   c_x, c_y;
	Unstructured_grid* ug = blockgrid->Give_unstructured_grid();

	int rank_send, rank_rec;  // for parallel

	Nxhex = blockgrid->Give_Nx_hexahedron ( id_hex );
	Nyhex = blockgrid->Give_Ny_hexahedron ( id_hex );
	Nzhex = blockgrid->Give_Nz_hexahedron ( id_hex );

	rank_rec = ug->Give_hexahedron ( id_hex )->Give_rank(); // for parallel

	// quadrangles
	for ( int fc=0; fc<6; ++fc ) {
		id_quad = ug->Give_hexahedron ( id_hex )->Give_id_quadrangle ( ( dir3D ) fc );
		rank_send = ug->Give_quadrangle ( id_quad )->Give_rank(); // for parallel
		if ( rank_send==my_rank || rank_rec==my_rank ) { // for parallel
			Nxquad = blockgrid->Give_Nx_quadrangle ( id_quad );
			Nyquad = blockgrid->Give_Ny_quadrangle ( id_quad );

			// a) coordinate points of quadrangle
			id_SW = ug->Give_quadrangle ( id_quad )->Give_id_corner ( SWdir2D );
			id_SE = ug->Give_quadrangle ( id_quad )->Give_id_corner ( SEdir2D );
			id_NW = ug->Give_quadrangle ( id_quad )->Give_id_corner ( NWdir2D );

			// b) Find corresponding points in hexahedron
			for ( int c=0; c<8; ++c ) {
				idsp = ug->Give_hexahedron ( id_hex )->Give_id_corner ( ( dir3D_sons ) c );
				if ( idsp==id_SW ) c_SW = ( dir3D_sons ) c;
				else if ( idsp==id_SE ) c_SE = ( dir3D_sons ) c;
				else if ( idsp==id_NW ) c_NW = ( dir3D_sons ) c;
			}

			c_x = direction_from_to ( c_SW,c_SE );
			c_y = direction_from_to ( c_SW,c_NW );

			// x - direction
			Nconst = Give_Nconst ( ( dir3D ) c_x, Nxhex, Nyhex, Nzhex );
			Ni     = Give_Ni (    ( dir3D ) c_x, Nxhex, Nyhex, Nzhex );
					
			// y - direction
			Nconst = Give_Nconst ( ( dir3D ) c_y, Nxhex, Nyhex, Nzhex ) + Nconst;
			Nj     = Give_Ni (    ( dir3D ) c_y, Nxhex, Nyhex, Nzhex );

			// empty - direction
			Nconst = Give_Nconst ( opposite3D ( ( dir3D ) fc ),
			                       Nxhex, Nyhex, Nzhex ) + Nconst;
			if ( rank_send == rank_rec ) {    // for parallel
			        // neu parallel
			        //bool do_parallel_here = UGBlocks::useOpenMP && (Ni==1);
			        //#pragma omp parallel for num_threads(UGBlocks::numThreadsToTake) if(do_parallel_here)			  
				for (int j=1; j<Nyquad; ++j )
					for (int i=1; i<Nxquad; ++i )
						data_hexahedra[id_hex][i*Ni+j*Nj+Nconst] =
						    data_quadrangles[id_quad][i+ ( Nxquad+1 ) *j];
			} else {
				if ( my_rank==rank_send )
					MPI_Send ( &data_quadrangles[id_quad][Nxquad+1], ( Nxquad+1 ) * ( Nyquad-1 ),
					           Give_MPI_data_typ<DTyp>::Give(), rank_rec, fc, comm );
				else {
					if ( developer_version )
						if ( my_rank!=rank_rec )
							cout << "error1 in Update_hexahedron! " << endl;
					MPI_Status stat;
					DTyp* data = Parallel_buffer<DTyp>::Give_receive ( ( Nxquad+1 ) * ( Nyquad-1 ) );
					MPI_Recv ( data, ( Nxquad+1 ) * ( Nyquad-1 ), Give_MPI_data_typ<DTyp>::Give(),
					           rank_send, fc, comm, &stat );
			                // neu parallel
			                bool do_parallel_here = UGBlocks::useOpenMP && (Ni==1);
			                #pragma omp parallel for num_threads(UGBlocks::numThreadsToTake) if(do_parallel_here)			  
					for(int j=1; j<Nyquad; ++j ) {
						for (int i=1; i<Nxquad; ++i ) {
							data_hexahedra[id_hex][i*Ni+j*Nj+Nconst] =
							    data[i+ ( Nxquad+1 ) * ( j-1 ) ];
						}
					}
				}
			}
		}
	}

	// edges
	for ( int ed=0; ed<12; ++ed ) {
		id_edge = ug->Give_hexahedron ( id_hex )->Give_id_edge ( ( Edges_cell ) ed );
		rank_send = ug->Give_edge ( id_edge )->Give_rank(); // for parallel
		if ( rank_send==my_rank || rank_rec==my_rank ) { // for parallel
			Nxedge = blockgrid->Give_Nx_edge ( id_edge );

			// a) coordinate points of edge
			id_W = ug->Give_edge ( id_edge )->Give_id_corner_W();
			id_E = ug->Give_edge ( id_edge )->Give_id_corner_E();

			// b) Find corresponding points in hexahedron
			for ( int c=0; c<8; ++c ) {
				idsp = ug->Give_hexahedron ( id_hex )->Give_id_corner ( ( dir3D_sons ) c );
				if ( idsp==id_W ) c_W = ( dir3D_sons ) c;
				else if ( idsp==id_E ) c_E = ( dir3D_sons ) c;
			}

			c_x = direction_from_to ( c_W,c_E );

			// x - direction
			Nconst = Give_Nconst ( ( dir3D ) c_x, Nxhex, Nyhex, Nzhex );
			Ni     = Give_Ni    ( ( dir3D ) c_x, Nxhex, Nyhex, Nzhex );

			// A. empty - direction
			Nconst = Give_Nconst ( opposite3D ( OrthoA ( ( Edges_cell ) ed ) ),
			                       Nxhex, Nyhex, Nzhex ) + Nconst;

			// B. empty - direction
			Nconst = Give_Nconst ( opposite3D ( OrthoB ( ( Edges_cell ) ed ) ),
			                       Nxhex, Nyhex, Nzhex ) + Nconst;

			if ( rank_send == rank_rec ) {    // for parallel

				for (int i=1; i<Nxedge; ++i )
					data_hexahedra[id_hex][i*Ni+Nconst]
					= data_edges[id_edge][i];
			} else {
				if ( my_rank==rank_send )
					MPI_Send ( &data_edges[id_edge][1], Nxedge-1,
					           Give_MPI_data_typ<DTyp>::Give(), rank_rec, 20+ed, comm );
				else {
					if ( developer_version )
						if ( my_rank!=rank_rec )
							cout << "error2 in Update_hexahedron! " << endl;
					MPI_Status stat;
					DTyp* data = Parallel_buffer<DTyp>::Give_receive ( Nxedge-1 );
					MPI_Recv ( data, Nxedge-1, Give_MPI_data_typ<DTyp>::Give(),
					           rank_send, 20+ed, comm, &stat );
					for(int i=1; i<Nxedge; ++i ) {
						data_hexahedra[id_hex][i*Ni+Nconst]
						= data[i-1];
					}
				}
			}
		}
	}

	// points
	for ( int c=0; c<8; ++c ) {
		id_point = ug->Give_hexahedron ( id_hex )->Give_id_corner ( ( dir3D_sons ) c );
		rank_send = ug->Give_point ( id_point )->Give_rank(); // for parallel
		if ( rank_send==my_rank || rank_rec==my_rank ) { // for parallel

			// x- direction
			Nconst = Give_Nconst ( opposite3D ( xCoord ( ( dir3D_sons ) c ) ),
			                       Nxhex, Nyhex, Nzhex );

			// y- direction
			Nconst = Give_Nconst ( opposite3D ( yCoord ( ( dir3D_sons ) c ) ),
			                       Nxhex, Nyhex, Nzhex ) + Nconst;

			// z- direction
			Nconst = Give_Nconst ( opposite3D ( zCoord ( ( dir3D_sons ) c ) ),
			                       Nxhex, Nyhex, Nzhex ) + Nconst;

			if ( rank_send == rank_rec )      // for parallel
				data_hexahedra[id_hex][Nconst]
				= data_points[id_point][0];
			else {
				if ( my_rank==rank_send )
					MPI_Send ( data_points[id_point], 1,
					           Give_MPI_data_typ<DTyp>::Give(), rank_rec, 40+c, comm );
				else {
					if ( developer_version )
						if ( my_rank!=rank_rec )
							cout << "error3 in Update_hexahedron! " << endl;
					MPI_Status stat;
					MPI_Recv ( &data_hexahedra[id_hex][Nconst], 1, Give_MPI_data_typ<DTyp>::Give(),
					           rank_send, 40+c, comm, &stat );
				}
			}
		}
	}
}

/*
template <>
template <>
void Variable<double>::Update<hexahedronEl>(int id_hex) const{
  Update_hexahedron<double>(id_hex, blockgrid, data_hexahedra, data_quadrangles,
                                               data_edges, data_points);
}

template <>
template <>
void Variable<std::complex<double> >::Update<hexahedronEl>(int id_hex) const{
  Update_hexahedron<std::complex<double> >(id_hex, blockgrid, data_hexahedra,
                                           data_quadrangles, data_edges, data_points);
}
*/
template <class DTyp>
void Update_quadrangle ( int id_quad, Blockgrid* blockgrid, DTyp** data_hexahedra,
                         DTyp** data_quadrangles, DTyp** data_edges, DTyp** data_points,
                         DTyp**  data_interior_quadrangles,
                         DTyp**  data_exterior_quadrangles, int my_rank, MPI_Comm comm )
{
	int i,j,num;
	int id_hex, id_face, id_edge, id_point;
	int Ni,Nj ( 0 ),Nconst;
	int Nxedge;
	int Nxquad,Nyquad,Nxface,Nyface;
	int Nxhex,Nyhex,Nzhex;
	int id_SW, id_SE, id_NW, id_NE, id_SoW, id_NoE, idsp;
	int id_SW_face, id_SE_face, id_NW_face, id_NE_face;
	int id_L, id_R;
	dir3D_sons c_SW ( ( dir3D_sons ) 0 ), c_SE ( ( dir3D_sons ) 0 ), c_NW ( ( dir3D_sons ) 0 );
	dir3D  c_x, c_y, quad;
	dir2D edge;
	dir2D_sons point;
	Unstructured_grid* ug = blockgrid->Give_unstructured_grid();

	int rank_send, rank_rec;  // for parallel

	Nxquad = blockgrid->Give_Nx_quadrangle ( id_quad );
	Nyquad = blockgrid->Give_Ny_quadrangle ( id_quad );

	rank_rec = ug->Give_quadrangle ( id_quad )->Give_rank(); // for parallel

	// edges Wdir2D, Edir2D, Sdir2D, Ndir2D

	// Wdir2D
	id_edge = ug->Give_quadrangle ( id_quad )->Give_id_edge ( Wdir2D );
	rank_send = ug->Give_edge ( id_edge )->Give_rank(); // for parallel
	if ( rank_send==my_rank || rank_rec==my_rank ) { // for parallel
		Nxedge = blockgrid->Give_Nx_edge ( id_edge );
		Nconst = 0;
		Ni     = Nxquad+1;
		if ( rank_send == rank_rec )    // for parallel

			for ( i=1; i<Nxedge; ++i )
				data_quadrangles[id_quad][i*Ni+Nconst]
				= data_edges[id_edge][i];
		else {
			if ( my_rank==rank_send )
				MPI_Send ( &data_edges[id_edge][1], Nxedge-1,
				           Give_MPI_data_typ<DTyp>::Give(), rank_rec, 10, comm );
			else {
				if ( developer_version )
					if ( my_rank!=rank_rec )
						cout << "error1 in Update_quadrangle! " << endl;
				MPI_Status stat;
				MPI_Datatype columntype;
				MPI_Type_vector ( Nxedge-1, 1, Ni, Give_MPI_data_typ<DTyp>::Give(), &columntype );
				MPI_Type_commit ( &columntype );
				MPI_Recv ( &data_quadrangles[id_quad][Ni], 1, columntype,
				           rank_send, 10, comm, &stat );
				MPI_Type_free ( &columntype );
			}
		}
	}

	// Edir2D
	id_edge = ug->Give_quadrangle ( id_quad )->Give_id_edge ( Edir2D );
	rank_send = ug->Give_edge ( id_edge )->Give_rank(); // for parallel
	if ( rank_send==my_rank || rank_rec==my_rank ) { // for parallel
		Nxedge = blockgrid->Give_Nx_edge ( id_edge );
		if ( ug->Give_quadrangle ( id_quad )->Give_id_corner ( SEdir2D )
		     < ug->Give_quadrangle ( id_quad )->Give_id_corner ( NEdir2D ) ) { //Assure that SEcorner has
			//smaller id than NEcorner
			Nconst = Nxquad;
			Ni     = Nxquad+1;
		} else {
			Nconst = Nxquad+ ( Nxquad+1 ) * ( Nyquad+1 );
			Ni     = - ( Nxquad+1 );
		}
		if ( rank_send == rank_rec )    // for parallel

			for ( i=1; i<Nxedge; ++i )
				data_quadrangles[id_quad][i*Ni+Nconst]
				= data_edges[id_edge][i];
		else {
			if ( my_rank==rank_send )
				MPI_Send ( &data_edges[id_edge][1], Nxedge-1,
				           Give_MPI_data_typ<DTyp>::Give(), rank_rec, 20, comm );
			else {
				if ( developer_version )
					if ( my_rank!=rank_rec )
						cout << "error2 in Update_quadrangle! " << endl;
				MPI_Status stat;
				DTyp* data = Parallel_buffer<DTyp>::Give_receive ( Nxedge-1 );
				MPI_Recv ( data, Nxedge-1, Give_MPI_data_typ<DTyp>::Give(),
				           rank_send, 20, comm, &stat );
				for ( i=1; i<Nxedge; ++i )
					data_quadrangles[id_quad][i*Ni+Nconst]
					= data[i-1];
			}
		}
	}

	// Sdir2D
	id_edge = ug->Give_quadrangle ( id_quad )->Give_id_edge ( Sdir2D );
	rank_send = ug->Give_edge ( id_edge )->Give_rank(); // for parallel
	if ( rank_send==my_rank || rank_rec==my_rank ) { // for parallel
		Nxedge = blockgrid->Give_Nx_edge ( id_edge );
		Nconst = 0;
		Ni     = 1;
		if ( rank_send == rank_rec )    // for parallel

			for ( i=1; i<Nxedge; ++i )
				data_quadrangles[id_quad][i*Ni+Nconst]
				= data_edges[id_edge][i];
		else {
			if ( my_rank==rank_send )
				MPI_Send ( &data_edges[id_edge][1], Nxedge-1,
				           Give_MPI_data_typ<DTyp>::Give(), rank_rec, 30, comm );
			else {
				if ( developer_version )
					if ( my_rank!=rank_rec )
						cout << "error3 in Update_quadrangle! " << endl;
				MPI_Status stat;
				MPI_Recv ( &data_quadrangles[id_quad][1], Nxedge-1, Give_MPI_data_typ<DTyp>::Give(),
				           rank_send, 30, comm, &stat );
			}
		}
	}

	// Ndir2D
	id_edge = ug->Give_quadrangle ( id_quad )->Give_id_edge ( Ndir2D );
	rank_send = ug->Give_edge ( id_edge )->Give_rank(); // for parallel
	if ( rank_send==my_rank || rank_rec==my_rank ) { // for parallel
		Nxedge = blockgrid->Give_Nx_edge ( id_edge );
		if ( ug->Give_quadrangle ( id_quad )->Give_id_corner ( NWdir2D )
		     < ug->Give_quadrangle ( id_quad )->Give_id_corner ( NEdir2D ) ) { //Assure that NWcorner has
			//smaller id than NEcorner
			Nconst = ( Nxquad+1 ) *Nyquad;
			Ni     = 1;
		} else {
			Nconst = ( Nxquad+1 ) * ( Nyquad+1 )-1;
			Ni     = -1;
		}
		if ( rank_send == rank_rec ) {   // for parallel

			for ( i=1; i<Nxedge; ++i ) {
				if ( developer_version && ( i*Ni+Nconst >= ( Nxquad+1 ) * ( Nyquad+1 ) ) ) {
					cout << "error4 in Update_quadrangle! " << endl;
				}
				data_quadrangles[id_quad][i*Ni+Nconst]
				= data_edges[id_edge][i];
			}
		} else {
			if ( my_rank==rank_send )
				MPI_Send ( &data_edges[id_edge][1], Nxedge-1,
				           Give_MPI_data_typ<DTyp>::Give(), rank_rec, 40, comm );
			else {
				if ( developer_version )
					if ( my_rank!=rank_rec )
						cout << "error5 in Update_quadrangle! " << endl;
				MPI_Status stat;
				DTyp* data = Parallel_buffer<DTyp>::Give_receive ( Nxedge-1 );
				MPI_Recv ( data, Nxedge-1, Give_MPI_data_typ<DTyp>::Give(),
				           rank_send, 40, comm, &stat );
				for ( i=1; i<Nxedge; ++i ) {
					if ( developer_version )
						if ( i*Ni+Nconst >= ( Nxquad+1 ) * ( Nyquad+1 ) )
							cout << "error6 in Update_quadrangle! " << endl;
					data_quadrangles[id_quad][i*Ni+Nconst]
					= data[i-1];
				}
			}
		}
	}

	// points SWdir2D, SEdir2D, NWdir2D, NEdir2D
	id_point = ug->Give_quadrangle ( id_quad )->Give_id_corner ( SWdir2D );
	rank_send = ug->Give_point ( id_point )->Give_rank(); // for parallel
	if ( rank_send==my_rank || rank_rec==my_rank ) { // for parallel
		if ( rank_send == rank_rec )    // for parallel
			data_quadrangles[id_quad][0] = data_points[id_point][0];
		else {
			if ( my_rank==rank_send )
				MPI_Send ( data_points[id_point], 1,
				           Give_MPI_data_typ<DTyp>::Give(), rank_rec, 11, comm );
			else {
				if ( developer_version )
					if ( my_rank!=rank_rec )
						cout << "error7 in Update_quadrangle! " << endl;
				MPI_Status stat;
				MPI_Recv ( data_quadrangles[id_quad], 1, Give_MPI_data_typ<DTyp>::Give(),
				           rank_send, 11, comm, &stat );
			}
		}
	}

	id_point = ug->Give_quadrangle ( id_quad )->Give_id_corner ( SEdir2D );
	rank_send = ug->Give_point ( id_point )->Give_rank(); // for parallel
	if ( rank_send==my_rank || rank_rec==my_rank ) { // for parallel
		if ( rank_send == rank_rec )    // for parallel
			data_quadrangles[id_quad][Nxquad] = data_points[id_point][0];
		else {
			if ( my_rank==rank_send )
				MPI_Send ( data_points[id_point], 1,
				           Give_MPI_data_typ<DTyp>::Give(), rank_rec, 21, comm );
			else {
				if ( developer_version )
					if ( my_rank!=rank_rec )
						cout << "error8 in Update_quadrangle! " << endl;
				MPI_Status stat;
				MPI_Recv ( &data_quadrangles[id_quad][Nxquad], 1, Give_MPI_data_typ<DTyp>::Give(),
				           rank_send, 21, comm, &stat );
			}
		}
	}

	id_point = ug->Give_quadrangle ( id_quad )->Give_id_corner ( NWdir2D );
	rank_send = ug->Give_point ( id_point )->Give_rank(); // for parallel
	if ( rank_send==my_rank || rank_rec==my_rank ) { // for parallel
		if ( rank_send == rank_rec )    // for parallel
			data_quadrangles[id_quad][ ( Nxquad+1 ) *Nyquad] = data_points[id_point][0];
		else {
			if ( my_rank==rank_send )
				MPI_Send ( data_points[id_point], 1,
				           Give_MPI_data_typ<DTyp>::Give(), rank_rec, 31, comm );
			else {
				if ( developer_version )
					if ( my_rank!=rank_rec )
						cout << "error9 in Update_quadrangle! " << endl;
				MPI_Status stat;
				MPI_Recv ( &data_quadrangles[id_quad][ ( Nxquad+1 ) *Nyquad], 1,
				           Give_MPI_data_typ<DTyp>::Give(),
				           rank_send, 31, comm, &stat );
			}
		}
	}

	id_point = ug->Give_quadrangle ( id_quad )->Give_id_corner ( NEdir2D );
	rank_send = ug->Give_point ( id_point )->Give_rank(); // for parallel
	if ( rank_send==my_rank || rank_rec==my_rank ) { // for parallel
		if ( rank_send == rank_rec )    // for parallel
			data_quadrangles[id_quad][ ( Nxquad+1 ) * ( Nyquad+1 )-1] = data_points[id_point][0];
		else {
			if ( my_rank==rank_send )
				MPI_Send ( data_points[id_point], 1,
				           Give_MPI_data_typ<DTyp>::Give(), rank_rec, 41, comm );
			else {
				if ( developer_version )
					if ( my_rank!=rank_rec )
						cout << "error10 in Update_quadrangle! " << endl;
				MPI_Status stat;
				MPI_Recv ( &data_quadrangles[id_quad][ ( Nxquad+1 ) * ( Nyquad+1 )-1], 1,
				           Give_MPI_data_typ<DTyp>::Give(),
				           rank_send, 41, comm, &stat );
			}
		}
	}

	// unstructured data

	// coordinate points of quadrangle
	id_SW = ug->Give_quadrangle ( id_quad )->Give_id_corner ( SWdir2D );
	id_SE = ug->Give_quadrangle ( id_quad )->Give_id_corner ( SEdir2D );
	id_NW = ug->Give_quadrangle ( id_quad )->Give_id_corner ( NWdir2D );
	id_NE = ug->Give_quadrangle ( id_quad )->Give_id_corner ( NEdir2D );

	// interior hexahedron
	id_hex = ug->Give_quadrangle ( id_quad )->Give_id_interior_vol();
	rank_send = ug->Give_hexahedron ( id_hex )->Give_rank(); // for parallel
	if ( rank_send==my_rank || rank_rec==my_rank ) { // for parallel
		Nxhex = blockgrid->Give_Nx_hexahedron ( id_hex );
		Nyhex = blockgrid->Give_Ny_hexahedron ( id_hex );
		Nzhex = blockgrid->Give_Nz_hexahedron ( id_hex );

		// Find corresponding points in hexahedron
		for ( int c=0; c<8; ++c ) {
			idsp = ug->Give_hexahedron ( id_hex )->Give_id_corner ( ( dir3D_sons ) c );
			if ( idsp==id_SW ) c_SW = ( dir3D_sons ) c;
			else if ( idsp==id_SE ) c_SE = ( dir3D_sons ) c;
			else if ( idsp==id_NW ) c_NW = ( dir3D_sons ) c;
		}

		c_x = direction_from_to ( c_SW,c_SE );
		c_y = direction_from_to ( c_SW,c_NW );

		for ( int fc=0; fc<6; ++fc )
			if ( id_quad == ug->Give_hexahedron ( id_hex )->Give_id_quadrangle ( ( dir3D ) fc ) ) {
				quad = ( dir3D ) fc;

				// x - direction
				Nconst = Give_Nconst ( ( dir3D ) c_x, Nxhex, Nyhex, Nzhex );
				Ni     = Give_Ni (    ( dir3D ) c_x, Nxhex, Nyhex, Nzhex );

				// y - direction
				Nconst = Give_Nconst ( ( dir3D ) c_y, Nxhex, Nyhex, Nzhex ) + Nconst;
				Nj     = Give_Ni (    ( dir3D ) c_y, Nxhex, Nyhex, Nzhex );

				// empty - direction
				Nconst = Give_Ni    ( opposite3D ( quad ),Nxhex, Nyhex, Nzhex )
				         + Give_Nconst ( opposite3D ( quad ),Nxhex, Nyhex, Nzhex )
				         + Nconst;
			}

		if ( rank_send == rank_rec ) {   // for parallel

			for ( j=1; j<Nyquad; ++j )
				for ( i=1; i<Nxquad; ++i )
					data_interior_quadrangles[id_quad][i+ ( Nxquad+1 ) *j] =
					    data_hexahedra[id_hex][i*Ni+j*Nj+Nconst];
		} else {
			if ( my_rank==rank_send ) {
				DTyp* data = Parallel_buffer<DTyp>::Give_send ( ( Nxquad-1 ) * ( Nyquad-1 ) );
				for ( j=1; j<Nyquad; ++j )
					for ( i=1; i<Nxquad; ++i )
						data[ ( i-1 ) + ( Nxquad-1 ) * ( j-1 ) ] =
						    data_hexahedra[id_hex][i*Ni+j*Nj+Nconst];
				MPI_Send ( data, ( Nxquad-1 ) * ( Nyquad-1 ),
				           Give_MPI_data_typ<DTyp>::Give(), rank_rec, 22, comm );
			} else {
				if ( developer_version )
					if ( my_rank!=rank_rec )
						cout << "error11 in Update_quadrangle! " << endl;
				MPI_Status stat;
				MPI_Datatype columntype;
				MPI_Type_vector ( Nyquad-1, Nxquad-1, Nxquad+1, Give_MPI_data_typ<DTyp>::Give(),
				                  &columntype );
				MPI_Type_commit ( &columntype );
				MPI_Recv ( &data_interior_quadrangles[id_quad][Nxquad+2], 1, columntype,
				           rank_send, 22, comm, &stat );
				MPI_Type_free ( &columntype );
			}
		}
	}

	// interior face Wdir2D, Edir2D, Sdir2D, Ndir2D
	for ( num=0; num<4; ++num ) { // Iterate over all edges of quad
		edge = ( dir2D ) num;
		id_face = ug->Give_quadrangle ( id_quad )->Give_id_interior_face ( edge ); // id of interior face
		rank_send = ug->Give_quadrangle ( id_face )->Give_rank(); // for parallel
		if ( rank_send==my_rank || rank_rec==my_rank ) { // for parallel

			Nxface = blockgrid->Give_Nx_quadrangle ( id_face );
			Nyface = blockgrid->Give_Ny_quadrangle ( id_face );

			// coordinate points of face
			id_SW_face = ug->Give_quadrangle ( id_face )->Give_id_corner ( SWdir2D );
			id_SE_face = ug->Give_quadrangle ( id_face )->Give_id_corner ( SEdir2D );
			id_NW_face = ug->Give_quadrangle ( id_face )->Give_id_corner ( NWdir2D );
			id_NE_face = ug->Give_quadrangle ( id_face )->Give_id_corner ( NEdir2D );

			// Set reference corners of quadrangle
			if ( edge == Wdir2D ) {
				id_SoW = id_SW;
				id_NoE = id_NW;
			} else if ( edge == Edir2D ) {
				id_SoW = id_SE;
				id_NoE = id_NE;
			} else if ( edge == Sdir2D ) {
				id_SoW = id_SW;
				id_NoE = id_SE;
			} else {
				if ( developer_version )
					if ( edge != Ndir2D ) cout << "error12 in Update_quadrangle! " << endl;
				id_SoW = id_NW;
				id_NoE = id_NE;
			}

			// Compare corners of quad with corners of face
			if ( id_SW_face == id_SoW ) {
				if ( id_NW_face == id_NoE ) {
					Ni = Nj = Nxface+1;
					Nconst = 1;
				} else {
					if ( developer_version )
						if ( id_SE_face != id_NoE ) cout << "error13 in Update_quadrangle! " << endl;
					Ni = Nj = 1;
					Nconst = Nxface+1;
				}
			} else if ( id_SE_face == id_SoW ) {
				if ( id_SW_face == id_NoE ) {
					Ni = Nj = -1;
					Nconst = 2*Nxface+1;
				} else {
					if ( developer_version )
						if ( id_NE_face != id_NoE ) cout << "error14 in Update_quadrangle! " << endl;
					Ni = Nj = Nxface+1;
					Nconst = Nxface-1;
				}
			} else if ( id_NW_face == id_SoW ) {
				if ( id_NE_face == id_NoE ) {
					Ni = Nj = 1;
					Nconst = ( Nyface-1 ) * ( Nxface+1 );
				} else {
					if ( developer_version )
						if ( id_SW_face != id_NoE ) cout << "error15 in Update_quadrangle! " << endl;
					Ni = Nj = - ( Nxface+1 );
					Nconst = Nyface* ( Nxface+1 ) +1;
				}
			} else {
				if ( developer_version )
					if ( id_NE_face != id_SoW ) cout << "error16 in Update_quadrangle! " << endl;
				if ( id_SE_face == id_NoE ) {
					Ni = Nj = - ( Nxface+1 );
					Nconst = ( Nyface+1 ) * ( Nxface+1 )-2;
				} else {
					if ( developer_version )
						if ( id_NW_face != id_NoE ) cout << "error17 in Update_quadrangle! " << endl;
					Ni = Nj = -1;
					Nconst = Nyface* ( Nxface+1 )-1;
				}
			}

			if ( rank_send == rank_rec ) {   // for parallel
				if ( edge == Wdir2D ) {
					i = 0;

					for ( j=1; j<Nyquad; ++j )
						data_interior_quadrangles[id_quad][i+ ( Nxquad+1 ) *j] =
						    data_quadrangles[id_face][j*Nj+Nconst];
				} else if ( edge == Edir2D ) {
					i = Nxquad;

					for ( j=1; j<Nyquad; ++j ) {
						data_interior_quadrangles[id_quad][i+ ( Nxquad+1 ) *j] =
						    data_quadrangles[id_face][j*Nj+Nconst];
					}
				} else if ( edge == Sdir2D ) {
					j = 0;

					for ( i=1; i<Nxquad; ++i )
						data_interior_quadrangles[id_quad][i+ ( Nxquad+1 ) *j] =
						    data_quadrangles[id_face][i*Ni+Nconst];
				} else {
					if ( developer_version )
						if ( edge != Ndir2D ) cout << "error18 in Update_quadrangle! " << endl;
					j = Nyquad;

					for ( i=1; i<Nxquad; ++i ) {
						data_interior_quadrangles[id_quad][i+ ( Nxquad+1 ) *j] =
						    data_quadrangles[id_face][i*Ni+Nconst];
					}
				}
			} else {
				if ( my_rank==rank_send ) {
					DTyp* data = Parallel_buffer<DTyp>::Give_send ( Nxquad-1 );
					for ( i=1; i<Nxquad; ++i )
						data[i-1] =
						    data_quadrangles[id_face][i*Ni+Nconst];
					MPI_Send ( data, Nxquad-1,
					           Give_MPI_data_typ<DTyp>::Give(), rank_rec, 33+num, comm );
				} else {
					if ( developer_version )
						if ( my_rank!=rank_rec )
							cout << "error19 in Update_quadrangle! " << endl;
					if ( edge == Wdir2D || edge == Edir2D ) {
						MPI_Status stat;
						MPI_Datatype columntype;
						MPI_Type_vector ( Nyquad-1, 1, Nxquad+1, Give_MPI_data_typ<DTyp>::Give(),
						                  &columntype );
						MPI_Type_commit ( &columntype );
						if ( edge == Wdir2D )
							i = 0;
						else {
							if ( developer_version )
								if ( edge != Edir2D )
									cout << "error20 in Update_quadrangle! " << endl;
							i = Nxquad;
						}
						MPI_Recv ( &data_interior_quadrangles[id_quad][i+Nxquad+1], 1, columntype,
						           rank_send, 33+num, comm, &stat );
						MPI_Type_free ( &columntype );
					} else {
						MPI_Status stat;
						if ( edge == Sdir2D )
							j = 0;
						else {
							if ( developer_version )
								if ( edge != Ndir2D ) cout << "error21 in Update_quadrangle! "
									                           << endl;
							j = Nyquad;
						}
						MPI_Recv ( &data_interior_quadrangles[id_quad][1+ ( Nxquad+1 ) *j], Nxquad-1,
						           Give_MPI_data_typ<DTyp>::Give(), rank_send, 33+num, comm, &stat );
					}
				}
			}
		}
	}

	// interior edge SWdir2D, SEdir2D, NWdir2D, NEdir2D
	for ( num=0; num<4; ++num ) { // Iterate over all corners of quad
		point = ( dir2D_sons ) num;
		id_edge = ug->Give_quadrangle ( id_quad )->Give_id_interior_edge ( point ); // id of interior edge
		rank_send = ug->Give_edge ( id_edge )->Give_rank(); // for parallel
		if ( rank_send==my_rank || rank_rec==my_rank ) { // for parallel

			id_point = ug->Give_quadrangle ( id_quad )->Give_id_corner ( point ); // id of corner

			Nxedge = blockgrid->Give_Nx_edge ( id_edge );

			// coordinate points of edge
			id_L = ug->Give_edge ( id_edge )->Give_id_corner ( Ldir1D );
			id_R = ug->Give_edge ( id_edge )->Give_id_corner ( Rdir1D );

			if ( rank_send == rank_rec ) {   // for parallel
				// Compare corners of quad with corners of edge
				if ( point == SWdir2D ) {
					if ( id_point == id_L )
						data_interior_quadrangles[id_quad][0] =
						    data_edges[id_edge][1];
					else {
						if ( developer_version )
							if ( id_point != id_R )
								cout << "error22 in Update_quadrangle! " << endl;
						data_interior_quadrangles[id_quad][0] =
						    data_edges[id_edge][Nxedge-1];
					}
				} else if ( point == SEdir2D ) {
					if ( id_point == id_L )
						data_interior_quadrangles[id_quad][Nxquad] =
						    data_edges[id_edge][1];
					else {
						if ( developer_version )
							if ( id_point != id_R )
								cout << "error23 in Update_quadrangle! " << endl;
						data_interior_quadrangles[id_quad][Nxquad] =
						    data_edges[id_edge][Nxedge-1];
					}
				} else if ( point == NWdir2D ) {
					if ( id_point == id_L )
						data_interior_quadrangles[id_quad][ ( Nxquad+1 ) *Nyquad] =
						    data_edges[id_edge][1];
					else {
						if ( developer_version )
							if ( id_point != id_R )
								cout << "error24 in Update_quadrangle! " << endl;
						data_interior_quadrangles[id_quad][ ( Nxquad+1 ) *Nyquad] =
						    data_edges[id_edge][Nxedge-1];
					}
				} else {
					if ( developer_version )
						if ( point != NEdir2D ) cout << "error25 in Update_quadrangle! " << endl;
					if ( id_point == id_L )
						data_interior_quadrangles[id_quad][ ( Nxquad+1 ) * ( Nyquad+1 )-1] =
						    data_edges[id_edge][1];
					else {
						if ( developer_version )
							if ( id_point != id_R )
								cout << "error26 in Update_quadrangle! " << endl;
						data_interior_quadrangles[id_quad][ ( Nxquad+1 ) * ( Nyquad+1 )-1] =
						    data_edges[id_edge][Nxedge-1];
					}
				}
			} else {
				if ( my_rank==rank_send ) {
					if ( id_point == id_L ) i=1;
					else {
						if ( developer_version )
							if ( id_point != id_R )
								cout << "error27 in Update_quadrangle! " << endl;
						i=Nxedge-1;
					}
					MPI_Send ( &data_edges[id_edge][i], 1,
					           Give_MPI_data_typ<DTyp>::Give(), rank_rec, 44+num, comm );
				} else {
					if ( developer_version )
						if ( my_rank!=rank_rec )
							cout << "error28 in Update_quadrangle! " << endl;
					if ( point == SWdir2D ) j=0;
					else if ( point == SEdir2D ) j=Nxquad;
					else if ( point == NWdir2D ) j= ( Nxquad+1 ) *Nyquad;
					else {
						if ( developer_version )
							if ( point != NEdir2D ) cout << "error29 in Update_quadrangle! " << endl;
						j= ( Nxquad+1 ) * ( Nyquad+1 )-1;
					}
					MPI_Status stat;
					MPI_Recv ( &data_interior_quadrangles[id_quad][j], 1,
					           Give_MPI_data_typ<DTyp>::Give(), rank_send, 44+num, comm, &stat );
				}
			}
		}
	}


	// exterior neighbor
	if ( ug->Give_quadrangle ( id_quad )->Give_exists_exterior() ) {

		// exterior hexahedron
		id_hex = ug->Give_quadrangle ( id_quad )->Give_id_exterior_vol();
		rank_send = ug->Give_hexahedron ( id_hex )->Give_rank(); // for parallel
		if ( rank_send==my_rank || rank_rec==my_rank ) { // for parallel
			Nxhex = blockgrid->Give_Nx_hexahedron ( id_hex );
			Nyhex = blockgrid->Give_Ny_hexahedron ( id_hex );
			Nzhex = blockgrid->Give_Nz_hexahedron ( id_hex );

			// Find corresponding points in hexahedron
			for ( int c=0; c<8; ++c ) {
				idsp = ug->Give_hexahedron ( id_hex )->Give_id_corner ( ( dir3D_sons ) c );
				if ( idsp==id_SW ) c_SW = ( dir3D_sons ) c;
				else if ( idsp==id_SE ) c_SE = ( dir3D_sons ) c;
				else if ( idsp==id_NW ) c_NW = ( dir3D_sons ) c;
			}

			c_x = direction_from_to ( c_SW,c_SE );
			c_y = direction_from_to ( c_SW,c_NW );

			for ( int fc=0; fc<6; ++fc )
				if ( id_quad == ug->Give_hexahedron ( id_hex )->Give_id_quadrangle ( ( dir3D ) fc ) ) {
					quad = ( dir3D ) fc;

					// x - direction
					Nconst = Give_Nconst ( ( dir3D ) c_x, Nxhex, Nyhex, Nzhex );
					Ni     = Give_Ni (    ( dir3D ) c_x, Nxhex, Nyhex, Nzhex );

					// y - direction
					Nconst = Give_Nconst ( ( dir3D ) c_y, Nxhex, Nyhex, Nzhex ) + Nconst;
					Nj     = Give_Ni (    ( dir3D ) c_y, Nxhex, Nyhex, Nzhex );

					// empty - direction
					Nconst = Give_Ni      ( opposite3D ( ( dir3D ) fc ),Nxhex, Nyhex, Nzhex )
					         + Give_Nconst ( opposite3D ( ( dir3D ) fc ),Nxhex, Nyhex, Nzhex )
					         + Nconst;
				}

			if ( rank_send == rank_rec )    // for parallel

				for ( j=1; j<Nyquad; ++j )
					for ( i=1; i<Nxquad; ++i )
						data_exterior_quadrangles[id_quad][i+ ( Nxquad+1 ) *j] =
						    data_hexahedra[id_hex][i*Ni+j*Nj+Nconst];
			else {
				if ( my_rank==rank_send ) {
					DTyp* data = Parallel_buffer<DTyp>::Give_send ( ( Nxquad-1 ) * ( Nyquad-1 ) );
					for ( j=1; j<Nyquad; ++j )
						for ( i=1; i<Nxquad; ++i )
							data[ ( i-1 ) + ( Nxquad-1 ) * ( j-1 ) ] =
							    data_hexahedra[id_hex][i*Ni+j*Nj+Nconst];
					MPI_Send ( data, ( Nxquad-1 ) * ( Nyquad-1 ),
					           Give_MPI_data_typ<DTyp>::Give(), rank_rec, 55, comm );
				} else {
					if ( developer_version )
						if ( my_rank!=rank_rec )
							cout << "error30 in Update_quadrangle! " << endl;
					MPI_Status stat;
					MPI_Datatype columntype;
					MPI_Type_vector ( Nyquad-1, Nxquad-1, Nxquad+1, Give_MPI_data_typ<DTyp>::Give(),
					                  &columntype );
					MPI_Type_commit ( &columntype );
					MPI_Recv ( &data_exterior_quadrangles[id_quad][Nxquad+2], 1, columntype,
					           rank_send, 55, comm, &stat );
					MPI_Type_free ( &columntype );
				}
			}
		}

		// exterior face Wdir2D, Edir2D, Sdir2D, Ndir2D
		for ( num=0; num<4; ++num ) { // Iterate over all edges of quad
			edge = ( dir2D ) num;
			id_face = ug->Give_quadrangle ( id_quad )->Give_id_exterior_face ( edge ); // id of exterior face
			rank_send = ug->Give_quadrangle ( id_face )->Give_rank(); // for parallel
			if ( rank_send==my_rank || rank_rec==my_rank ) { // for parallel
				Nxface = blockgrid->Give_Nx_quadrangle ( id_face );
				Nyface = blockgrid->Give_Ny_quadrangle ( id_face );

				// coordinate points of face
				id_SW_face = ug->Give_quadrangle ( id_face )->Give_id_corner ( SWdir2D );
				id_SE_face = ug->Give_quadrangle ( id_face )->Give_id_corner ( SEdir2D );
				id_NW_face = ug->Give_quadrangle ( id_face )->Give_id_corner ( NWdir2D );
				id_NE_face = ug->Give_quadrangle ( id_face )->Give_id_corner ( NEdir2D );

				// Set reference corners of quadrangle
				if ( edge == Wdir2D ) {
					id_SoW = id_SW;
					id_NoE = id_NW;
				} else if ( edge == Edir2D ) {
					id_SoW = id_SE;
					id_NoE = id_NE;
				} else if ( edge == Sdir2D ) {
					id_SoW = id_SW;
					id_NoE = id_SE;
				} else {
					if ( developer_version )
						if ( edge != Ndir2D ) cout << "error31 in Update_quadrangle! " << endl;
					id_SoW = id_NW;
					id_NoE = id_NE;
				}

				// Compare corners of quad with corners of face
				if ( id_SW_face == id_SoW ) {
					if ( id_NW_face == id_NoE ) {
						Ni = Nj = Nxface+1;
						Nconst = 1;
					} else {
						if ( developer_version )
							if ( id_SE_face != id_NoE ) cout << "error32 in Update_quadrangle! " << endl;
						Ni = Nj = 1;
						Nconst = Nxface+1;
					}
				} else if ( id_SE_face == id_SoW ) {
					if ( id_SW_face == id_NoE ) {
						Ni = Nj = -1;
						Nconst = 2*Nxface+1;
					} else {
						if ( developer_version )
							if ( id_NE_face != id_NoE ) cout << "error33 in Update_quadrangle! " << endl;
						Ni = Nj = Nxface+1;
						Nconst = Nxface-1;
					}
				} else if ( id_NW_face == id_SoW ) {
					if ( id_NE_face == id_NoE ) {
						Ni = Nj = 1;
						Nconst = ( Nyface-1 ) * ( Nxface+1 );
					} else {
						if ( developer_version )
							if ( id_SW_face != id_NoE ) cout << "error34 in Update_quadrangle! " << endl;
						Ni = Nj = - ( Nxface+1 );
						Nconst = Nyface* ( Nxface+1 ) +1;
					}
				} else {
					if ( developer_version )
						if ( id_NE_face != id_SoW ) cout << "error35 in Update_quadrangle! " << endl;
					if ( id_SE_face == id_NoE ) {
						Ni = Nj = - ( Nxface+1 );
						Nconst = ( Nyface+1 ) * ( Nxface+1 )-2;
					} else {
						if ( developer_version )
							if ( id_NW_face != id_NoE ) cout << "error36 in Update_quadrangle! " << endl;
						Ni = Nj = -1;
						Nconst = Nyface* ( Nxface+1 )-1;
					}
				}

				if ( rank_send == rank_rec ) {   // for parallel
					if ( edge == Wdir2D ) {
						i = 0;

						for ( j=1; j<Nyquad; ++j )
							data_exterior_quadrangles[id_quad][i+ ( Nxquad+1 ) *j] =
							    data_quadrangles[id_face][j*Nj+Nconst];
					} else if ( edge == Edir2D ) {
						i = Nxquad;

						for ( j=1; j<Nyquad; ++j ) {
							data_exterior_quadrangles[id_quad][i+ ( Nxquad+1 ) *j] =
							    data_quadrangles[id_face][j*Nj+Nconst];
						}
					} else if ( edge == Sdir2D ) {
						j = 0;

						for ( i=1; i<Nxquad; ++i )
							data_exterior_quadrangles[id_quad][i+ ( Nxquad+1 ) *j] =
							    data_quadrangles[id_face][i*Ni+Nconst];
					} else {
						if ( developer_version )
							if ( edge != Ndir2D ) cout << "error37 in Update_quadrangle! " << endl;
						j = Nyquad;

						for ( i=1; i<Nxquad; ++i ) {
							data_exterior_quadrangles[id_quad][i+ ( Nxquad+1 ) *j] =
							    data_quadrangles[id_face][i*Ni+Nconst];
						}
					}
				} else {
					if ( my_rank==rank_send ) {
						DTyp* data = Parallel_buffer<DTyp>::Give_send ( Nxquad-1 );
						for ( i=1; i<Nxquad; ++i )
							data[i-1] =
							    data_quadrangles[id_face][i*Ni+Nconst];
						MPI_Send ( data, Nxquad-1,
						           Give_MPI_data_typ<DTyp>::Give(), rank_rec, 66+num, comm );
					} else {
						if ( developer_version )
							if ( my_rank!=rank_rec )
								cout << "error38 in Update_quadrangle! " << endl;
						if ( edge == Wdir2D || edge == Edir2D ) {
							MPI_Status stat;
							MPI_Datatype columntype;
							MPI_Type_vector ( Nyquad-1, 1, Nxquad+1,
							                  Give_MPI_data_typ<DTyp>::Give(), &columntype );
							MPI_Type_commit ( &columntype );
							if ( edge == Wdir2D )
								i = 0;
							else {
								if ( developer_version )
									if ( edge != Edir2D )
										cout << "error39 in Update_quadrangle! "
										     << endl;
								i = Nxquad;
							}
							MPI_Recv ( &data_exterior_quadrangles[id_quad][i+Nxquad+1], 1,
							           columntype, rank_send, 66+num, comm, &stat );
							MPI_Type_free ( &columntype );
						} else {
							MPI_Status stat;
							if ( edge == Sdir2D )
								j = 0;
							else {
								if ( developer_version )
									if ( edge != Ndir2D ) cout << "error40 in Update_quadrangle! "
										                           << endl;
								j = Nyquad;
							}
							MPI_Recv ( &data_exterior_quadrangles[id_quad][1+ ( Nxquad+1 ) *j], Nxquad-1,
							           Give_MPI_data_typ<DTyp>::Give(), rank_send, 66+num, comm, &stat );
						}
					}
				}
			}
		}

		// exterior edge SWdir2D, SEdir2D, NWdir2D, NEdir2D
		for ( num=0; num<4; ++num ) { // Iterate over all corners of quad
			point = ( dir2D_sons ) num;
			id_edge = ug->Give_quadrangle ( id_quad )->Give_id_exterior_edge ( point ); // id of exterior edge
			rank_send = ug->Give_edge ( id_edge )->Give_rank(); // for parallel
			if ( rank_send==my_rank || rank_rec==my_rank ) { // for parallel
				id_point = ug->Give_quadrangle ( id_quad )->Give_id_corner ( point ); // id of corner

				Nxedge = blockgrid->Give_Nx_edge ( id_edge );

				// coordinate points of edge
				id_L = ug->Give_edge ( id_edge )->Give_id_corner ( Ldir1D );
				id_R = ug->Give_edge ( id_edge )->Give_id_corner ( Rdir1D );

				if ( rank_send == rank_rec ) {   // for parallel
					// Compare corners of quad with corners of edge
					if ( point == SWdir2D ) {
						if ( id_point == id_L )
							data_exterior_quadrangles[id_quad][0] =
							    data_edges[id_edge][1];
						else {
							if ( developer_version )
								if ( id_point != id_R ) cout << "error41 in Update_quadrangle! " << endl;
							data_exterior_quadrangles[id_quad][0] =
							    data_edges[id_edge][Nxedge-1];
						}
					} else if ( point == SEdir2D ) {
						if ( id_point == id_L )
							data_exterior_quadrangles[id_quad][Nxquad] =
							    data_edges[id_edge][1];
						else {
							if ( developer_version )
								if ( id_point != id_R ) cout << "error42 in Update_quadrangle! " << endl;
							data_exterior_quadrangles[id_quad][Nxquad] =
							    data_edges[id_edge][Nxedge-1];
						}
					} else if ( point == NWdir2D ) {
						if ( id_point == id_L )
							data_exterior_quadrangles[id_quad][ ( Nxquad+1 ) *Nyquad] =
							    data_edges[id_edge][1];
						else {
							if ( developer_version )
								if ( id_point != id_R ) cout << "error43 in Update_quadrangle! " << endl;
							data_exterior_quadrangles[id_quad][ ( Nxquad+1 ) *Nyquad] =
							    data_edges[id_edge][Nxedge-1];
						}
					} else {
						if ( developer_version )
							if ( point != NEdir2D ) cout << "error44 in Update_quadrangle! " << endl;
						if ( id_point == id_L )
							data_exterior_quadrangles[id_quad][ ( Nxquad+1 ) * ( Nyquad+1 )-1] =
							    data_edges[id_edge][1];
						else {
							if ( developer_version )
								if ( id_point != id_R ) cout << "error45 in Update_quadrangle! " << endl;
							data_exterior_quadrangles[id_quad][ ( Nxquad+1 ) * ( Nyquad+1 )-1] =
							    data_edges[id_edge][Nxedge-1];
						}
					}
				} else {
					if ( my_rank==rank_send ) {
						if ( id_point == id_L ) i=1;
						else {
							if ( developer_version )
								if ( id_point != id_R )
									cout << "error46 in Update_quadrangle! " << endl;
							i=Nxedge-1;
						}
						MPI_Send ( &data_edges[id_edge][i], 1,
						           Give_MPI_data_typ<DTyp>::Give(), rank_rec, 77+num, comm );
					} else {
						if ( developer_version )
							if ( my_rank!=rank_rec )
								cout << "error47 in Update_quadrangle! " << endl;
						if ( point == SWdir2D ) j=0;
						else if ( point == SEdir2D ) j=Nxquad;
						else if ( point == NWdir2D ) j= ( Nxquad+1 ) *Nyquad;
						else {
							if ( developer_version )
								if ( point != NEdir2D ) cout << "error48 in Update_quadrangle! " << endl;
							j= ( Nxquad+1 ) * ( Nyquad+1 )-1;
						}
						MPI_Status stat;
						MPI_Recv ( &data_exterior_quadrangles[id_quad][j], 1,
						           Give_MPI_data_typ<DTyp>::Give(), rank_send, 77+num, comm, &stat );
					}
				}
			}
		}
	}
}
/*
template <>
template <>
void Variable<double>::Update<quadrangleEl>(int id_quad) const{
  Update_quadrangle<double>(id_quad, blockgrid, data_hexahedra, data_quadrangles,
                            data_edges, data_points, data_interior_quadrangles,
                            data_exterior_quadrangles);
}

template <>
template <>
void Variable<std::complex<double> >::Update<quadrangleEl>(int id_quad) const{
  Update_quadrangle<std::complex<double> >(id_quad, blockgrid, data_hexahedra,
                                           data_quadrangles, data_edges, data_points,
                                           data_interior_quadrangles,
                                           data_exterior_quadrangles);
}
*/
template <class DTyp>
void Update_edge ( int id_edge, Blockgrid* blockgrid, DTyp** data_hexahedra,
                   DTyp** data_quadrangles, DTyp** data_edges, DTyp** data_points,
                   DTyp*** data_neighbor_edges, int my_rank, MPI_Comm comm )
{
	int i, j, ed, c;
	int id_point, id_ed, id_corner;
	int Nxedge, Nxhex, Nyhex, Nzhex, Nxquad, Nyquad, Nxed;
	int Nj, Nconst;
	int num_neighbor, neighbor_id;
	Unstructured_grid* ug = blockgrid->Give_unstructured_grid();
	elementTyp neighbor_typ;
	int W_corner, E_corner;
	dir3D_sons W_corner_ed_cell, E_corner_ed_cell;
	dir2D_sons W_corner_edge, E_corner_edge, corner;
	dir2D edge;
	Edges_cell ed_cell;

	int rank_send, rank_rec;  // for parallel

	Nxedge = blockgrid->Give_Nx_edge ( id_edge );

	rank_rec = ug->Give_edge ( id_edge )->Give_rank(); // for parallel

	// points corner_W, corner_E
	id_point = ug->Give_edge ( id_edge )->Give_id_corner_W();
	rank_send = ug->Give_point ( id_point )->Give_rank(); // for parallel
	if ( rank_send==my_rank || rank_rec==my_rank ) { // for parallel
		if ( rank_send == rank_rec )    // for parallel
			data_edges[id_edge][0] = data_points[id_point][0];
		else {
			if ( my_rank==rank_send )
				MPI_Send ( data_points[id_point], 1,
				           Give_MPI_data_typ<DTyp>::Give(), rank_rec, 10, comm );
			else {
				if ( developer_version )
					if ( my_rank!=rank_rec )
						cout << "error1 in Update_edge !" << endl;
				MPI_Status stat;
				MPI_Recv ( data_edges[id_edge], 1, Give_MPI_data_typ<DTyp>::Give(),
				           rank_send, 10, comm, &stat );
			}
		}
	}

	id_point = ug->Give_edge ( id_edge )->Give_id_corner_E();
	rank_send = ug->Give_point ( id_point )->Give_rank(); // for parallel
	if ( rank_send==my_rank || rank_rec==my_rank ) { // for parallel
		if ( rank_send == rank_rec )    // for parallel
			data_edges[id_edge][Nxedge] = data_points[id_point][0];
		else {
			if ( my_rank==rank_send )
				MPI_Send ( data_points[id_point], 1,
				           Give_MPI_data_typ<DTyp>::Give(), rank_rec, 20, comm );
			else {
				if ( developer_version )
					if ( my_rank!=rank_rec )
						cout << "error2 in Update_edge !" << endl;
				MPI_Status stat;
				MPI_Recv ( &data_edges[id_edge][Nxedge], 1, Give_MPI_data_typ<DTyp>::Give(),
				           rank_send, 20, comm, &stat );
			}
		}
	}

	// unstructured data
	num_neighbor = ug->Give_edge ( id_edge )->Give_number_neighbors();
	W_corner = ug->Give_edge ( id_edge )->Give_id_corner_W();
	E_corner = ug->Give_edge ( id_edge )->Give_id_corner_E();

	// data of "middle" neighbors
	for ( i=0; i<num_neighbor; ++i ) {
		neighbor_id = ug->Give_edge ( id_edge )->Give_id_M_el ( i );
		neighbor_typ = ug->Give_edge ( id_edge )->Give_typ_M_el ( i );
		if ( neighbor_typ == hexahedronEl ) {
			rank_send = ug->Give_hexahedron ( neighbor_id )->Give_rank(); // for parallel
			if ( rank_send==my_rank || rank_rec==my_rank ) { // for parallel
				Nxhex = blockgrid->Give_Nx_hexahedron ( neighbor_id );
				Nyhex = blockgrid->Give_Ny_hexahedron ( neighbor_id );
				Nzhex = blockgrid->Give_Nz_hexahedron ( neighbor_id );
				for ( ed=0; ed<12; ++ed ) {
					id_ed = ug->Give_hexahedron ( neighbor_id )->Give_id_edge ( ( Edges_cell ) ed );
					if ( id_ed == id_edge ) ed_cell = ( Edges_cell ) ed;
				}
				W_corner_ed_cell = Transform ( ed_cell, Ldir1D );
				E_corner_ed_cell = Transform ( ed_cell, Rdir1D );

				if ( ug->Give_hexahedron ( neighbor_id )->Give_id_corner ( W_corner_ed_cell ) <
				     ug->Give_hexahedron ( neighbor_id )->Give_id_corner ( E_corner_ed_cell ) ) {
					Nj = Give_Nj ( W_corner_ed_cell, E_corner_ed_cell, Nxhex, Nyhex, Nzhex );
					Nconst = Give_Nconst ( W_corner_ed_cell, E_corner_ed_cell, Nxhex, Nyhex, Nzhex );
				} else {
					Nj = Give_Nj ( E_corner_ed_cell, W_corner_ed_cell, Nxhex, Nyhex, Nzhex );
					Nconst = Give_Nconst ( E_corner_ed_cell, W_corner_ed_cell, Nxhex, Nyhex, Nzhex );
				}
				if ( rank_send == rank_rec )    // for parallel

					for ( j=1; j<Nxedge; ++j )
						data_neighbor_edges[id_edge][i][j] =
						    data_hexahedra[neighbor_id][ ( j-1 ) *Nj+Nconst];
				else {
					if ( my_rank==rank_send ) {
						DTyp* data = Parallel_buffer<DTyp>::Give_send ( Nxedge-1 );
						for ( j=1; j<Nxedge; ++j )
							data[j-1] =
							    data_hexahedra[neighbor_id][ ( j-1 ) *Nj+Nconst];
						MPI_Send ( data, Nxedge-1,
						           Give_MPI_data_typ<DTyp>::Give(), rank_rec, 11+num_neighbor, comm );
					} else {
						if ( developer_version )
							if ( my_rank!=rank_rec )
								cout << "error3 in Update_edge !" << endl;
						MPI_Status stat;
						MPI_Recv ( &data_neighbor_edges[id_edge][i][1], Nxedge-1,
						           Give_MPI_data_typ<DTyp>::Give(), rank_send, 11+num_neighbor, comm, &stat );
					}
				}
			}
		} else {
			if ( developer_version )
				if ( neighbor_typ != quadrangleEl ) cout << "error4 in Update_edge !" << endl;
			rank_send = ug->Give_quadrangle ( neighbor_id )->Give_rank(); // for parallel
			if ( rank_send==my_rank || rank_rec==my_rank ) { // for parallel
				Nxquad = blockgrid->Give_Nx_quadrangle ( neighbor_id );
				Nyquad = blockgrid->Give_Ny_quadrangle ( neighbor_id );
				for ( ed=0; ed<4; ++ed ) {
					id_ed = ug->Give_quadrangle ( neighbor_id )->Give_id_edge ( ( dir2D ) ed );
					if ( id_ed == id_edge ) edge = ( dir2D ) ed;
				}
				W_corner_edge = Transform ( edge, Ldir1D );
				E_corner_edge = Transform ( edge, Rdir1D );
				if ( ug->Give_quadrangle ( neighbor_id )->Give_id_corner ( W_corner_edge ) <
				     ug->Give_quadrangle ( neighbor_id )->Give_id_corner ( E_corner_edge ) ) {
					Nj = Give_Nj ( W_corner_edge, E_corner_edge, Nxquad, Nyquad );
					Nconst = Give_Nconst ( W_corner_edge, E_corner_edge, Nxquad, Nyquad );
				} else {
					Nj = Give_Nj ( E_corner_edge, W_corner_edge, Nxquad, Nyquad );
					Nconst = Give_Nconst ( E_corner_edge, W_corner_edge, Nxquad, Nyquad );
				}
				if ( rank_send == rank_rec )    // for parallel

					for ( j=1; j<Nxedge; ++j )
						data_neighbor_edges[id_edge][i][j] =
						    data_quadrangles[neighbor_id][ ( j-1 ) *Nj+Nconst];
				else {
					if ( my_rank==rank_send ) {
						DTyp* data = Parallel_buffer<DTyp>::Give_send ( Nxedge-1 );
						for ( j=1; j<Nxedge; ++j )
							data[j-1] =
							    data_quadrangles[neighbor_id][ ( j-1 ) *Nj+Nconst];
						MPI_Send ( data, Nxedge-1,
						           Give_MPI_data_typ<DTyp>::Give(), rank_rec, 11+num_neighbor, comm );
					} else {
						if ( developer_version )
							if ( my_rank!=rank_rec )
								cout << "error5 in Update_edge !" << endl;
						MPI_Status stat;
						MPI_Recv ( &data_neighbor_edges[id_edge][i][1], Nxedge-1,
						           Give_MPI_data_typ<DTyp>::Give(), rank_send, 11+num_neighbor, comm, &stat );
					}
				}
			}
		}
	}

	// data of "west" neighbors
	for ( i=0; i<num_neighbor; ++i ) {
		neighbor_id = ug->Give_edge ( id_edge )->Give_id_W_el ( i );
		neighbor_typ = ug->Give_edge ( id_edge )->Give_typ_W_el ( i );
		if ( neighbor_typ == quadrangleEl ) {
			rank_send = ug->Give_quadrangle ( neighbor_id )->Give_rank(); // for parallel
			if ( rank_send==my_rank || rank_rec==my_rank ) { // for parallel
				Nxquad = blockgrid->Give_Nx_quadrangle ( neighbor_id );
				Nyquad = blockgrid->Give_Ny_quadrangle ( neighbor_id );
				for ( c=0; c<4; ++c ) {
					id_corner = ug->Give_quadrangle ( neighbor_id )->Give_id_corner ( ( dir2D_sons ) c );
					if ( id_corner == W_corner ) corner = ( dir2D_sons ) c;
				}
				if ( rank_send == rank_rec ) {   // for parallel
					if ( corner == SWdir2D ) data_neighbor_edges[id_edge][i][0] =
						    data_quadrangles[neighbor_id][Nxquad+2];
					else if ( corner == SEdir2D ) data_neighbor_edges[id_edge][i][0] =
						    data_quadrangles[neighbor_id][2*Nxquad];
					else if ( corner == NWdir2D ) data_neighbor_edges[id_edge][i][0] =
						    data_quadrangles[neighbor_id][ ( Nyquad-1 ) * ( Nxquad+1 ) +1];
					else {
						if ( developer_version )
							if ( corner != NEdir2D ) cout << "error6 in Update_edge !" << endl;
						data_neighbor_edges[id_edge][i][0] =
						    data_quadrangles[neighbor_id][Nyquad* ( Nxquad+1 )-2];
					}
				} else {
					if ( my_rank==rank_send ) {
						if ( corner == SWdir2D ) j=Nxquad+2;
						else if ( corner == SEdir2D ) j=2*Nxquad;
						else if ( corner == NWdir2D ) j= ( Nyquad-1 ) * ( Nxquad+1 ) +1;
						else {
							if ( developer_version )
								if ( corner != NEdir2D ) cout << "error7 in Update_edge !" << endl;
							j=Nyquad* ( Nxquad+1 )-2;
						}
						MPI_Send ( &data_quadrangles[neighbor_id][j], 1,
						           Give_MPI_data_typ<DTyp>::Give(), rank_rec, 22+num_neighbor, comm );
					} else {
						if ( developer_version )
							if ( my_rank!=rank_rec )
								cout << "error8 in Update_edge !" << endl;
						MPI_Status stat;
						MPI_Recv ( data_neighbor_edges[id_edge][i], 1,
						           Give_MPI_data_typ<DTyp>::Give(), rank_send, 22+num_neighbor, comm, &stat );
					}
				}
			}
		} else {
			if ( developer_version )
				if ( neighbor_typ != edgeEl ) cout << "error9 in Update_edge !" << endl;
			rank_send = ug->Give_edge ( neighbor_id )->Give_rank(); // for parallel
			if ( rank_send==my_rank || rank_rec==my_rank ) { // for parallel
				Nxed = blockgrid->Give_Nx_edge ( neighbor_id );
				if ( rank_send == rank_rec ) {   // for parallel
					if ( ug->Give_edge ( neighbor_id )->Give_id_corner ( Ldir1D ) == W_corner )
						data_neighbor_edges[id_edge][i][0] =
						    data_edges[neighbor_id][1];
					else {
						if ( developer_version )
							if ( ug->Give_edge ( neighbor_id )->Give_id_corner ( Rdir1D ) != W_corner )
								cout << "error10 in Update_edge !" << endl;
						data_neighbor_edges[id_edge][i][0] =
						    data_edges[neighbor_id][Nxed-1];
					}
				} else {
					if ( my_rank==rank_send ) {
						if ( ug->Give_edge ( neighbor_id )->Give_id_corner ( Ldir1D ) == W_corner ) j=1;
						else {
							if ( developer_version )
								if ( ug->Give_edge ( neighbor_id )->Give_id_corner ( Rdir1D ) != W_corner )
									cout << "error11 in Update_edge !" << endl;
							j=Nxed-1;
						}
						MPI_Send ( &data_edges[neighbor_id][j], 1,
						           Give_MPI_data_typ<DTyp>::Give(), rank_rec, 22+num_neighbor, comm );
					} else {
						if ( developer_version )
							if ( my_rank!=rank_rec )
								cout << "error12 in Update_edge !" << endl;
						MPI_Status stat;
						MPI_Recv ( data_neighbor_edges[id_edge][i], 1,
						           Give_MPI_data_typ<DTyp>::Give(), rank_send, 22+num_neighbor, comm, &stat );
					}
				}
			}
		}
	}

	// data of "east" neighbors
	for ( i=0; i<num_neighbor; ++i ) {
		neighbor_id = ug->Give_edge ( id_edge )->Give_id_E_el ( i );
		neighbor_typ = ug->Give_edge ( id_edge )->Give_typ_E_el ( i );
		if ( neighbor_typ == quadrangleEl ) {
			rank_send = ug->Give_quadrangle ( neighbor_id )->Give_rank(); // for parallel
			if ( rank_send==my_rank || rank_rec==my_rank ) { // for parallel
				Nxquad = blockgrid->Give_Nx_quadrangle ( neighbor_id );
				Nyquad = blockgrid->Give_Ny_quadrangle ( neighbor_id );
				for ( c=0; c<4; ++c ) {
					id_corner = ug->Give_quadrangle ( neighbor_id )->Give_id_corner ( ( dir2D_sons ) c );
					if ( id_corner == E_corner ) corner = ( dir2D_sons ) c;
				}
				if ( rank_send == rank_rec ) {   // for parallel
					if ( corner == SWdir2D ) data_neighbor_edges[id_edge][i][Nxedge] =
						    data_quadrangles[neighbor_id][Nxquad+2];
					else if ( corner == SEdir2D ) data_neighbor_edges[id_edge][i][Nxedge] =
						    data_quadrangles[neighbor_id][2*Nxquad];
					else if ( corner == NWdir2D ) data_neighbor_edges[id_edge][i][Nxedge] =
						    data_quadrangles[neighbor_id][ ( Nyquad-1 ) * ( Nxquad+1 ) +1];
					else {
						if ( developer_version )
							if ( corner != NEdir2D ) cout << "error13 in Update_edge !" << endl;
						data_neighbor_edges[id_edge][i][Nxedge] =
						    data_quadrangles[neighbor_id][Nyquad* ( Nxquad+1 )-2];
					}
				} else {
					if ( my_rank==rank_send ) {
						if ( corner == SWdir2D ) j=Nxquad+2;
						else if ( corner == SEdir2D ) j=2*Nxquad;
						else if ( corner == NWdir2D ) j= ( Nyquad-1 ) * ( Nxquad+1 ) +1;
						else {
							if ( developer_version )
								if ( corner != NEdir2D ) cout << "error14 in Update_edge !" << endl;
							j=Nyquad* ( Nxquad+1 )-2;
						}
						MPI_Send ( &data_quadrangles[neighbor_id][j], 1,
						           Give_MPI_data_typ<DTyp>::Give(), rank_rec, 33+num_neighbor, comm );
					} else {
						if ( developer_version )
							if ( my_rank!=rank_rec )
								cout << "error15 in Update_edge !" << endl;
						MPI_Status stat;
						MPI_Recv ( &data_neighbor_edges[id_edge][i][Nxedge], 1,
						           Give_MPI_data_typ<DTyp>::Give(), rank_send, 33+num_neighbor, comm, &stat );
					}
				}
			}
		} else {
			if ( developer_version )
				if ( neighbor_typ != edgeEl ) cout << "error16 in Update_edge !" << endl;
			rank_send = ug->Give_edge ( neighbor_id )->Give_rank(); // for parallel
			if ( rank_send==my_rank || rank_rec==my_rank ) { // for parallel
				Nxed = blockgrid->Give_Nx_edge ( neighbor_id );
				if ( rank_send == rank_rec ) {   // for parallel
					if ( ug->Give_edge ( neighbor_id )->Give_id_corner ( Ldir1D ) == E_corner )
						data_neighbor_edges[id_edge][i][Nxedge] =
						    data_edges[neighbor_id][1];
					else {
						if ( developer_version )
							if ( ug->Give_edge ( neighbor_id )->Give_id_corner ( Rdir1D ) != E_corner )
								cout << "error17 in Update_edge !" << endl;
						data_neighbor_edges[id_edge][i][Nxedge] =
						    data_edges[neighbor_id][Nxed-1];
					}
				} else {
					if ( my_rank==rank_send ) {
						if ( ug->Give_edge ( neighbor_id )->Give_id_corner ( Ldir1D ) == E_corner ) j=1;
						else {
							if ( developer_version )
								if ( ug->Give_edge ( neighbor_id )->Give_id_corner ( Rdir1D ) != E_corner )
									cout << "error18 in Update_edge !" << endl;
							j=Nxed-1;
						}
						MPI_Send ( &data_edges[neighbor_id][j], 1,
						           Give_MPI_data_typ<DTyp>::Give(), rank_rec, 33+num_neighbor, comm );
					} else {
						if ( developer_version )
							if ( my_rank!=rank_rec )
								cout << "error19 in Update_edge !" << endl;
						MPI_Status stat;
						MPI_Recv ( &data_neighbor_edges[id_edge][i][Nxedge], 1,
						           Give_MPI_data_typ<DTyp>::Give(), rank_send, 33+num_neighbor, comm, &stat );
					}
				}
			}
		}
	}
}
/*
template <>
template <>
void Variable<double>::Update<edgeEl>(int id_edge) const{
  Update_edge<double>(id_edge, blockgrid, data_hexahedra,
                      data_quadrangles, data_edges, data_points, data_neighbor_edges);
}

template <>
template <>
void Variable<std::complex<double> >::Update<edgeEl>(int id_edge) const{
  Update_edge<std::complex<double> >(id_edge, blockgrid, data_hexahedra,
                                     data_quadrangles, data_edges, data_points,
                                     data_neighbor_edges);
}
*/
template <class DTyp>
void Update_point ( int id_point, Blockgrid* blockgrid, DTyp** data_hexahedra,
                    DTyp** data_quadrangles, DTyp** data_edges, DTyp** data_points,
                    DTyp**  data_neighbor_points, int my_rank, MPI_Comm comm )
{

	int id_corner, j;
	int Nxhex, Nyhex, Nzhex, Nxquad, Nyquad, Nxedge;
	int Nconst;
	int num_neighbor, neighbor_id;
	Unstructured_grid* ug = blockgrid->Give_unstructured_grid();
	elementTyp neighbor_typ;
	dir3D_sons corner;
	dir2D_sons corn;

	int rank_send, rank_rec;  // for parallel

	rank_rec = ug->Give_point ( id_point )->Give_rank(); // for parallel

	// unstructured data
	num_neighbor = ug->Give_point ( id_point )->Give_number_neighbors();
	for ( int i=0; i<num_neighbor; ++i ) {
		neighbor_id = ug->Give_point ( id_point )->Give_id_el_neighbor ( i );
		neighbor_typ = ug->Give_point ( id_point )->Give_typ_el_neighbor ( i );
		if ( neighbor_typ == hexahedronEl ) {
			rank_send = ug->Give_hexahedron ( neighbor_id )->Give_rank(); // for parallel
			if ( rank_send==my_rank || rank_rec==my_rank ) { // for parallel
				Nxhex = blockgrid->Give_Nx_hexahedron ( neighbor_id );
				Nyhex = blockgrid->Give_Ny_hexahedron ( neighbor_id );
				Nzhex = blockgrid->Give_Nz_hexahedron ( neighbor_id );
				for ( int c=0; c<8; ++c ) {
					id_corner = ug->Give_hexahedron ( neighbor_id )->Give_id_corner ( ( dir3D_sons ) c );
					if ( id_corner == id_point ) corner = ( dir3D_sons ) c;
				}
				Nconst = Give_Nconst ( corner, Nxhex, Nyhex, Nzhex );
				if ( rank_send == rank_rec )    // for parallel
					data_neighbor_points[id_point][i] = data_hexahedra[neighbor_id][Nconst];
				else {
					if ( my_rank==rank_send )
						MPI_Send ( &data_hexahedra[neighbor_id][Nconst], 1,
						           Give_MPI_data_typ<DTyp>::Give(), rank_rec, 10+num_neighbor, comm );
					else {
						if ( developer_version )
							if ( my_rank!=rank_rec )
								cout << "error1 in Update_point !" << endl;
						MPI_Status stat;
						MPI_Recv ( &data_neighbor_points[id_point][i], 1,
						           Give_MPI_data_typ<DTyp>::Give(), rank_send, 10+num_neighbor, comm, &stat );
					}
				}
			}
		} else if ( neighbor_typ == quadrangleEl ) {
			rank_send = ug->Give_quadrangle ( neighbor_id )->Give_rank(); // for parallel
			if ( rank_send==my_rank || rank_rec==my_rank ) { // for parallel
				Nxquad = blockgrid->Give_Nx_quadrangle ( neighbor_id );
				Nyquad = blockgrid->Give_Ny_quadrangle ( neighbor_id );
				for ( int c=0; c<4; ++c ) {
					id_corner = ug->Give_quadrangle ( neighbor_id )->Give_id_corner ( ( dir2D_sons ) c );
					if ( id_corner == id_point ) corn = ( dir2D_sons ) c;
				}
				Nconst = Give_Nconst ( corn, Nxquad, Nyquad );
				if ( rank_send == rank_rec )    // for parallel
					data_neighbor_points[id_point][i] = data_quadrangles[neighbor_id][Nconst];
				else {
					if ( my_rank==rank_send )
						MPI_Send ( &data_quadrangles[neighbor_id][Nconst], 1,
						           Give_MPI_data_typ<DTyp>::Give(), rank_rec, 10+num_neighbor, comm );
					else {
						if ( developer_version )
							if ( my_rank!=rank_rec )
								cout << "error2 in Update_point !" << endl;
						MPI_Status stat;
						MPI_Recv ( &data_neighbor_points[id_point][i], 1,
						           Give_MPI_data_typ<DTyp>::Give(), rank_send, 10+num_neighbor, comm, &stat );
					}
				}
			}
		} else {
			if ( developer_version )
				if ( neighbor_typ != edgeEl ) cout << "error3 in Update_point !" << endl;
			rank_send = ug->Give_edge ( neighbor_id )->Give_rank(); // for parallel
			if ( rank_send==my_rank || rank_rec==my_rank ) { // for parallel
				Nxedge = blockgrid->Give_Nx_edge ( neighbor_id );
				if ( rank_send == rank_rec ) {   // for parallel
					if ( ug->Give_edge ( neighbor_id )->Give_id_corner ( Ldir1D ) == id_point )
						data_neighbor_points[id_point][i] =
						    data_edges[neighbor_id][1];
					else {
						if ( developer_version )
							if ( ug->Give_edge ( neighbor_id )->Give_id_corner ( Rdir1D ) != id_point )
								cout << "error4 in Update_point !" << endl;
						data_neighbor_points[id_point][i] =
						    data_edges[neighbor_id][Nxedge-1];
					}
				} else {
					if ( my_rank==rank_send ) {
						if ( ug->Give_edge ( neighbor_id )->Give_id_corner ( Ldir1D ) == id_point ) j=1;
						else {
							if ( developer_version )
								if ( ug->Give_edge ( neighbor_id )->Give_id_corner ( Rdir1D ) != id_point )
									cout << "error5 in Update_point !" << endl;
							j=Nxedge-1;
						}
						MPI_Send ( &data_edges[neighbor_id][j], 1,
						           Give_MPI_data_typ<DTyp>::Give(), rank_rec, 10+num_neighbor, comm );
					} else {
						if ( developer_version )
							if ( my_rank!=rank_rec )
								cout << "error6 in Update_point !" << endl;
						MPI_Status stat;
						MPI_Recv ( &data_neighbor_points[id_point][i], 1,
						           Give_MPI_data_typ<DTyp>::Give(), rank_send, 10+num_neighbor, comm, &stat );
					}
				}
			}
		}
	}
}
/*
template <>
template <>
void Variable<double>::Update<pointEl>(int id_point) const{
  Update_point<double>(id_point, blockgrid, data_hexahedra,
                       data_quadrangles, data_edges, data_points, data_neighbor_points);
}

template <>
template <>
void Variable<std::complex<double> >::Update<pointEl>(int id_point) const{
  Update_point<std::complex<double> >(id_point, blockgrid, data_hexahedra,
                                      data_quadrangles, data_edges, data_points,
                                      data_neighbor_points);
}
*/

//----------------------------------------------------------------------


#endif // VARIABLE_CC_H