/**********************************************************************************
 * 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.
 **********************************************************************************/
 
// ------------------------------------------------------------
//
// cellvar.h
//
// ------------------------------------------------------------

#ifndef CE_VA_H
#define CE_VA_H

#ifdef _OPENMP
#include <omp.h>
#endif
#include "../grid/compose_grid.h"


//////////////////////////////////////////////////////////////
// 1. interpolation: variable -> cell variable
// 2. cell variable
// 2.1. definition
// 2.2. product_cell, L_infty_cell
// 3. Implementierungen
//////////////////////////////////////////////////////////////

//////////////////////////////////////////////////////////////
// 1. interpolation: variable -> cell variable
//////////////////////////////////////////////////////////////

template <class A>
class Expr_interpolant_point_to_cell : public Expr<Expr_interpolant_point_to_cell< A > >  {
 private:
  const A& a_;

 public:
  Expr_interpolant_point_to_cell(const Expr<A>& a) : a_(a) {}
    
  typedef typename A::Result Result;
                    
  bool totalCalcNotPossible()    const { return true; }	
                    
  Result Give_cell_hexahedra(params_in_cell) const {
    return 0.125 * (a_.template Give_data<hexahedronEl>(id_hex,i  ,j  ,k  ,Nx,Ny) +
		    a_.template Give_data<hexahedronEl>(id_hex,i+1,j  ,k  ,Nx,Ny) +
		    a_.template Give_data<hexahedronEl>(id_hex,i  ,j+1,k  ,Nx,Ny) +
		    a_.template Give_data<hexahedronEl>(id_hex,i+1,j+1,k  ,Nx,Ny) +
		    a_.template Give_data<hexahedronEl>(id_hex,i  ,j  ,k+1,Nx,Ny) +
		    a_.template Give_data<hexahedronEl>(id_hex,i+1,j  ,k+1,Nx,Ny) +
		    a_.template Give_data<hexahedronEl>(id_hex,i  ,j+1,k+1,Nx,Ny) +
		    a_.template Give_data<hexahedronEl>(id_hex,i+1,j+1,k+1,Nx,Ny));
  }
  
  Result Give_matrix_hexahedra(params_in_loc_matrix) const {
    return 0.125 * (a_.template Give_data<hexahedronEl>(id_hex,i  ,j  ,k  ,Nx,Ny) +
		    a_.template Give_data<hexahedronEl>(id_hex,i+1,j  ,k  ,Nx,Ny) +
		    a_.template Give_data<hexahedronEl>(id_hex,i  ,j+1,k  ,Nx,Ny) +
		    a_.template Give_data<hexahedronEl>(id_hex,i+1,j+1,k  ,Nx,Ny) +
		    a_.template Give_data<hexahedronEl>(id_hex,i  ,j  ,k+1,Nx,Ny) +
		    a_.template Give_data<hexahedronEl>(id_hex,i+1,j  ,k+1,Nx,Ny) +
		    a_.template Give_data<hexahedronEl>(id_hex,i  ,j+1,k+1,Nx,Ny) +
		    a_.template Give_data<hexahedronEl>(id_hex,i+1,j+1,k+1,Nx,Ny));
  }

  template <elementTyp TYP_EL>
  void Update(int id) const {
    a_.template Update<TYP_EL>(id);
  }
                  
  double Give_fromTotal(int i) const  { assert(false); return 0.0; }	
};

/** \addtogroup InterpolationOperators **/
/* @{ */  
/**
 * interpolates from point data to cell data
 **/
template <class A>
inline Expr_interpolant_point_to_cell<A> Cell_interpolation(const Expr<A>& ao) {
	const A& a ( ao );  
	Unstructured_grid* ug = a.Give_blockgrid()->Give_unstructured_grid();
  	for(int id = 0;id < ug->Give_number_hexahedra();++id ) {
		a.template Update<hexahedronEl> ( id );
	}
  return Expr_interpolant_point_to_cell<A>(ao);
}

/* @} */



/*
template <class DTyp>
inline Expr_interpolant_point_to_cell<Variable<DTyp> > Cell_interpolation(Variable<DTyp>& v) {
  v.UpdateHexahedra();
  return Expr_interpolant_point_to_cell<Variable<DTyp> >(v);
}
*/

//////////////////////////////////////////////////////
// 2. cell variable
//////////////////////////////////////////////////////
// 2.1. definition
//////////////////////////////////////////////////////////////

template <class DTyp, class Vari>
class VariableVector;


/** \addtogroup ExpressionTemplates **/
/* @{ */ 
template <class DTyp>
class Cell_variable : public Expr<Cell_variable< DTyp > > {

  friend  VariableVector<DTyp, Cell_variable<DTyp> >;

public:

  typedef DTyp Result;


  Cell_variable(Blockgrid& blockgrid_);
  Blockgrid* getBlockgrid() { return blockgrid; }   // WORK NOW  das sollte mal weg!!
		
  ~Cell_variable();

 // template <elementTyp TYP_EL>
 // void Update(int id) const {}

  Blockgrid* Give_blockgrid()  const { return blockgrid; };
  Unstructured_grid* Give_Ug() const { return ug; };
    
  void Print_VTK(std::ostream& Datei, double (*convert)(DTyp x), double stretch_z = 1., std::string title = "myData", Unstructured_grid_Marker * marker = NULL);
  void Print_VTK(std::ostream& Datei, double stretch_z = 1., std::string title = "myData",  Unstructured_grid_Marker * marker = NULL);

  void QPrint_VTK(QString DateiName, double (*convert)(DTyp x), double stretch_z = 1., QString title = "myData", Unstructured_grid_Marker * marker = NULL);
  void QPrint_VTK(QString DateiName, complex<double> (*convert)(DTyp x), double stretch_z = 1., QString title = "myData", Unstructured_grid_Marker * marker = NULL);
  void QPrint_VTK(QString DateiName, double stretch_z = 1., QString title = "myData",  Unstructured_grid_Marker * marker = NULL);

  inline DTyp Give_cell_hexahedra(params_in_cell)         const  { return data_cell[id_hex][ind_cell]; }
  inline DTyp Give_matrix_hexahedra(params_in_loc_matrix) const  { return data_cell[id_hex][ind_cell]; }
  inline DTyp Give_fromTotal(int i)                       const  { return dataTotal[i]; }

  /**
   * Man wird dies im Wesentlichen fuer die locale Steifigeitsmatrix von Helm anwenden
   * Dann wird das Volumen pro Zelle berechnet
   * @param localStiffnessMatrix das sollte locale Steifigeitsmatrix sein
   **/
  template <class LocalMat>
  void IntegrateOnCell(LocalMat& localStiffnessMatrix);
  
  template <class A>
  void operator=(const Expr<A>& a);
  void operator=(DTyp x);
  void operator=(const Cell_variable<DTyp>& varRight);

  void operator=(const DTyp_Restriction<DTyp>& a);

  template <class A>
  void operator=(const Expr_Restriction<A>& a);

  
  static DTyp product(Cell_variable& a, Cell_variable& b) { return product_cell(a,b,a.Give_blockgrid()->Give_unstructured_grid()->Give_all_points()); }
  
  // copy data from hex_a to hex_b (=)
  void Copy_invert_z(int hex_b,  Cell_variable<DTyp>& a, int hex_a);
  void Copy_invert_z(int hex_b_start, Cell_variable<DTyp>& a, int hex_a_start, int numBlocks);
  // copy and add data from hex_a to hex_b (+=)
  void Add_invert_z(int hex_b,  Cell_variable<DTyp>& a, int hex_a);
  void Add_invert_z(int hex_b_start, Cell_variable<DTyp>& a, int hex_a_start, int numBlocks);
  
  
  bool totalCalcNotPossible() const { return false; }
  int  getTotalNumberData()   const { return totalNumberData; }
  
  bool containsNaN();  
  
private:
  inline DTyp* give_startTotal() const  { return dataTotal; }  
  
  Blockgrid*        blockgrid;
  Unstructured_grid*       ug;

  // own data
  int totalNumberData;
  DTyp* dataTotal;
  
  int numberHex;  
  DTyp** data_cell;        // num_hexahedra

  // for parallel
  int my_rank;
  MPI_Comm comm;
};


//////////////////////////////////////////
// 2.2. product_cell, L_infty_cell
//////////////////////////////////////////////////////////////


template <class A, class B>
_TypeOf2_(A,B) product_cell(Expr<A>& a, Expr<B>& b, Marker& marker);

template <class A, class B>
_TypeOf2_(A,B) product_cell(Expr<A>& a, Expr<B>& b);


template <class A>
double L_infty_cell(const Expr<A>& a, Marker& marker);

template <class A>
double L_infty_cell(const Expr<A>& a);


template <class A>
double Maximum_cell(const Expr<A>& a);

template <class A>
double Minimum_cell(const Expr<A>& a);
/* @} */

//////////////////////////////////////////
// 3. Implementierungen
//////////////////////////////////////////////////////////////

template <class DTyp>
template <class LocalMat>
void Cell_variable<DTyp>::IntegrateOnCell(LocalMat& localStiffnessMatrix) {
     int Nx, Ny, Nz, N_total;
     
     assert(blockgrid == localStiffnessMatrix.Give_blockgrid());
     
     double** loc_m_hexahedra = localStiffnessMatrix.Give_loc_m_hexahedra();
     // hexahedra
     for(int id=0;id<numberHex;++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);
            N_total = Nx * Ny * Nz;

	    #pragma omp parallel for num_threads(UGBlocks::numThreadsToTake) if(UGBlocks::useOpenMP)
            for(int i=0;i<N_total;++i)
	        data_cell[id][i] = 0.0;
		
            #pragma omp parallel for num_threads(UGBlocks::numThreadsToTake) if(UGBlocks::useOpenMP)
            for(int i=0;i<N_total;++i) {
	        for(int s=0;s<64;++s) {
                    data_cell[id][i] = data_cell[id][i] + loc_m_hexahedra[id][i*64+s];
		}
	    }
	 }
     }
}


template <class DTyp>
Cell_variable<DTyp>::Cell_variable(Blockgrid& blockgrid_) {
  int Nx, Ny, Nz, N_total;

  blockgrid = &blockgrid_;
  ug = blockgrid->Give_unstructured_grid();

  // for parallel
  my_rank = ug->Give_my_rank();
  comm    = ug->Give_MPI_comm();

  numberHex = ug->Give_number_hexahedra();
  
  
  // hexahedra
  data_cell = new DTyp*[numberHex];

  totalNumberData = 0;
  for(int id=0;id<numberHex;++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);

         N_total = Nx * Ny * Nz;
         totalNumberData = totalNumberData + N_total;
      }
  }
  dataTotal = new DTyp[totalNumberData];
    
  int indexNow = 0;  
  for(int id=0;id<numberHex;++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);

         data_cell[id] = &(dataTotal[indexNow]);

         N_total = Nx * Ny * Nz;
         indexNow = indexNow + N_total;
      //data_cell[id] = new DTyp[N_total];

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

}

template <class DTyp>
Cell_variable<DTyp>::~Cell_variable() {
  if(dataTotal != NULL) {
     delete [] dataTotal;
     dataTotal = NULL;
  }
  if(data_cell != NULL) {
     /*
      for(int i=0; i<numberHex; ++i)
          if(data_cell[i] != NULL) 
              delete [] data_cell[i];
	  */
      delete [] data_cell;
  }
  data_cell = NULL;
}


template <class DTyp>
bool Cell_variable<DTyp>::containsNaN() {
    for(int i=0;i<totalNumberData;++i) {
        if(isNaN(dataTotal[i])) return true;
    }
    return false;
}


template <class DTyp>
void Cell_variable<DTyp>::operator=(DTyp x) {
#pragma omp parallel for num_threads(UGBlocks::numThreadsToTake) if(UGBlocks::useOpenMP)
    for(int i=0;i<totalNumberData;++i) {
        dataTotal[i] = x;
    }
  
  /*
  int Nx, Ny, Nz, N_total;

  // hexahedra
  for(int id=0;id<numberHex;++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);
      N_total = Nx * Ny * Nz;

#pragma omp parallel for num_threads(UGBlocks::numThreadsToTake) if(UGBlocks::useOpenMP)
      for(int i=0;i<N_total;++i)
        data_cell[id][i] = x;
    }
  }
  */
}

template <class DTyp>
void Cell_variable<DTyp>::operator=(const DTyp_Restriction<DTyp>& a) {
  int Nx, Ny, Nz, N_total;
  DTyp   x      = a.Give_x();

  // hexahedra
  for(int id=0;id<numberHex;++id) {

    if(ug->Give_hexahedron(id)->my_object(my_rank) &&
       a.template Give_marker<hexahedronEl>(id)==yes_mark) {

      Nx = blockgrid->Give_Nx_hexahedron(id);
      Ny = blockgrid->Give_Ny_hexahedron(id);
      Nz = blockgrid->Give_Nz_hexahedron(id);
      N_total = Nx * Ny * Nz;

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

template <class DTyp>
void Cell_variable<DTyp>::operator=(const Cell_variable<DTyp>& varRight) {
  int Nx, Ny, Nz;

  if(blockgrid->getId() != varRight.blockgrid->getId()) {
     assert(ug->isComposeGrid());
     ComposeUg* compUg = static_cast<ComposeUg*>(ug);

     int thisIdGrid = varRight.blockgrid->getId();
     int startHex = compUg->getStartHex(thisIdGrid);
     int endHex   = compUg->getEndHex(thisIdGrid);
  
     for(int id=startHex;id<endHex;++id) {
         int idRight = id - startHex;
//         varRight.template 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);

            int i,j,k;

// pragma laut gcc4.2 hier nicht erlaubt
#pragma omp parallel for private(i,j)  num_threads(UGBlocks::numThreadsToTake) if(UGBlocks::useOpenMP)
            for(k=0;k<Nz;++k) {
                for(j=0;j<Ny;++j)
                    for(i=0;i<Nx;++i){
                        data_cell[id][Ind_loc_matrix_hexahedra(i,j,k)] =
                           varRight.Give_cell_hexahedra(idRight,Ind_loc_matrix_hexahedra(i,j,k),i,j,k,Nx,Ny);
		    }
	    }
	 }
     }
     return;
  }
  
  // hexahedra
  for(int id=0;id<numberHex;++id) {
//    varRight.template 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);

      int i,j,k;

// pragma laut gcc4.2 hier nicht erlaubt
#pragma omp parallel for private(i,j) num_threads(UGBlocks::numThreadsToTake) if(UGBlocks::useOpenMP)
      for(k=0;k<Nz;++k) {
        for(j=0;j<Ny;++j)
          for(i=0;i<Nx;++i){
            data_cell[id][Ind_loc_matrix_hexahedra(i,j,k)] =
              varRight.Give_cell_hexahedra(id,Ind_loc_matrix_hexahedra(i,j,k),i,j,k,Nx,Ny);
	  }
      }
    }
  }
}


template <class DTyp>
template <class A>
void Cell_variable<DTyp>::operator=(const Expr<A>& a) {
   const A& ao ( a );

   if(ao.totalCalcNotPossible()==false) {
#pragma omp parallel for num_threads(UGBlocks::numThreadsToTake) if(UGBlocks::useOpenMP)     
      for(int i=0;i<totalNumberData;++i) {
          dataTotal[i] = ao.Give_fromTotal(i);
      }
   }
   else {  
      // hexahedra
      int Nx, Ny, Nz;  
      for(int id=0;id<numberHex;++id) {
          //ao.template 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);

#pragma omp parallel for num_threads(UGBlocks::numThreadsToTake) if(UGBlocks::useOpenMP)
             for(int k=0;k<Nz;++k) {
                 for(int j=0;j<Ny;++j)
                     for(int i=0;i<Nx;++i) {
                         data_cell[id][Ind_loc_matrix_hexahedra(i,j,k)] =
                             ao.Give_cell_hexahedra(id,Ind_loc_matrix_hexahedra(i,j,k),i,j,k,Nx,Ny);
		     }
	     }
	  }
      }
   }
}

template <class DTyp>
template <class A>
void Cell_variable<DTyp>::operator=(const Expr_Restriction<A>& a)  {
  int Nx, Ny, Nz;

  //  const A& ao ( a );

  // hexahedra
  for(int id=0;id<numberHex;++id) {
    //a.template Update<hexahedronEl>(id);

    if(ug->Give_hexahedron(id)->my_object(my_rank) &&
       a.template Give_marker<hexahedronEl>(id)==yes_mark) {
      
      Nx = blockgrid->Give_Nx_hexahedron(id);
      Ny = blockgrid->Give_Ny_hexahedron(id);
      Nz = blockgrid->Give_Nz_hexahedron(id);

      int i,j,k;
// pragma laut gcc4.2 hier nicht erlaubt
#pragma omp parallel for private(i,j) num_threads(UGBlocks::numThreadsToTake) if(UGBlocks::useOpenMP)
      for(k=0;k<Nz;++k) {
        for(j=0;j<Ny;++j)
          for(i=0;i<Nx;++i){
            data_cell[id][Ind_loc_matrix_hexahedra(i,j,k)] =
              a.Give_cell_hexahedra(id,Ind_loc_matrix_hexahedra(i,j,k),i,j,k,Nx,Ny);
	  }
      }
    }
  }
}

//////////////////////////////////////////
// 3.3. product_cell, L_infty_cell
//////////////////////////////////////////////////////////////

template<class DTyp>
struct ProductCalculator {
   template <class A, class B>
   static DTyp product_cell(Expr<A>& ao, Expr<B>& bo, Marker& marker);
   
   template <class A, class B>   
   static DTyp product_cell(Expr<A>& ao, Expr<B>& bo);   
};

template <class A, class B>
  _TypeOf2_(A,B) product_cell(Expr<A>& a, Expr<B>& b) {
  const A& ao ( a );
  typedef _TypeOf2_(A, B) DTyp;
  return ProductCalculator<DTyp>::product_cell(a,b);
}

template <class A, class B>
  _TypeOf2_(A,B) product_cell(Expr<A>& a, Expr<B>& b, Marker& marker) {
  const A& ao ( a );
  typedef _TypeOf2_(A, B) DTyp;
  return ProductCalculator<DTyp>::product_cell(a,b,marker);
}  
  

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


template<>
template <class A, class B>
double ProductCalculator<double>::product_cell(Expr<A>& ao, Expr<B>& bo) {
  int Nx, Ny, Nz;
  const A& a ( ao );
  const B& b ( bo );

  typedef double DTyp;
  
  double sum;
  double sum_total;

  sum = 0.0;

  if(developer_version) 
     if(a.Give_blockgrid() != b.Give_blockgrid()) std::cout << " error2 in product_cell!" << std::endl;

  Unstructured_grid* ug = a.Give_blockgrid()->Give_unstructured_grid();
  int my_rank   = ug->Give_my_rank();
  MPI_Comm comm = ug->Give_MPI_comm();

  
   if(a.totalCalcNotPossible() || b.totalCalcNotPossible()) {  
      for(int id=0;id<ug->Give_number_hexahedra();++id) {
          if(ug->Give_hexahedron(id)->my_object(my_rank)) {
             Nx = a.Give_blockgrid()->Give_Nx_hexahedron(id);
             Ny = a.Give_blockgrid()->Give_Ny_hexahedron(id);
             Nz = a.Give_blockgrid()->Give_Nz_hexahedron(id);

             double localSum = 0.0;
             #pragma omp parallel for reduction(+ : localSum) num_threads(UGBlocks::numThreadsToTake) if(UGBlocks::useOpenMP)     
             for(int k=0;k<Nz;++k) {
                 for(int j=0;j<Ny;++j)
                     for(int i=0;i<Nx;++i)
	                 localSum  +=   a.Give_cell_hexahedra(id,Ind_loc_matrix_hexahedra(i,j,k),i,j,k,Nx,Ny)
	                              * b.Give_cell_hexahedra(id,Ind_loc_matrix_hexahedra(i,j,k),i,j,k,Nx,Ny);
	     }
     
             sum = sum + localSum;
	  }
      }
  }
  else {          // WORK NOW muss man nich richtig machen mit getTotalNumberData
      #pragma omp parallel for reduction(+ : sum) num_threads(UGBlocks::numThreadsToTake) if(UGBlocks::useOpenMP)     
      for(int i=0;i<a.getTotalNumberData();++i) {
	  sum  = sum + a.Give_fromTotal(i) * b.Give_fromTotal(i); 
      }
  }

  sum_total = Make_MPI_All_Sum<DTyp>::Do(sum,comm);

  return sum_total;
}

template<>
template <class A, class B>
double ProductCalculator<double>::product_cell(Expr<A>& ao, Expr<B>& bo, Marker& marker) {
  int Nx, Ny, Nz;
  const A& a ( ao );
  const B& b ( bo );

  typedef double DTyp;
  
  double sum;
  double sum_total;

  sum = 0.0;

  if(developer_version) 
     if(a.Give_blockgrid() != b.Give_blockgrid()) std::cout << " error2 in product_cell!" << std::endl;

  Unstructured_grid* ug = a.Give_blockgrid()->Give_unstructured_grid();
  int my_rank   = ug->Give_my_rank();
  MPI_Comm comm = ug->Give_MPI_comm();

  
      for(int id=0;id<ug->Give_number_hexahedra();++id) {
          if(ug->Give_hexahedron(id)->my_object(my_rank)   &&
             marker.template Give_marker<hexahedronEl>(id)==yes_mark) {
             Nx = a.Give_blockgrid()->Give_Nx_hexahedron(id);
             Ny = a.Give_blockgrid()->Give_Ny_hexahedron(id);
             Nz = a.Give_blockgrid()->Give_Nz_hexahedron(id);

             double localSum = 0.0;
             #pragma omp parallel for reduction(+ : localSum) num_threads(UGBlocks::numThreadsToTake) if(UGBlocks::useOpenMP)     
             for(int k=0;k<Nz;++k) {
                 for(int j=0;j<Ny;++j)
                     for(int i=0;i<Nx;++i)
	                 localSum  +=   a.Give_cell_hexahedra(id,Ind_loc_matrix_hexahedra(i,j,k),i,j,k,Nx,Ny)
	                              * b.Give_cell_hexahedra(id,Ind_loc_matrix_hexahedra(i,j,k),i,j,k,Nx,Ny);
	     }
     
             sum = sum + localSum;
	  }
      }

  sum_total = Make_MPI_All_Sum<DTyp>::Do(sum,comm);

  return sum_total;
}

template <class DTyp>
template <class A, class B>
DTyp ProductCalculator<DTyp>::product_cell(Expr<A>& ao, Expr<B>& bo, Marker& marker) {
  int Nx, Ny, Nz;
  const A& a ( ao );
  const B& b ( bo );

  DTyp sum;
  DTyp sum_total;

  sum = (DTyp)0;

  if(developer_version) 
     if(a.Give_blockgrid() != b.Give_blockgrid()) cout << " error2 in product_cell!" << endl;

  Unstructured_grid* ug = a.Give_blockgrid()->Give_unstructured_grid();
  int my_rank   = ug->Give_my_rank();
  MPI_Comm comm = ug->Give_MPI_comm();

  // hexahedra
  for(int id=0;id<ug->Give_number_hexahedra();++id) {
     if(ug->Give_hexahedron(id)->my_object(my_rank)   &&
        marker.template Give_marker<hexahedronEl>(id)==yes_mark) {
      Nx = a.Give_blockgrid()->Give_Nx_hexahedron(id);
      Ny = a.Give_blockgrid()->Give_Ny_hexahedron(id);
      Nz = a.Give_blockgrid()->Give_Nz_hexahedron(id);

      for(int k=0;k<Nz;++k) {
        for(int j=0;j<Ny;++j)
          for(int i=0;i<Nx;++i)
	    sum+= a.Give_cell_hexahedra(id,
					Ind_loc_matrix_hexahedra(i,j,k),i,j,k,Nx,Ny)
	      * MyConj<DTyp>::Give(b.Give_cell_hexahedra(id,
							 Ind_loc_matrix_hexahedra(i,j,k),i,j,k,Nx,Ny));
      }
     }
  }

  sum_total = Make_MPI_All_Sum<DTyp>::Do(sum,comm);

  return sum_total;
}




template <class A>
double L_infty_cell(const Expr<A>& ao, Marker& marker) {
  int Nx, Ny, Nz;

  double infty;
  double infty_total;

  const A& a ( ao );

  infty = 0.0;

  Unstructured_grid* ug = a.Give_blockgrid()->Give_unstructured_grid();
  int my_rank   = ug->Give_my_rank();
  MPI_Comm comm = ug->Give_MPI_comm();

     for(int id=0; id<ug->Give_number_hexahedra(); ++id) {
         if(ug->Give_hexahedron(id)->my_object(my_rank) &&
                    marker.template Give_marker<hexahedronEl>(id)==yes_mark) {
            Nx = a.Give_blockgrid()->Give_Nx_hexahedron(id);
            Ny = a.Give_blockgrid()->Give_Ny_hexahedron(id);
            Nz = a.Give_blockgrid()->Give_Nz_hexahedron(id);

	    
	    double localInfy = 0.0;
#pragma omp parallel for reduction(max : localInfy) num_threads(UGBlocks::numThreadsToTake) if(UGBlocks::useOpenMP)    
            for(int k=0;k<Nz;++k) {
                for(int j=0;j<Ny;++j)
                    for(int i=0;i<Nx;++i)
	                if(localInfy < ABS(a.Give_cell_hexahedra(id,
						 Ind_loc_matrix_hexahedra(i,j,k),i,j,k,Nx,Ny)))
	                   localInfy = ABS(a.Give_cell_hexahedra(id,
						Ind_loc_matrix_hexahedra(i,j,k),i,j,k,Nx,Ny));
	    }
	    
	    if(infty < localInfy) infty = localInfy;
	 }
     }
  
  MPI_Allreduce(&infty,&infty_total,1,MPI_DOUBLE,MPI_MAX,comm);

  return infty_total;
}


template <class A>
double L_infty_cell(const Expr<A>& ao) {
  int Nx, Ny, Nz;

  double infty;
  double infty_total;

  const A& a ( ao );

  infty = 0.0;

  Unstructured_grid* ug = a.Give_blockgrid()->Give_unstructured_grid();
  int my_rank   = ug->Give_my_rank();
  MPI_Comm comm = ug->Give_MPI_comm();

  if(a.totalCalcNotPossible()) {  
     for(int id=0; id<ug->Give_number_hexahedra(); ++id) {
         if(ug->Give_hexahedron(id)->my_object(my_rank)) {
            Nx = a.Give_blockgrid()->Give_Nx_hexahedron(id);
            Ny = a.Give_blockgrid()->Give_Ny_hexahedron(id);
            Nz = a.Give_blockgrid()->Give_Nz_hexahedron(id);

	    
	    double localInfy = 0.0;
//#pragma omp parallel for reduction(max : localInfy) num_threads(UGBlocks::numThreadsToTake) if(UGBlocks::useOpenMP)    
            for(int k=0;k<Nz;++k) {
                for(int j=0;j<Ny;++j)
                    for(int i=0;i<Nx;++i)
	                if(localInfy < ABS(a.Give_cell_hexahedra(id,
						 Ind_loc_matrix_hexahedra(i,j,k),i,j,k,Nx,Ny)))
	                   localInfy = ABS(a.Give_cell_hexahedra(id,
						Ind_loc_matrix_hexahedra(i,j,k),i,j,k,Nx,Ny));
	    }
	    
	    if(infty < localInfy) infty = localInfy;
	 }
     }
  }
  else {
//#pragma omp parallel for reduction(max : infty) num_threads(UGBlocks::numThreadsToTake) if(UGBlocks::useOpenMP)     
      for(int i=0;i<a.getTotalNumberData();++i) {
	  if(infty < ABS(a.Give_fromTotal(i))) infty = ABS(a.Give_fromTotal(i)); 
      }
  }
  
  MPI_Allreduce(&infty,&infty_total,1,MPI_DOUBLE,MPI_MAX,comm);

  return infty_total;
}


template <class A>
double Maximum_cell(const Expr<A>& ao)  {
  int Nx, Ny, Nz;

  double maximum;
  double maximum_total;

  maximum = -std::numeric_limits<double>::max(); // früher: -1.0e50;
  
  const A& a ( ao );


  Unstructured_grid* ug = a.Give_blockgrid()->Give_unstructured_grid();
  int my_rank   = ug->Give_my_rank();
  MPI_Comm comm = ug->Give_MPI_comm();

  
  if(a.totalCalcNotPossible()) {  
     for(int id=0; id<ug->Give_number_hexahedra(); ++id) {
         if(ug->Give_hexahedron(id)->my_object(my_rank)) {
            Nx = a.Give_blockgrid()->Give_Nx_hexahedron(id);
            Ny = a.Give_blockgrid()->Give_Ny_hexahedron(id);
            Nz = a.Give_blockgrid()->Give_Nz_hexahedron(id);

	    double localMax = -std::numeric_limits<double>::max();
//#pragma omp parallel for reduction(max:localMax) num_threads(UGBlocks::numThreadsToTake) if(UGBlocks::useOpenMP)    	 
            for(int k=0;k<Nz;++k) {
                for(int j=0;j<Ny;++j)
                    for(int i=0;i<Nx;++i) {
	                if(maximum < a.Give_cell_hexahedra(id,
			   			 Ind_loc_matrix_hexahedra(i,j,k),i,j,k,Nx,Ny))
	                   maximum = a.Give_cell_hexahedra(id,
						Ind_loc_matrix_hexahedra(i,j,k),i,j,k,Nx,Ny);
		    }
	    } 
	    if(maximum < localMax) maximum = localMax;
	 }
     }
  }
  else {
//    #pragma omp parallel for reduction(max:maximum) num_threads(UGBlocks::numThreadsToTake) if(UGBlocks::useOpenMP)     
      for(int i=0;i<a.getTotalNumberData();++i) {
	  if(maximum < a.Give_fromTotal(i)) maximum = a.Give_fromTotal(i); 
      }
  }
  
  MPI_Allreduce(&maximum,&maximum_total,1,MPI_DOUBLE,MPI_MAX,comm);

  return maximum_total;
}

template <class A>
double Minimum_cell(const Expr<A>& ao)  {
  int Nx, Ny, Nz;

  double minimum;
  double minimum_total;

  minimum = std::numeric_limits<double>::max();
  
  const A& a ( ao );


  Unstructured_grid* ug = a.Give_blockgrid()->Give_unstructured_grid();
  int my_rank   = ug->Give_my_rank();
  MPI_Comm comm = ug->Give_MPI_comm();

  
  if(a.totalCalcNotPossible()) {  
     for(int id=0; id<ug->Give_number_hexahedra(); ++id) {
         if(ug->Give_hexahedron(id)->my_object(my_rank)) {
            Nx = a.Give_blockgrid()->Give_Nx_hexahedron(id);
            Ny = a.Give_blockgrid()->Give_Ny_hexahedron(id);
            Nz = a.Give_blockgrid()->Give_Nz_hexahedron(id);

	    double localMin = std::numeric_limits<double>::max();
//#pragma omp parallel for reduction(min:localMin) num_threads(UGBlocks::numThreadsToTake) if(UGBlocks::useOpenMP)    	 
            for(int k=0;k<Nz;++k) {
                for(int j=0;j<Ny;++j)
                    for(int i=0;i<Nx;++i) {
	                if(minimum > a.Give_cell_hexahedra(id,
			   			 Ind_loc_matrix_hexahedra(i,j,k),i,j,k,Nx,Ny))
	                   minimum = a.Give_cell_hexahedra(id,
						Ind_loc_matrix_hexahedra(i,j,k),i,j,k,Nx,Ny);
		    }
	    } 
	    if(minimum > localMin) minimum = localMin;
	 }
     }
  }
  else {
//    #pragma omp parallel for reduction(min:minimum) num_threads(UGBlocks::numThreadsToTake) if(UGBlocks::useOpenMP)     
      for(int i=0;i<a.getTotalNumberData();++i) {
	  if(minimum > a.Give_fromTotal(i)) minimum = a.Give_fromTotal(i); 
      }
  }
  
  MPI_Allreduce(&minimum,&minimum_total,1,MPI_DOUBLE,MPI_MIN,comm);

  return minimum_total;
}

/*
template <class A>
double Minimum_cell(const Expr<A>& ao) {
  int Nx, Ny, Nz;

  double minimum;
  double minimum_total;

  minimum = std::numeric_limits<double>::max();
  
  const A& a ( ao );


  Unstructured_grid* ug = a.Give_blockgrid()->Give_unstructured_grid();
  int my_rank   = ug->Give_my_rank();
  MPI_Comm comm = ug->Give_MPI_comm();

  
  
  for(int id=0; id<ug->Give_number_hexahedra(); ++id) {
    if(ug->Give_hexahedron(id)->my_object(my_rank)) {
      Nx = a.Give_blockgrid()->Give_Nx_hexahedron(id);
      Ny = a.Give_blockgrid()->Give_Ny_hexahedron(id);
      Nz = a.Give_blockgrid()->Give_Nz_hexahedron(id);

      /// todo hier sollte noch ein reduction rein 
      //#pragma omp parallel for num_threads(UGBlocks::numThreadsToTake) private(i,j)  if(UGBlocks::useOpenMP)
      for(int k=0;k<Nz;++k) {
        for(int j=0;j<Ny;++j)
          for(int i=0;i<Nx;++i) {
	    if(minimum > ABS(a.Give_cell_hexahedra(id,
						 Ind_loc_matrix_hexahedra(i,j,k),i,j,k,Nx,Ny)))
	      minimum = ABS(a.Give_cell_hexahedra(id,
						Ind_loc_matrix_hexahedra(i,j,k),i,j,k,Nx,Ny));
      }
    }
   }
  }

  MPI_Allreduce(&minimum,&minimum_total,1,MPI_DOUBLE,MPI_MAX,comm);

  return minimum_total;
}
*/



#endif // CE_VA_H