DenseVector3dField2d.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_DENSEVECTOR3DFIELD2D_HPP
#define SLIP_DENSEVECTOR3DFIELD2D_HPP

#include <iostream>
#include <fstream>
#include <iterator>
#include <cassert>
#include <numeric>
#include <cmath>
#include <string>
#include <cstddef>
#include "Matrix.hpp"
#include "stride_iterator.hpp"
#include "kstride_iterator.hpp"
#include "iterator2d_box.hpp"
#include "apply.hpp"
#include "Vector3d.hpp"
#include "GenericMultiComponent2d.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 iterator2d_box;

template <class T>
class DPoint2d;

template <class T>
class Point2d;

template <class T>
class Box2d;

template <typename T>
class DenseVector3dField2d;


template <typename T>
class DenseVector3dField2d:public GenericMultiComponent2d<slip::Vector3d<T> >
{
public :

  typedef slip::Vector3d<T> value_type;
  typedef DenseVector3dField2d<T> self;
  typedef GenericMultiComponent2d<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::GenericMultiComponent2d<Vector3d<T> >::iterator2d iterator2d;

  typedef T vector3d_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,2> vector3d_iterator;
  typedef slip::kstride_iterator<const_vector3d_pointer,2> const_vector3d_iterator;

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

  DenseVector3dField2d():
    base()
  {}
 
  DenseVector3dField2d(const size_type height,
                     const size_type width):
    base(height,width)
  {}
  
  DenseVector3dField2d(const size_type height,
                     const size_type width, 
                     const slip::Vector3d<T>& val):
    base(height,width,val)
  {}
 
  DenseVector3dField2d(const size_type height,
                     const size_type width, 
                     const T* val):
    base(height,width,val)
  {}
 
  DenseVector3dField2d(const size_type height,
                                 const size_type width, 
                     const slip::Vector3d<T>* val):
    base(height,width,val)
  {}
 

  template<typename InputIterator>
  DenseVector3dField2d(const size_type height,
                                       const size_type width, 
                                       InputIterator first,
                                       InputIterator last):
    base(height,width,first,last)
  {} 



  
  template<typename InputIterator>
  DenseVector3dField2d(const size_type height,
                                       const size_type width, 
                                       InputIterator first1,
                                       InputIterator last1,
                                       InputIterator first2,
                                       InputIterator first3):
    base(height,width)
  {
    
    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);
  
  } 


  DenseVector3dField2d(const self& rhs):
    base(rhs)
  {}

  ~DenseVector3dField2d()
  {}
  

  std::string name() const;


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


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




 

 

  self& operator=(const slip::Vector3d<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& Vx3(const size_type i,
                 const size_type j);
 
  const T& Vx3(const size_type i,
                       const size_type j) 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::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);
    

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

  typedef slip::DenseVector3dField2d<double> DenseVector3dField2d_d;
  typedef slip::DenseVector3dField2d<float> DenseVector3dField2d_f;
  typedef slip::DenseVector3dField2d<long> DenseVector3dField2d_l;
  typedef slip::DenseVector3dField2d<unsigned long> DenseVector3dField2d_ul;
  typedef slip::DenseVector3dField2d<short> DenseVector3dField2d_s;
  typedef slip::DenseVector3dField2d<unsigned short> DenseVector3dField2d_us;
  typedef slip::DenseVector3dField2d<int> DenseVector3dField2d_i;
  typedef slip::DenseVector3dField2d<unsigned int> DenseVector3dField2d_ui;
  typedef slip::DenseVector3dField2d<char> DenseVector3dField2d_c;
  typedef slip::DenseVector3dField2d<unsigned char> DenseVector3dField2d_uc;


}//slip::

namespace slip{
 

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

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


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

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



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

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


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


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

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

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

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

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

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

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


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

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


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

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

 


}//slip::

namespace slip
{


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

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


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

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

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


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


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

   template<typename T>
  inline
  const T& DenseVector3dField2d<T>::Vx3(const typename DenseVector3dField2d::size_type i,
                                                                                const typename DenseVector3dField2d::size_type j) const
  {
    return (*this)[i][j][2];
  } 
 
  template<typename T>
  inline
  DenseVector3dField2d<T>& DenseVector3dField2d<T>::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>
  inline
  DenseVector3dField2d<T>& DenseVector3dField2d<T>::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>
  inline
  DenseVector3dField2d<T>& DenseVector3dField2d<T>::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>
  inline
  DenseVector3dField2d<T>& DenseVector3dField2d<T>::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>
  inline
  DenseVector3dField2d<T>& DenseVector3dField2d<T>::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>
  inline
  DenseVector3dField2d<T>& DenseVector3dField2d<T>::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>
  inline
  DenseVector3dField2d<T>& DenseVector3dField2d<T>::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>
  inline
  DenseVector3dField2d<T>& DenseVector3dField2d<T>::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>
  inline
  DenseVector3dField2d<T> DenseVector3dField2d<T>::operator-() const
  {
    DenseVector3dField2d<T> tmp(*this);
    std::transform(this->begin(),this->end(),tmp.begin(),std::negate<slip::Vector3d<T> >());
    return tmp;
  }

  template<typename T>
  inline
  DenseVector3dField2d<T>& DenseVector3dField2d<T>::operator+=(const DenseVector3dField2d<T>& rhs)
  {
    assert(this->dim1() == rhs.dim1());
    assert(this->dim2() == rhs.dim2()); 
    std::transform(this->begin(),this->end(),rhs.begin(),this->begin(),std::plus<typename DenseVector3dField2d::value_type>());
    return *this;
  }

   template<typename T>
  inline
  DenseVector3dField2d<T>& DenseVector3dField2d<T>::operator-=(const DenseVector3dField2d<T>& rhs)
  {
    assert(this->dim1() == rhs.dim1());
    assert(this->dim2() == rhs.dim2());
    std::transform(this->begin(),this->end(),rhs.begin(),this->begin(),std::minus<typename DenseVector3dField2d::value_type>());
    return *this;
  }

  template<typename T>
  inline
  DenseVector3dField2d<T>& DenseVector3dField2d<T>::operator*=(const DenseVector3dField2d<T>& rhs)
  {
    assert(this->dim1() == rhs.dim1());
    assert(this->dim2() == rhs.dim2());
    std::transform(this->begin(),this->end(),rhs.begin(),this->begin(),std::multiplies<typename DenseVector3dField2d::value_type>());
    return *this;
  }

  template<typename T>
  inline
  DenseVector3dField2d<T>& DenseVector3dField2d<T>::operator/=(const DenseVector3dField2d<T>& rhs)
  {
    assert(this->dim1() == rhs.dim1());
    assert(this->dim2() == rhs.dim2());
    std::transform(this->begin(),this->end(),rhs.begin(),this->begin(),std::divides<typename DenseVector3dField2d::value_type>());
    return *this;
  }


  template<typename T>
  inline
  std::string 
  DenseVector3dField2d<T>::name() const {return "DenseVector3dField2d";} 
  
  template<typename T>
  inline
  void DenseVector3dField2d<T>::write_gnuplot(const std::string& file_path_name)
  {
    std::ofstream sortie(file_path_name.c_str(),std::ios::out);
    //Put an exception 
    if(!sortie)
    {
      std::cerr<<"Fail to open the file "<<file_path_name<<std::endl;
      sortie.close();
      exit(1);
    }
    
    typename DenseVector3dField2d::size_type dim1 = this->dim1();
    typename DenseVector3dField2d::size_type dim2 = this->dim2();
    
    sortie<<"#zone f=point, i="<<dim1<<", j="<<dim2<<"\n";
    for(typename DenseVector3dField2d::size_type i = 0; i < dim1; ++i)
      {
                for(typename DenseVector3dField2d::size_type j = 0; j < dim2; ++j)
                  {
                    sortie<<i<<" "<<j<<" "<<(*this)[i][j][0]<<" "<<(*this)[i][j][1]<<(*this)[i][j][2]<<"\n";
                  }
      }
  }

  template<typename T>
  inline
  void DenseVector3dField2d<T>::write_tecplot(const std::string& file_path_name,
                                                                                      const std::string& title,
                                                                                      const std::string& zone)
  {
     std::ofstream output(file_path_name.c_str(),std::ios::out);
    //Put an exception 
    if(!output)
    {
      std::cerr<<"Fail to open the file "<<file_path_name<<std::endl;
      output.close();
      exit(1);
    }
    
    typename DenseVector3dField2d::size_type dim1 = this->dim1();
    typename DenseVector3dField2d::size_type dim2 = this->dim2();

    //write file informations
    output<<"TITLE=\"" <<title<<"\""<<std::endl;
    output<<"VARIABLES=\"X\",\"Y\",\"U\",\"V\"" <<std::endl;
    output<<"ZONE T=\""<<zone<<"\", I="<<dim1<<", J="<<dim2<<std::endl;
    //write data
    for(typename DenseVector3dField2d::size_type i = 0; i < dim1; ++i)
      {
                for(typename DenseVector3dField2d::size_type j = 0; j < dim2; ++j)
                  {
                    output<<i<<" "<<j<<" "<<(*this)[i][j][0]<<" "<<(*this)[i][j][1]<<(*this)[i][j][2]<<"\n";
                  }
      }
  }

  /* @{ */
  template<typename T>
  inline
  DenseVector3dField2d<T> operator+(const DenseVector3dField2d<T>& M1, 
                                                                    const T& val)
  {
    DenseVector3dField2d<T>  tmp = M1;
    tmp += val;
    return tmp;
  }

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

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

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


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

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


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

  template<typename T>
  inline
  DenseVector3dField2d<T> operator+(const DenseVector3dField2d<T>& M1, 
                                                                    const DenseVector3dField2d<T>& M2)
  {
    assert(M1.dim1() == M2.dim1());
    assert(M1.dim2() == M2.dim2());
    DenseVector3dField2d<T> tmp(M1.dim1(),M1.dim2());
    std::transform(M1.begin(),M1.end(),M2.begin(),tmp.begin(),std::plus<slip::Vector3d<T> >());
    return tmp;
  }


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

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



 template<typename T>
  inline
  DenseVector3dField2d<T> operator-(const DenseVector3dField2d<T>& M1, 
                                                                    const DenseVector3dField2d<T>& M2)
  {
    assert(M1.dim1() == M2.dim1());
    assert(M1.dim2() == M2.dim2());
    DenseVector3dField2d<T> tmp(M1.dim1(),M1.dim2());
    std::transform(M1.begin(),M1.end(),M2.begin(),tmp.begin(),std::minus<slip::Vector3d<T> >());
    return tmp;
  }


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

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


  template<typename T>
  inline
  DenseVector3dField2d<T> operator*(const DenseVector3dField2d<T>& M1, 
                                                                    const DenseVector3dField2d<T>& M2)
  {
    assert(M1.dim1() == M2.dim1());
    assert(M1.dim2() == M2.dim2());
    DenseVector3dField2d<T> tmp(M1.dim1(),M1.dim2());
    std::transform(M1.begin(),M1.end(),M2.begin(),tmp.begin(),std::multiplies<slip::Vector3d<T> >());
    return tmp;
  }


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

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


 template<typename T>
  inline
  DenseVector3dField2d<T> operator/(const DenseVector3dField2d<T>& M1, 
                                                                    const DenseVector3dField2d<T>& M2)
  {
    assert(M1.dim1() == M2.dim1());
    assert(M1.dim2() == M2.dim2());
    DenseVector3dField2d<T> tmp(M1.dim1(),M1.dim2());
    std::transform(M1.begin(),M1.end(),M2.begin(),tmp.begin(),std::divides<slip::Vector3d<T> >());
    return tmp;
  }


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


/* @} */



}//slip::

#endif //SLIP_DENSEVECTOR3DFIELD2D_HPP