diff --git a/src/material_damage/material_brittle.cc b/src/material_damage/material_brittle.cc
index 2b9a63d98..8e3c3d964 100644
--- a/src/material_damage/material_brittle.cc
+++ b/src/material_damage/material_brittle.cc
@@ -1,122 +1,122 @@
 /**
  * @file   material_brittle.cc
  *
  * @author Aranda Ruiz Josue <josue.arandaruiz@epfl.ch>
  * @author Daniel Pino Muñoz <daniel.pinomunoz@epfl.ch>
  *
  *
  * @brief  Specialization of the material class for the brittle material
  *
  * @section LICENSE
  *
  * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  */
 
 /* -------------------------------------------------------------------------- */
 #include "material_brittle.hh"
 #include "solid_mechanics_model.hh"
 
 __BEGIN_AKANTU__
 
 /* -------------------------------------------------------------------------- */
 template<UInt spatial_dimension>
 MaterialBrittle<spatial_dimension>::MaterialBrittle(SolidMechanicsModel & model,
                                                     const ID & id):
   Material(model, id),
   MaterialDamage<spatial_dimension>(model, id),
   strain_rate_brittle("strain_rate_brittle", *this){
   AKANTU_DEBUG_IN();
 
   this->registerParam("S_0",        S_0, 157e6, ParamAccessType(_pat_parsable | _pat_modifiable));
   this->registerParam("E_0",        E_0,    27e3, _pat_parsable, "Strain rate threshold");
   this->registerParam("A",          A, 1.622e-5, _pat_parsable, "Polynome cubic constant");
   this->registerParam("B",          B, -1.3274, _pat_parsable, "Polynome quadratic constant");
   this->registerParam("C",          C, 3.6544e4, _pat_parsable, "Polynome linear constant");
   this->registerParam("D",          D, -181.38e6, _pat_parsable, "Polynome constant");
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template<UInt spatial_dimension>
 void MaterialBrittle<spatial_dimension>::initMaterial() {
   AKANTU_DEBUG_IN();
   MaterialDamage<spatial_dimension>::initMaterial();
 
   this->strain_rate_brittle.initialize(spatial_dimension * spatial_dimension);
 
   updateInternalParameters();
 
   //this->Yd.resize();
   // const Mesh & mesh = this->model->getFEEngine().getMesh();
 
   // Mesh::type_iterator it = mesh.firstType(spatial_dimension);
   // Mesh::type_iterator last_type = mesh.lastType(spatial_dimension);
 
   // for(; it != last_type; ++it) {
   //   UInt nb_element  = this->element_filter(*it).getSize();
   //   UInt nb_quad = this->model->getFEEngine().getNbQuadraturePoints(*it);
 
   //   Array <Real> & Yd_rand_vec = Yd_rand(*it);
   //   for(UInt e = 0; e < nb_element; ++e) {
   //     Real rand_part = (2 * drand48()-1) * Yd_randomness * Yd;
 
   //     for(UInt q = 0; q < nb_quad; ++q)
   //    Yd_rand_vec(nb_quad*e+q,0) = Yd + rand_part;
   //   }
   // }
 
   AKANTU_DEBUG_OUT();
 }
 
 
 /* -------------------------------------------------------------------------- */
 template<UInt spatial_dimension>
 void MaterialBrittle<spatial_dimension>::updateInternalParameters() {
   MaterialDamage<spatial_dimension>::updateInternalParameters();
 }
 
 /* -------------------------------------------------------------------------- */
 template<UInt spatial_dimension>
 void MaterialBrittle<spatial_dimension>::computeStress(ElementType el_type,
                                                       GhostType ghost_type) {
   AKANTU_DEBUG_IN();
 
   Real * dam = this->damage(el_type, ghost_type).storage();
 
   Array<Real> & velocity = this->model->getVelocity();
   Array<Real> & strain_rate_brittle = this->strain_rate_brittle(el_type, ghost_type);
   Array<UInt> & elem_filter = this->element_filter(el_type, ghost_type);
 
-  this->model->getFEEngine().gradientOnQuadraturePoints(velocity, strain_rate_brittle,
+  this->model->getFEEngine().gradientOnIntegrationPoints(velocity, strain_rate_brittle,
 						   spatial_dimension,
 						   el_type, ghost_type, elem_filter);
 
   Array<Real>::iterator< Matrix<Real> > strain_rate_it =
     this->strain_rate_brittle(el_type, ghost_type).begin(spatial_dimension, spatial_dimension);
 
 
   MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type);
 
   Real sigma_equivalent = 0;
   Real fracture_stress = 0;
   Matrix<Real> & grad_v = *strain_rate_it;
 
   computeStressOnQuad(grad_u, grad_v,  sigma, *dam, sigma_equivalent, fracture_stress);
 
   ++strain_rate_it;
   ++dam;
 
   MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END;
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 
 INSTANTIATE_MATERIAL(MaterialBrittle);
 
 
 __END_AKANTU__
diff --git a/src/material_damage/material_brittle_non_local_inline_impl.cc b/src/material_damage/material_brittle_non_local_inline_impl.cc
index f73deace6..216bc78d6 100644
--- a/src/material_damage/material_brittle_non_local_inline_impl.cc
+++ b/src/material_damage/material_brittle_non_local_inline_impl.cc
@@ -1,111 +1,111 @@
 /**
  * @file   material_brittle_non_local_inline_impl.cc
  *
  * @author Daniel Pino Muñoz <daniel.pinomunoz@epfl.ch>
  *
  *
  * @brief  MaterialBrittleNonLocal inline function implementation
  *
  * @section LICENSE
  *
  * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  */
 
 /* -------------------------------------------------------------------------- */
 __END_AKANTU__
 
 #if defined(AKANTU_DEBUG_TOOLS)
 #include "aka_debug_tools.hh"
 #include <string>
 #endif
 
 __BEGIN_AKANTU__
 
 /* -------------------------------------------------------------------------- */
 template<UInt spatial_dimension>
 MaterialBrittleNonLocal<spatial_dimension>::MaterialBrittleNonLocal(SolidMechanicsModel & model, const ID & id)  :
   Material(model, id),
   MaterialBrittleNonLocalParent(model, id),
   Sigma_max("Sigma max", *this),
   Sigma_maxnl("Sigma_max non local", *this),
   Sigma_fracture("Sigma_fracture", *this){
   AKANTU_DEBUG_IN();
   this->is_non_local = true;
   this->Sigma_max.initialize(1);
   this->Sigma_maxnl.initialize(1);
   this->Sigma_fracture.initialize(1);
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template<UInt spatial_dimension>
 void MaterialBrittleNonLocal<spatial_dimension>::initMaterial() {
   AKANTU_DEBUG_IN();
   this->model->getNonLocalManager().registerNonLocalVariable(this->Sigma_max.getName(), Sigma_maxnl.getName(),1);
   MaterialBrittleNonLocalParent::initMaterial();
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template<UInt spatial_dimension>
 void MaterialBrittleNonLocal<spatial_dimension>::computeStress(ElementType el_type, GhostType ghost_type) {
   AKANTU_DEBUG_IN();
 
   Real * dam = this->damage(el_type, ghost_type).storage();
   Real * Sigma_maxt  = this->Sigma_max(el_type, ghost_type).storage();
   Real * fracture_stress  = this->Sigma_fracture(el_type, ghost_type).storage();
 
   Array<Real> & velocity = this->model->getVelocity();
   Array<Real> & strain_rate_brittle = this->strain_rate_brittle(el_type, ghost_type);
   Array<UInt> & elem_filter = this->element_filter(el_type, ghost_type);
 
-  this->model->getFEEngine().gradientOnQuadraturePoints(velocity, strain_rate_brittle,
+  this->model->getFEEngine().gradientOnIntegrationPoints(velocity, strain_rate_brittle,
 						   spatial_dimension,
 						   el_type, ghost_type, elem_filter);
 
   Array<Real>::iterator< Matrix<Real> > strain_rate_it =
     this->strain_rate_brittle(el_type, ghost_type).begin(spatial_dimension, spatial_dimension);
 
   MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type);
 
   Matrix<Real> & grad_v = *strain_rate_it;
 
   MaterialBrittle<spatial_dimension>::computeStressOnQuad(grad_u, grad_v, sigma, *dam, *Sigma_maxt, *fracture_stress);
   ++dam;
   ++strain_rate_it;
   ++Sigma_maxt;
   ++fracture_stress;
   MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END;
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template<UInt spatial_dimension>
 void MaterialBrittleNonLocal<spatial_dimension>::computeNonLocalStress(ElementType type,
 								       GhostType ghost_type) {
   AKANTU_DEBUG_IN();
 
   Real * dam  = this->damage(type, ghost_type).storage();
   Real * Sigma_maxnlt = this->Sigma_maxnl(type, ghost_type).storage();
   Real * fracture_stress = this->Sigma_fracture(type, ghost_type).storage();
 
   MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(type, ghost_type);
   this->computeDamageAndStressOnQuad(sigma, *dam, *Sigma_maxnlt, *fracture_stress);
 
   ++dam;
   ++Sigma_maxnlt;
   ++fracture_stress;
   MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END;
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template<UInt spatial_dimension>
 void MaterialBrittleNonLocal<spatial_dimension>::nonLocalVariableToNeighborhood() {
   this->model->getNonLocalManager().nonLocalVariableToNeighborhood(Sigma_maxnl.getName(), this->name);
 }
 
diff --git a/src/material_damage/material_orthotropic_damage_iterative.hh b/src/material_damage/material_orthotropic_damage_iterative.hh
index fa367777e..0938a1faa 100644
--- a/src/material_damage/material_orthotropic_damage_iterative.hh
+++ b/src/material_damage/material_orthotropic_damage_iterative.hh
@@ -1,123 +1,123 @@
 /**
  * @file   material_orthtropic_damage_iterative.hh
  *
  * @author Aurelia Isabel Cuba Ramos <aurelia.cubaramos@epfl.ch>
  *
  * @date   Sun Mar  8 12:54:30 2015
  *
  * @brief Specialization of the class material orthotropic damage to
  * damage only one gauss point at a time and propagate damage in a
  * linear way. Max principal stress criterion is used as a failure
  * criterion.
  *
  * @section LICENSE
  *
  * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  */
 
 /* -------------------------------------------------------------------------- */
 #include "aka_common.hh"
 #include "material_orthotropic_damage.hh"
 #include "material.hh"
 /* -------------------------------------------------------------------------- */
 
 #ifndef __AKANTU_MATERIAL_ORTHOTROPIC_DAMAGE_ITERATIVE_HH__
 #define __AKANTU_MATERIAL_ORTHOTROPIC_DAMAGE_ITERATIVE_HH__
 
 __BEGIN_AKANTU__
 
 /**
  * Material damage iterative
  *
  * parameters in the material files :
  *   - Sc  
  */
 template<UInt spatial_dimension>
 class MaterialOrthotropicDamageIterative : public MaterialOrthotropicDamage<spatial_dimension> {
   /* ------------------------------------------------------------------------ */
   /* Constructors/Destructors                                                 */
   /* ------------------------------------------------------------------------ */
 public:
 
   MaterialOrthotropicDamageIterative(SolidMechanicsModel & model, const ID & id = "");
 
   virtual ~MaterialOrthotropicDamageIterative() {};
 
   /* ------------------------------------------------------------------------ */
   /* Methods                                                                  */
   /* ------------------------------------------------------------------------ */
 public:
 
   ///  virtual void updateInternalParameters();
 
   virtual void computeAllStresses(GhostType ghost_type = _not_ghost);
 
   /// update internal field damage
   UInt updateDamage();
 
   /// update energies after damage has been updated
   virtual void updateEnergiesAfterDamage(ElementType el_type, GhostType ghost_typ);
   
   virtual void onBeginningSolveStep(const AnalysisMethod & method) { };
 
   virtual void onEndSolveStep(const AnalysisMethod & method) { };
 protected:
   /// constitutive law for all element of a type
   virtual void computeStress(ElementType el_type, GhostType ghost_type = _not_ghost);
 
   ///compute the equivalent stress on each Gauss point (i.e. the max prinicpal stress) and normalize it by the tensile strength
   void computeNormalizedEquivalentStress(const Array<Real> & grad_u, ElementType el_type, GhostType ghost_type = _not_ghost);
 
   /// find max normalized equivalent stress
   void findMaxNormalizedEquivalentStress(ElementType el_type, GhostType ghost_type = _not_ghost);
 
   inline void computeDamageAndStressOnQuad(Matrix<Real> & sigma, Matrix<Real> & one_minus_D, Matrix<Real> & root_one_minus_D, Matrix<Real> & damage, Matrix<Real> & first_term, Matrix<Real> & third_term);
   /* ------------------------------------------------------------------------ */
   /* DataAccessor inherited members                                           */
   /* ------------------------------------------------------------------------ */
   /* ------------------------------------------------------------------------ */
   /* Accessors                                                                */
   /* ------------------------------------------------------------------------ */
 public:
   /// get max normalized equivalent stress
   AKANTU_GET_MACRO(NormMaxEquivalentStress, norm_max_equivalent_stress, Real);
 
   /* ------------------------------------------------------------------------ */
   /* Class Members                                                            */
   /* ------------------------------------------------------------------------ */
 protected:
 
   /// resistance to damage
   RandomInternalField<Real> Sc;
 
   /// internal field to store equivalent stress on each Gauss point
   InternalField<Real> equivalent_stress;
 
   /// internal field to store the direction of the current damage frame
   InternalField<Real> stress_dir;
 
   /// damage increment
   Real prescribed_dam;
 
   /// maximum equivalent stress
   Real norm_max_equivalent_stress;
 
   /// define damage threshold at which damage will be set to 1 
   Real dam_threshold; 
 
   /// quadrature point with highest equivalent Stress
-  QuadraturePoint q_max;
+  IntegrationPoint q_max;
 
 };
 
 /* -------------------------------------------------------------------------- */
 /* inline functions                                                           */
 /* -------------------------------------------------------------------------- */
 
 #include "material_orthotropic_damage_iterative_inline_impl.cc"
 
 __END_AKANTU__
 
 #endif /* __AKANTU_MATERIAL_ORTHOTROPIC_DAMAGE_ITERATIVE_HH__ */
diff --git a/src/material_damage/material_vreepeerlings_tmpl.hh b/src/material_damage/material_vreepeerlings_tmpl.hh
index 8f9bcbb37..f0406a03a 100644
--- a/src/material_damage/material_vreepeerlings_tmpl.hh
+++ b/src/material_damage/material_vreepeerlings_tmpl.hh
@@ -1,112 +1,112 @@
 /**
  * @file   material_vreepeerlings_tmpl.hh
  *
  * @author Cyprien Wolff <cyprien.wolff@epfl.ch>
  *
  *
  * @brief  Specialization of the material class for the VreePeerlings material
  *
  * @section LICENSE
  *
  * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  */
 
 /* -------------------------------------------------------------------------- */
 template<UInt spatial_dimension, template <UInt> class MatParent>
 MaterialVreePeerlings<spatial_dimension, MatParent>::MaterialVreePeerlings(SolidMechanicsModel & model,
                                              const ID & id)  :
   Material(model, id),
   MaterialVreePeerlingsParent(model, id),
   Kapa("Kapa", *this),
   strain_rate_vreepeerlings("strain-rate-vreepeerlings", *this),
   Full_dam_value_strain("fulldam-valstrain", *this),
   Full_dam_value_strain_rate("fulldam-valstrain-rate", *this),
   Number_damage("number-damage", *this),
   equi_strain("equi-strain", *this),
   equi_strain_rate("equi-strain-rate", *this) {
   AKANTU_DEBUG_IN();
 
   this->registerParam("Kapaoi"          , Kapaoi          , 0.0001, _pat_parsable);
   this->registerParam("Kapac"           , Kapac           , 0.0002, _pat_parsable);
   this->registerParam("Arate"           , Arate           , 0.    , _pat_parsable);
   this->registerParam("Brate"           , Brate           , 1.    , _pat_parsable);
   this->registerParam("Crate"           , Brate           , 1.    , _pat_parsable);
   this->registerParam("Kct"             , Kct             , 1.    , _pat_parsable);
   this->registerParam("Kapao_randomness", Kapao_randomness, 0.    , _pat_parsable);
 
   this->Kapa.initialize(1);
   this->equi_strain.initialize(1);
   this->equi_strain_rate.initialize(1);
   this->Full_dam_value_strain.initialize(1);
   this->Full_dam_value_strain_rate.initialize(1);
   this->Number_damage.initialize(1);
   this->strain_rate_vreepeerlings.initialize(spatial_dimension * spatial_dimension);
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template<UInt spatial_dimension, template <UInt> class MatParent>
 void MaterialVreePeerlings<spatial_dimension, MatParent>::initMaterial() {
   AKANTU_DEBUG_IN();
   MaterialVreePeerlingsParent::initMaterial();
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template<UInt spatial_dimension, template <UInt> class MatParent>
 void MaterialVreePeerlings<spatial_dimension, MatParent>::computeStress(ElementType el_type, GhostType ghost_type) {
   AKANTU_DEBUG_IN();
 
   MaterialVreePeerlingsParent::computeStress(el_type, ghost_type);
 
   Real * dam = this->damage(el_type, ghost_type).storage();
   Real * equi_straint = equi_strain(el_type, ghost_type).storage();
   Real * equi_straint_rate = equi_strain_rate(el_type, ghost_type).storage();
   Real * Kapaq = Kapa(el_type, ghost_type).storage();
   Real * FullDam_Valstrain = Full_dam_value_strain(el_type, ghost_type).storage();
   Real * FullDam_Valstrain_rate = Full_dam_value_strain_rate(el_type, ghost_type).storage();
   Real * Nb_damage = Number_damage(el_type, ghost_type).storage();
   Real dt = this->model->getTimeStep();
 
   Array<UInt> & elem_filter = this->element_filter(el_type, ghost_type);
   Array<Real> & velocity = this->model->getVelocity();
   Array<Real> & strain_rate_vrplgs = this->strain_rate_vreepeerlings(el_type, ghost_type);
 
-  this->model->getFEEngine().gradientOnQuadraturePoints(velocity, strain_rate_vrplgs,
+  this->model->getFEEngine().gradientOnIntegrationPoints(velocity, strain_rate_vrplgs,
 						   spatial_dimension,
 						   el_type, ghost_type, elem_filter);
 
   Array<Real>::matrix_iterator strain_rate_vrplgs_it =
         strain_rate_vrplgs.begin(spatial_dimension, spatial_dimension);
 
   MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type);
 
   Matrix<Real> & strain_rate = *strain_rate_vrplgs_it;
 
   computeStressOnQuad(grad_u, sigma, 
 		      *dam, 
 		      *equi_straint,
 		      *equi_straint_rate,
 		      *Kapaq,
 		      dt,
 		      strain_rate,
 		      *FullDam_Valstrain,
 		      *FullDam_Valstrain_rate,
 		      *Nb_damage);
   ++dam;
   ++equi_straint;
   ++equi_straint_rate;
   ++Kapaq;
   ++strain_rate_vrplgs_it;
   ++FullDam_Valstrain;
   ++FullDam_Valstrain_rate;
   ++Nb_damage;
 
   MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END;
 
   AKANTU_DEBUG_OUT();
 }
 
diff --git a/src/material_viscoelastic/material_stiffness_proportional.cc b/src/material_viscoelastic/material_stiffness_proportional.cc
index d07a40697..738cfbdcd 100644
--- a/src/material_viscoelastic/material_stiffness_proportional.cc
+++ b/src/material_viscoelastic/material_stiffness_proportional.cc
@@ -1,128 +1,128 @@
 /**
  * @file   material_stiffness_proportional.cc
  *
  * @author David Simon Kammer <david.kammer@epfl.ch>
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  *
  * @brief  Special. of the material class for the caughey viscoelastic material
  *
  * @section LICENSE
  *
  * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  */
 
 /* -------------------------------------------------------------------------- */
 #include "material_stiffness_proportional.hh"
 #include "solid_mechanics_model.hh"
 
 __BEGIN_AKANTU__
 
 /* -------------------------------------------------------------------------- */
 template<UInt spatial_dimension>
 MaterialStiffnessProportional<spatial_dimension>::MaterialStiffnessProportional(SolidMechanicsModel & model,
 										const ID & id)  :
   Material(model, id), MaterialElastic<spatial_dimension>(model, id),
   stress_viscosity("stress_viscosity", *this),
   stress_elastic("stress_elastic", *this) {
   AKANTU_DEBUG_IN();
 
   this->registerParam("Alpha", alpha, 0., ParamAccessType(_pat_parsable | _pat_modifiable), "Artificial viscous ratio");
 
   this->stress_viscosity.initialize(spatial_dimension * spatial_dimension);
   this->stress_elastic  .initialize(spatial_dimension * spatial_dimension);
 
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template<UInt spatial_dimension>
 void MaterialStiffnessProportional<spatial_dimension>::initMaterial() {
   AKANTU_DEBUG_IN();
   MaterialElastic<spatial_dimension>::initMaterial();
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template<UInt spatial_dimension>
 void MaterialStiffnessProportional<spatial_dimension>::computeStress(ElementType el_type,
 							      GhostType ghost_type) {
   AKANTU_DEBUG_IN();
   Array<UInt> & elem_filter = this->element_filter  (el_type, ghost_type);
   Array<Real> & stress_visc = stress_viscosity(el_type, ghost_type);
   Array<Real> & stress_el   = stress_elastic  (el_type, ghost_type);
 
   MaterialElastic<spatial_dimension>::computeStress(el_type, ghost_type);
 
   Array<Real> & velocity = this->model->getVelocity();
 
   Array<Real> strain_rate(0, spatial_dimension * spatial_dimension,
 			   "strain_rate");
 
-  this->model->getFEEngine().gradientOnQuadraturePoints(velocity, strain_rate,
+  this->model->getFEEngine().gradientOnIntegrationPoints(velocity, strain_rate,
 						   spatial_dimension,
 						   el_type, ghost_type, elem_filter);
 
   Array<Real>::matrix_iterator strain_rate_it =
     strain_rate.begin(spatial_dimension, spatial_dimension);
   Array<Real>::matrix_iterator stress_visc_it =
     stress_visc.begin(spatial_dimension, spatial_dimension);
   Array<Real>::matrix_iterator stress_el_it =
     stress_el.begin(spatial_dimension, spatial_dimension);
 
   MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type);
   Matrix<Real> & grad_v     = *strain_rate_it;
   Matrix<Real> & sigma_visc = *stress_visc_it;
   Matrix<Real> & sigma_el   = *stress_el_it;
 
   MaterialElastic<spatial_dimension>::computeStressOnQuad(grad_v, sigma_visc);
 
   sigma_visc *= alpha;
   sigma_el.copy(sigma);
   sigma += sigma_visc;
 
   ++strain_rate_it;
   ++stress_visc_it;
   ++stress_el_it;
 
   MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END;
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template<UInt spatial_dimension>
 void MaterialStiffnessProportional<spatial_dimension>::computePotentialEnergy(ElementType el_type,
 								       GhostType ghost_type) {
   AKANTU_DEBUG_IN();
 
   if(ghost_type != _not_ghost) return;
 
   Array<Real> & stress_el = stress_elastic(el_type, ghost_type);
   Array<Real>::matrix_iterator stress_el_it =
     stress_el.begin(spatial_dimension, spatial_dimension);
 
   Real * epot = this->potential_energy(el_type, ghost_type).storage();
 
   MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type);
   Matrix<Real> & sigma_el    = *stress_el_it;
   MaterialElastic<spatial_dimension>::computePotentialEnergyOnQuad(grad_u,
 								   sigma_el,
 								   *epot);
   epot++;
   ++stress_el_it;
   MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END;
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 
 
 INSTANTIATE_MATERIAL(MaterialStiffnessProportional);
 
 
 __END_AKANTU__