RegularVector3dField3d.hpp

Go to the documentation of this file.

/*
 * Copyright(c):
 * Signal Image and Communications (SIC) Department 
 * http://www.sic.sp2mi.univ-poitiers.fr/ 
 *  - University of Poitiers, France http://www.univ-poitiers.fr 
 *  - XLIM  Institute UMR CNRS 7252  http://www.xlim.fr/ 
 * 
 * and 
 * 
 * D2 Fluid, Thermic and Combustion 
 *  - University of Poitiers, France http://www.univ-poitiers.fr 
 *  - PPRIME Institute -  UPR CNRS 3346  http://www.pprime.fr
 *  - ISAE-ENSMA http://www.ensma.fr
 *
 * Contributor(s):
 * The SLIP team,
 * Benoit Tremblais <tremblais_AT_sic.univ-poitiers.fr>,
 * Laurent David <laurent.david_AT_lea.univ-poitiers.fr>, 
 * Ludovic Chatellier <ludovic.chatellier_AT_univ-poitiers.fr>, 
 * Lionel Thomas <lionel.thomas_AT_univ-poitiers.fr>, 
 * Denis Arrivault <arrivault_AT_sic.univ-poitiers.fr>, 
 * Julien Dombre <julien.dombre_AT_univ-poitiers.fr>.
 *
 * Description:
 * The Simple Library of Image Processing (SLIP) is a new image processing 
 * library. It is written in the C++ language following as much as possible 
 * the ISO/ANSI C++ standard. It is consequently compatible with any system  
 * satisfying the ANSI C++ complience. It works on different Unix , Linux ,  
 * Mircrosoft Windows and Mac OS X plateforms. SLIP is a research library that  
 * was created by the Signal, Image and Communications (SIC) departement of  
 * the XLIM, UMR 7252 CNRS Institute in collaboration with the Fluids, Thermic  
 * and Combustion departement of the P', UPR 3346 CNRS Institute of the  
 * University of Poitiers.
 * 
 * The SLIP Library source code has been registered to the APP (French Agency 
 * for the Protection of Programs) by the University of Poitiers and CNRS, 
 * under  registration number IDDN.FR.001.300034.000.S.P.2010.000.21000.

 * http://www.sic.sp2mi.univ-poitiers.fr/slip/
 *
 * This software is governed by the CeCILL-C license under French law and
 * abiding by the rules of distribution of free software.  You can  use, 
 * modify and/ or redistribute the software under the terms of the CeCILL-C 
 * license as circulated by CEA, CNRS and INRIA at the following URL 
 * http://www.cecill.info.
 * As a counterpart to the access to the source code and  rights to copy,
 * modify and redistribute granted by the license, users are provided only
 * with a limited warranty  and the software's author,  the holder of the
 * economic rights,  and the successive licensors  have only  limited
 * liability.
 * 
 * In this respect, the user's attention is drawn to the risks associated
 * with loading,  using,  modifying and/or developing or reproducing the
 * software by the user in light of its specific status of free software,
 * that may mean  that it is complicated to manipulate,  and  that  also
 * therefore means  that it is reserved for developers  and  experienced
 * professionals having in-depth computer knowledge. Users are therefore
 * encouraged to load and test the software's suitability as regards their
 * requirements in conditions enabling the security of their systems and/or
 * data to be ensured and,  more generally, to use and operate it in the 
 * same conditions as regards security. 
 *
 * The fact that you are presently reading this means that you have had
 * knowledge of the CeCILL-C license and that you accept its terms.
 */



#ifndef SLIP_REGULARVECTOR3DFIELD3D_HPP
#define SLIP_REGULARVECTOR3DFIELD3D_HPP

#include <iostream>
#include <fstream>
#include <iterator>
#include <cassert>
#include <numeric>
#include <cmath>
#include <string>
#include <cstddef>
#include "Matrix3d.hpp"
#include "stride_iterator.hpp"
#include "kstride_iterator.hpp"
#include "iterator3d_box.hpp"
#include "iterator3d_range.hpp"
#include "apply.hpp"
#include "Vector3d.hpp"
#include "GenericMultiComponent3d.hpp"
#include "derivatives.hpp"
#include "arithmetic_op.hpp"
#include "io_tools.hpp"
#include "linear_algebra_eigen.hpp"
#include "linear_algebra_traits.hpp"
#include "tecplot_binaries.hpp"

#include <boost/serialization/access.hpp>
#include <boost/serialization/split_member.hpp>
#include <boost/serialization/version.hpp>
#include <boost/serialization/base_object.hpp>

namespace slip
{

template<class T>
class stride_iterator;

template<class T>
class iterator3d_box;

template<typename T>
class iterator3d_range;

template<typename T>
class const_iterator3d_box;

template<typename T>
class const_iterator3d_range;

template <class T>
class DPoint3d;

template <class T>
class Point3d;

template <class T>
class Box3d;

  template <typename T, typename GridT>
class RegularVector3dField3d;

  template <typename T, typename GridT>
  std::ostream& operator<<(std::ostream & out, const slip::RegularVector3dField3d<T,GridT>& a);

  template <typename T, typename GridT = double>
class RegularVector3dField3d:public slip::GenericMultiComponent3d<slip::Vector3d<T> >
{
public :

  typedef slip::Vector3d<T> value_type;
  typedef RegularVector3dField3d<T,GridT> self;
  typedef const RegularVector3dField3d<T,GridT> const_self;
  typedef GenericMultiComponent3d<slip::Vector3d<T> > base;


  typedef value_type* pointer;
  typedef const value_type* const_pointer;
  typedef value_type& reference;
  typedef const value_type& const_reference;

  typedef ptrdiff_t difference_type;
  typedef std::size_t size_type;

  typedef pointer iterator;
  typedef const_pointer const_iterator;
  
  typedef std::reverse_iterator<iterator> reverse_iterator;
  typedef std::reverse_iterator<const_iterator> const_reverse_iterator;


  typedef pointer row_iterator;
  typedef const_pointer const_row_iterator;
  typedef slip::stride_iterator<pointer> col_iterator;
  typedef slip::stride_iterator<const_pointer> const_col_iterator;

  typedef std::reverse_iterator<iterator> reverse_row_iterator;
  typedef std::reverse_iterator<const_iterator> const_reverse_row_iterator;
  typedef std::reverse_iterator<col_iterator> reverse_col_iterator;
  typedef std::reverse_iterator<const_col_iterator> const_reverse_col_iterator;

  typedef typename slip::GenericMultiComponent3d<slip::Vector3d<T> >::iterator3d iterator3d;
 typedef typename slip::GenericMultiComponent3d<slip::Vector3d<T> >::const_iterator2d const_iterator3d;

  typedef T vector3d_value_type;
  typedef GridT grid_value_type;
  typedef vector3d_value_type* vector3d_pointer;
  typedef const vector3d_value_type* const_vector3d_pointer;
  typedef vector3d_value_type& vector3d_reference;
  typedef const vector3d_value_type const_vector3d_reference;
  typedef slip::kstride_iterator<vector3d_pointer,3> vector3d_iterator;
  typedef slip::kstride_iterator<const_vector3d_pointer,3> const_vector3d_iterator;

  typedef typename slip::lin_alg_traits<value_type>::value_type norm_type;
  static const std::size_t DIM = 3;
public:

  RegularVector3dField3d():
    base(),init_point_(slip::Point3d<GridT>()),grid_step_(slip::Point3d<GridT>(GridT(1),GridT(1),GridT(1)))
  {}
 
  RegularVector3dField3d(const size_type slices,
                                       const size_type rows,
                                       const size_type cols):
    base(slices,rows,cols),init_point_(slip::Point3d<GridT>()),grid_step_(slip::Point3d<GridT>(GridT(1),GridT(1),GridT(1)))
  {}
  
  RegularVector3dField3d(const size_type slices,
                                                 const size_type rows,
                                                 const size_type cols,
                                                 const slip::Point3d<GridT>& init_point,
                                                 const slip::Point3d<GridT>& grid_step):
    base(slices,rows,cols),init_point_(init_point),grid_step_(grid_step)
  {}
  
  
  RegularVector3dField3d(const size_type slices,
                                       const size_type rows,
                                       const size_type cols,
                                       const slip::Vector3d<T>& val):
    base(slices,rows,cols,val),init_point_(slip::Point3d<GridT>()),grid_step_(slip::Point3d<GridT>(GridT(1),GridT(1),GridT(1)))
  {}

  RegularVector3dField3d(const size_type slices,
                                                 const size_type rows,
                                                 const size_type cols,
                                                 const slip::Point3d<GridT>& init_point,
                                                 const slip::Point3d<GridT>& grid_step,
                                                 const slip::Vector3d<T>& val):
    base(slices,rows,cols,val),init_point_(init_point),grid_step_(grid_step)
  {}
 
  RegularVector3dField3d(const size_type slices,
                                       const size_type rows,
                                       const size_type cols,
                                       const T* val):
    base(slices,rows,cols,val),init_point_(slip::Point3d<GridT>()),grid_step_(slip::Point3d<GridT>(GridT(1),GridT(1),GridT(1)))
  {}

  RegularVector3dField3d(const size_type slices,
                                                 const size_type rows,
                                                 const size_type cols,
                                                 const slip::Point3d<GridT>& init_point,
                                                 const slip::Point3d<GridT>& grid_step,
                                                 const T* val):
    base(slices,rows,cols,val),init_point_(init_point),grid_step_(grid_step)
  {}
 
  RegularVector3dField3d(const size_type slices,
                                       const size_type rows,
                                       const size_type cols,
                                       const slip::Vector3d<T>* val):
    base(slices,rows,cols,val),init_point_(slip::Point3d<GridT>()),grid_step_(slip::Point3d<GridT>(GridT(1),GridT(1),GridT(1)))
  {}
 

  RegularVector3dField3d(const size_type slices,
                                       const size_type rows,
                                       const size_type cols,
                                                 const slip::Point3d<GridT>& init_point,
                                                 const slip::Point3d<GridT>& grid_step,
                                       const slip::Vector3d<T>* val):
    base(slices,rows,cols,val),init_point_(init_point),grid_step_(grid_step)
  {}
 

  template<typename InputIterator>
  RegularVector3dField3d(const size_type slices,
                                       const size_type rows,
                                       const size_type cols,
                                       InputIterator first,
                                       InputIterator last):
    base(slices,rows,cols,first,last),init_point_(slip::Point3d<GridT>()),grid_step_(slip::Point3d<GridT>(GridT(1),GridT(1),GridT(1)))
  {} 

  template<typename InputIterator>
  RegularVector3dField3d(const size_type slices,
                                       const size_type rows,
                                       const size_type cols,
                                                 const slip::Point3d<GridT>& init_point,
                                                 const slip::Point3d<GridT>& grid_step,
                                       InputIterator first,
                                       InputIterator last):
    base(slices,rows,cols,first,last),init_point_(init_point),grid_step_(grid_step)
  {}

  
  template<typename InputIterator>
  RegularVector3dField3d(const size_type slices,
                                       const size_type rows,
                                       const size_type cols,
                                       InputIterator first1,
                                       InputIterator last1,
                                       InputIterator first2,
                                       InputIterator first3):
    base(slices,rows,cols),init_point_(slip::Point3d<GridT>()),grid_step_(slip::Point3d<GridT>(GridT(1),GridT(1),GridT(1)))
  {
    
    std::vector<InputIterator> first_iterators_list(3);
    first_iterators_list[0] = first1;
    first_iterators_list[1] = first2;
    first_iterators_list[2] = first3;
    this->fill(first_iterators_list,last1);
  
  } 

  template<typename InputIterator>
  RegularVector3dField3d(const size_type slices,
                                                 const size_type rows,
                                                 const size_type cols,
                                                 const slip::Point3d<GridT>& init_point,
                                                 const slip::Point3d<GridT>& grid_step,
                                                 InputIterator first1,
                                                 InputIterator last1,
                                                 InputIterator first2,
                                                 InputIterator first3):
    base(slices,rows,cols),init_point_(init_point),grid_step_(grid_step)
  {
    
    std::vector<InputIterator> first_iterators_list(3);
    first_iterators_list[0] = first1;
    first_iterators_list[1] = first2;
    first_iterators_list[2] = first3;
    this->fill(first_iterators_list,last1);
  
  } 

  RegularVector3dField3d(const self& rhs):
    base(rhs),init_point_(rhs.init_point_),grid_step_(rhs.grid_step_)
  {}

  ~RegularVector3dField3d()
  {}
  

  std::string name() const;



   
  friend std::ostream& operator<< <>(std::ostream & out, const self& a);
 
 

  void write_tecplot(const std::string& file_path_name,
                                     const std::string& title,
                                     const std::string& zone) const;


  void write_tecplot_bin(const std::string& file_path_name,
                                                 const std::string& title,
                                                 const std::string& zone) const;


 

  void read_tecplot(const std::string& file_path_name);

  void read_tecplot_bin(const std::string& file_path_name,
                                                const int zone_loaded=1);



 
  self& operator=(const slip::Vector3d<T>& val);


  self& operator=(const T& val);



  
 
  T& Vx1(const size_type k,
                 const size_type i,
                 const size_type j);

  const T& Vx1(const size_type k,
                       const size_type i,
                       const size_type j) const;

  T& u(const size_type k,
       const size_type i,
       const size_type j);

  const T& u(const size_type k,
                     const size_type i,
                     const size_type j) const;

 

 
  T& Vx2(const size_type k,
                 const size_type i,
                 const size_type j);

  const T& Vx2(const size_type k,
                       const size_type i,
                       const size_type j) const;

  T& v(const size_type k,
       const size_type i,
       const size_type j);

  const T& v(const size_type k,
                     const size_type i,
                     const size_type j) const;

 
  T& Vx3(const size_type k,
                 const size_type i,
                 const size_type j);

  const T& Vx3(const size_type k,
                       const size_type i,
                       const size_type j) const;

  T& w(const size_type k,
                 const size_type i,
                 const size_type j);

  const T& w(const size_type k,
                     const size_type i,
                     const size_type j) const;

  norm_type norm(const size_type k,
                                 const size_type i,
                                 const size_type j) const;
 

  const GridT x1(const size_type k,
                     const size_type i,
                     const size_type j) const;

  const GridT x2(const size_type k,
                     const size_type i,
                     const size_type j) const;

  const GridT x3(const size_type k,
                     const size_type i,
                     const size_type j) const;

  const GridT x(const size_type k,
                     const size_type i,
                     const size_type j) const;

  const GridT y(const size_type k,
                     const size_type i,
                     const size_type j) const;

  const GridT z(const size_type k,
                     const size_type i,
                     const size_type j) const;

  void set_init_point(const slip::Point3d<GridT>& init_point);

  const slip::Point3d<GridT>& get_init_point() const;
 
  void set_grid_step(const slip::Point3d<GridT>& grid_step);

  const slip::Point3d<GridT>& get_grid_step() const;
  self& operator+=(const T& val);
  self& operator-=(const T& val);
  self& operator*=(const T& val);
  self& operator/=(const T& val);



   self  operator-() const;

  self& operator+=(const slip::Vector3d<T>& val);
  self& operator-=(const slip::Vector3d<T>& val);
  self& operator*=(const slip::Vector3d<T>& val);
  self& operator/=(const slip::Vector3d<T>& val);

  

  self& operator+=(const self& rhs);
  self& operator-=(const self& rhs);
  self& operator*=(const self& rhs);
  self& operator/=(const self& rhs);
    


  template<typename Container3D>
  void derivative(const std::size_t component,
                                  const slip::SPATIAL_DIRECTION der_dir,
                                  const std::size_t der_order, 
                                  const std::size_t sch_order,
                                  Container3D& result) const
  {
    assert(component < 3);
    assert(sch_order >= der_order);
    assert(   ((der_dir == X_DIRECTION) && (this->rows() > sch_order))
                    ||((der_dir == Y_DIRECTION) && (this->cols() > sch_order)) 
                    ||((der_dir == Z_DIRECTION) && (this->cols() > sch_order)));
    assert(sch_order / 2 <= sch_order);
    const std::size_t sch_shift = sch_order / 2;
    //computes all kernels
    slip::Matrix<double> kernels(sch_order + 1, sch_order + 1);
    std::vector<slip::Matrix<double>::iterator> kernels_iterators(sch_order + 1);   
    for(std::size_t i = 0; i < (sch_order + 1); ++i)
      {
                kernels_iterators[i] = kernels.row_begin(i);
      }
   
    if(der_dir == X_DIRECTION)
      {
                slip::finite_diff_kernels(der_order,sch_order,sch_shift,
                                                                  grid_step_[0],
                                                                  kernels_iterators);
                const size_type rows = this->rows();
                const size_type slices = this->slices();
                for(size_type k = 0; k < slices; ++k)
                  {
                    for(size_type i = 0; i < rows; ++i)
                      {
                                slip::derivative(this->row_begin(component,k,i),
                                                                 this->row_end(component,k,i), 
                                                                 der_order,
                                                                 sch_order,
                                                                 sch_shift,
                                                                 kernels_iterators,
                                                                 result.row_begin(k,i));
                      }
                  }
      }
    if(der_dir == Y_DIRECTION)
      {
                slip::finite_diff_kernels(der_order,sch_order,sch_shift,
                                                                  grid_step_[1],
                                                                  kernels_iterators);
                const size_type cols = this->cols();
                const size_type slices = this->slices();
                for(size_type k = 0; k < slices; ++k)
                  {
                    for(size_type j = 0; j < cols; ++j)
                      {
                                slip::derivative(this->col_rbegin(component,k,j),
                                                                 this->col_rend(component,k,j), 
                                                                 der_order,
                                                                 sch_order,
                                                                 sch_shift,
                                                                 kernels_iterators,
                                                                 result.col_rbegin(k,j));
                      }
                  }
      }
    if(der_dir == Z_DIRECTION)
      {
                slip::finite_diff_kernels(der_order,sch_order,sch_shift,
                                                                  grid_step_[2],
                                                                  kernels_iterators);
                const size_type rows = this->rows();
                const size_type cols = this->cols();
                for(size_type i = 0; i < rows; ++i)
                  {
                    for(size_type j = 0; j < cols; ++j)
                      {
                                slip::derivative(this->slice_rbegin(component,i,j),
                                                                 this->slice_rend(component,i,j), 
                                                                 der_order,
                                                                 sch_order,
                                                                 sch_shift,
                                                                 kernels_iterators,
                                                                 result.slice_rbegin(i,j));
                      }
                  }
      }
  }


  template<typename Container3D>
  void derivative(const std::size_t component,
      const slip::SPATIAL_DIRECTION der_dir,
      const std::size_t der_order, 
      const std::size_t sch_order,
      const slip::Box3d<int> box,
      Container3D& result) const
  {
    assert(component < 3);
    assert(sch_order >= der_order);
    assert(   ((der_dir == X_DIRECTION) && (this->rows() > sch_order))
                      ||((der_dir == Y_DIRECTION) && (this->cols() > sch_order)) 
                      ||((der_dir == Z_DIRECTION) && (this->cols() > sch_order)));
    assert(sch_order / 2 <= sch_order);
    const std::size_t sch_shift = sch_order / 2;
    //computes all kernels
    slip::Matrix<double> kernels(sch_order + 1, sch_order + 1);
    std::vector<slip::Matrix<double>::iterator> kernels_iterators(sch_order + 1);   
    for(std::size_t i = 0; i < (sch_order + 1); ++i)
      {
                kernels_iterators[i] = kernels.row_begin(i);
      }
   
    const size_type slices=box.depth();
    const size_type rows=box.height();
    const size_type cols=box.width();

    if(der_dir == X_DIRECTION)
      {
                slip::finite_diff_kernels(der_order,sch_order,sch_shift,
                                                                  grid_step_[0],
                                                                  kernels_iterators);
                const slip::Range<int> range_row(box.front_upper_left()[2],box.back_bottom_right()[2],1);
                for(int k = box.front_upper_left()[0]; k <= box.back_bottom_right()[0]; ++k)
                  {
                    for(int i = box.front_upper_left()[1]; i <= box.back_bottom_right()[1]; ++i)
                      {
                                slip::derivative(this->row_begin(component, k,i,range_row),
                                                                 this->row_end(component, k,i, range_row), 
                                                                 der_order,
                                                                 sch_order,
                                                                 sch_shift,
                                                                 kernels_iterators,
                                                                 result.row_begin(k,i, range_row));
                      }
                  }
      }
    if(der_dir == Y_DIRECTION)
      {
                slip::finite_diff_kernels(der_order,sch_order,sch_shift,
                                                                  grid_step_[1],
                                                                  kernels_iterators);
                const slip::Range<int> range_col(box.front_upper_left()[1],box.back_bottom_right()[1],1);
                for(int k = box.front_upper_left()[0]; k <= box.back_bottom_right()[0]; ++k)
                  {
                    for(int j = box.front_upper_left()[2]; j <= box.back_bottom_right()[2]; ++j)
                      {
                                slip::derivative(this->col_rbegin(component, k,j,range_col),
                                                                 this->col_rend(component, k,j, range_col), 
                                                                 der_order,
                                                                 sch_order,
                                                                 sch_shift,
                                                                 kernels_iterators,
                                                                 result.col_rbegin(k,j, range_col));
                      }
                  }
      }
    if(der_dir == Z_DIRECTION)
      {
                slip::finite_diff_kernels(der_order,sch_order,sch_shift,
                                                                  grid_step_[2],
                                                                  kernels_iterators);
                const slip::Range<int> range_slice(box.front_upper_left()[0],box.back_bottom_right()[0],1);
                for(int i = box.front_upper_left()[1]; i < box.back_bottom_right()[1]; ++i)
                  {
                    for(int j = box.front_upper_left()[2]; j < box.back_bottom_right()[2]; ++j)
                      {
                                slip::derivative(this->slice_rbegin(component, i,j,range_slice),
                                                                 this->slice_rend(component, i,j, range_slice),
                                                                 der_order,
                                                                 sch_order,
                                                                 sch_shift,
                                                                 kernels_iterators,
                                                                 result.slice_rbegin(i,j, range_slice));
                      }
                  }
      }
  }

  template<typename Container3D>
  void divergence(const std::size_t sch_order,
                                  Container3D& result) const
  {
    assert(sch_order >= 1);
    assert(sch_order / 2 <= sch_order);
    const std::size_t der_order = 1;
  
    //x-derivative of the x component (first component)
    slip::Matrix3d<double> Mtmp(this->slices(),this->rows(),this->cols()); 
    std::size_t component = 0;
    this->derivative(component,slip::X_DIRECTION,der_order,sch_order,Mtmp);
 
    //y-derivative of the y component (second component) 
    component = 1;
    this->derivative(component,slip::Y_DIRECTION,der_order,sch_order,result);
 
    //pre-compute the divergence
    slip::plus(Mtmp.begin(),Mtmp.end(),result.begin(),result.begin());
    //z-derivative of the 2 component (third component) 
    component = 2;
    this->derivative(component,slip::Z_DIRECTION,der_order,sch_order,Mtmp);
     //compute the divergence
    slip::plus(result.begin(),result.end(),Mtmp.begin(),result.begin());
                
    
  }
 

  template<typename Container3D>
  void vorticity(const std::size_t sch_order,
                                  Container3D& result) const
  {
    assert(sch_order >= 1);
    assert(sch_order / 2 <= sch_order);
    const std::size_t sch_shift = sch_order / 2;
    const size_type slices = this->slices();
    const size_type rows = this->rows();
    const size_type cols = this->cols();
    //computes all kernels
    const std::size_t der_order = 1;
    //x kernels
    slip::Matrix<double> kernelsx(sch_order + 1, sch_order + 1);
    std::vector<slip::Matrix<double>::iterator> kernelsx_iterators(sch_order + 1);   
    for(std::size_t i = 0; i < (sch_order + 1); ++i)
      {
                kernelsx_iterators[i] = kernelsx.row_begin(i);
      }
    slip::finite_diff_kernels(der_order,sch_order,sch_shift,
                                                      grid_step_[0],
                                                      kernelsx_iterators);
     //y kernels
     slip::Matrix<double> kernelsy(sch_order + 1, sch_order + 1);
    std::vector<slip::Matrix<double>::iterator> kernelsy_iterators(sch_order + 1);   
    for(std::size_t i = 0; i < (sch_order + 1); ++i)
      {
                kernelsy_iterators[i] = kernelsy.row_begin(i);
      }
    slip::finite_diff_kernels(der_order,sch_order,sch_shift,
                                                      grid_step_[1],
                                                      kernelsy_iterators);
     //z kernels
 slip::Matrix<double> kernelsz(sch_order + 1, sch_order + 1);
    std::vector<slip::Matrix<double>::iterator> kernelsz_iterators(sch_order + 1);   
    for(std::size_t i = 0; i < (sch_order + 1); ++i)
      {
                kernelsz_iterators[i] = kernelsz.row_begin(i);
      }
    slip::finite_diff_kernels(der_order,sch_order,sch_shift,
                                                      grid_step_[2],
                                                      kernelsz_iterators);
    // compute z component
    std::size_t vort_component = 2;
    //x-derivative of the y component (second component)
    slip::Matrix3d<double> Mtmp(this->slices(),this->rows(),this->cols()); 
    std::size_t component = 1;

    for(size_type k = 0; k < slices; ++k)
      {
                for(size_type i = 0; i < rows; ++i)
                  {
                    slip::derivative(this->row_begin(component,k,i),
                                                     this->row_end(component,k,i), 
                                                     der_order,
                                                     sch_order,
                                                     sch_shift,
                                                     kernelsx_iterators,
                                                     Mtmp.row_begin(k,i));
                  }
      }
  
    
    //y-derivative of the x component (first component) 
    component = 0;
    for(size_type k = 0; k < slices; ++k)
      {
                for(size_type j = 0; j < cols; ++j)
                  {
                    slip::derivative(this->col_rbegin(component,k,j),
                                                     this->col_rend(component,k,j), 
                                                     der_order,
                                                     sch_order,
                                                     sch_shift,
                                                     kernelsy_iterators,
                                                     result.col_rbegin(vort_component,k,j));
                  }
      } 
    // z component = dVy/dx-dVx/dy
    slip::minus(Mtmp.begin(),Mtmp.end(),result.begin(vort_component),result.begin(vort_component));

    // compute y component
    vort_component = 1;
    //z-derivative of the x component (first component) 
    component = 0;
    for(size_type i = 0; i < rows; ++i)
      {
                for(size_type j = 0; j < cols; ++j)
                  {
                    slip::derivative(this->slice_rbegin(component,i,j),
                                                     this->slice_rend(component,i,j), 
                                                     der_order,
                                                     sch_order,
                                                     sch_shift,
                                                     kernelsz_iterators,
                                                     Mtmp.slice_rbegin(i,j));
                  }
      }
    //x-derivative of the z component (third component)
    component = 2;
    for(size_type k = 0; k < slices; ++k)
      {
                for(size_type i = 0; i < rows; ++i)
                  {
                    slip::derivative(this->row_begin(component,k,i),
                                                     this->row_end(component,k,i), 
                                                     der_order,
                                                     sch_order,
                                                     sch_shift,
                                                     kernelsx_iterators,
                                                     result.row_begin(vort_component,k,i));
                  }
      }
    // y component = dVx/dz - dVz/dx
    slip::minus(Mtmp.begin(),Mtmp.end(),result.begin(vort_component),result.begin(vort_component));

    // compute x component
    vort_component = 0;
    //z-derivative of the y component (third component) 
    component = 1;
    for(size_type i = 0; i < rows; ++i)
      {
                for(size_type j = 0; j < cols; ++j)
                  {
                    slip::derivative(this->slice_rbegin(component,i,j),
                                                     this->slice_rend(component,i,j), 
                                                     der_order,
                                                     sch_order,
                                                     sch_shift,
                                                     kernelsz_iterators,
                                                     Mtmp.slice_rbegin(i,j));
                  }
      }
    //y-derivative of the z component (third component) 
    component = 2;
    for(size_type k = 0; k < slices; ++k)
      {
                for(size_type j = 0; j < cols; ++j)
                  {
                    slip::derivative(this->col_rbegin(component,k,j),
                                                     this->col_rend(component,k,j), 
                                                     der_order,
                                                     sch_order,
                                                     sch_shift,
                                                     kernelsy_iterators,
                                                     result.col_rbegin(vort_component,k,j));
                  }
      } 
    // x component = dVz/dy - dVy/dz
    slip::minus(result.begin(vort_component),result.end(vort_component),Mtmp.begin(),result.begin(vort_component));
  }



  template<typename Container3D, typename Container3D1, typename Container3D2, typename Container3D3>
  void lambda(const std::size_t sch_order,
                      Container3D& Lambda, Container3D1& LambdaV1, Container3D2& LambdaV2, Container3D3& LambdaV3) const
  {
    //D=d(this)/d(x,y,z)
    //O=(D-Dt)/2
    //S=(D+Dt)/2
    //lambda2=eig(S^2+O^2,2)
    assert(sch_order >= 1);
    const std::size_t der_order = 1;
    const size_type slices = this->slices();
    const size_type rows = this->rows();
    const size_type cols = this->cols();
 
    // gradient tensor
    //derivatives of x component
    std::size_t component = 0;
    //x-derivative
    slip::Matrix3d<double> D11(slices,rows,cols); 
    this->derivative(component,slip::X_DIRECTION,der_order,sch_order,D11);
    //y-derivative
    slip::Matrix3d<double> D21(slices,rows,cols); 
    this->derivative(component,slip::Y_DIRECTION,der_order,sch_order,D21);
    //z-derivative
    slip::Matrix3d<double> D31(slices,rows,cols); 
    this->derivative(component,slip::Z_DIRECTION,der_order,sch_order,D31);
    //derivatives of y component
    component = 1;
    //x-derivative
    slip::Matrix3d<double> D12(slices,rows,cols); 
    this->derivative(component,slip::X_DIRECTION,der_order,sch_order,D12);
    //y-derivative
    slip::Matrix3d<double> D22(slices,rows,cols); 
    this->derivative(component,slip::Y_DIRECTION,der_order,sch_order,D22);
    //z-derivative
    slip::Matrix3d<double> D32(slices,rows,cols); 
    this->derivative(component,slip::Z_DIRECTION,der_order,sch_order,D32);
    //derivatives of z component
    component = 2;
    //x-derivative
    slip::Matrix3d<double> D13(slices,rows,cols); 
    this->derivative(component,slip::X_DIRECTION,der_order,sch_order,D13);
    //y-derivative
    slip::Matrix3d<double> D23(slices,rows,cols); 
    this->derivative(component,slip::Y_DIRECTION,der_order,sch_order,D23);
    //z-derivative
    slip::Matrix3d<double> D33(slices,rows,cols); 
    this->derivative(component,slip::Z_DIRECTION,der_order,sch_order,D33);

    //symetric tensor
    slip::Matrix<double> S(3,3);
    //antisymetric tensor
    slip::Matrix<double> O(3,3);
    // tensor to diagonalize S2O2=S^2+O^2
    slip::Matrix<double> S2(3,3),O2(3,3),S2O2(3,3); 
   //  // eigenvalues
//     slip::Vector<double> lambda(6); // 2*3 components <=> 3 complex eigenvalues
//     // eigenvectors
//     slip::Matrix<double> eigenvectors(6,3);// 3 vectors * 2*3 components
//     // sort eigenvalues ? no <- sorting is based on absolute value, not suitable with real eigens
//     bool sort=false;
//     // normalize eigenvectors ? yes
//     bool normalize=true;
//     // error flag for slip::eigen
//     std::size_t ier;

      // eigenvalues
    slip::Vector<double> lambda(3);
    // eigenvectors
    slip::Matrix<double> eigenvectors(3,3);
   

    Lambda.resize(slices,rows,cols);
    LambdaV1.resize(slices,rows,cols);
    LambdaV2.resize(slices,rows,cols);
    LambdaV3.resize(slices,rows,cols);

    //ça serait mieux avec un iterateur...
    for(size_type k = 0; k < slices; ++k)
      {
                for(size_type i = 0; i < rows; ++i) 
                  {
                    for(size_type j = 0; j < cols; ++j)
                      {
                                S(0,0)=D11[k][i][j];
                                S(1,1)=D22[k][i][j];
                                S(2,2)=D33[k][i][j];
                                S(0,1)=0.5*(D12[k][i][j]+D21[k][i][j]);
                                S(0,2)=0.5*(D13[k][i][j]+D31[k][i][j]);
                                S(1,0)=S(0,1);
                                S(1,2)=0.5*(D23[k][i][j]+D32[k][i][j]);
                                S(2,0)=S(0,2);
                                S(2,1)=S(1,2);
                                O(0,1)=0.5*(D12[k][i][j]-D21[k][i][j]);
                                O(0,2)=0.5*(D13[k][i][j]-D31[k][i][j]);
                                O(1,0)=-O(0,1);
                                O(1,2)=0.5*(D23[k][i][j]-D32[k][i][j]);
                                O(2,0)=-O(0,2);
                                O(2,1)=-O(1,2);
                                slip::matrix_matrix_multiplies(S,S,S2);
                                slip::matrix_matrix_multiplies(O,O,O2);
                                S2O2=0.0;
                                S2O2+=S2;
                                S2O2+=O2;
                //              //slip::eigen(S2O2, lambda, ier,sort);
//                              //slip::eigen(S2O2, lambda, eigenvectors, ier, sort,normalize);
//                              //std::cout<<k<<" "<<i<<" "<<j<<std::endl;
//                              eigen(S2O2, lambda, eigenvectors, ier, sort,normalize);
//                              // lambda2 is the 2nd highest eigenvalue of symetric tensor S^2+O^2 (should be real)
//                              // sorting is made upon absolute value of the a priori complex eigenvalues -> no sorting here
                                slip::hermitian_eigen(S2O2,lambda,eigenvectors);
                //              for(std::size_t kc=0; kc<2; ++kc)
//                                {
//                                  Lambda[k][i][j][kc]=lambda[2*kc];
                                    
//                                  LambdaV1[k][i][j][kc]=eigenvectors[2*kc][0];
//                                  LambdaV2[k][i][j][kc]=eigenvectors[2*kc][1];
//                                  LambdaV3[k][i][j][kc]=eigenvectors[2*kc][2];
                                    
//                                }
                                for(std::size_t kc=0; kc<2; ++kc)
                                  {
                                    Lambda[k][i][j][kc]=lambda[kc];
                                    
                                    LambdaV1[k][i][j][kc]=eigenvectors[kc][0];
                                    LambdaV2[k][i][j][kc]=eigenvectors[kc][1];
                                    LambdaV3[k][i][j][kc]=eigenvectors[kc][2];
                                    
                                  }
                      }
                  }
      }
  }



  template<typename Container3D>
  void lambda(const std::size_t sch_order,
                                  Container3D& Lambda) const
  {
    //D=d(this)/d(x,y,z)
    //O=(D-Dt)/2
    //S=(D+Dt)/2
    //lambda2=eig(S^2+O^2,2)
    assert(sch_order >= 1);
    const std::size_t der_order = 1;
    const size_type slices = this->slices();
    const size_type rows = this->rows();
    const size_type cols = this->cols();
 
    // gradient tensor
    //derivatives of x component
    std::size_t component = 0;
    //x-derivative
    slip::Matrix3d<double> D11(slices,rows,cols); 
    this->derivative(component,slip::X_DIRECTION,der_order,sch_order,D11);
    //y-derivative
    slip::Matrix3d<double> D21(slices,rows,cols); 
    this->derivative(component,slip::Y_DIRECTION,der_order,sch_order,D21);
    //z-derivative
    slip::Matrix3d<double> D31(slices,rows,cols); 
    this->derivative(component,slip::Z_DIRECTION,der_order,sch_order,D31);
    //derivatives of y component
    component = 1;
    //x-derivative
    slip::Matrix3d<double> D12(slices,rows,cols); 
    this->derivative(component,slip::X_DIRECTION,der_order,sch_order,D12);
    //y-derivative
    slip::Matrix3d<double> D22(slices,rows,cols); 
    this->derivative(component,slip::Y_DIRECTION,der_order,sch_order,D22);
    //z-derivative
    slip::Matrix3d<double> D32(slices,rows,cols); 
    this->derivative(component,slip::Z_DIRECTION,der_order,sch_order,D32);
    //derivatives of z component
    component = 2;
    //x-derivative
    slip::Matrix3d<double> D13(slices,rows,cols); 
    this->derivative(component,slip::X_DIRECTION,der_order,sch_order,D13);
    //y-derivative
    slip::Matrix3d<double> D23(slices,rows,cols); 
    this->derivative(component,slip::Y_DIRECTION,der_order,sch_order,D23);
    //z-derivative
    slip::Matrix3d<double> D33(slices,rows,cols); 
    this->derivative(component,slip::Z_DIRECTION,der_order,sch_order,D33);

    //symetric tensor
    slip::Matrix<double> S(3,3);
    //antisymetric tensor
    slip::Matrix<double> O(3,3);
    // tensor to diagonalize S2O2=S^2+O^2
    slip::Matrix<double> S2(3,3),O2(3,3),S2O2(3,3); 
   //  // eigenvalues
//     slip::Vector<double> lambda(6); // 2*3 components <=> 3 complex eigenvalues
//     // eigenvectors
//     slip::Matrix<double> eigenvectors(6,3);// 3 vectors * 2*3 components
//     // sort eigenvalues ? no
//     bool sort=false;
//     // error flag for slip::eigen
//     std::size_t ier;
 
    // eigenvalues
    slip::Vector<double> lambda(3); 
    // eigenvectors
    slip::Matrix<double> eigenvectors(3,3);
   
    Lambda.resize(slices,rows,cols);

    //ça serait mieux avec un iterateur...
    for(size_type k = 0; k < slices; ++k)
      {
                for(size_type i = 0; i < rows; ++i) 
                  {
                    for(size_type j = 0; j < cols; ++j)
                      {
                                S(0,0)=D11[k][i][j];
                                S(1,1)=D22[k][i][j];
                                S(2,2)=D33[k][i][j];
                                S(0,1)=0.5*(D12[k][i][j]+D21[k][i][j]);
                                S(0,2)=0.5*(D13[k][i][j]+D31[k][i][j]);
                                S(1,0)=S(0,1);
                                S(1,2)=0.5*(D23[k][i][j]+D32[k][i][j]);
                                S(2,0)=S(0,2);
                                S(2,1)=S(1,2);
                                O(0,1)=0.5*(D12[k][i][j]-D21[k][i][j]);
                                O(0,2)=0.5*(D13[k][i][j]-D31[k][i][j]);
                                O(1,0)=-O(0,1);
                                O(1,2)=0.5*(D23[k][i][j]-D32[k][i][j]);
                                O(2,0)=-O(0,2);
                                O(2,1)=-O(1,2);
                                slip::matrix_matrix_multiplies(S,S,S2);
                                slip::matrix_matrix_multiplies(O,O,O2);
                                S2O2=0.0;
                                S2O2+=S2;
                                S2O2+=O2;
                                
                //              slip::eigen(S2O2, lambda, ier,sort);
//                              // lambda2 is the 2nd highest eigenvalue of symetric tensor S^2+O^2 (should be real)
//                              // sorting is made upon absolute value of the a priori complex eigenvalues -> no sorting here
                                slip::hermitian_eigen(S2O2,lambda,eigenvectors);
                                // for(std::size_t kc=0; kc<2; ++kc)
//                                {
//                                Lambda[k][i][j][kc]=lambda[2*kc];
//                                }
                                for(std::size_t kc = 0; kc < 2; ++kc)
                                  {
                                    Lambda[k][i][j][kc]=lambda[kc];
                                  }
                      }
                  }
      }

}

  template<typename Container3D>
  void norm(Container3D& result) const
  {
    //typedef slip::Vector3d<double> vector;
    //std::transform(this->begin(),this->end(),result.begin(),std::mem_fun_ref(&vector::Euclidean_norm));
    std::transform(this->begin(),this->end(),result.begin(),std::mem_fun_ref(&value_type::Euclidean_norm));

  }

private:
  slip::Point3d<GridT> init_point_;
  slip::Point3d<GridT> grid_step_;

  

private :
   friend class boost::serialization::access;
 template<class Archive>
    void save(Archive & ar, const unsigned int version) const
    {
      ar & init_point_ & grid_step_;
      ar & boost::serialization::base_object<slip::GenericMultiComponent3d<slip::Vector3d<T> > >(*this);
    }
   template<class Archive>
    void load(Archive & ar, const unsigned int version)
    {
      ar & init_point_ & grid_step_;
      ar & boost::serialization::base_object<slip::GenericMultiComponent3d<slip::Vector3d<T> > >(*this);
    }
  BOOST_SERIALIZATION_SPLIT_MEMBER()


 };

  typedef RegularVector3dField3d<double,double> Regular3D3Cd;
  typedef RegularVector3dField3d<float,double> Regular3D3Cf;
  typedef RegularVector3dField3d<double,double> Regular3D3C_d;
  typedef RegularVector3dField3d<float,double> Regular3D3C_f;
  
  typedef slip::RegularVector3dField3d<double,double> RegularVector3dField3d_d;
  typedef slip::RegularVector3dField3d<float,float> RegularVector3dField3d_f;
  typedef slip::RegularVector3dField3d<long,double> RegularVector3dField3d_l;
  typedef slip::RegularVector3dField3d<unsigned long,double> RegularVector3dField3d_ul;
  typedef slip::RegularVector3dField3d<short,double> RegularVector3dField3d_s;
  typedef slip::RegularVector3dField3d<unsigned short,double> RegularVector3dField3d_us;
  typedef slip::RegularVector3dField3d<int,double> RegularVector3dField3d_i;
  typedef slip::RegularVector3dField3d<unsigned int,double> RegularVector3dField3d_ui;
  typedef slip::RegularVector3dField3d<char,double> RegularVector3dField3d_c;
  typedef slip::RegularVector3dField3d<unsigned char,double> RegularVector3dField3d_uc;


}//slip::

namespace slip{
 
  
template<typename T, typename GridT>
RegularVector3dField3d<T,GridT> operator+(const RegularVector3dField3d<T,GridT>& M1, 
                                                                  const T& val);

  template<typename T, typename GridT>
  RegularVector3dField3d<T,GridT> operator+(const T& val, 
                                                                    const RegularVector3dField3d<T,GridT>& M1);


template<typename T, typename GridT>
RegularVector3dField3d<T,GridT> operator-(const RegularVector3dField3d<T,GridT>& M1, 
                                                                  const T& val);

  template<typename T, typename GridT>
  RegularVector3dField3d<T,GridT> operator-(const T& val, 
                                                                    const RegularVector3dField3d<T,GridT>& M1);



  template<typename T, typename GridT>
  RegularVector3dField3d<T,GridT> operator*(const RegularVector3dField3d<T,GridT>& M1, 
                                                                    const T& val);

  template<typename T, typename GridT>
  RegularVector3dField3d<T,GridT> operator*(const T& val, 
                                                                    const RegularVector3dField3d<T,GridT>& M1);
  


  template<typename T, typename GridT>
  RegularVector3dField3d<T,GridT> operator/(const RegularVector3dField3d<T,GridT>& M1, 
                                                                    const T& val);


template<typename T, typename GridT>
RegularVector3dField3d<T,GridT> operator+(const RegularVector3dField3d<T,GridT>& M1, 
                                                                  const RegularVector3dField3d<T,GridT>& M2);

template<typename T, typename GridT>
RegularVector3dField3d<T,GridT> operator+(const RegularVector3dField3d<T,GridT>& M1, 
                                                                  const slip::Vector3d<T>& val);

  template<typename T, typename GridT>
  RegularVector3dField3d<T,GridT> operator+(const slip::Vector3d<T>& val, 
                                                                    const RegularVector3dField3d<T,GridT>& M1);
  

 
template<typename T, typename GridT>
RegularVector3dField3d<T,GridT> operator-(const RegularVector3dField3d<T,GridT>& M1, 
                                                                  const RegularVector3dField3d<T,GridT>& M2);

 
template<typename T, typename GridT>
RegularVector3dField3d<T,GridT> operator-(const RegularVector3dField3d<T,GridT>& M1, 
                                                                  const slip::Vector3d<T>& val);

  template<typename T, typename GridT>
  RegularVector3dField3d<T,GridT> operator-(const slip::Vector3d<T>& val, 
                                                                    const RegularVector3dField3d<T,GridT>& M1);
  

  
  template<typename T, typename GridT>
  RegularVector3dField3d<T,GridT> operator*(const RegularVector3dField3d<T,GridT>& M1, 
                                      const RegularVector3dField3d<T,GridT>& M2);


template<typename T, typename GridT>
RegularVector3dField3d<T,GridT> operator*(const RegularVector3dField3d<T,GridT>& M1, 
                                                                  const slip::Vector3d<T>& val);

  template<typename T, typename GridT>
  RegularVector3dField3d<T,GridT> operator*(const slip::Vector3d<T>& val, 
                                                                    const RegularVector3dField3d<T,GridT>& M1);
  


  template<typename T, typename GridT>
  RegularVector3dField3d<T,GridT> operator/(const RegularVector3dField3d<T,GridT>& M1, 
                                      const RegularVector3dField3d<T,GridT>& M2);

template<typename T, typename GridT>
RegularVector3dField3d<T,GridT> operator/(const RegularVector3dField3d<T,GridT>& M1, 
                                                                  const slip::Vector3d<T>& val);

 


}//slip::

namespace slip
{


//   template<typename T, typename GridT>
//   inline
//   RegularVector3dField3d<T,GridT>& RegularVector3dField3d<T,GridT>::operator=(const RegularVector3dField3d<T,GridT> & rhs)
//   {
//     if(this != &rhs)
//       {
//              *matrix_ = *(rhs.matrix_);
//       }
//     return *this;
//   }

  template<typename T, typename GridT>
  inline
  RegularVector3dField3d<T,GridT>& RegularVector3dField3d<T,GridT>::operator=(const slip::Vector3d<T>& val)
  {
    std::fill_n(this->begin(),this->size(),val);
    return *this;
  }
  
 
 template<typename T, typename GridT>
  inline
  RegularVector3dField3d<T,GridT>& RegularVector3dField3d<T,GridT>::operator=(const T& val)
  {
    std::fill_n(this->begin(),this->size(),val);
    return *this;
  }


  template<typename T, typename GridT>
  inline
  T& RegularVector3dField3d<T,GridT>::Vx1(const typename RegularVector3dField3d::size_type k,
                                                                  const typename RegularVector3dField3d::size_type i,
                                                                  const typename RegularVector3dField3d::size_type j)
  {
    return (*this)[k][i][j][0];
  }

  template<typename T, typename GridT>
  inline
  const T& RegularVector3dField3d<T,GridT>::Vx1(const typename RegularVector3dField3d::size_type k,
                                                                                const typename RegularVector3dField3d::size_type i,
                                                                                const typename RegularVector3dField3d::size_type j) const
  {
    return (*this)[k][i][j][0];
  }

   template<typename T, typename GridT>
  inline
  T& RegularVector3dField3d<T,GridT>::u(const typename RegularVector3dField3d::size_type k,
                                                                  const typename RegularVector3dField3d::size_type i,
                                                                  const typename RegularVector3dField3d::size_type j)
  {
    return this->Vx1(k,i,j);
  }

  template<typename T, typename GridT>
  inline
  const T& RegularVector3dField3d<T,GridT>::u(const typename RegularVector3dField3d::size_type k,
                                                                                const typename RegularVector3dField3d::size_type i,
                                                                                const typename RegularVector3dField3d::size_type j) const
  {
    return this->Vx1(k,i,j);
  }


   template<typename T, typename GridT>
  inline
  T& RegularVector3dField3d<T,GridT>::Vx2(const typename RegularVector3dField3d::size_type k,
                                                                  const typename RegularVector3dField3d::size_type i,
                                                                  const typename RegularVector3dField3d::size_type j)
  {
    return (*this)[k][i][j][1];
  }


  template<typename T, typename GridT>
  inline
  const T& RegularVector3dField3d<T,GridT>::Vx2(const typename RegularVector3dField3d::size_type k,
                                                                                const typename RegularVector3dField3d::size_type i,
                                                                                const typename RegularVector3dField3d::size_type j) const
  {
    return (*this)[k][i][j][1];
  }
 template<typename T, typename GridT>
 inline
  T& RegularVector3dField3d<T,GridT>::v(const typename RegularVector3dField3d::size_type k,
                                                                  const typename RegularVector3dField3d::size_type i,
                                                                  const typename RegularVector3dField3d::size_type j)
  {
    return this->Vx2(k,i,j);
  }


  template<typename T, typename GridT>
  inline
  const T& RegularVector3dField3d<T,GridT>::v(const typename RegularVector3dField3d::size_type k,
                                                                                const typename RegularVector3dField3d::size_type i,
                                                                                const typename RegularVector3dField3d::size_type j) const
  {
    return this->Vx2(k,i,j);
  }


  template<typename T, typename GridT>
  inline
  T& RegularVector3dField3d<T,GridT>::Vx3(const typename RegularVector3dField3d::size_type k,
                                                                  const typename RegularVector3dField3d::size_type i,
                                                                  const typename RegularVector3dField3d::size_type j)
  {
    return (*this)[k][i][j][2];
  }

   template<typename T, typename GridT>
  inline
  const T& RegularVector3dField3d<T,GridT>::Vx3(const typename RegularVector3dField3d::size_type k,
                                                                                const typename RegularVector3dField3d::size_type i,
                                                                                const typename RegularVector3dField3d::size_type j) const
  {
    return (*this)[k][i][j][2];
  } 
 
   template<typename T, typename GridT>
  inline
  T& RegularVector3dField3d<T,GridT>::w(const typename RegularVector3dField3d::size_type k,
                                                                  const typename RegularVector3dField3d::size_type i,
                                                                  const typename RegularVector3dField3d::size_type j)
  {
    return this->Vx3(k,i,j);
  }

   template<typename T, typename GridT>
  inline
  const T& RegularVector3dField3d<T,GridT>::w(const typename RegularVector3dField3d::size_type k,
                                                                                const typename RegularVector3dField3d::size_type i,
                                                                                const typename RegularVector3dField3d::size_type j) const
  {
    return this->Vx3(k,i,j);
  } 
 

  template<typename T, typename GridT>
  inline
  typename RegularVector3dField3d<T,GridT>::norm_type 
  RegularVector3dField3d<T,GridT>::norm(const typename RegularVector3dField3d::size_type k,
                                                                const typename RegularVector3dField3d::size_type i,
                                                                const typename RegularVector3dField3d::size_type j) const
  {
    return (*this)[k][i][j].Euclidean_norm();
  }

   template<typename T, typename GridT>
  inline
   const GridT RegularVector3dField3d<T,GridT>::x1(const typename RegularVector3dField3d::size_type k,
                                                                                const typename RegularVector3dField3d::size_type i,
                                                                                const typename RegularVector3dField3d::size_type j) const
  {
    return this->init_point_[0] + GridT(j) * this->grid_step_[0];
  }

    template<typename T, typename GridT>
  inline
   const GridT RegularVector3dField3d<T,GridT>::x2(const typename RegularVector3dField3d::size_type k,
                                                                                const typename RegularVector3dField3d::size_type i,
                                                                                const typename RegularVector3dField3d::size_type j) const
  {
    return this->init_point_[1] + GridT(this->rows() - 1 - i) * this->grid_step_[1];
  }

  template<typename T, typename GridT>
  inline
   const GridT RegularVector3dField3d<T,GridT>::x3(const typename RegularVector3dField3d::size_type k,
                                                                                const typename RegularVector3dField3d::size_type i,
                                                                                const typename RegularVector3dField3d::size_type j) const
  {
    return this->init_point_[2] + GridT(this->slices() - 1 - k) * this->grid_step_[2];
  }

 template<typename T, typename GridT>
  inline
   const GridT RegularVector3dField3d<T,GridT>::x(const typename RegularVector3dField3d::size_type k,
                                                                                const typename RegularVector3dField3d::size_type i,
                                                                                const typename RegularVector3dField3d::size_type j) const
  {
    return this->x1(k,i,j);
  }

    template<typename T, typename GridT>
  inline
   const GridT RegularVector3dField3d<T,GridT>::y(const typename RegularVector3dField3d::size_type k,
                                                                                const typename RegularVector3dField3d::size_type i,
                                                                                const typename RegularVector3dField3d::size_type j) const
  {
    return this->x2(k,i,j);
  }

  template<typename T, typename GridT>
  inline
   const GridT RegularVector3dField3d<T,GridT>::z(const typename RegularVector3dField3d::size_type k,
                                                                                const typename RegularVector3dField3d::size_type i,
                                                                                const typename RegularVector3dField3d::size_type j) const
  {
    return this->x3(k,i,j);
  }

 template<typename T, typename GridT>
  inline
  const slip::Point3d<GridT>& RegularVector3dField3d<T,GridT>::get_init_point() const
  {
    return this->init_point_;
  }

  template<typename T, typename GridT>
  inline
  void RegularVector3dField3d<T,GridT>::set_init_point(const slip::Point3d<GridT>& init_point)
  {
    this->init_point_ = init_point;
  }

  template<typename T, typename GridT>
  inline
  const slip::Point3d<GridT>& RegularVector3dField3d<T,GridT>::get_grid_step() const
  {
    return this->grid_step_;
  }

  template<typename T, typename GridT>
  inline
  void RegularVector3dField3d<T,GridT>::set_grid_step(const slip::Point3d<GridT>& grid_step)
  {
    this->grid_step_ = grid_step;
  }

  template<typename T, typename GridT>
  inline
  RegularVector3dField3d<T,GridT>& RegularVector3dField3d<T,GridT>::operator+=(const T&val)
  {
     std::transform(this->begin(),this->end(),this->begin(),std::bind2nd(std::plus<slip::Vector3d<T> >(),val));
    return *this;
  }

  template<typename T, typename GridT>
  inline
  RegularVector3dField3d<T,GridT>& RegularVector3dField3d<T,GridT>::operator-=(const T&val)
  {
     std::transform(this->begin(),this->end(),this->begin(),std::bind2nd(std::minus<slip::Vector3d<T> >(),val));
    return *this;
  }

  template<typename T, typename GridT>
  inline
  RegularVector3dField3d<T,GridT>& RegularVector3dField3d<T,GridT>::operator*=(const T&val)
  {
     std::transform(this->begin(),this->end(),this->begin(),std::bind2nd(std::multiplies<slip::Vector3d<T> >(),val));
    return *this;
  }

  template<typename T, typename GridT>
  inline
  RegularVector3dField3d<T,GridT>& RegularVector3dField3d<T,GridT>::operator/=(const T&val)
  {
     std::transform(this->begin(),this->end(),this->begin(),std::bind2nd(std::divides<slip::Vector3d<T> >(),val));
    return *this;
  }



 template<typename T, typename GridT>
  inline
  RegularVector3dField3d<T,GridT>& RegularVector3dField3d<T,GridT>::operator+=(const slip::Vector3d<T>&val)
  {
     std::transform(this->begin(),this->end(),this->begin(),std::bind2nd(std::plus<slip::Vector3d<T> >(),val));
    return *this;
  }

  template<typename T, typename GridT>
  inline
  RegularVector3dField3d<T,GridT>& RegularVector3dField3d<T,GridT>::operator-=(const slip::Vector3d<T>&val)
  {
     std::transform(this->begin(),this->end(),this->begin(),std::bind2nd(std::minus<slip::Vector3d<T> >(),val));
    return *this;
  }

  template<typename T, typename GridT>
  inline
  RegularVector3dField3d<T,GridT>& RegularVector3dField3d<T,GridT>::operator*=(const slip::Vector3d<T>&val)
  {
     std::transform(this->begin(),this->end(),this->begin(),std::bind2nd(std::multiplies<slip::Vector3d<T> >(),val));
    return *this;
  }

  template<typename T, typename GridT>
  inline
  RegularVector3dField3d<T,GridT>& RegularVector3dField3d<T,GridT>::operator/=(const slip::Vector3d<T>&val)
  {
     std::transform(this->begin(),this->end(),this->begin(),std::bind2nd(std::divides<slip::Vector3d<T> >(),val));
    return *this;
  }
  

  
  template<typename T, typename GridT>
  inline
  RegularVector3dField3d<T,GridT> RegularVector3dField3d<T,GridT>::operator-() const
  {
    RegularVector3dField3d<T,GridT> tmp(*this);
    std::transform(this->begin(),this->end(),tmp.begin(),std::negate<slip::Vector3d<T> >());
    return tmp;
  }

  template<typename T, typename GridT>
  inline
  RegularVector3dField3d<T,GridT>& RegularVector3dField3d<T,GridT>::operator+=(const RegularVector3dField3d<T,GridT>& rhs)
  {
    assert(this->dim1() == rhs.dim1());
    assert(this->dim2() == rhs.dim2()); 
    assert(this->dim3() == rhs.dim3()); 
     if(this->get_init_point() != rhs.get_init_point())
      {
                throw std::runtime_error(slip::GRID_INIT_POINT_ERROR);
      }
     if(this->get_grid_step() != rhs.get_grid_step())
      {
                throw std::runtime_error(slip::GRID_STEP_ERROR);
      }
    std::transform(this->begin(),this->end(),rhs.begin(),this->begin(),std::plus<typename RegularVector3dField3d::value_type>());
    return *this;
  }

   template<typename T, typename GridT>
  inline
  RegularVector3dField3d<T,GridT>& RegularVector3dField3d<T,GridT>::operator-=(const RegularVector3dField3d<T,GridT>& rhs)
  {
    assert(this->dim1() == rhs.dim1());
    assert(this->dim2() == rhs.dim2());
    assert(this->dim3() == rhs.dim3()); 
     if(this->get_init_point() != rhs.get_init_point())
      {
                throw std::runtime_error(slip::GRID_INIT_POINT_ERROR);
      }
     if(this->get_grid_step() != rhs.get_grid_step())
      {
                throw std::runtime_error(slip::GRID_STEP_ERROR);
      }
    std::transform(this->begin(),this->end(),rhs.begin(),this->begin(),std::minus<typename RegularVector3dField3d::value_type>());
    return *this;
  }

  template<typename T, typename GridT>
  inline
  RegularVector3dField3d<T,GridT>& RegularVector3dField3d<T,GridT>::operator*=(const RegularVector3dField3d<T,GridT>& rhs)
  {
    assert(this->dim1() == rhs.dim1());
    assert(this->dim2() == rhs.dim2());
    assert(this->dim3() == rhs.dim3()); 
     if(this->get_init_point() != rhs.get_init_point())
      {
                throw std::runtime_error(slip::GRID_INIT_POINT_ERROR);
      }
     if(this->get_grid_step() != rhs.get_grid_step())
      {
                throw std::runtime_error(slip::GRID_STEP_ERROR);
      }
    std::transform(this->begin(),this->end(),rhs.begin(),this->begin(),std::multiplies<typename RegularVector3dField3d::value_type>());
    return *this;
  }

  template<typename T, typename GridT>
  inline
  RegularVector3dField3d<T,GridT>& RegularVector3dField3d<T,GridT>::operator/=(const RegularVector3dField3d<T,GridT>& rhs)
  {
    assert(this->dim1() == rhs.dim1());
    assert(this->dim2() == rhs.dim2());
    assert(this->dim3() == rhs.dim3()); 
     if(this->get_init_point() != rhs.get_init_point())
      {
                throw std::runtime_error(slip::GRID_INIT_POINT_ERROR);
      }
     if(this->get_grid_step() != rhs.get_grid_step())
      {
                throw std::runtime_error(slip::GRID_STEP_ERROR);
      }
    std::transform(this->begin(),this->end(),rhs.begin(),this->begin(),std::divides<typename RegularVector3dField3d::value_type>());
    return *this;
  }


  template <typename T, typename GridT>
  inline
  std::ostream& operator<<(std::ostream & out, 
                                                   const RegularVector3dField3d<T,GridT>& a)
  {
    
    out<<"init point: "<<a.init_point_<<std::endl;
    out<<"grid step: "<<a.grid_step_<<std::endl;
    out<<"data: "<<std::endl;
    typedef typename RegularVector3dField3d<T,GridT>::size_type size_type;
    const size_type slices  = a.slices();
    const size_type rows  = a.rows();
    const size_type cols = a.cols();
    for(size_type k = 0; k < slices; ++k)
      {
                for(size_type i = 0; i < rows; ++i)
                  {
                    for(size_type j = 0; j < cols; ++j)
                      {
                                out<<a[k][i][j]<<" ";
                      }
                    out<<std::endl;
                  }             
                out<<std::endl;
      }
    return out;
  }

  template<typename T, typename GridT>
  inline
  std::string 
  RegularVector3dField3d<T,GridT>::name() const {return "RegularVector3dField3d";} 
 

 //  template<typename T, typename GridT>
//   inline
//   void RegularVector3dField3d<T,GridT>::write_gnuplot(const std::string& file_path_name)
//   {
//     slip::write_gnuplot_vect3d<RegularVector3dField3d<T,GridT> >(*this,file_path_name);
//   }

  template<typename T, typename GridT>
  inline
  void RegularVector3dField3d<T,GridT>::write_tecplot(const std::string& file_path_name,
                                                                                      const std::string& title,
                                                                                      const std::string& zone) const
  {
    slip::write_tecplot_vect3d<RegularVector3dField3d<T,GridT> >(*this,
                                                                                                                   this->init_point_,
      this->grid_step_,
      file_path_name,
      title,zone);
 
  }

  template<typename T, typename GridT>
  inline
  void RegularVector3dField3d<T,GridT>::write_tecplot_bin(const std::string& file_path_name,
                                                                                                                  const std::string& title,
                                                                                                                  const std::string& zone) const
  {
    slip::RegularVector3dField3d_to_plt(file_path_name,
                                                                                *this,
                                                                                title,
                                                                                zone);
 
  }

  //  template<typename T, typename GridT>
//   inline
//   void RegularVector3dField3d<T,GridT>::read_gnuplot(const std::string& file_path_name)
//   {
//     slip::read_gnuplot_vect3d<RegularVector3dField3d<T,GridT> >(file_path_name,*this);
//   }

  template<typename T, typename GridT>
  inline
  void RegularVector3dField3d<T,GridT>::read_tecplot(const std::string& file_path_name)
  {
    slip::read_tecplot_vect3d<RegularVector3dField3d<T,GridT> >(file_path_name,
                                                                                                                  *this,
                                                                                                                  this->init_point_,
      this->grid_step_);
  }

   template<typename T, typename GridT>
  inline
   void RegularVector3dField3d<T,GridT>::read_tecplot_bin(const std::string& file_path_name,
                                                                                                                  const int zone_loaded)
  {
    std::string _useless_string;
    std::vector<std::string> _useless_varnames;
    slip::plt_to_RegularVector3dField3d(file_path_name, 
                                                                                *this, 
                                                                                _useless_string, 
                                                                                _useless_string, 
                                                                                _useless_varnames,
                                                                                zone_loaded);
  }

  template<typename T, typename GridT>
  inline
  RegularVector3dField3d<T,GridT> operator+(const RegularVector3dField3d<T,GridT>& M1, 
                                                                    const T& val)
  {
    RegularVector3dField3d<T,GridT>  tmp = M1;
    tmp += val;
    return tmp;
  }

  template<typename T, typename GridT>
  inline
  RegularVector3dField3d<T,GridT> operator+(const T& val,
                                                                    const RegularVector3dField3d<T,GridT>& M1)
  {
    return M1 + val;
  }

  template<typename T, typename GridT>
  inline
  RegularVector3dField3d<T,GridT> operator-(const RegularVector3dField3d<T,GridT>& M1, 
                                                                    const T& val)
  {
    RegularVector3dField3d<T,GridT>  tmp = M1;
    tmp -= val;
    return tmp;
  }

  template<typename T, typename GridT>
  inline
  RegularVector3dField3d<T,GridT> operator-(const T& val,
                                                                    const RegularVector3dField3d<T,GridT>& M1)
  {
    return -(M1 - val);
  }


  template<typename T, typename GridT>
  inline
  RegularVector3dField3d<T,GridT> operator*(const RegularVector3dField3d<T,GridT>& M1, 
                                                                    const T& val)
  {
    RegularVector3dField3d<T,GridT>  tmp = M1;
    tmp *= val;
    return tmp;
  }

  template<typename T, typename GridT>
  inline
  RegularVector3dField3d<T,GridT> operator*(const T& val,
                                                                    const RegularVector3dField3d<T,GridT>& M1)
  {
    return M1 * val;
  }


  template<typename T, typename GridT>
  inline
  RegularVector3dField3d<T,GridT> operator/(const RegularVector3dField3d<T,GridT>& M1, 
                                                                    const T& val)
  {
    RegularVector3dField3d<T,GridT>  tmp = M1;
    tmp /= val;
    return tmp;
  }

  template<typename T, typename GridT>
  inline
  RegularVector3dField3d<T,GridT> operator+(const RegularVector3dField3d<T,GridT>& M1, 
                                                                    const RegularVector3dField3d<T,GridT>& M2)
  {

    RegularVector3dField3d<T,GridT> tmp(M1);
    tmp+=M2;
    return tmp;
  }


  template<typename T, typename GridT>
  inline
  RegularVector3dField3d<T,GridT> operator+(const RegularVector3dField3d<T,GridT>& M1, 
                                                      const slip::Vector3d<T>&val)
  {
    RegularVector3dField3d<T,GridT>  tmp = M1;
    tmp += val;
    return tmp;
  }

  template<typename T, typename GridT>
  inline
  RegularVector3dField3d<T,GridT> operator+(const slip::Vector3d<T>&val,
                                                      const RegularVector3dField3d<T,GridT>& M1)
  {
    return M1 + val;
  }



 template<typename T, typename GridT>
  inline
  RegularVector3dField3d<T,GridT> operator-(const RegularVector3dField3d<T,GridT>& M1, 
                                                                    const RegularVector3dField3d<T,GridT>& M2)
  {
    
    RegularVector3dField3d<T,GridT> tmp(M1);
    tmp-=M2;
    return tmp;
  }


  template<typename T, typename GridT>
  inline
  RegularVector3dField3d<T,GridT> operator-(const RegularVector3dField3d<T,GridT>& M1, 
                                                      const slip::Vector3d<T>&val)
  {
    RegularVector3dField3d<T,GridT>  tmp = M1;
    tmp -= val;
    return tmp;
  }

  template<typename T, typename GridT>
  inline
  RegularVector3dField3d<T,GridT> operator-(const slip::Vector3d<T>&val,
                                                      const RegularVector3dField3d<T,GridT>& M1)
  {
    return -(M1 - val);
  }


  template<typename T, typename GridT>
  inline
  RegularVector3dField3d<T,GridT> operator*(const RegularVector3dField3d<T,GridT>& M1, 
                                                                    const RegularVector3dField3d<T,GridT>& M2)
  {
   
    RegularVector3dField3d<T,GridT> tmp(M1);
    tmp*=M2;
    return tmp;
  }


  template<typename T, typename GridT>
  inline
  RegularVector3dField3d<T,GridT> operator*(const RegularVector3dField3d<T,GridT>& M1, 
                                                                    const slip::Vector3d<T>&val)
  {
    RegularVector3dField3d<T,GridT>  tmp = M1;
    tmp *= val;
    return tmp;
  }

  template<typename T, typename GridT>
  inline
  RegularVector3dField3d<T,GridT> operator*(const slip::Vector3d<T>&val,
                                                                    const RegularVector3dField3d<T,GridT>& M1)
  {
    return M1 * val;
  }


 template<typename T, typename GridT>
  inline
  RegularVector3dField3d<T,GridT> operator/(const RegularVector3dField3d<T,GridT>& M1, 
                                                                    const RegularVector3dField3d<T,GridT>& M2)
  {
   
    RegularVector3dField3d<T,GridT> tmp(M1);
    tmp/=M2;
    return tmp;
  }


  template<typename T, typename GridT>
  inline
  RegularVector3dField3d<T,GridT> operator/(const RegularVector3dField3d<T,GridT>& M1, 
                                                                    const slip::Vector3d<T>&val)
  {
    RegularVector3dField3d<T,GridT>  tmp = M1;
    tmp /= val;
    return tmp;
  }




}//slip::

#endif //SLIP_REGULARVECTOR3DFIELD3D_HPP