RegularVector2dField2d.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_REGULARVECTOR2DFIELD2D_HPP
#define SLIP_REGULARVECTOR2DFIELD2D_HPP

#include <iostream>
#include <fstream>
#include <iterator>
#include <cassert>
#include <numeric>
#include <algorithm>
#include <cmath>
#include <string>
#include <cstddef>
#include "Matrix.hpp"
#include "stride_iterator.hpp"
#include "kstride_iterator.hpp"
#include "iterator2d_box.hpp"
#include "iterator2d_range.hpp"
#include "apply.hpp"
#include "Vector2d.hpp"
#include "Point2d.hpp"
#include "GenericMultiComponent2d.hpp"
#include "derivatives.hpp"
#include "arithmetic_op.hpp"
#include "io_tools.hpp"
#include "linear_algebra_traits.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<typename T>
class stride_iterator;

template<typename T>
class iterator2d_box;

template<typename T>
class iterator2d_range;

template<typename T>
class const_iterator2d_box;

template<typename T>
class const_iterator2d_range;


template <typename T>
class DPoint2d;

template <typename T>
class Point2d;

template <typename T>
class Box2d;

  template <typename T, typename GridT>
class RegularVector2dField2d;

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

  template <typename T, typename GridT = double>
class RegularVector2dField2d:public slip::GenericMultiComponent2d<slip::Vector2d<T> >
{
public :

  typedef slip::Vector2d<T> value_type;
  typedef RegularVector2dField2d<T,GridT> self;
  typedef const RegularVector2dField2d<T,GridT> const_self;
  typedef GenericMultiComponent2d<slip::Vector2d<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::GenericMultiComponent2d<slip::Vector2d<T> >::iterator2d iterator2d;
  typedef typename slip::GenericMultiComponent2d<slip::Vector2d<T> >::const_iterator2d const_iterator2d;

  typedef T vector2d_value_type;
  typedef GridT grid_value_type;
  typedef vector2d_value_type* vector2d_pointer;
  typedef const vector2d_value_type* const_vector2d_pointer;
  typedef vector2d_value_type& vector2d_reference;
  typedef const vector2d_value_type const_vector2d_reference;
  typedef slip::kstride_iterator<vector2d_pointer,2> vector2d_iterator;
  typedef slip::kstride_iterator<const_vector2d_pointer,2> const_vector2d_iterator;

  typedef typename slip::lin_alg_traits<value_type>::value_type norm_type;

  static const std::size_t DIM = 2;
public:

  RegularVector2dField2d():
    base(),init_point_(slip::Point2d<GridT>()),grid_step_(slip::Point2d<GridT>(GridT(1),GridT(1)))
  {}
 
  RegularVector2dField2d(const size_type height,
                                                 const size_type width):
    base(height,width),init_point_(slip::Point2d<GridT>()),grid_step_(slip::Point2d<GridT>(GridT(1),GridT(1)))
  {}

  RegularVector2dField2d(const size_type height,
                                                 const size_type width,
                                                 const slip::Point2d<GridT>& init_point,
                                                 const slip::Point2d<GridT>& grid_step):
    base(height,width),init_point_(init_point),grid_step_(grid_step)
  {}
  
  RegularVector2dField2d(const size_type height,
                     const size_type width, 
                     const slip::Vector2d<T>& val):
    base(height,width,val),init_point_(slip::Point2d<GridT>()),grid_step_(slip::Point2d<GridT>(GridT(1),GridT(1)))
  {}
  
  RegularVector2dField2d(const size_type height,
                                                 const size_type width, 
                                                 const slip::Point2d<GridT>& init_point,
                                                 const slip::Point2d<GridT>& grid_step,
                                                 const slip::Vector2d<T>& val):
    base(height,width,val),init_point_(init_point),grid_step_(grid_step)
  {}
 
  RegularVector2dField2d(const size_type height,
                                                 const size_type width, 
                                                 const T* val):
    base(height,width,val),init_point_(slip::Point2d<GridT>()),grid_step_(slip::Point2d<GridT>(GridT(1),GridT(1)))
  {}

 


  RegularVector2dField2d(const size_type height,
                                                 const size_type width, 
                                                 const slip::Point2d<GridT>& init_point,
                                                 const slip::Point2d<GridT>& grid_step,
                                                 const T* val):
    base(height,width,val),init_point_(init_point),grid_step_(grid_step)
  {}
 
  RegularVector2dField2d(const size_type height,
                                                 const size_type width, 
                                                 const slip::Vector2d<T>* val):
    base(height,width,val),init_point_(slip::Point2d<GridT>()),grid_step_(slip::Point2d<GridT>(GridT(1),GridT(1)))
  {}
 
  RegularVector2dField2d(const size_type height,
                                                 const size_type width, 
                                                 const slip::Point2d<GridT>& init_point,
                                                 const slip::Point2d<GridT>& grid_step,
                                                 const slip::Vector2d<T>* val):
    base(height,width,val),init_point_(init_point),grid_step_(grid_step)
  {}
 

  template<typename InputIterator>
  RegularVector2dField2d(const size_type height,
                                       const size_type width, 
                                       InputIterator first,
                                       InputIterator last):
    base(height,width,first,last),init_point_(slip::Point2d<GridT>()),grid_step_(slip::Point2d<GridT>(GridT(1),GridT(1)))
  {
    // this->init_point_ = slip::Point2d<GridT>();
//     this->grid_step_  = slip::Point2d<GridT>(GridT(1),GridT(1));
  } 


  template<typename InputIterator>
  RegularVector2dField2d(const size_type height,
                                                 const size_type width, 
                                                 const slip::Point2d<GridT>& init_point,
                                                 const slip::Point2d<GridT>& grid_step,
                                                 InputIterator first,
                                                 InputIterator last):
    base(height,width,first,last),init_point_(init_point),grid_step_(grid_step)
  {
  //   this->init_point_ = init_point;
//     this->grid_step_ = grid_step;

  } 



  
  template<typename InputIterator>
  RegularVector2dField2d(const size_type height,
                                       const size_type width, 
                                       InputIterator first1,
                                       InputIterator last1,
                                       InputIterator first2):
    base(height,width),init_point_(slip::Point2d<GridT>()),grid_step_(slip::Point2d<GridT>(GridT(1),GridT(1)))
  {
    
    std::vector<InputIterator> first_iterators_list(2);
    first_iterators_list[0] = first1;
    first_iterators_list[1] = first2;
    this->fill(first_iterators_list,last1);
  
  } 

  template<typename InputIterator>
  RegularVector2dField2d(const size_type height,
                                                 const size_type width, 
                                                 const slip::Point2d<GridT>& init_point,
                                                 const slip::Point2d<GridT>& grid_step,
                                                 InputIterator first1,
                                                 InputIterator last1,
                                                 InputIterator first2):
    base(height,width),init_point_(init_point),grid_step_(grid_step)
  {
    
    std::vector<InputIterator> first_iterators_list(2);
    first_iterators_list[0] = first1;
    first_iterators_list[1] = first2;
    this->fill(first_iterators_list,last1);
  
  } 


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

  ~RegularVector2dField2d()
  {}
  

  std::string name() const;


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

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


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


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


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


 

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


  self& operator=(const T& val);



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

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

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

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

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

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


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


  norm_type angle(const size_type i,
                                   const size_type j) const;
  const GridT x1(const size_type i,
                                 const size_type j) const;

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

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

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

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


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

  const slip::Point2d<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 T& val);
//   self& operator^=(const T& val);
//   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;

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

  

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


  template<typename Container2D>
  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,
                                  Container2D& result) const
  {
    assert(component < 2);
    assert(sch_order >= der_order);
    assert(   ((der_dir == X_DIRECTION) && (this->rows() > sch_order))
                    ||((der_dir == Y_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();
                for(size_type i = 0; i < rows; ++i)
                  {
                    slip::derivative(this->row_begin(component,i),
                                                     this->row_end(component,i), 
                                                     der_order,
                                                     sch_order,
                                                     sch_shift,
                                                     kernels_iterators,
                                                     result.row_begin(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();
                for(size_type j = 0; j < cols; ++j)
                  {
                    slip::derivative(this->col_rbegin(component,j),
                                                     this->col_rend(component,j), 
                                                     der_order,
                                                     sch_order,
                                                     sch_shift,
                                                     kernels_iterators,
                                                     result.col_rbegin(j));
                  }
      }
    // slip::derivative_2d(this->upper_left(component),
//                                              this->bottom_right(component),
//                                              this->grid_step_,
//                                              der_dir,
//                                              der_order,
//                                              sch_order,
//                                              result.upper_left(),
//                                              result.bottom_right());
                                                
  }


  template<typename Container2D>
  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::Box2d<int> box,
      Container2D& result) const
  {
    assert(component < 2);
    assert(sch_order >= der_order);
    assert(   ((der_dir == X_DIRECTION) && (this->rows() > sch_order))
                      ||((der_dir == Y_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 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.upper_left()[1],box.bottom_right()[1],1);
                for(int i = box.upper_left()[0]; i <= box.bottom_right()[0]; ++i)
                  {
                    slip::derivative(this->row_begin(component,i, range_row),
                                                     this->row_end(component,i, range_row), 
                                                     der_order,
                                                     sch_order,
                                                     sch_shift,
                                                     kernels_iterators,
                                                     result.row_begin(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.upper_left()[0],box.bottom_right()[0],1);
                for(int j = box.upper_left()[1]; j <= box.bottom_right()[1]; ++j)
                  {
                    slip::derivative(this->col_rbegin(component,j, range_col),
                                                     this->col_rend(component,j, range_col), 
                                                     der_order,
                                                     sch_order,
                                                     sch_shift,
                                                     kernels_iterators,
                                                     result.col_rbegin(j, range_col));
                  }
      }
  }


  template<typename Container2D>
  void divergence(const std::size_t sch_order,
                                  Container2D& 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::Matrix<double> Mtmp(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);

    //compute the divergence
    slip::plus(Mtmp.begin(),Mtmp.end(),result.begin(),result.begin());
                
    
  }


  template<typename Container2D>
  void vorticity(const std::size_t sch_order,
                                 Container2D& result) const
  {
    assert(sch_order >= 1);
    assert((sch_order / 2) <= sch_order);
        
    const std::size_t der_order = 1;
    
    //derivate the y component (second component)
    slip::Matrix<double> Mtmp(this->rows(),this->cols()); 
    std::size_t component = 1;
    this->derivative(component,slip::X_DIRECTION,der_order,sch_order,Mtmp);

    //derivate the x component (first component) 
     component = 0;
    this->derivative(component,slip::Y_DIRECTION,der_order,sch_order,result);
    
    //compute the divergence
    slip::minus(Mtmp.begin(),Mtmp.end(),result.begin(),result.begin());
    
  }


  template<typename Container2D>
  void lambda2(const std::size_t sch_order,
                       Container2D& result)
  {
    assert(sch_order >= 1);
    assert((sch_order / 2) <= sch_order);
    
    const std::size_t der_order = 1;
    const size_type rows = this->rows();
    const size_type cols = this->cols();

    //vectorial x derivative of the RegularVector2dField2d
    self Mtmpx(rows,cols); 
    slip::derivative_2d(this->upper_left(),
                                                this->bottom_right(),
                                                this->grid_step_,
                                                slip::X_DIRECTION,
                                                der_order,
                                                sch_order,
                                                Mtmpx.upper_left(),
                                                Mtmpx.bottom_right());
    
    //vectorial y derivative of the RegularVector2dField2d
    self Mtmpy(rows,cols); 
    slip::derivative_2d(this->upper_left(),
                                                this->bottom_right(),
                                                this->grid_step_,
                                                slip::Y_DIRECTION,
                                                der_order,
                                                sch_order,
                                                Mtmpy.upper_left(),
                                                Mtmpy.bottom_right());
    //compute the lamda2
    for(size_type i = 0; i < rows; ++i)
      {
                for(size_type j = 0; j < cols; ++j)
                  {
                    result[i][j] = Mtmpx[i][j][1] * Mtmpy[i][j][0] - Mtmpx[i][j][0] * Mtmpy[i][j][1];
                  }
      }
                
    
  }


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

  }

  template<typename Container2D>
  void angle(Container2D& result)
  {
    
        std::transform(this->begin(),this->end(),result.begin(),std::mem_fun_ref(&slip::Vector2d<T>::angle));

  }

private:
  slip::Point2d<GridT> init_point_;
  slip::Point2d<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::GenericMultiComponent2d<slip::Vector2d<T> > >(*this);
    }
   template<class Archive>
    void load(Archive & ar, const unsigned int version)
    {
      ar & init_point_ & grid_step_;
      ar & boost::serialization::base_object<slip::GenericMultiComponent2d<slip::Vector2d<T> > >(*this);
    }
  BOOST_SERIALIZATION_SPLIT_MEMBER()
 };

  typedef RegularVector2dField2d<double,double> Regular2D3Cd;
  typedef RegularVector2dField2d<float,double> Regular2D3Cf;
  typedef RegularVector2dField2d<double,double> Regular2D3C_d;
  typedef RegularVector2dField2d<float,double> Regular2D3C_f;

  typedef slip::RegularVector2dField2d<double,double> RegularVector2dField2d_d;
  typedef slip::RegularVector2dField2d<float,double> RegularVector2dField2d_f;
  typedef slip::RegularVector2dField2d<long,double> RegularVector2dField2d_l;
  typedef slip::RegularVector2dField2d<unsigned long,double> RegularVector2dField2d_ul;
  typedef slip::RegularVector2dField2d<short,double> RegularVector2dField2d_s;
  typedef slip::RegularVector2dField2d<unsigned short,double> RegularVector2dField2d_us;
  typedef slip::RegularVector2dField2d<int,double> RegularVector2dField2d_i;
  typedef slip::RegularVector2dField2d<unsigned int,double> RegularVector2dField2d_ui;
  typedef slip::RegularVector2dField2d<char,double> RegularVector2dField2d_c;
  typedef slip::RegularVector2dField2d<unsigned char,double> RegularVector2dField2d_uc;

  

}//slip::

namespace slip{
 

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

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


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

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



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

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


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


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

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

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

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

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

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

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


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

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


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

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

 


}//slip::

namespace slip
{


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

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


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

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

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

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

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

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

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

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

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

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

  template<typename T, typename GridT>
  inline
  const GridT RegularVector2dField2d<T,GridT>::x2(const typename RegularVector2dField2d::size_type i,
                                                                                const typename RegularVector2dField2d::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 RegularVector2dField2d<T,GridT>::x(const typename RegularVector2dField2d::size_type i,
                                                                                const typename RegularVector2dField2d::size_type j) const
  {
    return this->x1(i,j);
  }

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

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

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

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

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

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

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

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

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



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

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

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

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

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

  template<typename T, typename GridT>
  inline
  RegularVector2dField2d<T,GridT>& RegularVector2dField2d<T,GridT>::operator+=(const RegularVector2dField2d<T,GridT>& rhs)
  {
    assert(this->dim1() == rhs.dim1());
    assert(this->dim2() == rhs.dim2()); 
     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 RegularVector2dField2d::value_type>());
    return *this;
  }

   template<typename T, typename GridT>
  inline
  RegularVector2dField2d<T,GridT>& RegularVector2dField2d<T,GridT>::operator-=(const RegularVector2dField2d<T,GridT>& rhs)
  {
    assert(this->dim1() == rhs.dim1());
    assert(this->dim2() == rhs.dim2());
    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 RegularVector2dField2d::value_type>());
    return *this;
  }

  template<typename T, typename GridT>
  inline
  RegularVector2dField2d<T,GridT>& RegularVector2dField2d<T,GridT>::operator*=(const RegularVector2dField2d<T,GridT>& rhs)
  {
    assert(this->dim1() == rhs.dim1());
    assert(this->dim2() == rhs.dim2());
    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 RegularVector2dField2d::value_type>());
    return *this;
  }

  template<typename T, typename GridT>
  inline
  RegularVector2dField2d<T,GridT>& RegularVector2dField2d<T,GridT>::operator/=(const RegularVector2dField2d<T,GridT>& rhs)
  {
    assert(this->dim1() == rhs.dim1());
    assert(this->dim2() == rhs.dim2());
    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 RegularVector2dField2d::value_type>());
    return *this;
  }

  template <typename T, typename GridT>
  inline
  std::ostream& operator<<(std::ostream & out, 
                                                   const RegularVector2dField2d<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 RegularVector2dField2d<T,GridT>::size_type size_type;
    size_type rows  = a.rows();
    size_type cols = a.cols();
    
    for(size_type i = 0; i < rows; ++i)
      {
                for(size_type j = 0; j < cols; ++j)
                  {
                    out<<a[i][j]<<" ";
                  }
                out<<std::endl;
      }         
    return out;
  }


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

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

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

  template<typename T, typename GridT>
  inline
  void RegularVector2dField2d<T,GridT>::read_tecplot(const std::string& file_path_name)

  {
    slip::read_tecplot_vect2d<slip::RegularVector2dField2d<T,GridT> >(file_path_name,*this,this->init_point_,this->grid_step_);
  }


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

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

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

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


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

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


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

  template<typename T, typename GridT>
  inline
  RegularVector2dField2d<T,GridT> operator+(const RegularVector2dField2d<T,GridT>& M1, 
                                                                    const RegularVector2dField2d<T,GridT>& M2)
  {
    RegularVector2dField2d<T,GridT> tmp(M1);
    tmp+=M2;

    return tmp;
  }


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

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



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


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

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


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

    return tmp;
  }


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

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


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

    return tmp;
  }


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



}//slip::

#endif //SLIP_REGULARVECTOR2DFIELD2D_HPP