diff --git a/src/common/aka_types.hh b/src/common/aka_types.hh
index d73dcb9e2..96f4b5a0b 100644
--- a/src/common/aka_types.hh
+++ b/src/common/aka_types.hh
@@ -1,977 +1,987 @@
 /**
  * @file   aka_types.hh
  *
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Thu Feb 17 2011
  * @date last modification: Tue Aug 19 2014
  *
  * @brief  description of the "simple" types
  *
  * @section LICENSE
  *
  * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  * Akantu is free  software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as  published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Akantu is  distributed in the  hope that it  will be useful, but  WITHOUT ANY
  * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
  * A  PARTICULAR PURPOSE. See  the GNU  Lesser General  Public License  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 #include "aka_error.hh"
 #include "aka_fwd.hh"
 #include "aka_math.hh"
 #include "aka_array.hh"
 
 /* -------------------------------------------------------------------------- */
 #include <iomanip>
 
 #ifndef __INTEL_COMPILER
 #include <tr1/unordered_map>
 #else
 #include <map>
 #endif
 
 /* -------------------------------------------------------------------------- */
 #ifndef __AKANTU_AKA_TYPES_HH__
 #define __AKANTU_AKA_TYPES_HH__
 
 __BEGIN_AKANTU__
 
 /* -------------------------------------------------------------------------- */
 /* maps                                                                       */
 /* -------------------------------------------------------------------------- */
 #ifndef __INTEL_COMPILER
 template<class Key, class Ty>
 struct unordered_map { typedef typename std::tr1::unordered_map<Key, Ty> type; };
 #else
 template<class Key, class Ty>
 struct unordered_map { typedef typename std::map<Key, Ty> type; };
 #endif
 
 enum NormType {
   L_1 = 1,
   L_2 = 2,
   L_inf = UInt(-1)
 };
 
 
 template<UInt dim>
 struct DimHelper {
   static inline void setDims(UInt m, UInt n, UInt p, UInt dims[dim]);
 };
 
 template<>
 struct DimHelper<1> {
   static inline void setDims(UInt m,
 			     __attribute__((unused)) UInt n,
 			     __attribute__((unused)) UInt p,
 			     UInt dims[1]) {
     dims[0] = m;
   }
 };
 
 template<>
 struct DimHelper<2> {
   static inline void setDims(UInt m,
 			     UInt n,
 			     __attribute__((unused)) UInt p,
 			     UInt dims[2]) {
     dims[0] = m;
     dims[1] = n;
   }
 };
 
 template<>
 struct DimHelper<3> {
   static inline void setDims(UInt m,
 			     UInt n,
 			     UInt p,
 			     UInt dims[3]) {
     dims[0] = m;
     dims[1] = n;
     dims[2] = p;
   }
 };
 
 
 /* -------------------------------------------------------------------------- */
 template<typename T, UInt ndim>
 class TensorProxy {
 protected:
   TensorProxy(T * data, UInt m, UInt n, UInt p) {
     DimHelper<ndim>::setDims(m, n, p, this->n);
     this->values = data;
   }
 
   TensorProxy(const TensorProxy & other) {
     this->values = other.storage();
     for (UInt i = 0; i < ndim; ++i)
       this->n[i] = other.n[i];
   }
 
 public:
   UInt size(UInt i) const {
     AKANTU_DEBUG_ASSERT(i < ndim,
 			"This tensor has only "
 			<< ndim << " dimensions, not " << (i + 1));
     return n[i];
   }
 
   inline UInt size() const{
     UInt _size = 1;
     for (UInt d = 0; d < ndim; ++d) _size *= this->n[d];
     return _size;
   }
 
   T * storage() const { return values; }
 
 //   TensorProxy operator=(const TensorProxy & other) {
 //     UInt _size = this->size();
 // #ifndef AKANTU_NDEBUG
 //     UInt _size_other = other.size();
 //     AKANTU_DEBUG_ASSERT(_size == _size_other, "The two tensor are not compatible in size");
 // #endif
 //     memcpy(this->values, other.storage(), _size * sizeof(T));
 //     return *this;
 //   }
 
 protected:
   T * values;
   UInt n[ndim];
 };
 
 
 /* -------------------------------------------------------------------------- */
 template<typename T>
 class VectorProxy : public TensorProxy<T, 1> {
   typedef TensorProxy<T, 1> parent;
 public:
   VectorProxy(T * data = NULL, UInt n = 0) : parent(data, n, 0, 0) { }
   VectorProxy(const VectorProxy & src) : parent(src) {  }
   // VectorProxy & operator=(const VectorProxy & other) {
   //   parent::operator=(other);
   //   return *this;
   // }
 };
 
 template<typename T>
 class MatrixProxy : public TensorProxy<T, 2> {
   typedef TensorProxy<T, 2> parent;
 public:
   MatrixProxy(T * data = NULL, UInt m = 0, UInt n = 0) : parent(data, m, n, 0) { }
   MatrixProxy(const MatrixProxy & src) : parent(src) {  }
   // MatrixProxy & operator=(const MatrixProxy & other) {
   //   parent::operator=(other);
   //   return *this;
   // }
 
 };
 
 template<typename T>
 class Tensor3Proxy : public TensorProxy<T, 3> {
   typedef TensorProxy<T, 3> parent;
 public:
   Tensor3Proxy(T * data = NULL, UInt m = 0, UInt n = 0, UInt k = 0) :
     parent(data, m, n, k) { }
   Tensor3Proxy(const Tensor3Proxy & src) : parent(src) {  }
   // Tensor3Proxy & operator=(const Tensor3Proxy & other) {
   //   parent::operator=(other);
   //   return *this;
   // }
 };
 
 /* -------------------------------------------------------------------------- */
 /* Tensor base class                                                          */
 /* -------------------------------------------------------------------------- */
 template<typename T, UInt ndim, class RetType>
 class TensorStorage {
 public:
   typedef T value_type;
 protected:
   template<class TensorType>
   void copySize(const TensorType & src) {
     for (UInt d = 0; d < ndim; ++d) this->n[d] = src.size(d);
     this->_size = src.size();
   }
 
   TensorStorage() :
     values(NULL), wrapped(false) {
     for (UInt d = 0; d < ndim; ++d) this->n[d] = 0;
     _size = 0;
   }
 
   TensorStorage(const TensorProxy<T, ndim> & proxy) {
     this->copySize(proxy);
     this->values = proxy.storage();
     this->wrapped = true;
   }
 
 protected:
   TensorStorage(const TensorStorage & src) { }
 
 public:
 
   TensorStorage(const TensorStorage & src, bool deep_copy) :
     values(NULL), wrapped(false) {
     if(deep_copy) this->deepCopy(src);
     else this->shallowCopy(src);
   }
 
 protected:
   TensorStorage(UInt m, UInt n, UInt p, const T & def) {
     DimHelper<ndim>::setDims(m, n, p, this->n);
 
     this->computeSize();
     this->values = new T[this->_size];
     this->set(def);
     this->wrapped = false;
   }
 
   TensorStorage(T * data, UInt m, UInt n, UInt p) {
     DimHelper<ndim>::setDims(m, n, p, this->n);
 
     this->computeSize();
     this->values = data;
     this->wrapped = true;
   }
 
 public:
   /* ------------------------------------------------------------------------ */
   template<class TensorType>
   inline void shallowCopy(const TensorType & src) {
     this->copySize(src);
     if(!this->wrapped) delete[] this->values;
     this->values = src.storage();
     this->wrapped = true;
   }
 
   /* ------------------------------------------------------------------------ */
   template<class TensorType>
   inline void deepCopy(const TensorType & src) {
     this->copySize(src);
     if(!this->wrapped) delete [] this->values;
     this->values = new T[this->_size];
     memcpy(this->values, src.storage(), this->_size * sizeof(T));
     this->wrapped = false;
   }
 
   virtual ~TensorStorage() {
     if(!this->wrapped)
       delete [] this->values;
   }
 
   inline TensorStorage & operator=(const RetType & src) {
     if(this != &src) {
       if (this->wrapped) {
 	AKANTU_DEBUG_ASSERT(this->_size == src.size(), "vectors of different size");
 	memcpy(this->values, src.storage(), this->_size * sizeof(T));
       } else {
 	deepCopy(src);
       }
     }
     return *this;
   }
 
   /* ------------------------------------------------------------------------ */
   template<class R>
   inline RetType & operator+=(const TensorStorage<T, ndim, R> & other) {
     T * a = this->storage();
     T * b = other.storage();
     AKANTU_DEBUG_ASSERT(_size == other.size(),
 			"The two tensors do not have the same size, they cannot be subtracted");
     for (UInt i = 0; i < _size; ++i) *(a++) += *(b++);
     return *(static_cast<RetType *>(this));
   }
 
   /* ------------------------------------------------------------------------ */
   template<class R>
   inline RetType & operator-=(const TensorStorage<T, ndim, R> & other) {
     T * a = this->storage();
     T * b = other.storage();
     AKANTU_DEBUG_ASSERT(_size == other.size(),
 			"The two tensors do not have the same size, they cannot be subtracted");
     for (UInt i = 0; i < _size; ++i) *(a++) -= *(b++);
     return *(static_cast<RetType *>(this));
   }
 
   /* ------------------------------------------------------------------------ */
   inline RetType & operator+=(const T & x) {
     T * a = this->values;
     for (UInt i = 0; i < _size; ++i) *(a++) += x;
     return *(static_cast<RetType *>(this));
   }
 
   /* ------------------------------------------------------------------------ */
   inline RetType & operator-=(const T & x) {
     T * a = this->values;
     for (UInt i = 0; i < _size; ++i) *(a++) -= x;
     return *(static_cast<RetType *>(this));
   }
 
   /* ------------------------------------------------------------------------ */
   inline RetType & operator*=(const T & x) {
     T * a = this->storage();
     for (UInt i = 0; i < _size; ++i) *(a++) *= x;
     return *(static_cast<RetType *>(this));
   }
 
   /* ---------------------------------------------------------------------- */
   inline RetType & operator/=(const T & x) {
     T * a = this->values;
     for (UInt i = 0; i < _size; ++i) *(a++) /= x;
     return *(static_cast<RetType *>(this));
   }
 
   /* ------------------------------------------------------------------------ */
   T * storage() const { return values; }
   UInt size() const { return _size; }
   UInt size(UInt i) const {
     AKANTU_DEBUG_ASSERT(i < ndim,
 			"This tensor has only "
 			<< ndim << " dimensions, not " << (i + 1));
     return n[i];
   };
   /* ------------------------------------------------------------------------ */
   inline void clear() { memset(values, 0, _size * sizeof(T)); };
   inline void set(const T & t) { std::fill_n(values, _size, t); };
 
 
   template<class TensorType>
   inline void copy(const TensorType & other) {
     AKANTU_DEBUG_ASSERT(_size == other.size(),
 			"The two tensors do not have the same size, they cannot be copied");
     memcpy(values, other.storage(), _size * sizeof(T));
   }
 
 protected:
   friend class Array<T>;
 
   inline void computeSize() {
     _size = 1;
     for (UInt d = 0; d < ndim; ++d) _size *= this->n[d];
   }
 
 protected:
   template<typename R, NormType norm_type>
   struct NormHelper {
     template<class Ten>
     static R norm(const Ten & ten) {
       R _norm = 0.;
       R * it = ten.storage();
       R * end = ten.storage() + ten.size();
       for (; it < end; ++it) _norm += std::pow(std::abs(*it), norm_type);
       return std::pow(_norm, 1./norm_type);
     }
   };
 
   template<typename R>
   struct NormHelper<R, L_1> {
     template<class Ten>
     static R norm(const Ten & ten) {
       R _norm = 0.;
       R * it = ten.storage();
       R * end = ten.storage() + ten.size();
       for (; it < end; ++it) _norm += std::abs(*it);
       return _norm;
     }
   };
 
   template<typename R>
   struct NormHelper<R, L_2> {
     template<class Ten>
     static R norm(const Ten & ten) {
       R _norm = 0.;
       R * it = ten.storage();
       R * end = ten.storage() + ten.size();
       for (; it < end; ++it) _norm += *it * *it;
       return sqrt(_norm);
     }
   };
 
   template<typename R>
   struct NormHelper<R, L_inf> {
     template<class Ten>
     static R norm(const Ten & ten) {
       R _norm = 0.;
       R * it = ten.storage();
       R * end = ten.storage() + ten.size();
       for (; it < end; ++it) _norm = std::max(std::abs(*it), _norm);
       return _norm;
     }
   };
 
 public:
 
   /*----------------------------------------------------------------------- */
   /// "Entrywise" norm norm<L_p> @f[ \|\boldsymbol{T}\|_p = \left(
   /// \sum_i^{n[0]}\sum_j^{n[1]}\sum_k^{n[2]} |T_{ijk}|^p \right)^{\frac{1}{p}}
   /// @f]
   template<NormType norm_type>
   inline T norm() const { return NormHelper<T, norm_type>::norm(*this); }
 
 protected:
   UInt n[ndim];
   UInt _size;
   T * values;
   bool wrapped;
 };
 
 
 /* -------------------------------------------------------------------------- */
 /* Vector                                                                     */
 /* -------------------------------------------------------------------------- */
 template<typename T>
 class Vector : public TensorStorage< T, 1, Vector<T> > {
   typedef TensorStorage< T, 1, Vector<T> > parent;
 public:
   typedef typename parent::value_type value_type;
 public:
   Vector() : parent() {}
   Vector(UInt n, const T & def = T()) : parent(n, 0, 0, def) { }
   Vector(T * data, UInt n) : parent(data, n, 0, 0) { }
   Vector(const Vector & src, bool deep_copy = true) : parent(src, deep_copy) { }
   Vector(const VectorProxy<T> & src) : parent(src) { }
 
 public:
   virtual ~Vector() { };
 
   /* ------------------------------------------------------------------------ */
   inline Vector & operator=(const Vector & src) {
     parent::operator=(src);
     return *this;
   }
 
   /* ------------------------------------------------------------------------ */
   inline T& operator()(UInt i) { return *(this->values + i); };
   inline const T& operator()(UInt i) const { return *(this->values + i); };
   inline T& operator[](UInt i) { return *(this->values + i); };
   inline const T& operator[](UInt i) const { return *(this->values + i); };
 
   /* ------------------------------------------------------------------------ */
   inline Vector<T> & operator*=(Real x) { return parent::operator*=(x); }
   inline Vector<T> & operator/=(Real x) { return parent::operator/=(x); }
   /* ------------------------------------------------------------------------ */
   inline Vector<T> & operator*=(const Vector<T> & vect) {
     T * a = this->storage();
     T * b = vect.storage();
     for (UInt i = 0; i < this->_size; ++i) *(a++) *= *(b++);
     return *this;
   }
 
   /* ------------------------------------------------------------------------ */
   inline Real dot(const Vector<T> & vect) const {
     return Math::vectorDot(this->values, vect.storage(), this->_size);
   }
 
   /* ------------------------------------------------------------------------ */
   inline Vector & crossProduct(const Vector<T> & v1,
 			       const Vector<T> & v2) {
     AKANTU_DEBUG_ASSERT(this->size() == 3,
 			"crossProduct is only defined in 3D (n=" << this->size() << ")");
     AKANTU_DEBUG_ASSERT(this->size() == v1.size() && this->size() == v2.size(),
 			"crossProduct is not a valid operation non matching size vectors");
     Math::vectorProduct3(v1.storage(), v2.storage(), this->values);
     return *this;
   }
 
   /* ------------------------------------------------------------------------ */
   inline void solve(Matrix<T> & A, const Vector<T> & b) {
     AKANTU_DEBUG_ASSERT(this->size() == A.rows() && this->_size = A.cols(),
 			"The solution vector as a mismatch in size with the matrix");
     AKANTU_DEBUG_ASSERT(this->_size == b._size, "The rhs vector as a mismatch in size with the matrix");
     Math::solve(this->_size, A.storage(), this->values, b.storage());
   }
 
   /* ------------------------------------------------------------------------ */
   template<bool tr_A>
   inline void mul(const Matrix<T> & A,
 		  const Vector<T> & x,
 		  Real alpha = 1.0);
   /* ------------------------------------------------------------------------ */
   inline Real norm() const {
     return parent::template norm<L_2>();
   }
 
   template<NormType nt>
   inline Real norm() const {
     return parent::template norm<nt>();
   }
 
   /* ------------------------------------------------------------------------ */
   inline void normalize() {
     Real n = norm();
     operator/=(n);
   }
 
   /* ------------------------------------------------------------------------ */
   /// norm of (*this - x)
   inline Real distance(const Vector<T> & y) const {
     Real * vx = this->values; Real * vy = y.storage();
     Real sum_2 = 0;
     for (UInt i = 0; i < this->_size; ++i, ++vx, ++vy) sum_2 += (*vx - *vy)*(*vx - *vy);
     return sqrt(sum_2);
   }
 
   /* ------------------------------------------------------------------------ */
   inline bool equal(const Vector<T> & v, Real tolerance = Math::getTolerance()) const {
     T * a = this->storage();
     T * b = v.storage();
     UInt i = 0;
     while (i < this->_size && (std::abs(*(a++) - *(b++)) < tolerance)) ++i;
     return i == this->_size;
   }
 
   /* ------------------------------------------------------------------------ */
   inline short compare(const Vector<T> & v, Real tolerance = Math::getTolerance()) const {
     T * a = this->storage();
     T * b = v.storage();
     for (UInt i(0); i < this->_size; ++i, ++a, ++b) {
       if(std::abs(*a - *b) > tolerance)
         return (((*a - *b) > tolerance) ? 1 : -1);
     }
     return 0;
   }
 
   /* ------------------------------------------------------------------------ */
   inline bool operator==(const Vector<T> & v) const { return equal(v); }
   inline bool operator<(const Vector<T> & v) const { return compare(v) == -1; }
   inline bool operator>(const Vector<T> & v) const { return compare(v) == 1; }
 
   /* ------------------------------------------------------------------------ */
   /// function to print the containt of the class
   virtual void printself(std::ostream & stream, int indent = 0) const {
     std::string space;
     for(Int i = 0; i < indent; i++, space += AKANTU_INDENT);
 
     stream << space << "Vector<" << debug::demangle(typeid(T).name()) << ">(" << this->_size <<") : [";
     for (UInt i = 0; i < this->_size; ++i) {
       if(i != 0) stream << ", ";
       stream << this->values[i];
     }
     stream << "]";
   }
 
   friend class ::akantu::Array<T>;
 };
 
 typedef Vector<Real> RVector;
 
 /* -------------------------------------------------------------------------- */
 // support operations for the creation of other vectors
 template <typename T> Vector<T> operator*(T scalar, const Vector<T> & a);
 template <typename T> Vector<T> operator+(const Vector<T> & a, const Vector<T> & b);
 template <typename T> Vector<T> operator-(const Vector<T> & a, const Vector<T> & b);
 
 /* -------------------------------------------------------------------------- */
 template <typename T>
 Vector<T> operator*(T scalar, const Vector<T> & a) {
   Vector<T> r(a.size());
   r = a;
   r *= scalar;
   return r;
 }
 
 template <typename T>
 Vector<T> operator+(const Vector<T> & a, const Vector<T> & b) {
   Vector<T> r(a.size());
   r = a;
   r += b;
   return r;
 }
 
 template <typename T>
 Vector<T> operator-(const Vector<T>& a, const Vector<T>& b) {
   Vector<T> r(a.size());
   r = a;
   r -= b;
   return r;
 }
 
 template <typename T>
 Matrix<T> operator*(T scalar, const Matrix<T>& a) {
   Matrix<T> r(a.rows(), a.cols());
   r = a;
   r *= scalar;
   return r;
 }
 
 template <typename T>
 Matrix<T> operator+(const Matrix<T>& a, const Matrix<T>& b) {
   Matrix<T> r(a.rows(), a.cols());
   r = a;
   r += b;
   return r;
 }
 
 template <typename T>
 Matrix<T> operator-(const Matrix<T>& a, const Matrix<T>& b) {
   Matrix<T> r(a.rows(), a.cols());
   r = a;
   r -= b;
   return r;
 }
 
 /* ------------------------------------------------------------------------ */
 template<>
 inline bool Vector<UInt>::equal(const Vector<UInt> & v, __attribute__((unused)) Real tolerance) const {
   UInt * a = this->storage();
   UInt * b = v.storage();
   UInt i = 0;
   while (i < this->_size && (*(a++) == *(b++))) ++i;
   return i == this->_size;
 }
 
 
 /* ------------------------------------------------------------------------ */
 /* Matrix                                                                   */
 /* ------------------------------------------------------------------------ */
 template<typename T>
 class Matrix : public TensorStorage< T, 2, Matrix<T> > {
   typedef TensorStorage< T, 2, Matrix<T> > parent;
 public:
   typedef typename parent::value_type value_type;
 public:
   Matrix() : parent() {}
   Matrix(UInt m, UInt n, const T & def = T()) : parent(m, n, 0, def) { }
   Matrix(T * data, UInt m, UInt n) : parent(data, m, n, 0) { }
   Matrix(const Matrix & src, bool deep_copy = true) : parent(src, deep_copy) { }
   Matrix(const MatrixProxy<T> & src) : parent(src) { }
 
   virtual ~Matrix() { }
   /* ------------------------------------------------------------------------ */
   inline Matrix & operator=(const Matrix & src) {
     parent::operator=(src);
     return *this;
   }
 
 public:
   /* ---------------------------------------------------------------------- */
   UInt rows() const { return this->n[0]; }
   UInt cols() const { return this->n[1]; }
 
   /* ---------------------------------------------------------------------- */
   inline T& operator()(UInt i, UInt j) {
+    return this->at(i,j);
+  }
+
+  inline T& at(UInt i, UInt j) {
     return *(this->values + i + j*this->n[0]);
   }
+  
   inline const T& operator()(UInt i, UInt j) const {
+    return this->at(i,j);
+  }
+
+  inline T& at(UInt i, UInt j) const {
     return *(this->values + i + j*this->n[0]);
   }
 
+  
   /// give a line vector wrapped on the column i
   inline VectorProxy<T> operator()(UInt j) {
     return VectorProxy<T>(this->values + j*this->n[0], this->n[0]);
   }
   inline const VectorProxy<T> operator()(UInt j) const {
     return VectorProxy<T>(this->values + j*this->n[0], this->n[0]);
   }
 
   inline T& operator[](UInt idx) { return *(this->values + idx); };
   inline const T& operator[](UInt idx) const { return *(this->values + idx); };
 
   /* ---------------------------------------------------------------------- */
   inline Matrix operator* (const Matrix & B) {
     Matrix C(this->rows(), B.cols());
     C.mul<false, false>(*this, B);
     return C;
   }
 
   /* ----------------------------------------------------------------------- */
   inline Matrix & operator*= (const T & x) {
     return parent::operator*= (x);
   }
 
   inline Matrix & operator*= (const Matrix & B) {
     Matrix C(*this);
     this->mul<false, false>(C, B);
     return *this;
   }
 
   /* ---------------------------------------------------------------------- */
   template<bool tr_A, bool tr_B>
   inline void mul(const Matrix & A, const Matrix & B, T alpha = 1.0) {
     UInt k = A.cols();
     if(tr_A) k = A.rows();
 
 #ifndef AKANTU_NDEBUG
     if (tr_B){
       AKANTU_DEBUG_ASSERT(k == B.cols(),
 			  "matrices to multiply have no fit dimensions");
       AKANTU_DEBUG_ASSERT(this->cols() == B.rows(),
 			  "matrices to multiply have no fit dimensions");
     }
     else {
       AKANTU_DEBUG_ASSERT(k == B.rows(),
 			  "matrices to multiply have no fit dimensions");
       AKANTU_DEBUG_ASSERT(this->cols() == B.cols(),
 			  "matrices to multiply have no fit dimensions");
     }
     if (tr_A){
       AKANTU_DEBUG_ASSERT(this->rows() == A.cols(),
 			  "matrices to multiply have no fit dimensions");
     }
     else{
       AKANTU_DEBUG_ASSERT(this->rows() == A.rows(),
 			  "matrices to multiply have no fit dimensions");
     }
 #endif //AKANTU_NDEBUG
 
     Math::matMul<tr_A, tr_B>(this->rows(), this->cols(), k,
 			     alpha, A.storage(), B.storage(),
 			     0., this->storage());
   }
 
   /* ---------------------------------------------------------------------- */
   inline void outerProduct(const Vector<T> & A,
 			   const Vector<T> & B) {
     AKANTU_DEBUG_ASSERT(A.size() == this->rows() && B.size() == this->cols(),
 			"A and B are not compatible with the size of the matrix");
     for (UInt i = 0; i < this->rows(); ++i) {
       for (UInt j = 0; j < this->cols(); ++j) {
 	this->values[i + j * this->rows()] += A[i] * B[j];
       }
     }
   }
 
 private:
   class EigenSorter {
   public:
     EigenSorter(const Vector<T> & eigs) : eigs(eigs) { }
 
     bool operator()(const UInt & a, const UInt & b) const {
       return (eigs(a) > eigs(b));
     }
 
   private:
     const Vector<T> & eigs;
   };
 
 public:
   /* ---------------------------------------------------------------------- */
   inline void eig(Vector<T> & eigenvalues, Matrix<T> & eigenvectors) const {
     AKANTU_DEBUG_ASSERT(this->cols() == this->rows(),
 			"eig is not a valid operation on a rectangular matrix");
     AKANTU_DEBUG_ASSERT(eigenvalues.size() == this->cols(),
 			"eigenvalues should be of size " << this->cols() << ".");
 #ifndef AKANTU_NDEBUG
     if(eigenvectors.storage() != NULL)
       AKANTU_DEBUG_ASSERT((eigenvectors.rows() == eigenvectors.cols()) &&
 			  (eigenvectors.rows() == this->cols()),
 			  "Eigenvectors needs to be a square matrix of size "
 			  << this->cols() << " x " << this->cols() << ".");
 #endif
 
     Matrix<T> tmp = *this;
     Vector<T> tmp_eigs(eigenvalues.size());
     Matrix<T> tmp_eig_vects(eigenvectors.rows(), eigenvectors.cols());
 
     if(tmp_eig_vects.rows() == 0 || tmp_eig_vects.cols() == 0)
       Math::matrixEig(tmp.cols(), tmp.storage(), tmp_eigs.storage());
     else
       Math::matrixEig(tmp.cols(), tmp.storage(), tmp_eigs.storage(), tmp_eig_vects.storage());
 
     Vector<UInt> perm(eigenvalues.size());
     for (UInt i = 0; i < perm.size(); ++i) perm(i) = i;
 
     std::sort(perm.storage(), perm.storage() + perm.size(), EigenSorter(tmp_eigs));
 
     for (UInt i = 0; i < perm.size(); ++i) eigenvalues(i) = tmp_eigs(perm(i));
 
     if(tmp_eig_vects.rows() != 0 && tmp_eig_vects.cols() != 0)
       for (UInt i = 0; i < perm.size(); ++i) {
 	for (UInt j = 0; j < eigenvectors.rows(); ++j) {
 	  eigenvectors(j, i) = tmp_eig_vects(j, perm(i));
 	}
       }
   }
 
   /* ---------------------------------------------------------------------- */
   inline void eig(Vector<T> & eigenvalues) const {
     Matrix<T> empty;
     eig(eigenvalues, empty);
   }
 
   /* ---------------------------------------------------------------------- */
   inline void eye(T alpha = 1.) {
     AKANTU_DEBUG_ASSERT(this->cols() == this->rows(), "eye is not a valid operation on a rectangular matrix");
     this->clear();
     for (UInt i = 0; i < this->cols(); ++i) {
       this->values[i + i * this->rows()] = alpha;
     }
   }
 
   /* ---------------------------------------------------------------------- */
   static inline Matrix<T> eye(UInt m, T alpha = 1.) {
     Matrix<T> tmp(m, m);
     tmp.eye(alpha);
     return tmp;
   }
 
   /* ---------------------------------------------------------------------- */
   inline T trace() const {
     AKANTU_DEBUG_ASSERT(this->cols() == this->rows(), "trace is not a valid operation on a rectangular matrix");
     T trace = 0.;
     for (UInt i = 0; i < this->rows(); ++i) {
       trace += this->values[i + i * this->rows()];
     }
     return trace;
   }
 
   /* ---------------------------------------------------------------------- */
   inline Matrix transpose() const {
     Matrix tmp(this->cols(), this->rows());
     for (UInt i = 0; i < this->rows(); ++i) {
       for (UInt j = 0; j < this->cols(); ++j) {
 	tmp(j,i) = operator()(i, j);
       }
     }
     return tmp;
   }
 
   /* ---------------------------------------------------------------------- */
   inline void inverse(const Matrix & A) {
     AKANTU_DEBUG_ASSERT(A.cols() == A.rows(),
 			"inv is not a valid operation on a rectangular matrix");
     AKANTU_DEBUG_ASSERT(this->cols() == A.cols(),
 			"the matrix should have the same size as its inverse");
 
     if(this->cols() == 1) *this->values = 1./ *A.storage();
     else if(this->cols() == 2) Math::inv2(A.storage(), this->values);
     else if(this->cols() == 3) Math::inv3(A.storage(), this->values);
     else Math::inv(this->cols(), A.storage(), this->values);
   }
 
   /* --------------------------------------------------------------------- */
   inline T det() const {
     AKANTU_DEBUG_ASSERT(this->cols() == this->rows(),
 			"inv is not a valid operation on a rectangular matrix");
     if(this->cols() == 1) return *(this->values);
     else if(this->cols() == 2) return Math::det2(this->values);
     else if(this->cols() == 3) return Math::det3(this->values);
     else return Math::det(this->cols(), this->values);
   }
 
   /* --------------------------------------------------------------------- */
   inline T doubleDot(const Matrix<T> & other) const {
      AKANTU_DEBUG_ASSERT(this->cols() == this->rows(),
 			 "doubleDot is not a valid operation on a rectangular matrix");
      if(this->cols() == 1) return *(this->values) * *(other.storage());
      else if(this->cols() == 2) return Math::matrixDoubleDot22(this->values, other.storage());
      else if(this->cols() == 3) return Math::matrixDoubleDot33(this->values, other.storage());
      else AKANTU_DEBUG_ERROR("doubleDot is not defined for other spatial dimensions"
                              << " than 1, 2 or 3.");
      return T();
   }
 
   /* ---------------------------------------------------------------------- */
   /// function to print the containt of the class
   virtual void printself(std::ostream & stream, int indent = 0) const {
     std::string space;
     for(Int i = 0; i < indent; i++, space += AKANTU_INDENT);
 
     stream << space << "Matrix<" << debug::demangle(typeid(T).name())
 	   << ">(" << this->n[0] << "," << this->n[1] <<") :" << "[";
     for (UInt i = 0; i < this->n[0]; ++i) {
       if(i != 0) stream << ", ";
       stream << "[";
       for (UInt j = 0; j < this->n[1]; ++j) {
 	if(j != 0) stream << ", ";
 	stream << operator()(i, j);
       }
       stream << "]";
     }
     stream << "]";
   };
 };
 
 /* ------------------------------------------------------------------------ */
 template<typename T>
 template<bool tr_A>
 inline void Vector<T>::mul(const Matrix<T> & A,
 			  const Vector<T> & x,
 			  Real alpha) {
 #ifndef AKANTU_NDEBUG
   UInt n = x.size();
   if (tr_A){
     AKANTU_DEBUG_ASSERT(n == A.rows(), "matrix and vector to multiply have no fit dimensions");
     AKANTU_DEBUG_ASSERT(this->size() == A.cols(), "matrix and vector to multiply have no fit dimensions");
   } else {
     AKANTU_DEBUG_ASSERT(n == A.cols(), "matrix and vector to multiply have no fit dimensions");
     AKANTU_DEBUG_ASSERT(this->size() == A.rows(), "matrix and vector to multiply have no fit dimensions");
   }
 #endif
   Math::matVectMul<tr_A>(A.rows(), A.cols(), alpha, A.storage(), x.storage(), 0., this->storage());
 }
 
 
 /* -------------------------------------------------------------------------- */
 template<typename T>
 inline std::ostream & operator<<(std::ostream & stream, const Matrix<T> & _this)
 {
   _this.printself(stream);
   return stream;
 }
 
 /* -------------------------------------------------------------------------- */
 template<typename T>
 inline std::ostream & operator<<(std::ostream & stream, const Vector<T> & _this)
 {
   _this.printself(stream);
   return stream;
 }
 
 /* ------------------------------------------------------------------------ */
 /* Tensor3                                                                  */
 /* ------------------------------------------------------------------------ */
 template<typename T>
 class Tensor3 : public TensorStorage< T, 3, Tensor3<T> > {
   typedef TensorStorage< T, 3, Tensor3<T> > parent;
 public:
   typedef typename parent::value_type value_type;
 public:
   Tensor3() : parent() { };
 
   Tensor3(UInt m, UInt n, UInt p, const T & def = T()) : parent(m, n, p, def) {  }
 
   Tensor3(T * data, UInt m, UInt n, UInt p) : parent(data, m, n, p) {  }
 
   Tensor3(const Tensor3 & src, bool deep_copy = true) : parent(src, deep_copy) { }
 
 public:
   /* ------------------------------------------------------------------------ */
   inline Tensor3 & operator=(const Tensor3 & src) {
     parent::operator=(src);
     return *this;
   }
 
   /* ---------------------------------------------------------------------- */
   inline T& operator()(UInt i, UInt j, UInt k)
   { return *(this->values + (k*this->n[0] + i)*this->n[1] + j); };
   inline const T& operator()(UInt i, UInt j, UInt k) const
   { return *(this->values + (k*this->n[0] + i)*this->n[1] + j); };
 
   inline MatrixProxy<T> operator()(UInt k)
   { return MatrixProxy<T>(this->values + k*this->n[0]*this->n[1], this->n[0], this->n[1]); }
   inline const MatrixProxy<T> operator()(UInt k) const
   { return MatrixProxy<T>(this->values + k*this->n[0]*this->n[1], this->n[0], this->n[1]); }
 
   inline MatrixProxy<T> operator[](UInt k)
   { return Matrix<T>(this->values + k*this->n[0]*this->n[1], this->n[0], this->n[1]); }
   inline const MatrixProxy<T> operator[](UInt k) const
   { return MatrixProxy<T>(this->values + k*this->n[0]*this->n[1], this->n[0], this->n[1]); }
 };
 
 __END_AKANTU__
 
 #endif /* __AKANTU_AKA_TYPES_HH__ */
diff --git a/src/model/solid_mechanics/materials/material_elastic_linear_anisotropic.cc b/src/model/solid_mechanics/materials/material_elastic_linear_anisotropic.cc
index ddd7ecfff..a9ce71605 100644
--- a/src/model/solid_mechanics/materials/material_elastic_linear_anisotropic.cc
+++ b/src/model/solid_mechanics/materials/material_elastic_linear_anisotropic.cc
@@ -1,331 +1,332 @@
 /**
  * @file   material_elastic_linear_anisotropic.cc
  *
  * @author Till Junge <till.junge@epfl.ch>
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Wed Sep 25 2013
  * @date last modification: Fri Sep 19 2014
  *
  * @brief  Anisotropic elastic material
  *
  * @section LICENSE
  *
  * Copyright (©) 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  * Akantu is free  software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as  published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Akantu is  distributed in the  hope that it  will be useful, but  WITHOUT ANY
  * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
  * A  PARTICULAR PURPOSE. See  the GNU  Lesser General  Public License  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 
 /* -------------------------------------------------------------------------- */
 #include "material_elastic_linear_anisotropic.hh"
 #include "solid_mechanics_model.hh"
 #include <algorithm>
 #include <sstream>
 
 
 __BEGIN_AKANTU__
 
 /* -------------------------------------------------------------------------- */
 template<UInt spatial_dimension>
 MaterialElasticLinearAnisotropic<spatial_dimension>::
 MaterialElasticLinearAnisotropic(SolidMechanicsModel & model,
                                  const ID & id,
                                  bool symmetric)  :
   Material(model, id),
   rot_mat(spatial_dimension, spatial_dimension),
   Cprime(spatial_dimension*spatial_dimension,
          spatial_dimension*spatial_dimension),
   C(this->voigt_h.size, this->voigt_h.size),
   eigC(this->voigt_h.size),
   symmetric(symmetric),
   alpha(0) {
   AKANTU_DEBUG_IN();
 
   this->dir_vecs.push_back(new Vector<Real>(spatial_dimension));
   (*this->dir_vecs.back())[0] = 1.;
   this->registerParam("n1", *(this->dir_vecs.back()), _pat_parsmod,
                       "Direction of main material axis");
   this->dir_vecs.push_back(new Vector<Real>(spatial_dimension));
   (*this->dir_vecs.back())[1] = 1.;
   this->registerParam("n2", *(this->dir_vecs.back()), _pat_parsmod,
                       "Direction of secondary material axis");
 
   if (spatial_dimension > 2) {
     this->dir_vecs.push_back(new Vector<Real>(spatial_dimension));
     (*this->dir_vecs.back())[2] = 1.;
     this->registerParam("n3", *(this->dir_vecs.back()), _pat_parsmod,
                         "Direction of tertiary material axis");
   }
 
   for (UInt i = 0 ;  i < this->voigt_h.size ; ++i) {
     UInt start = 0;
     if (this->symmetric) {
       start = i;
     }
     for (UInt j = start ;  j < this->voigt_h.size ; ++j) {
       std::stringstream param("C");
       param << "C" << i+1 << j+1;
       this->registerParam(param.str() , this->Cprime(i,j), 0., _pat_parsmod,
                           "Coefficient " + param.str());
     }
   }
 
   this->registerParam("alpha", this->alpha, _pat_parsmod,
                       "Proportion of viscous stress");
 
   AKANTU_DEBUG_OUT();
   }
 
 /* -------------------------------------------------------------------------- */
 template<UInt spatial_dimension>
 MaterialElasticLinearAnisotropic<spatial_dimension>::~MaterialElasticLinearAnisotropic() {
   for (UInt i = 0 ;  i < spatial_dimension ; ++i) {
     delete this->dir_vecs[i];
   }
 }
 
 /* -------------------------------------------------------------------------- */
 template<UInt spatial_dimension>
 void MaterialElasticLinearAnisotropic<spatial_dimension>::initMaterial() {
   AKANTU_DEBUG_IN();
   Material::initMaterial();
 
   updateInternalParameters();
 
   AKANTU_DEBUG_OUT();
 }
 
 
 /* -------------------------------------------------------------------------- */
 template<UInt spatial_dimension>
 void MaterialElasticLinearAnisotropic<spatial_dimension>::updateInternalParameters() {
   Material::updateInternalParameters();
   if (this->symmetric) {
     for (UInt i = 0 ;  i < this->voigt_h.size ; ++i) {
       for (UInt j = i+1 ;  j < this->voigt_h.size ; ++j) {
         this->Cprime(j, i) = this->Cprime(i, j);
       }
     }
   }
   this->rotateCprime();
   this->C.eig(this->eigC);
 }
 
 /* -------------------------------------------------------------------------- */
 template<UInt Dim>
 void MaterialElasticLinearAnisotropic<Dim>::rotateCprime() {
   // start by filling the empty parts fo Cprime
   UInt diff = Dim*Dim - this->voigt_h.size;
   for (UInt i = this->voigt_h.size ;  i < Dim*Dim ; ++i) {
     for (UInt j = 0 ;  j < Dim*Dim ; ++j) {
       this->Cprime(i, j) = this->Cprime(i-diff, j);
     }
   }
   for (UInt i = 0 ;  i < Dim*Dim ; ++i) {
     for (UInt j = this->voigt_h.size ;  j < Dim*Dim ; ++j) {
       this->Cprime(i, j) = this->Cprime(i, j-diff);
     }
   }
   // construction of rotator tensor
   // normalise rotation matrix
   for (UInt j = 0 ;  j < Dim ; ++j) {
     Vector<Real> rot_vec = this->rot_mat(j);
     rot_vec = *this->dir_vecs[j];
     rot_vec.normalize();
   }
 
   // make sure the vectors form a right-handed base
   Vector<Real> test_axis(3);
   Vector<Real> v1(3), v2(3), v3(3);
 
   if (Dim == 2){
     for (UInt i = 0; i < Dim; ++i) {
       v1[i] = this->rot_mat(0, i);
       v2[i] = this->rot_mat(1, i);
       v3[i] = 0.;
     }
     v3[2] = 1.;
     v1[2] = 0.;
     v2[2] = 0.;
   }
   else if (Dim == 3){
     v1 = this->rot_mat(0);
     v2 = this->rot_mat(1);
     v3 = this->rot_mat(2);
   }
 
   test_axis.crossProduct(v1,v2);
   test_axis -= v3;
   if (test_axis.norm() > 8*std::numeric_limits<Real>::epsilon()) {
     AKANTU_DEBUG_ERROR("The axis vectors do not form a right-handed coordinate "
                        << "system. I. e., ||n1 x n2 - n3|| should be zero, but "
                        << "it is " << test_axis.norm() << ".");
   }
 
   // create the rotator and the reverse rotator
   Matrix<Real> rotator(Dim * Dim, Dim * Dim);
   Matrix<Real> revrotor(Dim * Dim, Dim * Dim);
   for (UInt i = 0 ;  i < Dim ; ++i) {
     for (UInt j = 0 ;  j < Dim ; ++j) {
       for (UInt k = 0 ;  k < Dim ; ++k) {
         for (UInt l = 0 ;  l < Dim ; ++l) {
           UInt I = this->voigt_h.mat[i][j];
           UInt J = this->voigt_h.mat[k][l];
           rotator (I, J) = this->rot_mat(k, i) * this->rot_mat(l, j);
           revrotor(I, J) = this->rot_mat(i, k) * this->rot_mat(j, l);
         }
       }
     }
   }
 
   // create the full rotated matrix
   Matrix<Real> Cfull(Dim*Dim, Dim*Dim);
   Cfull = rotator*Cprime*revrotor;
 
   for (UInt i = 0 ;  i < this->voigt_h.size ; ++i) {
     for (UInt j = 0 ;  j < this->voigt_h.size ; ++j) {
       this->C(i, j) = Cfull(i, j);
     }
   }
 
 }
 
 /* -------------------------------------------------------------------------- */
 template<UInt spatial_dimension>
 void MaterialElasticLinearAnisotropic<spatial_dimension>::computeStress(ElementType el_type,
 									GhostType ghost_type) {
   // Wikipedia convention:
   // 2*eps_ij (i!=j) = voigt_eps_I
   // http://en.wikipedia.org/wiki/Voigt_notation
   AKANTU_DEBUG_IN();
 
   Array<Real>::iterator< Matrix<Real> > gradu_it =
     this->gradu(el_type, ghost_type).begin(spatial_dimension,
 					   spatial_dimension);
   Array<Real>::iterator< Matrix<Real> > gradu_end =
     this->gradu(el_type, ghost_type).end(spatial_dimension,
 					 spatial_dimension);
 
   UInt nb_quad_pts = gradu_end - gradu_it;
 
   // create array for strains and stresses of all dof of all gauss points
   // for efficient computation of stress
   Matrix<Real> voigt_strains(this->voigt_h.size, nb_quad_pts);
   Matrix<Real> voigt_stresses(this->voigt_h.size, nb_quad_pts);
 
   // copy strains
   Matrix<Real> strain(spatial_dimension, spatial_dimension);
 
   for (UInt q = 0; gradu_it != gradu_end ; ++gradu_it, ++q) {
     Matrix<Real> & grad_u = *gradu_it;
 
     for(UInt I = 0; I < this->voigt_h.size; ++I) {
       Real voigt_factor = this->voigt_h.factors[I];
       UInt i = this->voigt_h.vec[I][0];
       UInt j = this->voigt_h.vec[I][1];
 
       voigt_strains(I, q) = voigt_factor * (grad_u(i, j) + grad_u(j, i)) / 2.;
     }
   }
 
 
   // compute the strain rate proportional part if needed
-  bool viscous = this->alpha == 0.; // only works id default value
+  //bool viscous = this->alpha == 0.; // only works if default value
+  bool viscous = false;
   if(viscous) {
     Array<Real> strain_rate(0, spatial_dimension * spatial_dimension,
 			    "strain_rate");
 
     Array<Real> & velocity = this->model->getVelocity();
     const Array<UInt> & elem_filter = this->element_filter(el_type, ghost_type);
 
     this->model->getFEEngine().gradientOnQuadraturePoints(velocity, strain_rate,
 							  spatial_dimension, el_type,
 							  ghost_type, elem_filter);
 
     Array<Real>::matrix_iterator gradu_dot_it  = strain_rate.begin(spatial_dimension,
 								   spatial_dimension);
     Array<Real>::matrix_iterator gradu_dot_end = strain_rate.end(spatial_dimension,
 								 spatial_dimension);
 
     Matrix<Real> strain_dot(spatial_dimension, spatial_dimension);
     for (UInt q = 0; gradu_dot_it != gradu_dot_end ; ++gradu_dot_it, ++q) {
       Matrix<Real> & grad_u_dot = *gradu_dot_it;
 
       for(UInt I = 0; I < this->voigt_h.size; ++I) {
 	Real voigt_factor = this->voigt_h.factors[I];
 	UInt i = this->voigt_h.vec[I][0];
 	UInt j = this->voigt_h.vec[I][1];
 
 	voigt_strains(I, q) = this->alpha * voigt_factor * (grad_u_dot(i, j) + grad_u_dot(j, i)) / 2.;
       }
     }
   }
 
   // compute stresses
   voigt_stresses = this->C * voigt_strains;
 
   // copy stresses back
   Array<Real>::iterator< Matrix<Real> > stress_it =
     this->stress(el_type, ghost_type).begin(spatial_dimension,
 					    spatial_dimension);
 
   Array<Real>::iterator< Matrix<Real> > stress_end =
     this->stress(el_type, ghost_type).end(spatial_dimension,
 					  spatial_dimension);
 
   for (UInt q = 0 ; stress_it != stress_end; ++stress_it, ++q) {
     Matrix<Real> & stress = *stress_it;
 
     for(UInt I = 0; I < this->voigt_h.size; ++I) {
       UInt i = this->voigt_h.vec[I][0];
       UInt j = this->voigt_h.vec[I][1];
       stress(i, j) = stress(j, i) = voigt_stresses(I, q);
     }
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template<UInt spatial_dimension>
 void MaterialElasticLinearAnisotropic<spatial_dimension>::computeTangentModuli(const ElementType & el_type,
                                                                                Array<Real> & tangent_matrix,
                                                                                GhostType ghost_type) {
   AKANTU_DEBUG_IN();
 
   MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_BEGIN(tangent_matrix);
   tangent.copy(this->C);
   MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_END;
 
   AKANTU_DEBUG_OUT();
 }
 
 
 /* -------------------------------------------------------------------------- */
 template<UInt spatial_dimension>
 Real MaterialElasticLinearAnisotropic<spatial_dimension>::getCelerity(__attribute__((unused)) const Element & element) const {
   return std::sqrt( this->eigC(0) / rho);
 }
 
 
 /* -------------------------------------------------------------------------- */
 
 INSTANSIATE_MATERIAL(MaterialElasticLinearAnisotropic);
 
 
 __END_AKANTU__