/***************************************************************************
                          tnlGrid.h  -  description
                             -------------------
    begin                : Dec 12, 2010
    copyright            : (C) 2010 by Tomas Oberhuber
    email                : tomas.oberhuber@fjfi.cvut.cz
 ***************************************************************************/

/***************************************************************************
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 ***************************************************************************/

#ifndef TNLGRID_H_
#define TNLGRID_H_

#include <iomanip>
#include <fstream>
#include <core/tnlAssert.h>
#include <core/tnlArray.h>
#include <core/tnlVector.h>

using namespace std;

template< int Dimensions, typename Real = double, tnlDevice Device = tnlHost, typename Index = int >
class tnlGrid : public tnlArray< Dimensions, Real, Device, Index >
{
   //! We do not allow constructor without parameters.
   tnlGrid();

   //! We do not allow copy constructor without object name.
   tnlGrid( const tnlGrid< Dimensions, Real, Device, Index >& a );

   public:

   tnlGrid( const tnlString& name );

   tnlGrid( const tnlString& name,
            const tnlGrid< Dimensions, Real, tnlHost, Index >& grid );

   tnlGrid( const tnlString& name,
            const tnlGrid< Dimensions, Real, tnlCuda, Index >& grid );

   const tnlTuple< Dimensions, Index >& getDimensions() const;

   //! Sets the dimensions
   /***
    * This method also must recompute space steps. It is save to call setDimensions and
    * setDomain in any order. Both recompute the space steps.
    */
   bool setDimensions( const tnlTuple< Dimensions, Index >& dimensions );

   //! Sets the computation domain in form of "rectangle".
   /***
    * This method also must recompute space steps. It is save to call setDimensions and
    * setDomain in any order. Both recompute the space steps.
    */

   bool setDomain( const tnlTuple< Dimensions, Real >& lowerCorner,
                   const tnlTuple< Dimensions, Real >& upperCorner );

   //! Set dimensions and domain of the grid using another grid as a template
   bool setLike( const tnlGrid< Dimensions, Real, tnlHost, Index >& v );

   //! Set dimensions and domain of the grid using another grid as a template
   bool setLike( const tnlGrid< Dimensions, Real, tnlCuda, Index >& v );

   const tnlTuple< Dimensions, Real >& getDomainLowerCorner() const;

   const tnlTuple< Dimensions, Real >& getDomainUpperCorner() const;

   const tnlTuple< Dimensions, Real >& getSpaceSteps() const;

   tnlString getType() const;

   bool operator == ( const tnlGrid< Dimensions, Real, Device, Index >& array ) const;

   bool operator != ( const tnlGrid< Dimensions, Real, Device, Index >& array ) const;

   template< typename Real2, tnlDevice Device2, typename Index2 >
   tnlGrid< Dimensions, Real, Device, Index >& operator = ( const tnlGrid< Dimensions, Real2, Device2, Index2 >& array );

   //! This method interpolates value at given point.
   Real getValue( const tnlTuple< Dimensions, Real >& point ) const;

   //! Interpolation for 1D grid.
   Real getValue( const Real& x ) const;

   //! Interpolation for 2D grid.
   Real getValue( const Real& x,
                  const Real& y ) const;

   //! Interpolation for 3D grid.
   Real getValue( const Real& x,
                  const Real& y,
                  const Real& z ) const;

   //! Forward difference w.r.t x
   Real Partial_x_f( const Index i1 ) const;

   //! Forward difference w.r.t x in two dimensions
   Real Partial_x_f( const Index i1,
                     const Index i2 ) const;

   //! Forward difference w.r.t x in three dimensions
   Real Partial_x_f( const Index i1,
                     const Index i2,
                     const Index i3 ) const;

   //! Backward difference w.r.t x
   Real Partial_x_b( const Index i1 ) const;

   //! Backward difference w.r.t x in two dimensions
   Real Partial_x_b( const Index i1,
                     const Index i2 ) const;

   //! Backward difference w.r.t x in three dimensions
   Real Partial_x_b( const Index i1,
                     const Index i2,
                     const Index i3 ) const;

   //! Central difference w.r.t. x
   Real Partial_x( const Index i1 ) const;

   //! Central difference w.r.t. x in two dimensions
   Real Partial_x( const Index i1,
                   const Index i2 ) const;

   //! Central difference w.r.t. x
   Real Partial_x( const Index i1,
                   const Index i2,
                   const Index i3 ) const;

   //! Second order difference w.r.t. x
   Real Partial_xx( const Index i1 ) const;

   //! Second order difference w.r.t. x in two dimensions
   Real Partial_xx( const Index i1,
                    const Index i2 ) const;

   //! Second order difference w.r.t. x
   Real Partial_xx( const Index i1,
                    const Index i2,
                    const Index i3 ) const;

   //! Forward difference w.r.t y
   Real Partial_y_f( const Index i1,
                     const Index i2 ) const;

   //! Forward difference w.r.t y in three dimensions
   Real Partial_y_f( const Index i1,
                     const Index i2,
                     const Index i3 ) const;

   //! Backward difference w.r.t y
   Real Partial_y_b( const Index i1,
                     const Index i2 ) const;

   //! Backward difference w.r.t y in three dimensions
   Real Partial_y_b( const Index i1,
                     const Index i2,
                     const Index i3 ) const;

   //! Central difference w.r.t y
   Real Partial_y( const Index i1,
                   const Index i2 ) const;

   //! Central difference w.r.t y
   Real Partial_y( const Index i1,
                   const Index i2,
                   const Index i3 ) const;

   //! Second order difference w.r.t. y
   Real Partial_yy( const Index i1,
                    const Index i2 ) const;

   //! Second order difference w.r.t. y in three dimensions
   Real Partial_yy( const Index i1,
                    const Index i2,
                    const Index i3 ) const;

   //! Forward difference w.r.t z
   Real Partial_z_f( const Index i1,
                     const Index i2,
                     const Index i3 ) const;

   //! Backward difference w.r.t z
   Real Partial_z_b( const Index i1,
                     const Index i2,
                     const Index i3 ) const;

   //! Central difference w.r.t z
   Real Partial_z( const Index i1,
                   const Index i2,
                   const Index i3 ) const;

   //! Second order difference w.r.t. z
   Real Partial_zz( const Index i1,
                    const Index i2,
                    const Index i3 ) const;

   //! Set space dependent Dirichlet boundary conditions
   void setDirichletBC( const tnlGrid< Dimensions, Real, Device, Index >&bc,
                        const tnlTuple< Dimensions, bool >& lowerBC,
                        const tnlTuple< Dimensions, bool >& upperBC );

   //! Set constant Dirichlet boundary conditions
   void setDirichletBC( const Real& bc,
                        const tnlTuple< Dimensions, bool >& lowerBC,
                        const tnlTuple< Dimensions, bool >& upperBC );

   //! Set space dependent Neumann boundary conditions
   void setNeumannBC( const tnlGrid< Dimensions, Real, Device, Index >&bc,
                      const tnlTuple< Dimensions, bool >& lowerBC,
                      const tnlTuple< Dimensions, bool >& upperBC );

   //! Set constant Neumann boundary conditions
   void setNeumannBC( const Real& bc,
                      const tnlTuple< Dimensions, bool >& lowerBC,
                      const tnlTuple< Dimensions, bool >& upperBC );

   Real getMax() const;

   Real getMin() const;

   Real getAbsMax() const;

   Real getAbsMin() const;

   Real getLpNorm() const;

   Real getSum() const;

   Real getDifferenceMax( const tnlVector< Real, tnlHost, Index >& v ) const;

   Real getDifferenceMin( const tnlVector< Real, tnlHost, Index >& v ) const;

   Real getDifferenceAbsMax( const tnlVector< Real, tnlHost, Index >& v ) const;

   Real getDifferenceAbsMin( const tnlVector< Real, tnlHost, Index >& v ) const;

   Real getDifferenceLpNorm( const tnlVector< Real, tnlHost, Index >& v, const Real& p ) const;

   Real getDifferenceSum( const tnlVector< Real, tnlHost, Index >& v ) const;

   //! Method for saving the object to a file as a binary data
   bool save( tnlFile& file ) const;

   //! Method for restoring the object from a file
   bool load( tnlFile& file );

   bool save( const tnlString& fileName ) const;

   bool load( const tnlString& fileName );

   //! This method writes the grid in some format suitable for some other preprocessing.
   /*! Possible formats are:
    *  1) Gnuplot format (gnuplot)
    *  2) VTI format (vti)
    *  3) Povray format (povray) - only for 3D.
    */
   bool draw( const tnlString& fileName,
              const tnlString& format,
              const tnlTuple< Dimensions, Index > steps = ( tnlTuple< Dimensions, Index > ) 1 ) const;

   protected:
   tnlTuple< Dimensions, Real > domainLowerCorner, domainUpperCorner, spaceSteps;
};

#ifdef HAVE_CUDA
template< int Dimensions, typename Real, typename Index >
__global__ void setConstantDirichletBC( const Index xSize,
                                        const Index ySize,
                                        const Real bc,
                                        Real* u );

template< int Dimensions, typename Real, typename Index >
__global__ void setDirichletBC( const Index xSize,
                                const Index ySize,
                                const Real* bc,
                                Real* u );

template< int Dimensions, typename Real, typename Index >
__global__ void setConstantNeumannBC( const Index xSize,
                                      const Index ySize,
                                      const Real hx,
                                      const Real hy,
                                      const Real bc,
                                      Real* u );

template< int Dimensions, typename Real, typename Index >
__global__ void setNeumannBC( const Index xSize,
                              const Index ySize,
                              const Real hx,
                              const Real hy,
                              const Real* bc,
                              Real* u );
#endif

template< int Dimensions, typename Real, tnlDevice Device, typename Index >
tnlGrid< Dimensions, Real, Device, Index > :: tnlGrid( const tnlString& name )
: tnlArray< Dimensions, Real, Device, Index >( name )
  {
  }

template< int Dimensions, typename Real, tnlDevice Device, typename Index >
tnlGrid< Dimensions, Real, Device, Index > :: tnlGrid( const tnlString& name,
                                                       const tnlGrid< Dimensions, Real, tnlHost, Index >& grid )
: tnlArray< Dimensions, Real, Device, Index >( name, grid )
{
   this -> setDomain( grid. getDomainLowerCorner(),
                      grid. getDomainUpperCorner() );
}

template< int Dimensions, typename Real, tnlDevice Device, typename Index >
tnlGrid< Dimensions, Real, Device, Index > :: tnlGrid( const tnlString& name,
                                                       const tnlGrid< Dimensions, Real, tnlCuda, Index >& grid )
: tnlArray< Dimensions, Real, Device, Index >( name, grid )
{
   this -> setDomain( grid. getDomainLowerCorner(),
                      grid. getDomainUpperCorner() );
}

template< int Dimensions, typename Real, tnlDevice Device, typename Index >
const tnlTuple< Dimensions, Index >& tnlGrid< Dimensions, Real, Device, Index > :: getDimensions() const
{
   return tnlArray< Dimensions, Real, Device, Index > :: getDimensions();
}

template< int Dimensions, typename Real, tnlDevice Device, typename Index >
bool tnlGrid< Dimensions, Real, Device, Index > :: setDimensions( const tnlTuple< Dimensions, Index >& dimensions )
{
   if( ! tnlArray< Dimensions, Real, Device, Index > :: setDimensions( dimensions ) )
      return false;
   for( int i = 0; i < Dimensions; i ++ )
      spaceSteps[ i ] = ( domainUpperCorner[ i ] - domainLowerCorner[ i ] ) / ( Real ) ( this -> getDimensions()[ i ] - 1 );
   return true;
}

template< int Dimensions, typename Real, tnlDevice Device, typename Index >
bool tnlGrid< Dimensions, Real, Device, Index > :: setDomain( const tnlTuple< Dimensions, Real >& lowerCorner,
                                                              const tnlTuple< Dimensions, Real >& upperCorner )
{
   if( lowerCorner >= upperCorner )
   {
      cerr << "Wrong parameters for the grid domain of " << this -> getName() << ". The low corner must by smaller than the high corner." << endl
               << "lowerCorner = " << lowerCorner << endl << "upperCorner = " << upperCorner << endl;
      return false;
   }
   domainLowerCorner = lowerCorner;
   domainUpperCorner = upperCorner;
   for( int i = 0; i < Dimensions; i ++ )
      spaceSteps[ i ] = ( domainUpperCorner[ i ] - domainLowerCorner[ i ] ) / ( Real ) ( this -> getDimensions()[ i ] - 1 );
   return true;
}

template< int Dimensions, typename Real, tnlDevice Device, typename Index >
bool tnlGrid< Dimensions, Real, Device, Index > :: setLike( const tnlGrid< Dimensions, Real, tnlHost, Index >& v )
{
   return tnlArray< Dimensions, Real, Device, Index > :: setLike( v ) &&
          this -> setDomain( v. getDomainLowerCorner(), v. getDomainUpperCorner() );
}

template< int Dimensions, typename Real, tnlDevice Device, typename Index >
bool tnlGrid< Dimensions, Real, Device, Index > :: setLike( const tnlGrid< Dimensions, Real, tnlCuda, Index >& v )
{
   return tnlArray< Dimensions, Real, Device, Index > :: setLike( v ) &&
          this -> setDomain( v. getDomainLowerCorner(), v. getDomainUpperCorner() );
}

template< int Dimensions, typename Real, tnlDevice Device, typename Index >
const tnlTuple< Dimensions, Real >& tnlGrid< Dimensions, Real, Device, Index > :: getDomainLowerCorner() const
{
   return this -> domainLowerCorner;
}

template< int Dimensions, typename Real, tnlDevice Device, typename Index >
const tnlTuple< Dimensions, Real >& tnlGrid< Dimensions, Real, Device, Index > :: getDomainUpperCorner() const
{
   return this -> domainUpperCorner;
}

template< int Dimensions, typename Real, tnlDevice Device, typename Index >
const tnlTuple< Dimensions, Real >& tnlGrid< Dimensions, Real, Device, Index > :: getSpaceSteps() const
{
   return spaceSteps;
}

template< int Dimensions, typename Real, tnlDevice Device, typename Index >
tnlString tnlGrid< Dimensions, Real, Device, Index > :: getType() const
{
   return tnlString( "tnlGrid< ") +
          tnlString( Dimensions ) +
          tnlString( ", " ) +
          tnlString( getParameterType< Real >() ) +
          tnlString( ", " ) +
          getDeviceType( Device ) +
          tnlString( ", " ) +
          tnlString( getParameterType< Index >() ) +
          tnlString( " >" );
}

template< int Dimensions, typename Real, tnlDevice Device, typename Index >
bool tnlGrid< Dimensions, Real, Device, Index > :: operator == ( const tnlGrid< Dimensions, Real, Device, Index >& grid ) const
{
   tnlAssert( this -> getDomainLowerCorner() == grid. getDomainLowerCorner() &&
              this -> getDomainUpperCorner() == grid. getDomainUpperCorner(),
              cerr << "You are attempting to compare two grids with different domains." << endl
                   << "First grid name is " << this -> getName()
                   << " domain is ( " << this -> getDomainLowerCorner() << " )- ("
                                      << this -> getDomainUpperCorner() << ")" << endl
                   << "Second grid is " << grid. getName()
                   << " domain is ( " << grid. getDomainLowerCorner() << " )- ("
                                      << grid. getDomainUpperCorner() << ")" << endl; );
   return tnlArray< Dimensions, Real, Device, Index > :: operator == ( grid );
}

template< int Dimensions, typename Real, tnlDevice Device, typename Index >
bool tnlGrid< Dimensions, Real, Device, Index > :: operator != ( const tnlGrid< Dimensions, Real, Device, Index >& grid ) const
{
   return ! ( (* this ) == grid );
}

template< int Dimensions, typename Real, tnlDevice Device, typename Index >
  template< typename Real2, tnlDevice Device2, typename Index2 >
tnlGrid< Dimensions, Real, Device, Index >&
tnlGrid< Dimensions, Real, Device, Index >
:: operator = ( const tnlGrid< Dimensions, Real2, Device2, Index2 >& grid )
{
   tnlAssert( this -> getDimensions() == grid. getDimensions(),
              cerr << "You are attempting to assign two arrays with different dimensions." << endl
                   << "First array name is " << this -> getName()
                   << " dimensions are ( " << this -> getDimensions() << " )" << endl
                   << "Second array is " << grid. getName()
                   << " dimensions are ( " << grid. getDimensions() << " )" << endl; );
   tnlVector< Real, Device, Index > :: operator = ( grid );
   return ( *this );
}

template< int Dimensions, typename Real, tnlDevice Device, typename Index >
Real tnlGrid< Dimensions, Real, Device, Index > :: getValue( const tnlTuple< Dimensions, Real >& point ) const
{
   if( Dimensions == 1)
   {
      tnlAssert( 0, cerr << "Interpolation is not implemented yet.");
   }
   if( Dimensions == 2 )
   {
      Real x = ( point[ 0 ] - domainLowerCorner[ 0 ] ) / spaceSteps[ 0 ];
      Real y = ( point[ 1 ] - domainLowerCorner[ 1 ] ) / spaceSteps[ 1 ];
      Index ix = ( Index ) ( x );
      Index iy = ( Index ) ( y );
      Real dx = x - ( Real ) ix;
      Real dy = y - ( Real ) iy;
      if( iy >= this -> getDimensions()[ tnlY ] - 1 )
      {
         if( ix >= this -> getDimensions()[ tnlX ] - 1 )
            return  this -> getElement( this -> getDimensions()[ tnlX ] - 1,
                                        this -> getDimensions()[ tnlY ] - 1 );
         return ( Real( 1.0 ) - dx ) * this -> getElement( ix,
                                                           this -> getDimensions()[ tnlY ] - 1 ) +
                                                           dx * this -> getElement( ix + 1,
                                                           this -> getDimensions()[ tnlY ] - 1 );
      }
      if( ix >= this -> getDimensions()[ tnlX ] - 1 )
         return ( Real( 1.0 ) - dy ) * this -> getElement( this -> getDimensions()[ tnlX ] - 1,
                                                           iy ) +
                                                           dy * this -> getElement( this -> getDimensions()[ tnlX ] - 1,
                                                           iy + 1 );
      Real a1, a2;
      a1 = ( Real( 1.0 ) - dx ) * this -> getElement( ix, iy ) +
                             dx * this -> getElement( ix + 1, iy );

      a2 = ( Real( 1.0 ) - dx ) * this -> getElement( ix, iy + 1 ) +
                             dx * this -> getElement( ix + 1, iy + 1 );
      return ( Real( 1.0 ) - dy ) * a1 + dy * a2;
   }
   if( Dimensions == 3 )
   {
      tnlAssert( 0, cerr << "Interpolation is not implemented yet.");
   }
   if( Dimensions > 3 )
   {
      cerr << "Interpolation for 4D grids is not implemented yet." << endl;
      abort();
   }
}

template< int Dimensions, typename Real, tnlDevice Device, typename Index >
Real tnlGrid< Dimensions, Real, Device, Index > :: getValue( const Real& x ) const
{
   return this -> getValue( tnlTuple< 1, Real >( x ) );
}

template< int Dimensions, typename Real, tnlDevice Device, typename Index >
Real tnlGrid< Dimensions, Real, Device, Index > :: getValue( const Real& x,
                                                             const Real& y ) const
{
   return this -> getValue( tnlTuple< 2, Real >( x, y ) );
}

template< int Dimensions, typename Real, tnlDevice Device, typename Index >
Real tnlGrid< Dimensions, Real, Device, Index > :: getValue( const Real& x,
                                                             const Real& y,
                                                             const Real& z ) const
{
   return this -> getValue( tnlTuple< 3, Real >( x, y, z ) );
}

template< int Dimensions, typename Real, tnlDevice Device, typename Index >
Real tnlGrid< Dimensions, Real, Device, Index > :: Partial_x_f( const Index i1 ) const
{
   tnlAssert( Dimensions == 1,
              cerr << "The array " << this -> getName()
                   << " has " << Dimensions << " but 1 is expected." << endl; );
   tnlAssert( i1 >= 0 &&
              i1 < ( this -> getDimensions(). x() - 1 ),
              cerr << " i1 = " << i1 << " and it should be in [ 0, " <<
                   ( this -> getDimensions(). y() - 1  ) << " )" << endl; );

   const Real& Hx = spaceSteps[ 0 ];
   tnlAssert( Hx > 0, cerr << "Hx = " << Hx << endl; );
   return ( tnlArray< 1, Real, tnlHost, Index > :: getElement( i1 + 1 ) -
            tnlArray< 1, Real, tnlHost, Index > ::  getElement( i1 ) ) / Hx;
};

template< int Dimensions, typename Real, tnlDevice Device, typename Index >
Real tnlGrid< Dimensions, Real, Device, Index > :: Partial_x_f( const Index i1,
                                                                const Index i2 ) const
{
   tnlAssert( Dimensions == 2,
              cerr << "The array " << this -> getName()
                   << " has " << Dimensions << " but 2 are expected." << endl; );
   tnlAssert( i1 >= 0 && i2 >= 0 &&
              i1 < this -> getDimensions(). x() - 1 && i2 <= this -> getDimensions(). y() - 1,
              cerr << " i1 = " << i1 << " and it should be in [ 0, " <<
                      ( this -> getDimensions() .x() - 1  ) << " ) " << endl
                   << " i2 = " << i2 << " and it should be in [ 0, " <<
                      ( this -> getDimensions(). y() - 1  ) << " ] " << endl; );

   const Real& Hx = spaceSteps[ 0 ];
   tnlAssert( Hx > 0, cerr << "Hx = " << Hx << endl; );
   return ( this -> getElement( i1 + 1, i2 ) -
            this ->  getElement( i1, i2 ) ) / Hx;
};

template< int Dimensions, typename Real, tnlDevice Device, typename Index >
Real tnlGrid< Dimensions, Real, Device, Index > :: Partial_x_f( const Index i1,
                                                                const Index i2,
                                                                const Index i3 ) const
{
   tnlAssert( Dimensions == 3,
              cerr << "The array " << this -> getName()
                   << " has " << Dimensions << " but 3 are expected." << endl; );
   tnlAssert( i1 >= 0 && i2 >= 0 && i3 >= 0 &&
              i1 < this -> getDimensions()[ tnlX ] - 1 &&
              i2 <= this -> getDimensions()[ tnlY ] - 1 &&
              i3 <= this -> getDimensions()[ tnlZ ] - 1,
              cerr << " i1 = " << i1 << " and it should be in [ 0, " <<
                      ( this -> getDimensions()[ tnlX ] - 1  ) << " ) " << endl
                   << " i2 = " << i2 << " and it should be in [ 0, " <<
                      ( this -> getDimensions()[ tnlY ] - 1  ) << " ] " << endl
                   << " i3 = " << i3 << " and it should be in [ 0, " <<
                      ( this -> getDimensions()[ tnlZ ] - 1  ) << " ] " << endl; );

   const Real& Hx = spaceSteps[ 0 ];
   tnlAssert( Hx > 0, cerr << "Hx = " << Hx << endl; );
   return ( this -> getElement( i1 + 1, i2, i3 ) -
            this ->  getElement( i1, i2, i3 ) ) / Hx;
};

template< int Dimensions, typename Real, tnlDevice Device, typename Index >
Real tnlGrid< Dimensions, Real, Device, Index > :: Partial_x_b( const Index i1 ) const
{
   tnlAssert( Dimensions == 1,
              cerr << "The array " << this -> getName()
                   << " has " << Dimensions << " but 1 is expected." << endl; );
   tnlAssert( i1 > 0 &&
              i1 <= ( tnlArray< 1, Real, tnlHost, Index > :: getDimensions()[ tnlX ] - 1 ),
              cerr << " i1 = " << i1 << " and it should be in ( 0, " <<
                   ( tnlArray< 1, Real, tnlHost, Index > :: getDimensions()[ tnlX ] - 1  ) << " ] " << endl; );

   const Real& Hx = spaceSteps[ 0 ];
   tnlAssert( Hx > 0, cerr << "Hx = " << Hx << endl; );
   return ( tnlArray< 1, Real, tnlHost, Index > :: getElement( i1 ) -
            tnlArray< 1, Real, tnlHost, Index > :: getElement( i1 - 1 ) ) / Hx;
};

template< int Dimensions, typename Real, tnlDevice Device, typename Index >
Real tnlGrid< Dimensions, Real, Device, Index > :: Partial_x_b( const Index i1,
                                                                const Index i2 ) const
{
   tnlAssert( Dimensions == 2,
              cerr << "The array " << this -> getName()
                   << " has " << Dimensions << " but 2 are expected." << endl; );
   tnlAssert( i1 > 0 && i2 >= 0 &&
              i1 <= this -> getDimensions()[ tnlX ] - 1 && i2 <= this -> getDimensions()[ tnlY ] - 1,
              cerr << " i1 = " << i1 << " and it should be in ( 0, " <<
                      ( this -> getDimensions()[ tnlX ] - 1  ) << " ] " << endl
                   << " i2 = " << i2 << " and it should be in [ 0, " <<
                      ( this -> getDimensions()[ tnlY ] - 1  ) << " ] " << endl; );

   const Real& Hx = spaceSteps[ 0 ];
   tnlAssert( Hx > 0, cerr << "Hx = " << Hx << endl; );
   return ( this -> getElement( i1, i2 ) -
            this -> getElement( i1 - 1, i2 ) ) / Hx;
};

template< int Dimensions, typename Real, tnlDevice Device, typename Index >
Real tnlGrid< Dimensions, Real, Device, Index > :: Partial_x_b( const Index i1,
                                                                const Index i2,
                                                                const Index i3 ) const
{
   tnlAssert( Dimensions == 3,
              cerr << "The array " << this -> getName()
                   << " has " << Dimensions << " but 3 are expected." << endl; );
   tnlAssert( i1 > 0 && i2 >= 0 && i3 >= 0 &&
              i1 <= this -> getDimensions()[ tnlX ] - 1 &&
              i2 <= this -> getDimensions()[ tnlY ] - 1 &&
              i3 <= this -> getDimensions()[ tnlZ ] - 1,
              cerr << " i1 = " << i1 << " and it should be in ( 0, " <<
                      ( this -> getDimensions()[ tnlX ] - 1  ) << " ] " << endl
                   << " i2 = " << i2 << " and it should be in [ 0, " <<
                      ( this -> getDimensions()[ tnlY ] - 1  ) << " ] " << endl
                   << " i3 = " << i3 << " and it should be in [ 0, " <<
                      ( this -> getDimensions()[ tnlZ ] - 1  ) << " ] " << endl; );

   const Real& Hx = spaceSteps[ 0 ];
   tnlAssert( Hx > 0, cerr << "Hx = " << Hx << endl; );
   return ( this -> getElement( i1, i2, i3 ) -
            this ->  getElement( i1 - 1, i2, i3 ) ) / Hx;
};

template< int Dimensions, typename Real, tnlDevice Device, typename Index >
Real tnlGrid< Dimensions, Real, Device, Index > :: Partial_x( const Index i1 ) const
{
   tnlAssert( Dimensions == 1,
              cerr << "The array " << this -> getName()
                   << " has " << Dimensions << " but 1 is expected." << endl; );
   tnlAssert( i1 > 0 &&
              i1 < ( this -> getDimensions()[ tnlX ] - 1 ),
              cerr << " i1 = " << i1 << " and it should be in ( 0, " <<
                   ( this -> getDimensions()[ tnlX ] - 1  ) << " ) " << endl; );

   const Real& Hx = spaceSteps[ 0 ];
   tnlAssert( Hx > 0, cerr << "Hx = " << Hx << endl; );
   return ( tnlArray< 1, Real, tnlHost, Index > :: getElement( i1 + 1 ) -
            tnlArray< 1, Real, tnlHost, Index > :: getElement( i1 - 1 ) ) / ( 2.0 * Hx );
};

template< int Dimensions, typename Real, tnlDevice Device, typename Index >
Real tnlGrid< Dimensions, Real, Device, Index > :: Partial_x( const Index i1,
                                                              const Index i2 ) const
{
   tnlAssert( Dimensions == 2,
              cerr << "The array " << this -> getName()
                   << " has " << Dimensions << " but 2 are expected." << endl; );
   tnlAssert( i1 > 0 && i2 >= 0 &&
              i1 < this -> getDimensions()[ tnlX ] - 1 && i2 <= this -> getDimensions()[ tnlY ] - 1,
              cerr << " i1 = " << i1 << " and it should be in ( 0, " <<
                      ( this -> getDimensions()[ tnlX ] - 1  ) << " ) " << endl
                   << " i2 = " << i2 << " and it should be in [ 0, " <<
                      ( this -> getDimensions()[ tnlY ] - 1  ) << " ] " << endl; );

   const Real& Hx = spaceSteps[ 0 ];
   tnlAssert( Hx > 0, cerr << "Hx = " << Hx << endl; );
   return ( this -> getElement( i1 + 1, i2 ) -
            this -> getElement( i1 - 1, i2 ) ) / ( 2.0 * Hx );
};

template< int Dimensions, typename Real, tnlDevice Device, typename Index >
Real tnlGrid< Dimensions, Real, Device, Index > :: Partial_x( const Index i1,
                                                              const Index i2,
                                                              const Index i3 ) const
{
   tnlAssert( Dimensions == 3,
              cerr << "The array " << this -> getName()
                   << " has " << Dimensions << " but 3 are expected." << endl; );
   tnlAssert( i1 > 0 && i2 >= 0 && i3 >= 0 &&
              i1 < this -> getDimensions()[ tnlX ] - 1 &&
              i2 <= this -> getDimensions()[ tnlY ] - 1 &&
              i3 <= this -> getDimensions()[ tnlZ ] - 1,
              cerr << " i1 = " << i1 << " and it should be in ( 0, " <<
                      ( this -> getDimensions()[ tnlX ] - 1  ) << " ) " << endl
                   << " i2 = " << i2 << " and it should be in [ 0, " <<
                      ( this -> getDimensions()[ tnlY ] - 1  ) << " ] " << endl
                   << " i3 = " << i3 << " and it should be in [ 0, " <<
                      ( this -> getDimensions()[ tnlZ ] - 1  ) << " ] " << endl; );

   const Real& Hx = spaceSteps[ 0 ];
   tnlAssert( Hx > 0, cerr << "Hx = " << Hx << endl; );
   return ( this -> getElement( i1 + 1, i2, i3 ) -
            this -> getElement( i1 - 1, i2, i3 ) ) / ( 2.0 * Hx );
};

template< int Dimensions, typename Real, tnlDevice Device, typename Index >
Real tnlGrid< Dimensions, Real, Device, Index > :: Partial_xx( const Index i1 ) const
{
   tnlAssert( Dimensions == 1,
              cerr << "The array " << this -> getName()
                   << " has " << Dimensions << " but 1 is expected." << endl; );
   tnlAssert( i1 > 0 &&
              i1 < ( tnlArray< 1, Real, tnlHost, Index > :: getDimensions()[ tnlX ] - 1 ),
              cerr << " i1 = " << i1 << " and it should be in ( 0, " <<
                   ( tnlArray< 1, Real, tnlHost, Index > :: getDimensions()[ tnlX ] - 1  ) << " ) " << endl; );

   const Real& Hx = spaceSteps[ 0 ];
   tnlAssert( Hx > 0, cerr << "Hx = " << Hx << endl; );
   return ( tnlArray< 1, Real, tnlHost, Index > :: getElement( i1 + 1 ) -
            2.0 * tnlArray< 1, Real, tnlHost, Index > :: getElement( i1 ) +
            tnlArray< 1, Real, tnlHost, Index > :: getElement( i1 - 1 ) ) / ( Hx * Hx );
};

template< int Dimensions, typename Real, tnlDevice Device, typename Index >
Real tnlGrid< Dimensions, Real, Device, Index > :: Partial_xx( const Index i1,
                                                               const Index i2 ) const
{
   tnlAssert( Dimensions == 2,
              cerr << "The array " << this -> getName()
                   << " has " << Dimensions << " but 2 are expected." << endl; );
   tnlAssert( i1 > 0 && i2 >= 0 &&
              i1 < this -> getDimensions()[ tnlX ] - 1 && i2 <= this -> getDimensions()[ tnlY ] - 1,
              cerr << " i1 = " << i1 << " and it should be in ( 0, " <<
                      ( this -> getDimensions()[ tnlX ] - 1  ) << " ) " << endl
                   << " i2 = " << i2 << " and it should be in [ 0, " <<
                      ( this -> getDimensions()[ tnlY ] - 1  ) << " ] " << endl; );

   const Real& Hx = spaceSteps[ 0 ];
   tnlAssert( Hx > 0, cerr << "Hx = " << Hx << endl; );
   return ( this -> getElement( i1 + 1, i2 ) -
            2.0 * this -> getElement( i1, i2 ) +
            this -> getElement( i1 - 1, i2 ) ) / ( Hx * Hx );
};

template< int Dimensions, typename Real, tnlDevice Device, typename Index >
Real tnlGrid< Dimensions, Real, Device, Index > :: Partial_xx( const Index i1,
                                                               const Index i2,
                                                               const Index i3 ) const
{
   tnlAssert( Dimensions == 3,
              cerr << "The array " << this -> getName()
                   << " has " << Dimensions << " but 3 are expected." << endl; );
   tnlAssert( i1 > 0 && i2 >= 0 && i3 >= 0 &&
              i1 < this -> getDimensions()[ tnlX ] - 1 &&
              i2 <= this -> getDimensions()[ tnlY ] - 1 &&
              i3 <= this -> getDimensions()[ tnlZ ] - 1,
              cerr << " i1 = " << i1 << " and it should be in ( 0, " <<
                      ( this -> getDimensions()[ tnlX ] - 1  ) << " ) " << endl
                   << " i2 = " << i2 << " and it should be in [ 0, " <<
                      ( this -> getDimensions()[ tnlY ] - 1  ) << " ] " << endl
                   << " i3 = " << i3 << " and it should be in [ 0, " <<
                      ( this -> getDimensions()[ tnlZ ] - 1  ) << " ] " << endl; );

   const Real& Hx = spaceSteps[ 0 ];
   tnlAssert( Hx > 0, cerr << "Hx = " << Hx << endl; );
   return ( this -> getElement( i1 + 1, i2, i3 ) -
            2.0 * this -> getElement( i1, i2, i3 ) +
            this -> getElement( i1 - 1, i2, i3 ) ) / ( Hx * Hx );
};

template< int Dimensions, typename Real, tnlDevice Device, typename Index >
Real tnlGrid< Dimensions, Real, Device, Index > :: Partial_y_f( const Index i1,
                                                                const Index i2 ) const
{
   tnlAssert( Dimensions == 2,
              cerr << "The array " << this -> getName()
                   << " has " << Dimensions << " but 2 are expected." << endl; );
   tnlAssert( i1 >= 0 && i2 >= 0 &&
              i1 <= this -> getDimensions()[ tnlX ] - 1 && i2 < this -> getDimensions()[ tnlY ] - 1,
              cerr << " i1 = " << i1 << " and it should be in [ 0, " <<
                      ( this -> getDimensions()[ tnlX ] - 1  ) << " ] " << endl
                   << " i2 = " << i2 << " and it should be in [ 0, " <<
                      ( this -> getDimensions()[ tnlY ] - 1  ) << " ) " << endl; );

   const Real& Hy = spaceSteps[ 1 ];
   tnlAssert( Hy > 0, cerr << "Hy = " << Hy << endl; );
   return ( this -> getElement( i1, i2 + 1 ) -
            this -> getElement( i1, i2 ) ) / Hy;
};

template< int Dimensions, typename Real, tnlDevice Device, typename Index >
Real tnlGrid< Dimensions, Real, Device, Index > :: Partial_y_f( const Index i1,
                                                                const Index i2,
                                                                const Index i3 ) const
{
   tnlAssert( Dimensions == 3,
              cerr << "The array " << this -> getName()
                   << " has " << Dimensions << " but 3 are expected." << endl; );
   tnlAssert( i1 >= 0 && i2 >= 0 && i3 >= 0 &&
              i1 <= this -> getDimensions()[ tnlX ] - 1 &&
              i2 < this -> getDimensions()[ tnlY ] - 1 &&
              i3 <= this -> getDimensions()[ tnlZ ] - 1,
              cerr << " i1 = " << i1 << " and it should be in [ 0, " <<
                      ( this -> getDimensions()[ tnlX ] - 1  ) << " ] " << endl
                   << " i2 = " << i2 << " and it should be in [ 0, " <<
                      ( this -> getDimensions()[ tnlY ] - 1  ) << " ) " << endl
                   << " i3 = " << i3 << " and it should be in [ 0, " <<
                      ( this -> getDimensions()[ tnlZ ] - 1  ) << " ] " << endl; );

   const Real& Hy = spaceSteps[ 1 ];
   tnlAssert( Hy > 0, cerr << "Hy = " << Hy << endl; );
   return ( this -> getElement( i1, i2 + 1, i3 ) -
            this -> getElement( i1, i2, i3 ) ) / Hy;
};

template< int Dimensions, typename Real, tnlDevice Device, typename Index >
Real tnlGrid< Dimensions, Real, Device, Index > :: Partial_y_b( const Index i1,
                                                                const Index i2 ) const
{
   tnlAssert( Dimensions == 2,
              cerr << "The array " << this -> getName()
                   << " has " << Dimensions << " but 2 are expected." << endl; );
   tnlAssert( i1 >= 0 && i2 > 0 &&
              i1 <= this -> getDimensions()[ tnlX ] - 1 && i2 <= this -> getDimensions()[ tnlY ] - 1,
              cerr << " i1 = " << i1 << " and it should be in [ 0, " <<
                      ( this -> getDimensions()[ tnlX ] - 1  ) << " ] " << endl
                   << " i2 = " << i2 << " and it should be in ( 0, " <<
                      ( this -> getDimensions()[ tnlY ] - 1  ) << " ] " << endl; );

   const Real& Hy = spaceSteps[ 1 ];
   tnlAssert( Hy > 0, cerr << "Hy = " << Hy << endl; );
   return ( this -> getElement( i1, i2 ) -
            this -> getElement( i1, i2 - 1 ) ) / Hy;
};

template< int Dimensions, typename Real, tnlDevice Device, typename Index >
Real tnlGrid< Dimensions, Real, Device, Index > :: Partial_y_b( const Index i1,
                                                                const Index i2,
                                                                const Index i3 ) const
{
   tnlAssert( Dimensions == 3,
              cerr << "The array " << this -> getName()
                   << " has " << Dimensions << " but 3 are expected." << endl; );
   tnlAssert( i1 >= 0 && i2 > 0 && i3 >= 0 &&
              i1 <= this -> getDimensions()[ tnlX ] - 1 &&
              i2 <= this -> getDimensions()[ tnlY ] - 1 &&
              i3 <= this -> getDimensions()[ tnlZ ] - 1,
              cerr << " i1 = " << i1 << " and it should be in [ 0, " <<
                      ( this -> getDimensions()[ tnlX ] - 1  ) << " ] " << endl
                   << " i2 = " << i2 << " and it should be in ( 0, " <<
                      ( this -> getDimensions()[ tnlY ] - 1  ) << " ] " << endl
                   << " i3 = " << i3 << " and it should be in [ 0, " <<
                      ( this -> getDimensions()[ tnlZ ] - 1  ) << " ] " << endl; );

   const Real& Hy = spaceSteps[ 1 ];
   tnlAssert( Hy > 0, cerr << "Hy = " << Hy << endl; );
   return ( this -> getElement( i1, i2, i3 ) -
            this -> getElement( i1, i2 - 1, i3 ) ) / Hy;
};

template< int Dimensions, typename Real, tnlDevice Device, typename Index >
Real tnlGrid< Dimensions, Real, Device, Index > :: Partial_y( const Index i1,
                                                              const Index i2 ) const
{
   tnlAssert( Dimensions == 2,
              cerr << "The array " << this -> getName()
                   << " has " << Dimensions << " but 2 are expected." << endl; );
   tnlAssert( i1 >= 0 && i2 > 0 &&
              i1 <= this -> getDimensions()[ tnlX ] - 1 && i2 < this -> getDimensions()[ tnlY ] - 1,
              cerr << " i1 = " << i1 << " and it should be in [ 0, " <<
                      ( this -> getDimensions()[ tnlX ] - 1  ) << " ] " << endl
                   << " i2 = " << i2 << " and it should be in ( 0, " <<
                      ( this -> getDimensions()[ tnlY ] - 1  ) << " ) " << endl; );

   const Real& Hy = spaceSteps[ 1 ];
   tnlAssert( Hy > 0, cerr << "Hy = " << Hy << endl; );
   return ( this -> getElement( i1, i2 + 1 ) -
            this -> getElement( i1, i2 - 1 ) ) / ( 2.0 * Hy );
};

template< int Dimensions, typename Real, tnlDevice Device, typename Index >
Real tnlGrid< Dimensions, Real, Device, Index > :: Partial_y( const Index i1,
                                                              const Index i2,
                                                              const Index i3 ) const
{
   tnlAssert( Dimensions == 3,
              cerr << "The array " << this -> getName()
                   << " has " << Dimensions << " but 3 are expected." << endl; );
   tnlAssert( i1 >= 0 && i2 > 0 && i3 >= 0 &&
              i1 <= this -> getDimensions()[ tnlX ] - 1 &&
              i2 < this -> getDimensions()[ tnlY ] - 1 &&
              i3 <= this -> getDimensions()[ tnlZ ] - 1,
              cerr << " i1 = " << i1 << " and it should be in [ 0, " <<
                      ( this -> getDimensions()[ tnlX ] - 1  ) << " ] " << endl
                   << " i2 = " << i2 << " and it should be in ( 0, " <<
                      ( this -> getDimensions()[ tnlY ] - 1  ) << " ) " << endl
                   << " i3 = " << i3 << " and it should be in [ 0, " <<
                      ( this -> getDimensions()[ tnlZ ] - 1  ) << " ] " << endl; );

   const Real& Hy = spaceSteps[ 1 ];
   tnlAssert( Hy > 0, cerr << "Hy = " << Hy << endl; );
   return ( this -> getElement( i1, i2 + 1, i3 ) -
            this -> getElement( i1, i2 - 1, i3 ) ) / ( 2.0 * Hy );
};

template< int Dimensions, typename Real, tnlDevice Device, typename Index >
Real tnlGrid< Dimensions, Real, Device, Index > :: Partial_yy( const Index i1,
                                                               const Index i2 ) const
{
   tnlAssert( Dimensions == 2,
              cerr << "The array " << this -> getName()
                   << " has " << Dimensions << " but 2 are expected." << endl; );
   tnlAssert( i1 >= 0 && i2 > 0 &&
              i1 <= this -> getDimensions()[ tnlX ] - 1 && i2 < this -> getDimensions()[ tnlY ] - 1,
              cerr << " i1 = " << i1 << " and it should be in [ 0, " <<
                      ( this -> getDimensions()[ tnlX ] - 1  ) << " ] " << endl
                   << " i2 = " << i2 << " and it should be in ( 0, " <<
                      ( this -> getDimensions()[ tnlY ] - 1  ) << " ) " << endl; );

   const Real& Hy = spaceSteps[ 1 ];
   tnlAssert( Hy > 0, cerr << "Hy = " << Hy << endl; );
   return ( this -> getElement( i1, i2 + 1 ) -
            2.0 * this -> getElement( i1, i2 ) +
            this -> getElement( i1, i2 - 1 ) ) / ( Hy * Hy );
};

template< int Dimensions, typename Real, tnlDevice Device, typename Index >
Real tnlGrid< Dimensions, Real, Device, Index > :: Partial_yy( const Index i1,
                                                               const Index i2,
                                                               const Index i3 ) const
{
   tnlAssert( Dimensions == 3,
              cerr << "The array " << this -> getName()
                   << " has " << Dimensions << " but 3 are expected." << endl; );
   tnlAssert( i1 >= 0 && i2 > 0 && i3 >= 0 &&
              i1 <= this -> getDimensions()[ tnlX ] - 1 &&
              i2 < this -> getDimensions()[ tnlY ] - 1 &&
              i3 <= this -> getDimensions()[ tnlZ ] - 1,
              cerr << " i1 = " << i1 << " and it should be in [ 0, " <<
                      ( this -> getDimensions()[ tnlX ] - 1  ) << " ] " << endl
                   << " i2 = " << i2 << " and it should be in ( 0, " <<
                      ( this -> getDimensions()[ tnlY ] - 1  ) << " ) " << endl
                   << " i3 = " << i3 << " and it should be in [ 0, " <<
                      ( this -> getDimensions()[ tnlZ ] - 1  ) << " ] " << endl; );

   const Real& Hy = spaceSteps[ 1 ];
   tnlAssert( Hy > 0, cerr << "Hy = " << Hy << endl; );
   return ( this -> getElement( i1, i2 + 1, i3 ) -
            2.0 * this -> getElement( i1, i2, i3 ) +
            this -> getElement( i1, i2 - 1, i3 ) ) / ( Hy * Hy );
};

template< int Dimensions, typename Real, tnlDevice Device, typename Index >
Real tnlGrid< Dimensions, Real, Device, Index > :: Partial_z_f( const Index i1,
                                                                const Index i2,
                                                                const Index i3 ) const
{
   tnlAssert( Dimensions == 3,
              cerr << "The array " << this -> getName()
                   << " has " << Dimensions << " but 3 are expected." << endl; );
   tnlAssert( i1 >= 0 && i2 >= 0 && i3 >= 0 &&
              i1 <= this -> getDimensions()[ tnlX ] - 1 &&
              i2 <= this -> getDimensions()[ tnlY ] - 1 &&
              i3 < this -> getDimensions()[ tnlZ ] - 1,
              cerr << " i1 = " << i1 << " and it should be in [ 0, " <<
                      ( this -> getDimensions()[ tnlX ] - 1  ) << " ] " << endl
                   << " i2 = " << i2 << " and it should be in [ 0, " <<
                      ( this -> getDimensions()[ tnlY ] - 1  ) << " ] " << endl
                   << " i3 = " << i3 << " and it should be in [ 0, " <<
                      ( this -> getDimensions()[ tnlZ ] - 1  ) << " ) " << endl; );

   const Real& Hz = spaceSteps[ 2 ];
   tnlAssert( Hz > 0, cerr << "Hz = " << Hz << endl; );
   return ( this -> getElement( i1, i2, i3 + 1 ) -
            this -> getElement( i1, i2, i3 ) ) / Hz;
};

template< int Dimensions, typename Real, tnlDevice Device, typename Index >
Real tnlGrid< Dimensions, Real, Device, Index > :: Partial_z_b( const Index i1,
                                                                const Index i2,
                                                                const Index i3 ) const
{
   tnlAssert( Dimensions == 3,
              cerr << "The array " << this -> getName()
                   << " has " << Dimensions << " but 3 are expected." << endl; );
   tnlAssert( i1 >= 0 && i2 >= 0 && i3 > 0 &&
              i1 <= this -> getDimensions()[ tnlX ] - 1 &&
              i2 <= this -> getDimensions()[ tnlY ] - 1 &&
              i3 <= this -> getDimensions()[ tnlZ ] - 1,
              cerr << " i1 = " << i1 << " and it should be in [ 0, " <<
                      ( this -> getDimensions()[ tnlX ] - 1  ) << " ] " << endl
                   << " i2 = " << i2 << " and it should be in [ 0, " <<
                      ( this -> getDimensions()[ tnlY ] - 1  ) << " ] " << endl
                   << " i3 = " << i3 << " and it should be in ( 0, " <<
                      ( this -> getDimensions()[ tnlZ ] - 1  ) << " ] " << endl; );

   const Real& Hz = spaceSteps[ 2 ];
   tnlAssert( Hz > 0, cerr << "Hz = " << Hz << endl; );
   return ( this -> getElement( i1, i2, i3 ) -
            this -> getElement( i1, i2, i3 - 1 ) ) / Hz;
};

template< int Dimensions, typename Real, tnlDevice Device, typename Index >
Real tnlGrid< Dimensions, Real, Device, Index > :: Partial_z( const Index i1,
                                                              const Index i2,
                                                              const Index i3 ) const
{
   tnlAssert( Dimensions == 3,
              cerr << "The array " << this -> getName()
                   << " has " << Dimensions << " but 3 are expected." << endl; );
   tnlAssert( i1 >= 0 && i2 >= 0 && i3 > 0 &&
              i1 <= this -> getDimensions()[ tnlX ] - 1 &&
              i2 <= this -> getDimensions()[ tnlY ] - 1 &&
              i3 < this -> getDimensions()[ tnlZ ] - 1,
              cerr << " i1 = " << i1 << " and it should be in [ 0, " <<
                      ( this -> getDimensions()[ tnlX ] - 1  ) << " ] " << endl
                   << " i2 = " << i2 << " and it should be in [ 0, " <<
                      ( this -> getDimensions()[ tnlY ] - 1  ) << " ] " << endl
                   << " i3 = " << i3 << " and it should be in ( 0, " <<
                      ( this -> getDimensions()[ tnlZ ] - 1  ) << " ) " << endl; );

   const Real& Hz = spaceSteps[ 2 ];
   tnlAssert( Hz > 0, cerr << "Hz = " << Hz << endl; );
   return ( this -> getElement( i1, i2, i3 + 1 ) -
            this -> getElement( i1, i2, i3 - 1 ) ) / ( 2.0 * Hz );
};

template< int Dimensions, typename Real, tnlDevice Device, typename Index >
Real tnlGrid< Dimensions, Real, Device, Index > :: Partial_zz( const Index i1,
                                                               const Index i2,
                                                               const Index i3 ) const
{
   tnlAssert( Dimensions == 3,
              cerr << "The array " << this -> getName()
                   << " has " << Dimensions << " but 3 are expected." << endl; );
   tnlAssert( i1 >= 0 && i2 >= 0 && i3 > 0 &&
              i1 <= this -> getDimensions()[ tnlX ] - 1 &&
              i2 <= this -> getDimensions()[ tnlY ] - 1 &&
              i3 < this -> getDimensions()[ tnlZ ] - 1,
              cerr << " i1 = " << i1 << " and it should be in [ 0, " <<
                      ( this -> getDimensions()[ tnlX ] - 1  ) << " ] " << endl
                   << " i2 = " << i2 << " and it should be in [ 0, " <<
                      ( this -> getDimensions()[ tnlY ] - 1  ) << " ] " << endl
                   << " i3 = " << i3 << " and it should be in ( 0, " <<
                      ( this -> getDimensions()[ tnlZ ] - 1  ) << " ) " << endl; );

   const Real& Hz = spaceSteps[ 2 ];
   tnlAssert( Hz > 0, cerr << "Hz = " << Hz << endl; );
   return ( this -> getElement( i1, i2, i3 + 1 ) -
            2.0 * this -> getElement( i1, i2, i3 ) +
            this -> getElement( i1, i2, i3 - 1 ) ) / ( Hz * Hz );
};

template< int Dimensions, typename Real, tnlDevice Device, typename Index >
void tnlGrid< Dimensions, Real, Device, Index > :: setDirichletBC( const tnlGrid< Dimensions, Real, Device, Index >&bc,
                                                                   const tnlTuple< Dimensions, bool >& lowerBC,
                                                                   const tnlTuple< Dimensions, bool >& upperBC )
{
   if( Device == tnlHost )
   {
      if( Dimensions == 1 )
      {
         const Index xSize = this -> getDimensions(). x();
         if( lowerBC. x() )
            ( *this )( 0 ) = bc( 0 );
         if( upperBC. x() )
            ( *this )( xSize - 1 ) = bc( xSize - 1 );
      }
      if( Dimensions == 2 )
      {
         const Index xSize = this -> getDimensions(). x();
         const Index ySize = this -> getDimensions(). y();
         Index i, j;

         if( lowerBC. x() )
            for( i = 0; i < xSize; i ++ )
               ( *this )( i, 0 ) = bc( i, 0 );
         if( upperBC. x() )
            for( i = 0; i < xSize; i ++ )
               ( *this )( i, ySize - 1 ) = bc( i, ySize - 1 );

         if( lowerBC. y() )
            for( j = 0; j < ySize; j ++ )
               ( *this )( 0, j ) = bc( 0, j );
         if( upperBC. y() )
            for( j = 0; j < ySize; j ++ )
               ( *this )( xSize - 1, j ) = bc( xSize - 1, j );
      }
      if( Dimensions == 3 )
      {
         const Index xSize = this -> getDimensions(). x();
         const Index ySize = this -> getDimensions(). y();
         const Index zSize = this -> getDimensions(). z();
         Index i, j, k;

         if( lowerBC. x() )
            for( i = 0; i < xSize; i ++ )
               for( k = 0; k < zSize; k ++ )
                  ( *this )( i, 0, k ) = bc( i, 0, k );

         if( upperBC. x() )
            for( i = 0; i < xSize; i ++ )
               for( k = 0; k < zSize; k ++ )
                  ( *this )( i, ySize - 1, k ) = ( i, ySize - 1, k );

         if( lowerBC. y() )
            for( j = 0; j < ySize; j ++ )
               for( k = 0; k < zSize; k ++ )
                  ( *this )( 0, j, k ) = bc( 0, j, k );

         if( upperBC. y() )
            for( j = 0; j < ySize; j ++ )
               for( k = 0; k < zSize; k ++ )
                  ( *this )( xSize - 1, j, k ) = bc( xSize - 1, j, k );

         if( lowerBC. z() )
            for( i = 0; i < xSize; i ++ )
               for( j = 0; j < ySize; j ++ )
                  ( *this )( i, j, 0 ) = bc( i, j, 0 );

         if( upperBC. z() )
            for( i = 0; i < xSize; i ++ )
               for( j = 0; j < ySize; j ++ )
                  ( *this )( i, j, zSize - 1 ) = bc( i, j, zSize - 1 );
      }
   }
   if( Device == tnlCuda )
   {
#ifdef HAVE_CUDA
      if( Dimensions == 2 )
      {
         const Index size = this -> getSize();
         const Index desBlockSize = 64;
         const Index gridSize = size / desBlockSize + ( size % desBlockSize != 0 );
         dim3 gridDim( 0 ), blockDim( 0 );
         gridDim. x = gridSize;
         blockDim. x = desBlockSize;
         :: setDirichletBC< 2, Real, Index ><<< gridDim, blockDim >>>( this -> getDimensions(). x(),
                                                                       this -> getDimensions(). y(),
                                                                       bc. getData(),
                                                                       this -> getData() );
      }
#endif
   }
}

template< int Dimensions, typename Real, tnlDevice Device, typename Index >
void tnlGrid< Dimensions, Real, Device, Index > :: setDirichletBC( const Real& bcValue,
                                                                   const tnlTuple< Dimensions, bool >& lowerBC,
                                                                   const tnlTuple< Dimensions, bool >& upperBC )
{
   if( Device == tnlHost )
   {
      if( Dimensions == 1 )
      {
         const Index xSize = this -> getDimensions(). x();
         if( lowerBC. x() )
            ( *this )( 0 ) = bcValue;
         if( upperBC. x() )
            ( *this )( xSize - 1 ) = bcValue;
      }
      if( Dimensions == 2 )
      {
         const Index xSize = this -> getDimensions(). x();
         const Index ySize = this -> getDimensions(). y();
         Index i, j;

         if( lowerBC. x() )
            for( i = 0; i < xSize; i ++ )
               ( *this )( i, 0 ) = bcValue;
         if( upperBC. x() )
            for( i = 0; i < xSize; i ++ )
               ( *this )( i, ySize - 1 ) = bcValue;

         if( lowerBC. y() )
            for( j = 0; j < ySize; j ++ )
               ( *this )( 0, j ) = bcValue;
         if( upperBC. y() )
            for( j = 0; j < ySize; j ++ )
               ( *this )( xSize - 1, j ) = bcValue;
      }
      if( Dimensions == 3 )
      {
         const Index xSize = this -> getDimensions(). x();
         const Index ySize = this -> getDimensions(). y();
         const Index zSize = this -> getDimensions(). z();
         Index i, j, k;

         if( lowerBC. x() )
            for( i = 0; i < xSize; i ++ )
               for( k = 0; k < zSize; k ++ )
                  ( *this )( i, 0, k ) = bcValue;

         if( upperBC. x() )
            for( i = 0; i < xSize; i ++ )
               for( k = 0; k < zSize; k ++ )
                  ( *this )( i, ySize - 1, k ) = bcValue;

         if( lowerBC. y() )
            for( j = 0; j < ySize; j ++ )
               for( k = 0; k < zSize; k ++ )
                  ( *this )( 0, j, k ) = bcValue;

         if( upperBC. y() )
            for( j = 0; j < ySize; j ++ )
               for( k = 0; k < zSize; k ++ )
                  ( *this )( xSize - 1, j, k ) = bcValue;

         if( lowerBC. z() )
            for( i = 0; i < xSize; i ++ )
               for( j = 0; j < ySize; j ++ )
                  ( *this )( i, j, 0 ) = bcValue;

         if( upperBC. z() )
            for( i = 0; i < xSize; i ++ )
               for( j = 0; j < ySize; j ++ )
                  ( *this )( i, j, zSize - 1 ) =  bcValue;
      }
   }
   if( Device == tnlCuda )
   {
#ifdef HAVE_CUDA
      if( Dimensions == 2 )
      {
         const Index size = this -> getSize();
         const Index desBlockSize = 64;
         const Index gridSize = size / desBlockSize + ( size % desBlockSize != 0 );
         dim3 gridDim( 0 ), blockDim( 0 );
         gridDim. x = gridSize;
         blockDim. x = desBlockSize;
         :: setConstantDirichletBC< 2, Real, Index ><<< gridDim, blockDim >>>( this -> getDimensions(). x(),
                                                                               this -> getDimensions(). y(),
                                                                               bcValue,
                                                                               this -> getData() );
      }
#endif
   }
}

template< int Dimensions, typename Real, tnlDevice Device, typename Index >
void tnlGrid< Dimensions, Real, Device, Index > :: setNeumannBC( const tnlGrid< Dimensions, Real, Device, Index >&bc,
                                                                 const tnlTuple< Dimensions, bool >& lowerBC,
                                                                 const tnlTuple< Dimensions, bool >& upperBC )
{
   if( Device == tnlHost )
   {
      if( Dimensions == 1 )
      {
         const Real& hx = this -> getSpaceSteps(). x();

         const Index xSize = this -> getDimensions(). x();
         if( lowerBC. x() )
            ( *this )( 0 ) = ( *this )( 1 ) + hx * bc( 0 );
         if( upperBC. x() )
            ( *this )( xSize - 1 ) = ( *this )( xSize - 1 ) + hx * bc( xSize -1 );
      }
      if( Dimensions == 2 )
      {
         const Real& hx = this -> getSpaceSteps(). x();
         const Real& hy = this -> getSpaceSteps(). y();

         const Index xSize = this -> getDimensions(). x();
         const Index ySize = this -> getDimensions(). y();
         Index i, j;

         if( lowerBC. x() )
            for( i = 0; i < xSize; i ++ )
                ( *this )( i, 0 ) = ( *this )( i, 1 ) + hy * bc( i, 0 );
         if( upperBC. x() )
            for( i = 0; i < xSize; i ++ )
               ( *this )( i, ySize - 1 ) = ( *this )( i, ySize - 2 ) + hy * bc( i, ySize - 1 );

         if( lowerBC. y() )
            for( j = 0; j < ySize; j ++ )
               ( *this )( 0, j ) = ( *this )( 1, j ) + hx * bc( 0, j );
         if( upperBC. y() )
            for( j = 0; j < ySize; j ++ )
               ( *this )( xSize - 1, j ) = ( *this )( xSize - 2, j ) + hx * bc( xSize - 1, j );
      }
      if( Dimensions == 3 )
      {
         const Real& hx = this ->  getSpaceSteps(). x();
         const Real& hy = this ->  getSpaceSteps(). y();
         const Real& hz = this ->  getSpaceSteps(). z();

         const Index xSize = this -> getDimensions(). x();
         const Index ySize = this -> getDimensions(). y();
         const Index zSize = this -> getDimensions(). z();
         Index i, j, k;

         if( lowerBC. x() )
            for( i = 0; i < xSize; i ++ )
               for( k = 0; k < zSize; k ++ )
                  ( *this )( i, 0, k ) = ( *this )( i, 1, k ) + hy * bc. getElement( i, 0, k );
         if( upperBC. x() )
            for( i = 0; i < xSize; i ++ )
               for( k = 0; k < zSize; k ++ )
                  ( *this )( i, ySize - 1, k ) = ( *this )( i, ySize - 2, k ) + hy * bc( i, ySize - 1, k );

         if( lowerBC. y() )
            for( j = 0; j < ySize; j ++ )
               for( k = 0; k < zSize; k ++ )
                  ( *this )( 0, j, k ) =  ( *this )( 1, j, k ) + hx * bc( 0, j, k );

         if( upperBC. y() )
            for( j = 0; j < ySize; j ++ )
               for( k = 0; k < zSize; k ++ )
                  ( *this )( xSize - 1, j, k ) = ( *this )( xSize - 2, j, k ) + hx * bc( xSize - 1, j, k );

         if( lowerBC. z() )
            for( i = 0; i < xSize; i ++ )
               for( j = 0; j < ySize; j ++ )
                  ( *this )( i, j, 0 ) = ( *this )( i, j, 1 ) + hz * bc( i, j, 0 );

         if( upperBC. z() )
            for( i = 0; i < xSize; i ++ )
               for( j = 0; j < ySize; j ++ )
                  ( *this )( i, j, zSize - 1 ) = ( *this )( i, j, zSize - 2 ) + hz * bc( i, j, zSize - 1 );
      }
   }
   if( Device == tnlCuda )
   {
#ifdef HAVE_CUDA
      if( Dimensions == 2 )
      {
         const Index size = this -> getSize();
         const Index desBlockSize = 64;
         const Index gridSize = size / desBlockSize + ( size % desBlockSize != 0 );
         dim3 gridDim( 0 ), blockDim( 0 );
         gridDim. x = gridSize;
         blockDim. x = desBlockSize;
         ::setNeumannBC< 2, Real, Index ><<< gridDim, blockDim >>>( this -> getDimensions(). x(),
                                                                    this -> getDimensions(). y(),
                                                                    this -> getSpaceSteps(). x(),
                                                                    this -> getSpaceSteps(). y(),
                                                                    bc. getData(),
                                                                    this -> getData() );
      }
#endif
   }
}

template< int Dimensions, typename Real, tnlDevice Device, typename Index >
void tnlGrid< Dimensions, Real, Device, Index > :: setNeumannBC( const Real& bcValue,
                                                                 const tnlTuple< Dimensions, bool >& lowerBC,
                                                                 const tnlTuple< Dimensions, bool >& upperBC )
{
   if( Device == tnlHost )
   {
      if( Dimensions == 1 )
      {
         const Real& hx = this -> getSpaceSteps(). x();

         const Index xSize = this -> getDimensions(). x();
         if( lowerBC. x() )
            ( *this )( 0 ) = ( *this )( 1 ) + hx * bcValue;
         if( upperBC. x() )
            ( *this )( xSize - 1 ) = ( *this )( xSize - 1 ) + hx * bcValue;
      }
      if( Dimensions == 2 )
      {
         const Real& hx = this -> getSpaceSteps(). x();
         const Real& hy = this -> getSpaceSteps(). y();

         const Index xSize = this -> getDimensions(). x();
         const Index ySize = this -> getDimensions(). y();
         Index i, j;

         if( lowerBC. x() )
            for( i = 0; i < xSize; i ++ )
               ( *this )( i, 0 ) = ( *this )( i, 1 ) + hy * bcValue;
         if( upperBC. x() )
            for( i = 0; i < xSize; i ++ )
               ( *this )( i, ySize - 1 ) = ( *this )( i, ySize - 2 ) + hy * bcValue;

         if( lowerBC. y() )
            for( j = 0; j < ySize; j ++ )
               ( *this )( 0, j ) = ( *this )( 1, j ) + hx * bcValue;
         if( upperBC. y() )
            for( j = 0; j < ySize; j ++ )
               ( *this )( xSize - 1, j ) = ( *this )( xSize - 2, j ) + hx * bcValue;
      }
      if( Dimensions == 3 )
      {
         const Real& hx = this -> getSpaceSteps(). x();
         const Real& hy = this -> getSpaceSteps(). y();
         const Real& hz = this -> getSpaceSteps(). z();

         const Index xSize = this -> getDimensions(). x();
         const Index ySize = this -> getDimensions(). y();
         const Index zSize = this -> getDimensions(). z();
         Index i, j, k;

         if( lowerBC. x() )
            for( i = 0; i < xSize; i ++ )
               for( k = 0; k < zSize; k ++ )
                  ( *this )( i, 0, k ) = ( *this )( i, 1, k ) + hy * bcValue;
         if( upperBC. x() )
            for( i = 0; i < xSize; i ++ )
               for( k = 0; k < zSize; k ++ )
                  ( *this )( i, ySize - 1, k ) = ( *this )( i, ySize - 2, k ) + hy * bcValue;

         if( lowerBC. y() )
            for( j = 0; j < ySize; j ++ )
               for( k = 0; k < zSize; k ++ )
                  ( *this )( 0, j, k ) = ( *this )( 1, j, k ) + hx * bcValue;

         if( upperBC. y() )
            for( j = 0; j < ySize; j ++ )
               for( k = 0; k < zSize; k ++ )
                  ( *this )( xSize - 1, j, k ) = ( *this )( xSize - 2, j, k ) + hx * bcValue;

         if( lowerBC. z() )
            for( i = 0; i < xSize; i ++ )
               for( j = 0; j < ySize; j ++ )
                  ( *this )( i, j, 0 ) = ( *this )( i, j, 1 ) + hz * bcValue;

         if( upperBC. z() )
            for( i = 0; i < xSize; i ++ )
               for( j = 0; j < ySize; j ++ )
                  ( *this )( i, j, zSize - 1 ) = ( *this )( i, j, zSize - 2 ) + hz * bcValue;
      }
   }
   if( Device == tnlCuda )
   {
#ifdef HAVE_CUDA
      if( Dimensions == 2 )
      {
         const Index size = this -> getSize();
         const Index desBlockSize = 64;
         const Index gridSize = size / desBlockSize + ( size % desBlockSize != 0 );
         dim3 gridDim( 0 ), blockDim( 0 );
         gridDim. x = gridSize;
         blockDim. x = desBlockSize;
         :: setConstantNeumannBC< 2, Real, Index ><<< gridDim, blockDim >>>( this -> getDimensions(). x(),
                                                                             this -> getDimensions(). y(),
                                                                             this -> getSpaceSteps(). x(),
                                                                             this -> getSpaceSteps(). y(),
                                                                             bcValue,
                                                                             this -> getData() );
      }
#endif
   }
}

template< int Dimensions, typename Real, tnlDevice Device, typename Index >
Real tnlGrid< Dimensions, Real, Device, Index > :: getMax() const
{
   return tnlMax( * this );
}

template< int Dimensions, typename Real, tnlDevice Device, typename Index >
Real tnlGrid< Dimensions, Real, Device, Index > :: getMin() const
{
   return tnlMin( * this );
}

template< int Dimensions, typename Real, tnlDevice Device, typename Index >
Real tnlGrid< Dimensions, Real, Device, Index > :: getAbsMax() const
{
   return tnlAbsMax( * this );
}

template< int Dimensions, typename Real, tnlDevice Device, typename Index >
Real tnlGrid< Dimensions, Real, Device, Index > :: getAbsMin() const
{
   return tnlAbsMin( * this );
}

template< int Dimensions, typename Real, tnlDevice Device, typename Index >
Real tnlGrid< Dimensions, Real, Device, Index > :: getLpNorm() const
{
   Real result = tnlLpNorm( * this );
   if( Dimensions == 1 )
      return result * getSpaceSteps(). x();
   if( Dimensions == 2 )
      return result * getSpaceSteps(). x()
                    * getSpaceSteps(). y();
   if( Dimensions == 3 )
      return result * getSpaceSteps(). x()
                    * getSpaceSteps(). y()
                    * getSpaceSteps(). z();
   for( int i = 0; i < Dimensions; i ++ )
      result *= getSpaceSteps()[ i ];
   return result;
}

template< int Dimensions, typename Real, tnlDevice Device, typename Index >
Real tnlGrid< Dimensions, Real, Device, Index > :: getSum() const
{
   return tnlSum( * this );
}

template< int Dimensions, typename Real, tnlDevice Device, typename Index >
Real tnlGrid< Dimensions, Real, Device, Index > :: getDifferenceMax( const tnlVector< Real, tnlHost, Index >& v ) const
{
   tnlAssert( this -> getDimensions() == v. getDimensions(),
              cerr << "The grid names are " << this -> getName()
                   << " and " << v. getName()
                   << "To get grids with the same parameters use the method setLike." << endl; );
   tnlAssert( this -> getDomainLoweCorner() == v. getDomainLowerCorner() &&
              this -> getDomainUpperCorner() == v. getDomainUpperCorner(),
              cerr << "The grid names are " << this -> getName()
                   << " and " << v. getName()
                   << "To get grids with the same parameters use the method setLike." << endl; );
   return tnlDifferenceMax( *this, v );
}

template< int Dimensions, typename Real, tnlDevice Device, typename Index >
Real tnlGrid< Dimensions, Real, Device, Index > :: getDifferenceMin( const tnlVector< Real, tnlHost, Index >& v ) const
{
   tnlAssert( this -> getDimensions() == v. getDimensions(),
              cerr << "The grid names are " << this -> getName()
                   << " and " << v. getName()
                   << "To get grids with the same parameters use the method setLike." << endl; );
   tnlAssert( this -> getDomainLoweCorner() == v. getDomainLowerCorner() &&
              this -> getDomainUpperCorner() == v. getDomainUpperCorner(),
              cerr << "The grid names are " << this -> getName()
                   << " and " << v. getName()
                   << "To get grids with the same parameters use the method setLike." << endl; );
   return tnlDifferenceMin( *this, v );
}

template< int Dimensions, typename Real, tnlDevice Device, typename Index >
Real tnlGrid< Dimensions, Real, Device, Index > :: getDifferenceAbsMax( const tnlVector< Real, tnlHost, Index >& v ) const
{
   tnlAssert( this -> getDimensions() == v. getDimensions(),
              cerr << "The grid names are " << this -> getName()
                   << " and " << v. getName()
                   << "To get grids with the same parameters use the method setLike." << endl; );
   tnlAssert( this -> getDomainLoweCorner() == v. getDomainLowerCorner() &&
              this -> getDomainUpperCorner() == v. getDomainUpperCorner(),
              cerr << "The grid names are " << this -> getName()
                   << " and " << v. getName()
                   << "To get grids with the same parameters use the method setLike." << endl; );
   return tnlDifferenceAbsMax( *this, v );
}

template< int Dimensions, typename Real, tnlDevice Device, typename Index >
Real tnlGrid< Dimensions, Real, Device, Index > :: getDifferenceAbsMin( const tnlVector< Real, tnlHost, Index >& v ) const
{
   tnlAssert( this -> getDimensions() == v. getDimensions(),
              cerr << "The grid names are " << this -> getName()
                   << " and " << v. getName()
                   << "To get grids with the same parameters use the method setLike." << endl; );
   tnlAssert( this -> getDomainLoweCorner() == v. getDomainLowerCorner() &&
              this -> getDomainUpperCorner() == v. getDomainUpperCorner(),
              cerr << "The grid names are " << this -> getName()
                   << " and " << v. getName()
                   << "To get grids with the same parameters use the method setLike." << endl; );
   return tnlDifferenceAbsMin( *this, v );
}

template< int Dimensions, typename Real, tnlDevice Device, typename Index >
Real tnlGrid< Dimensions, Real, Device, Index > :: getDifferenceLpNorm( const tnlVector< Real, tnlHost, Index >& v, const Real& p ) const
{
   tnlAssert( this -> getDimensions() == v. getDimensions(),
              cerr << "The grid names are " << this -> getName()
                   << " and " << v. getName()
                   << "To get grids with the same parameters use the method setLike." << endl; );
   tnlAssert( this -> getDomainLoweCorner() == v. getDomainLowerCorner() &&
              this -> getDomainUpperCorner() == v. getDomainUpperCorner(),
              cerr << "The grid names are " << this -> getName()
                   << " and " << v. getName()
                   << "To get grids with the same parameters use the method setLike." << endl; );

   Real result = tnlDifferenceLpNorm( * this, v, p );
   if( Dimensions == 1 )
      return result * getSpaceSteps(). x();
   if( Dimensions == 2 )
      return result * getSpaceSteps(). x()
                    * getSpaceSteps(). y();
   if( Dimensions == 3 )
      return result * getSpaceSteps(). x()
                    * getSpaceSteps(). y()
                    * getSpaceSteps(). z();
   for( int i = 0; i < Dimensions; i ++ )
      result *= getSpaceSteps()[ i ];
   return result;
}

template< int Dimensions, typename Real, tnlDevice Device, typename Index >
Real tnlGrid< Dimensions, Real, Device, Index > :: getDifferenceSum( const tnlVector< Real, tnlHost, Index >& v ) const
{
   tnlAssert( this -> getDimensions() == v. getDimensions(),
              cerr << "The grid names are " << this -> getName()
                   << " and " << v. getName()
                   << "To get grids with the same parameters use the method setLike." << endl; );
   tnlAssert( this -> getDomainLoweCorner() == v. getDomainLowerCorner() &&
              this -> getDomainUpperCorner() == v. getDomainUpperCorner(),
              cerr << "The grid names are " << this -> getName()
                   << " and " << v. getName()
                   << "To get grids with the same parameters use the method setLike." << endl; );
   return tnlDifferenceSum( *this, v );
}


template< int Dimensions, typename Real, tnlDevice Device, typename Index >
bool tnlGrid< Dimensions, Real, Device, Index > :: save( tnlFile& file ) const
{
   if( ! tnlArray< Dimensions, Real, Device, Index > :: save( file ) )
   {
      cerr << "I was not able to write the tnlArray of the tnlGrid "
           << this -> getName() << endl;
      return false;
   }
   if( ! domainLowerCorner. save( file ) ||
       ! domainUpperCorner. save( file ) ||
       ! spaceSteps. save( file ) )
   {
      cerr << "I was not able to write the domain decsription of the tnlGrid "
           << this -> getName() << endl;
      return false;
   }
   return true;
}

template< int Dimensions, typename Real, tnlDevice Device, typename Index >
bool tnlGrid< Dimensions, Real, Device, Index > :: load( tnlFile& file )
{
   if( ! tnlArray< Dimensions, Real, Device, Index > :: load( file ) )
   {
      cerr << "I was not able to read the tnlArray of the tnlGrid "
           << this -> getName() << endl;
      return false;
   }
   if( ! domainLowerCorner. load( file ) ||
       ! domainUpperCorner. load( file ) ||
       ! spaceSteps. load( file ) )
   {
      cerr << "I was not able to read the domain description of the tnlGrid "
           << this -> getName() << endl;
      return false;
   }
   return true;
}

template< int Dimensions, typename Real, tnlDevice Device, typename Index >
bool tnlGrid< Dimensions, Real, Device, Index > :: save( const tnlString& fileName ) const
{
   return tnlObject :: save( fileName );
}

template< int Dimensions, typename Real, tnlDevice Device, typename Index >
bool tnlGrid< Dimensions, Real, Device, Index > :: load( const tnlString& fileName )
{
   return tnlObject :: load( fileName );
}

template< int Dimensions, typename Real, tnlDevice Device, typename Index >
bool tnlGrid< Dimensions, Real, Device, Index > :: draw( const tnlString& fileName,
                                                         const tnlString& format,
                                                         const tnlTuple< Dimensions, Index > steps ) const
{
   tnlAssert( steps > ( tnlTuple< Dimensions, Index >( 0 ) ),
              cerr << "Wrong steps of increment ( " << steps << " )"
                   << " for drawing the tnlGrid " << this -> getName() << "." << endl; );
   if( format == "tnl" )
      return this -> save( fileName );
   fstream file;
   file. open( fileName. getString(), ios :: out );
   if( ! file )
   {
      cerr << " I am not able to open the file " << fileName
           << " for drawing the tnlGrid " << this -> getName() << "." << endl;
      return false;
   }
   if( Dimensions == 1 )
   {
      if( format == "gnuplot" )
      {
         const Index xSize = this -> getDimensions()[ tnlX ];
         const Real& ax = this -> getDomainLowerCorner()[ tnlX ];
         const Real& hx = this -> getSpaceSteps()[ tnlX ];
         for( Index i = 0; i < xSize; i += steps[ tnlX ] )
                  file << setprecision( 12 )
                       << ax + i * hx * steps[ tnlX ]
                       << " "
                       << this -> getElement( i )
                       << endl;
         return true;
      }
   }
   if( Dimensions == 2 )
   {
      const Index xSize = this -> getDimensions()[ tnlX ];
      const Index ySize = this -> getDimensions()[ tnlY ];
      const Real& ax = this -> getDomainLowerCorner()[ tnlX ];
      const Real& ay = this -> getDomainLowerCorner()[ tnlY ];
      const Real& hx = this -> getSpaceSteps()[ tnlX ];
      const Real& hy = this -> getSpaceSteps()[ tnlY ];
      if( format == "gnuplot" )
      {
         for( Index i = 0; i < xSize; i += steps[ tnlX ] )
         {
            for( Index j = 0; j < ySize; j += steps[ tnlY ] )
            {
               file << setprecision( 12 )
                    << ax + Real( i ) * hx * steps[ tnlX ]
                    << " "
                    << ay + Real( j ) * hy * steps[ tnlY ]
                    << " "
                    << this -> getElement( i, j )
                    << endl;
            }
            file << endl;
         }
         return true;
      }
      if( format == "vti" )
      {
         file << "<VTKFile type=\"ImagegetString\" version=\"0.1\" byte_order=\"LittleEndian\">" << endl;
         file << "<ImagegetString WholeExtent=\""
              << 0 << " " << xSize - 1 << " " << 0 << " " << ySize - 1
              << " 0 0\" Origin=\"0 0 0\" Spacing=\""
              << hx * steps[ tnlX ] << " " << hy * steps[ tnlY ] << " 0\">" << endl;
         file << "<Piece Extent=\"0 " << xSize - 1 << " 0 " << ySize - 1 <<" 0 0\">" << endl;
         file << "<PointgetString Scalars=\"order_parameter\">" << endl;
         file << "<getStringArray Name=\"order_parameter\" type=\"Float32\" format=\"ascii\">" << endl;
         file. flags( ios_base::scientific );
         Index iStep = steps[ tnlX ];
         Index jStep = steps[ tnlY ];
         for( Index j = 0; j <= ySize - jStep; j += jStep )
            for( Index i = 0; i <= xSize - iStep; i += iStep )
                 file << this -> getElement( i, j ) << " ";
         file << endl;
         file << "</getStringArray>" << endl;
         file << "</PointgetString>" << endl;
         file << "</Piece>" << endl;
         file << "</ImagegetString>" << endl;
         file << "</VTKFile>" << endl;
         return true;
      }
   }
   if( Dimensions == 3 )
   {
      const Index xSize = this -> getDimensions()[ tnlX ];
      const Index ySize = this -> getDimensions()[ tnlY ];
      const Index zSize = this -> getDimensions()[ tnlZ ];
      const Real& ax = this -> getDomainLowerCorner()[ tnlX ];
      const Real& ay = this -> getDomainLowerCorner()[ tnlY ];
      const Real& az = this -> getDomainLowerCorner()[ tnlZ ];
      const Real& hx = this -> getSpaceSteps()[ tnlX ];
      const Real& hy = this -> getSpaceSteps()[ tnlY ];
      const Real& hz = this -> getSpaceSteps()[ tnlZ ];
      if( format == "gnuplot" )
      {
         cout << "GNUPLOT is not supported for tnlGrid3D." << endl;
         return false;
      }
      if( format == "vti" )
      {
         file << "<VTKFile type=\"ImagegetString\" version=\"0.1\" byte_order=\"LittleEndian\">" << endl;
         file << "<ImagegetString WholeExtent=\""
              << "0 " << xSize - 1
              << " 0 " << ySize - 1
              << " 0 " << zSize - 1 << "\" Origin=\"0 0 0\" Spacing=\""
              << hx * steps[ tnlX ] << " " << hy * steps[ tnlY ] << " " << hz * steps[ tnlZ ] << "\">" << endl;
         file << "<Piece Extent=\"0 "
              << xSize - 1 << " 0 "
              << ySize - 1 << " 0 "
              << zSize - 1 << "\">" << endl;
         file << "<PointgetString Scalars=\"order_parameter\">" << endl;
         file << "<getStringArray Name=\"order_parameter\" type=\"Float32\" format=\"ascii\">" << endl;
         file. flags( ios_base::scientific );
         Index iStep = steps[ tnlX ];
         Index jStep = steps[ tnlY ];
         Index kStep = steps[ tnlZ ];
         for( Index k = 0; k <= zSize - kStep; k += kStep )
            for( Index j = 0; j <= ySize - jStep; j += jStep )
               for( Index i = 0; i <= xSize - iStep; i += iStep )
               {
                 file << this -> getElement( i, j, k ) << " ";
               }
         file << endl;
         file << "</getStringArray>" << endl;
         file << "</PointgetString>" << endl;
         file << "</Piece>" << endl;
         file << "</ImagegetString>" << endl;
         file << "</VTKFile>" << endl;
         return true;
      }
      if( format == "povray" )
      {
         file. put( ( char ) ( xSize >> 8 ) );
         file. put( ( char ) ( xSize & 0xff ) );
         file. put( ( char ) ( ySize >> 8 ) );
         file. put( ( char ) ( ySize & 0xff ) );
         file. put( ( char ) ( zSize >> 8 ) );
         file. put( ( char ) ( zSize & 0xff ) );
         Real min( this -> getElement( 0, 0, 0 ) ), max( this -> getElement( 0, 0, 0 ) );
         for( Index k = 0; k < zSize; k ++ )
            for( Index j = 0; j < ySize; j ++ )
               for( Index i = 0; i < xSize; i ++ )
               {
                  min = Min( min, this -> getElement( i, j, k ) );
                  max = Max( max, this -> getElement( i, j, k ) );
               }

         for( Index k = 0; k < zSize; k ++ )
            for( Index j = 0; j < ySize; j ++ )
               for( Index i = 0; i < xSize; i ++ )
               {
                  int v = Real( 255.0 ) * ( this -> getElement( i, j, k ) - min ) / ( max - min );
                  file. write( ( char* ) &v, sizeof( int ) );
               }
         return true;
      }

   }
   cerr << endl << "I do not know a format " << format << " for tnlGrid with " << Dimensions << " dimensions.";
   return false;
}

#ifdef HAVE_CUDA
template< int Dimensions, typename Real, typename Index >
__global__ void  setConstantDirichletBC( const Index xSize,
                                         const Index ySize,
                                         const Real bc,
                                         Real* u )
{
   const Index ij = blockIdx. x * blockDim. x + threadIdx. x;
   const Index i = ij / ySize;
   const Index j = ij % ySize;

   if( ij < xSize * ySize &&
       ( i == 0 || j == 0 || i == xSize - 1 || j == ySize - 1 ) )
   {
      u[ ij ] = bc;
   }

}

template< int Dimensions, typename Real, typename Index >
__global__ void setDirichletBC( const Index xSize,
                                const Index ySize,
                                const Real* bc,
                                Real* u )
{
   const Index ij = blockIdx. x * blockDim. x + threadIdx. x;
   const Index i = ij / ySize;
   const Index j = ij % ySize;

   if( ij < xSize * ySize &&
       ( i == 0 || j == 0 || i == xSize - 1 || j == ySize - 1 ) )
   {
      u[ ij ] = bc[ ij ];
   }
}

template< int Dimensions, typename Real, typename Index >
__global__ void setConstantNeumannBC( const Index xSize,
                                      const Index ySize,
                                      const Real hx,
                                      const Real hy,
                                      const Real bc,
                                      Real* u )
{
   const Index ij = blockIdx. x * blockDim. x + threadIdx. x;
   const Index i = ij / ySize;
   const Index j = ij % ySize;

   if( ij < xSize * ySize &&
       ( i == 0 || j == 0 || i == xSize - 1 || j == ySize - 1 ) )
   {
     if( i == 0 )
        u[ ij ] = u[ ij + ySize ] + hx * bc;
     if( i == xSize - 1 )
        u[ ij ] = u[ ij - ySize ] + hx * bc;

     __syncthreads();

     if( j == 0 )
        u[ ij ] = u[ ij + 1 ] + hy * bc;
     if( j == ySize - 1 )
        u[ ij ] = u[ ij - 1 ] + hy * bc;
   }
}

template< int Dimensions, typename Real, typename Index >
__global__ void setNeumannBC( const Index xSize,
                              const Index ySize,
                              const Real hx,
                              const Real hy,
                              const Real* bc,
                              Real* u )
{
   const Index ij = blockIdx. x * blockDim. x + threadIdx. x;
   const Index i = ij / ySize;
   const Index j = ij % ySize;

   if( ij < xSize * ySize &&
       ( i == 0 || j == 0 || i == xSize - 1 || j == ySize - 1 ) )
   {
     if( i == 0 )
        u[ ij ] = u[ ij + ySize ] + hx * bc[ ij ];
     if( i == xSize - 1 )
        u[ ij ] = u[ ij - ySize ] + hx * bc[ ij ];

     __syncthreads();

     if( j == 0 )
        u[ ij ] = u[ ij + 1 ] + hy * bc[ ij ];
     if( j == ySize - 1 )
        u[ ij ] = u[ ij - 1 ] + hy * bc[ ij ];
   }
}

#endif


#endif /* TNLGRID_H_ */
