diff --git a/src/model/solid_mechanics/material.cc b/src/model/solid_mechanics/material.cc
index 788575218..c194b02bd 100644
--- a/src/model/solid_mechanics/material.cc
+++ b/src/model/solid_mechanics/material.cc
@@ -1,1772 +1,1772 @@
 /**
  * @file   material.cc
  *
  * @author Aurelia Isabel Cuba Ramos <aurelia.cubaramos@epfl.ch>
  * @author Marco Vocialta <marco.vocialta@epfl.ch>
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  * @author Daniel Pino Muñoz <daniel.pinomunoz@epfl.ch>
  *
  * @date creation: Tue Jul 27 2010
  * @date last modification: Tue Sep 16 2014
  *
  * @brief  Implementation of the common part of the material class
  *
  * @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 "material.hh"
 #include "solid_mechanics_model.hh"
 #include "sparse_matrix.hh"
 #include "dof_synchronizer.hh"
 /* -------------------------------------------------------------------------- */
 
 
 __BEGIN_AKANTU__
 
 /* -------------------------------------------------------------------------- */
 Material::Material(SolidMechanicsModel & model, const ID & id) :
   Memory(id, model.getMemoryID()),
   Parsable(_st_material, id),
   is_init(false),
   finite_deformation(false),
   name(""),
   model(&model),
   spatial_dimension(this->model->getSpatialDimension()),
   element_filter("element_filter", id, this->memory_id),
   stress("stress", *this),
   eigenstrain("eigenstrain", *this),
   gradu("grad_u", *this),
   green_strain("green_strain",*this),
   piola_kirchhoff_2("piola_kirchhoff_2", *this),
   //  potential_energy_vector(false),
   potential_energy("potential_energy", *this),
   is_non_local(false),
   use_previous_stress(false),
   use_previous_gradu(false),
   interpolation_inverse_coordinates("interpolation inverse coordinates", *this),
   interpolation_points_matrices("interpolation points matrices", *this) {
   AKANTU_DEBUG_IN();
 
   /// for each connectivity types allocate the element filer array of the material
   model.getMesh().initElementTypeMapArray(element_filter,
 					 1,
 					 spatial_dimension,
 					 false,
 					 _ek_regular);
 
 
   this->initialize();
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 Material::Material(SolidMechanicsModel & model,
                    UInt dim,
                    const Mesh & mesh,
                    FEEngine & fe_engine,
                    const ID & id) :
   Memory(id, model.getMemoryID()),
   Parsable(_st_material, id),
   is_init(false),
   finite_deformation(false),
   name(""),
   model(&model),
   spatial_dimension(dim),
   element_filter("element_filter", id, this->memory_id),
   stress("stress", *this, dim, fe_engine, this->element_filter),
   eigenstrain("eignestrain", *this, dim, fe_engine, this->element_filter),
   gradu("gradu", *this, dim, fe_engine, this->element_filter),
   green_strain("green_strain", *this, dim, fe_engine, this->element_filter),
   piola_kirchhoff_2("poila_kirchhoff_2", *this, dim, fe_engine, this->element_filter),
   potential_energy("potential_energy", *this, dim, fe_engine, this->element_filter),
   is_non_local(false),
   use_previous_stress(false),
   use_previous_gradu(false),
   interpolation_inverse_coordinates("interpolation inverse_coordinates", *this,
                                     dim, fe_engine, this->element_filter),
   interpolation_points_matrices("interpolation points matrices", *this,
                                 dim, fe_engine, this->element_filter) {
 
   AKANTU_DEBUG_IN();
   mesh.initElementTypeMapArray(element_filter,
                                1,
                                spatial_dimension,
                                false,
                                _ek_regular);
 
   this->initialize();
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 Material::~Material() {
   AKANTU_DEBUG_IN();
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 void Material::initialize() {
   registerParam("rho"                  , rho                  , 0.           , _pat_parsable | _pat_modifiable, "Density");
   registerParam("name"                 , name                 , std::string(), _pat_parsable | _pat_readable);
   registerParam("finite_deformation"   , finite_deformation   , false        , _pat_parsable | _pat_readable, "Is finite deformation");
   registerParam("inelastic_deformation", inelastic_deformation, false        ,  _pat_internal, "Is inelastic deformation");
 
   /// allocate gradu stress for local elements
   eigenstrain.initialize(spatial_dimension * spatial_dimension);
   gradu.initialize(spatial_dimension * spatial_dimension);
   stress.initialize(spatial_dimension * spatial_dimension);
 
   this->model->registerEventHandler(*this);
 }
 
 /* -------------------------------------------------------------------------- */
 void Material::initMaterial() {
   AKANTU_DEBUG_IN();
 
   if(finite_deformation) {
     this->piola_kirchhoff_2.initialize(spatial_dimension * spatial_dimension);
     if(use_previous_stress) this->piola_kirchhoff_2.initializeHistory();
     this->green_strain.initialize(spatial_dimension * spatial_dimension);
   }
 
   if(use_previous_stress) this->stress.initializeHistory();
   if(use_previous_gradu) this->gradu.initializeHistory();
 
   for (std::map<ID, InternalField<Real> *>::iterator it = internal_vectors_real.begin();
        it != internal_vectors_real.end();
        ++it) it->second->resize();
 
   for (std::map<ID, InternalField<UInt> *>::iterator it = internal_vectors_uint.begin();
        it != internal_vectors_uint.end();
        ++it) it->second->resize();
 
   for (std::map<ID, InternalField<bool> *>::iterator it = internal_vectors_bool.begin();
        it != internal_vectors_bool.end();
        ++it) it->second->resize();
 
   is_init = true;
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 void Material::savePreviousState() {
   AKANTU_DEBUG_IN();
 
   for (std::map<ID, InternalField<Real> *>::iterator it = internal_vectors_real.begin();
        it != internal_vectors_real.end();
        ++it) {
     if(it->second->hasHistory()) it->second->saveCurrentValues();
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 
 /* -------------------------------------------------------------------------- */
 /**
  * Compute  the  residual  by  assembling  @f$\int_{e}  \sigma_e  \frac{\partial
  * \varphi}{\partial X} dX @f$
  *
  * @param[in] displacements nodes displacements
  * @param[in] ghost_type compute the residual for _ghost or _not_ghost element
  */
 void Material::updateResidual(GhostType ghost_type) {
   AKANTU_DEBUG_IN();
 
   computeAllStresses(ghost_type);
   assembleResidual(ghost_type);
 
   AKANTU_DEBUG_OUT();
 }
 
 
 /* -------------------------------------------------------------------------- */
 void Material::assembleResidual(GhostType ghost_type) {
   AKANTU_DEBUG_IN();
 
   UInt spatial_dimension = model->getSpatialDimension();
 
   if(!finite_deformation){
 
     Array<Real> & residual = const_cast<Array<Real> &>(model->getResidual());
 
     Mesh & mesh = model->getFEEngine().getMesh();
     Mesh::type_iterator it = element_filter.firstType(spatial_dimension, ghost_type);
     Mesh::type_iterator last_type = element_filter.lastType(spatial_dimension, ghost_type);
     for(; it != last_type; ++it) {
       const Array<Real> & shapes_derivatives = model->getFEEngine().getShapesDerivatives(*it, ghost_type);
       Array<UInt> & elem_filter = element_filter(*it, ghost_type);
 
       UInt size_of_shapes_derivatives = shapes_derivatives.getNbComponent();
       UInt nb_nodes_per_element       = Mesh::getNbNodesPerElement(*it);
       UInt nb_quadrature_points       = model->getFEEngine().getNbQuadraturePoints(*it, ghost_type);
 
       UInt nb_element = elem_filter.getSize();
       /// compute @f$\sigma \frac{\partial \varphi}{\partial X}@f$ by @f$\mathbf{B}^t \mathbf{\sigma}_q@f$
       Array<Real> * sigma_dphi_dx =
 	new Array<Real>(nb_element*nb_quadrature_points,
 			size_of_shapes_derivatives, "sigma_x_dphi_/_dX");
 
       Array<Real> * shapesd_filtered =
 	new Array<Real>(0, size_of_shapes_derivatives, "filtered shapesd");
 
       FEEngine::filterElementalData(mesh, shapes_derivatives, *shapesd_filtered,
 			       *it, ghost_type, elem_filter);
 
       Array<Real> & stress_vect = stress(*it, ghost_type);
 
       Array<Real>::matrix_iterator sigma =
 	stress_vect.begin(spatial_dimension, spatial_dimension);
       Array<Real>::matrix_iterator B =
 	shapesd_filtered->begin(spatial_dimension, nb_nodes_per_element);
       Array<Real>::matrix_iterator Bt_sigma_it =
 	sigma_dphi_dx->begin(spatial_dimension, nb_nodes_per_element);
 
       for (UInt q = 0; q < nb_element*nb_quadrature_points; ++q, ++sigma, ++B, ++Bt_sigma_it)
 	Bt_sigma_it->mul<false,false>(*sigma, *B);
 
 
       delete shapesd_filtered;
 
       /**
        * compute @f$\int \sigma  * \frac{\partial \varphi}{\partial X}dX@f$ by  @f$ \sum_q \mathbf{B}^t
        * \mathbf{\sigma}_q \overline w_q J_q@f$
        */
       Array<Real> * int_sigma_dphi_dx = new Array<Real>(nb_element,
 							nb_nodes_per_element * spatial_dimension,
 							"int_sigma_x_dphi_/_dX");
 
       model->getFEEngine().integrate(*sigma_dphi_dx, *int_sigma_dphi_dx,
 				size_of_shapes_derivatives,
 				*it, ghost_type,
 				elem_filter);
       delete sigma_dphi_dx;
 
 
       /// assemble
       model->getFEEngine().assembleArray(*int_sigma_dphi_dx, residual,
 				    model->getDOFSynchronizer().getLocalDOFEquationNumbers(),
 				    residual.getNbComponent(),
 				    *it, ghost_type, elem_filter, -1);
       delete int_sigma_dphi_dx;
     }
 
   }
   else{
     switch (spatial_dimension){
     case 1:
       this->assembleResidual<1>(ghost_type);
       break;
     case 2:
       this->assembleResidual<2>(ghost_type);
       break;
     case 3:
       this->assembleResidual<3>(ghost_type);
       break;
     }
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 /**
  * Compute  the  stress from the gradu
  *
  * @param[in] current_position nodes postition + displacements
  * @param[in] ghost_type compute the residual for _ghost or _not_ghost element
  */
 void Material::computeAllStresses(GhostType ghost_type) {
   AKANTU_DEBUG_IN();
 
   UInt spatial_dimension = model->getSpatialDimension();
 
   Mesh::type_iterator it = model->getFEEngine().getMesh().firstType(spatial_dimension, ghost_type);
   Mesh::type_iterator last_type = model->getFEEngine().getMesh().lastType(spatial_dimension, ghost_type);
 
   for(; it != last_type; ++it) {
     Array<UInt> & elem_filter = element_filter(*it, ghost_type);
     Array<Real> & gradu_vect = gradu(*it, ghost_type);
 
     /// compute @f$\nabla u@f$
     model->getFEEngine().gradientOnQuadraturePoints(model->getDisplacement(), gradu_vect,
 						    spatial_dimension,
 						    *it, ghost_type, elem_filter);
 
     gradu_vect -= eigenstrain(*it, ghost_type);
 
     /// compute @f$\mathbf{\sigma}_q@f$ from @f$\nabla u@f$
     computeStress(*it, ghost_type);
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 void Material::computeAllCauchyStresses(GhostType ghost_type) {
   AKANTU_DEBUG_IN();
 
   AKANTU_DEBUG_ASSERT(finite_deformation,"The Cauchy stress can only be computed if you are working in finite deformation.");
 
   //resizeInternalArray(stress);
 
   Mesh::type_iterator it = model->getFEEngine().getMesh().firstType(spatial_dimension, ghost_type);
   Mesh::type_iterator last_type = model->getFEEngine().getMesh().lastType(spatial_dimension, ghost_type);
 
   for(; it != last_type; ++it)
     switch (spatial_dimension){
     case 1:
       this->computeCauchyStress<1>(*it, ghost_type);
       break;
     case 2:
       this->computeCauchyStress<2>(*it, ghost_type);
       break;
     case 3:
       this->computeCauchyStress<3>(*it, ghost_type);
       break;
     }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template<UInt dim>
 void Material::computeCauchyStress(ElementType el_type, GhostType ghost_type) {
   AKANTU_DEBUG_IN();
 
   Array<Real>::matrix_iterator gradu_it =
     this->gradu(el_type, ghost_type).begin(dim, dim);
 
   Array<Real>::matrix_iterator gradu_end =
     this->gradu(el_type, ghost_type).end(dim, dim);
 
   Array<Real>::matrix_iterator piola_it =
     this->piola_kirchhoff_2(el_type, ghost_type).begin(dim, dim);
 
   Array<Real>::matrix_iterator stress_it =
     this->stress(el_type, ghost_type).begin(dim, dim);
 
   Matrix<Real> F_tensor(dim, dim);
 
   for (; gradu_it != gradu_end; ++gradu_it, ++piola_it, ++stress_it) {
     Matrix<Real> & grad_u = *gradu_it;
     Matrix<Real> & piola  = *piola_it;
     Matrix<Real> & sigma  = *stress_it;
 
     gradUToF<dim > (grad_u, F_tensor);
     this->computeCauchyStressOnQuad<dim >(F_tensor, piola, sigma);
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 void Material::setToSteadyState(GhostType ghost_type) {
   AKANTU_DEBUG_IN();
 
   const Array<Real> & displacement = model->getDisplacement();
 
   //resizeInternalArray(gradu);
 
   UInt spatial_dimension = model->getSpatialDimension();
 
   Mesh::type_iterator it = model->getFEEngine().getMesh().firstType(spatial_dimension, ghost_type);
   Mesh::type_iterator last_type = model->getFEEngine().getMesh().lastType(spatial_dimension, ghost_type);
 
   for(; it != last_type; ++it) {
     Array<UInt> & elem_filter = element_filter(*it, ghost_type);
     Array<Real> & gradu_vect = gradu(*it, ghost_type);
 
     /// compute @f$\nabla u@f$
     model->getFEEngine().gradientOnQuadraturePoints(displacement, gradu_vect,
 					       spatial_dimension,
 					       *it, ghost_type, elem_filter);
 
     setToSteadyState(*it, ghost_type);
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 
 /* -------------------------------------------------------------------------- */
 /**
  * Compute  the stiffness  matrix by  assembling @f$\int_{\omega}  B^t  \times D
  * \times B d\omega @f$
  *
  * @param[in] current_position nodes postition + displacements
  * @param[in] ghost_type compute the residual for _ghost or _not_ghost element
  */
 void Material::assembleStiffnessMatrix(GhostType ghost_type) {
   AKANTU_DEBUG_IN();
 
   UInt spatial_dimension = model->getSpatialDimension();
 
   Mesh::type_iterator it = element_filter.firstType(spatial_dimension, ghost_type);
   Mesh::type_iterator last_type = element_filter.lastType(spatial_dimension, ghost_type);
   for(; it != last_type; ++it) {
     if(finite_deformation){
       switch (spatial_dimension) {
       case 1:
 	{
 	  assembleStiffnessMatrixNL < 1 > (*it, ghost_type);
 	  assembleStiffnessMatrixL2 < 1 > (*it, ghost_type);
 	  break;
 	}
       case 2:
 	{
 	  assembleStiffnessMatrixNL < 2 > (*it, ghost_type);
 	  assembleStiffnessMatrixL2 < 2 > (*it, ghost_type);
 	  break;
 	}
       case 3:
 	{
 	  assembleStiffnessMatrixNL < 3 > (*it, ghost_type);
 	  assembleStiffnessMatrixL2 < 3 > (*it, ghost_type);
 	  break;
 	}
       }
     } else {
       switch(spatial_dimension) {
       case 1: { assembleStiffnessMatrix<1>(*it, ghost_type); break; }
       case 2: { assembleStiffnessMatrix<2>(*it, ghost_type); break; }
       case 3: { assembleStiffnessMatrix<3>(*it, ghost_type); break; }
       }
     }
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template<UInt dim>
 void Material::assembleStiffnessMatrix(const ElementType & type,
 				       GhostType ghost_type) {
   AKANTU_DEBUG_IN();
 
   SparseMatrix & K = const_cast<SparseMatrix &>(model->getStiffnessMatrix());
 
   const Array<Real> & shapes_derivatives = model->getFEEngine().getShapesDerivatives(type,
 										ghost_type);
 
   Array<UInt> & elem_filter = element_filter(type, ghost_type);
   Array<Real> & gradu_vect = gradu(type, ghost_type);
 
   UInt nb_element           = elem_filter.getSize();
   UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
   UInt nb_quadrature_points = model->getFEEngine().getNbQuadraturePoints(type,
 								    ghost_type);
 
   gradu_vect.resize(nb_quadrature_points * nb_element);
 
   model->getFEEngine().gradientOnQuadraturePoints(model->getDisplacement(),
 					     gradu_vect, dim, type, ghost_type,
 					     elem_filter);
 
   UInt tangent_size = getTangentStiffnessVoigtSize(dim);
 
   Array<Real> * tangent_stiffness_matrix =
     new Array<Real>(nb_element*nb_quadrature_points, tangent_size * tangent_size,
 		    "tangent_stiffness_matrix");
 
   tangent_stiffness_matrix->clear();
 
   computeTangentModuli(type, *tangent_stiffness_matrix, ghost_type);
 
   Array<Real> * shapesd_filtered =
     new Array<Real>(0, dim * nb_nodes_per_element, "filtered shapesd");
 
   FEEngine::filterElementalData(model->getFEEngine().getMesh(), shapes_derivatives,
 			   *shapesd_filtered, type, ghost_type, elem_filter);
 
 
   /// compute @f$\mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$
   UInt bt_d_b_size = dim * nb_nodes_per_element;
 
   Array<Real> * bt_d_b = new Array<Real>(nb_element * nb_quadrature_points,
 					 bt_d_b_size * bt_d_b_size,
 					 "B^t*D*B");
 
   Matrix<Real> B(tangent_size, dim * nb_nodes_per_element);
   Matrix<Real> Bt_D(dim * nb_nodes_per_element, tangent_size);
 
   Array<Real>::matrix_iterator shapes_derivatives_filtered_it =
     shapesd_filtered->begin(dim, nb_nodes_per_element);
 
   Array<Real>::matrix_iterator Bt_D_B_it
     = bt_d_b->begin(dim*nb_nodes_per_element, dim*nb_nodes_per_element);
 
   Array<Real>::matrix_iterator D_it
     = tangent_stiffness_matrix->begin(tangent_size, tangent_size);
 
   Array<Real>::matrix_iterator D_end
     = tangent_stiffness_matrix->end  (tangent_size, tangent_size);
 
 
   for(; D_it != D_end; ++D_it, ++Bt_D_B_it, ++shapes_derivatives_filtered_it) {
     Matrix<Real> & D = *D_it;
     Matrix<Real> & Bt_D_B = *Bt_D_B_it;
 
     VoigtHelper<dim>::transferBMatrixToSymVoigtBMatrix(
 						       *shapes_derivatives_filtered_it, B, nb_nodes_per_element);
     Bt_D.mul<true, false>(B, D);
     Bt_D_B.mul<false, false>(Bt_D, B);
   }
 
   delete tangent_stiffness_matrix;
   delete shapesd_filtered;
 
   /// compute @f$ k_e = \int_e \mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$
   Array<Real> * K_e = new Array<Real>(nb_element,
 				      bt_d_b_size * bt_d_b_size,
 				      "K_e");
 
   model->getFEEngine().integrate(*bt_d_b, *K_e,
 			    bt_d_b_size * bt_d_b_size,
 			    type, ghost_type,
 			    elem_filter);
 
   delete bt_d_b;
 
   model->getFEEngine().assembleMatrix(*K_e, K, spatial_dimension, type, ghost_type,
 				 elem_filter);
   delete K_e;
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template<UInt dim>
 void Material::assembleStiffnessMatrixNL(const ElementType & type,
 					 GhostType ghost_type) {
   AKANTU_DEBUG_IN();
 
   SparseMatrix & K = const_cast<SparseMatrix &> (model->getStiffnessMatrix());
 
   const Array<Real> & shapes_derivatives = model->getFEEngine().getShapesDerivatives(type, ghost_type);
 
   Array<UInt> & elem_filter = element_filter(type, ghost_type);
   //Array<Real> & gradu_vect = delta_gradu(type, ghost_type);
 
   UInt nb_element = elem_filter.getSize();
   UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
   UInt nb_quadrature_points = model->getFEEngine().getNbQuadraturePoints(type, ghost_type);
 
   //gradu_vect.resize(nb_quadrature_points * nb_element);
 
   // model->getFEEngine().gradientOnQuadraturePoints(model->getIncrement(), gradu_vect,
   //        dim, type, ghost_type, &elem_filter);
 
   Array<Real> * shapes_derivatives_filtered = new Array<Real > (nb_element * nb_quadrature_points,
 								dim * nb_nodes_per_element,
 								"shapes derivatives filtered");
 
 
   Array<Real>::const_matrix_iterator shapes_derivatives_it = shapes_derivatives.begin(spatial_dimension,
 											       nb_nodes_per_element);
 
   Array<Real>::matrix_iterator shapes_derivatives_filtered_it = shapes_derivatives_filtered->begin(spatial_dimension,
 													    nb_nodes_per_element);
   UInt * elem_filter_val = elem_filter.storage();
   for (UInt e = 0; e < nb_element; ++e, ++elem_filter_val)
     for (UInt q = 0; q < nb_quadrature_points; ++q, ++shapes_derivatives_filtered_it)
       *shapes_derivatives_filtered_it = shapes_derivatives_it[*elem_filter_val * nb_quadrature_points + q];
 
   /// compute @f$\mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$
   UInt bt_s_b_size = dim * nb_nodes_per_element;
 
   Array<Real> * bt_s_b = new Array<Real > (nb_element * nb_quadrature_points,
 					   bt_s_b_size * bt_s_b_size,
 					   "B^t*D*B");
 
   UInt piola_matrix_size = getCauchyStressMatrixSize(dim);
 
   Matrix<Real> B(piola_matrix_size, bt_s_b_size);
   Matrix<Real> Bt_S(bt_s_b_size, piola_matrix_size);
   Matrix<Real> S(piola_matrix_size, piola_matrix_size);
 
   shapes_derivatives_filtered_it = shapes_derivatives_filtered->begin(spatial_dimension, nb_nodes_per_element);
 
   Array<Real>::matrix_iterator Bt_S_B_it = bt_s_b->begin(bt_s_b_size,
 								  bt_s_b_size);
   Array<Real>::matrix_iterator Bt_S_B_end = bt_s_b->end(bt_s_b_size,
 								 bt_s_b_size);
 
   Array<Real>::matrix_iterator piola_it = piola_kirchhoff_2(type, ghost_type).begin(dim, dim);
 
   for (; Bt_S_B_it != Bt_S_B_end; ++Bt_S_B_it, ++shapes_derivatives_filtered_it, ++piola_it) {
     Matrix<Real> & Bt_S_B = *Bt_S_B_it;
     Matrix<Real> & Piola_kirchhoff_matrix = *piola_it;
 
     setCauchyStressMatrix< dim >(Piola_kirchhoff_matrix, S);
     VoigtHelper<dim>::transferBMatrixToBNL(*shapes_derivatives_filtered_it, B, nb_nodes_per_element);
     Bt_S.mul < true, false > (B, S);
     Bt_S_B.mul < false, false > (Bt_S, B);
   }
 
   delete shapes_derivatives_filtered;
 
   /// compute @f$ k_e = \int_e \mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$
   Array<Real> * K_e = new Array<Real > (nb_element,
 					bt_s_b_size * bt_s_b_size,
 					"K_e");
 
   model->getFEEngine().integrate(*bt_s_b, *K_e,
 			    bt_s_b_size * bt_s_b_size,
 			    type, ghost_type,
 			    elem_filter);
 
   delete bt_s_b;
 
   model->getFEEngine().assembleMatrix(*K_e, K, spatial_dimension, type, ghost_type, elem_filter);
   delete K_e;
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template<UInt dim>
 void Material::assembleStiffnessMatrixL2(const ElementType & type,
 					 GhostType ghost_type) {
   AKANTU_DEBUG_IN();
 
   SparseMatrix & K = const_cast<SparseMatrix &> (model->getStiffnessMatrix());
 
   const Array<Real> & shapes_derivatives = model->getFEEngine().getShapesDerivatives(type, ghost_type);
 
   Array<UInt> & elem_filter = element_filter(type, ghost_type);
   Array<Real> & gradu_vect = gradu(type, ghost_type);
 
   UInt nb_element = elem_filter.getSize();
   UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
   UInt nb_quadrature_points = model->getFEEngine().getNbQuadraturePoints(type, ghost_type);
 
   gradu_vect.resize(nb_quadrature_points * nb_element);
 
   model->getFEEngine().gradientOnQuadraturePoints(model->getDisplacement(), gradu_vect,
 					     dim, type, ghost_type, elem_filter);
 
   UInt tangent_size = getTangentStiffnessVoigtSize(dim);
 
   Array<Real> * tangent_stiffness_matrix =
     new Array<Real > (nb_element*nb_quadrature_points, tangent_size * tangent_size,
 		      "tangent_stiffness_matrix");
 
   tangent_stiffness_matrix->clear();
 
   computeTangentModuli(type, *tangent_stiffness_matrix, ghost_type);
 
 
   Array<Real> * shapes_derivatives_filtered = new Array<Real > (nb_element * nb_quadrature_points,
 								dim * nb_nodes_per_element,
 								"shapes derivatives filtered");
 
 
   Array<Real>::const_matrix_iterator shapes_derivatives_it = shapes_derivatives.begin(spatial_dimension,
 											       nb_nodes_per_element);
 
   Array<Real>::matrix_iterator shapes_derivatives_filtered_it = shapes_derivatives_filtered->begin(spatial_dimension,
 													    nb_nodes_per_element);
   UInt * elem_filter_val = elem_filter.storage();
   for (UInt e = 0; e < nb_element; ++e, ++elem_filter_val)
     for (UInt q = 0; q < nb_quadrature_points; ++q, ++shapes_derivatives_filtered_it)
       *shapes_derivatives_filtered_it = shapes_derivatives_it[*elem_filter_val * nb_quadrature_points + q];
 
   /// compute @f$\mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$
   UInt bt_d_b_size = dim * nb_nodes_per_element;
 
   Array<Real> * bt_d_b = new Array<Real > (nb_element * nb_quadrature_points,
 					   bt_d_b_size * bt_d_b_size,
 					   "B^t*D*B");
 
   Matrix<Real> B(tangent_size, dim * nb_nodes_per_element);
   Matrix<Real> B2(tangent_size, dim * nb_nodes_per_element);
   Matrix<Real> Bt_D(dim * nb_nodes_per_element, tangent_size);
 
   shapes_derivatives_filtered_it = shapes_derivatives_filtered->begin(spatial_dimension, nb_nodes_per_element);
 
   Array<Real>::matrix_iterator Bt_D_B_it = bt_d_b->begin(dim*nb_nodes_per_element,
 								  dim * nb_nodes_per_element);
 
   Array<Real>::matrix_iterator grad_u_it = gradu_vect.begin(dim, dim);
 
   Array<Real>::matrix_iterator D_it = tangent_stiffness_matrix->begin(tangent_size,
 									       tangent_size);
   Array<Real>::matrix_iterator D_end = tangent_stiffness_matrix->end(tangent_size,
 									      tangent_size);
 
 
   for (; D_it != D_end; ++D_it, ++Bt_D_B_it, ++shapes_derivatives_filtered_it, ++grad_u_it) {
     Matrix<Real> & grad_u = *grad_u_it;
     Matrix<Real> & D = *D_it;
     Matrix<Real> & Bt_D_B = *Bt_D_B_it;
 
     //transferBMatrixToBL1<dim > (*shapes_derivatives_filtered_it, B, nb_nodes_per_element);
     VoigtHelper<dim>::transferBMatrixToSymVoigtBMatrix(*shapes_derivatives_filtered_it, B, nb_nodes_per_element);
     VoigtHelper<dim>::transferBMatrixToBL2(*shapes_derivatives_filtered_it, grad_u, B2, nb_nodes_per_element);
     B += B2;
     Bt_D.mul < true, false > (B, D);
     Bt_D_B.mul < false, false > (Bt_D, B);
   }
 
   delete tangent_stiffness_matrix;
   delete shapes_derivatives_filtered;
 
   /// compute @f$ k_e = \int_e \mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$
   Array<Real> * K_e = new Array<Real > (nb_element,
 					bt_d_b_size * bt_d_b_size,
 					"K_e");
 
   model->getFEEngine().integrate(*bt_d_b, *K_e,
 			    bt_d_b_size * bt_d_b_size,
 			    type, ghost_type,
 			    elem_filter);
 
   delete bt_d_b;
 
   model->getFEEngine().assembleMatrix(*K_e, K, spatial_dimension, type, ghost_type, elem_filter);
   delete K_e;
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template<UInt dim>
 void Material::assembleResidual(GhostType ghost_type){
 
   AKANTU_DEBUG_IN();
 
   Array<Real> & residual = const_cast<Array<Real> &> (model->getResidual());
 
   Mesh & mesh = model->getFEEngine().getMesh();
   Mesh::type_iterator it = element_filter.firstType(dim, ghost_type);
   Mesh::type_iterator last_type = element_filter.lastType(dim, ghost_type);
   for (; it != last_type; ++it) {
     const Array<Real> & shapes_derivatives = model->getFEEngine().getShapesDerivatives(*it, ghost_type);
 
     Array<UInt> & elem_filter = element_filter(*it, ghost_type);
 
     UInt size_of_shapes_derivatives = shapes_derivatives.getNbComponent();
     UInt nb_element = elem_filter.getSize();
     UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(*it);
     UInt nb_quadrature_points = model->getFEEngine().getNbQuadraturePoints(*it, ghost_type);
 
     Array<Real> * shapesd_filtered =
       new Array<Real>(0, size_of_shapes_derivatives, "filtered shapesd");
 
     FEEngine::filterElementalData(mesh, shapes_derivatives, *shapesd_filtered,
 			     *it, ghost_type, elem_filter);
 
     Array<Real>::matrix_iterator shapes_derivatives_filtered_it = shapesd_filtered->begin(dim,
 												   nb_nodes_per_element);
 
     //Set stress vectors
 
     UInt stress_size = getTangentStiffnessVoigtSize(dim);
 
     //Set matrices B and BNL*
     UInt bt_s_size = dim * nb_nodes_per_element;
 
     Array<Real> * bt_s = new Array<Real > (nb_element * nb_quadrature_points,
 					   bt_s_size,
 					   "B^t*S");
 
     Array<Real>::matrix_iterator grad_u_it = this->gradu(*it, ghost_type).begin(dim, dim);
 
     Array<Real>::matrix_iterator grad_u_end = this->gradu(*it, ghost_type).end(dim, dim);
 
     Array<Real>::matrix_iterator stress_it = this->piola_kirchhoff_2(*it, ghost_type).begin(dim, dim);
 
     shapes_derivatives_filtered_it = shapesd_filtered->begin(dim, nb_nodes_per_element);
 
     Array<Real>::matrix_iterator bt_s_it = bt_s->begin(bt_s_size, 1);
 
     Matrix<Real> S_vect(stress_size, 1);
     Matrix<Real> B_tensor(stress_size, bt_s_size);
     Matrix<Real> B2_tensor(stress_size, bt_s_size);
 
     for (; grad_u_it != grad_u_end; ++grad_u_it, ++stress_it, ++shapes_derivatives_filtered_it, ++bt_s_it) {
       Matrix<Real> & grad_u = *grad_u_it;
       Matrix<Real> & r_it = *bt_s_it;
       Matrix<Real> & S_it = *stress_it;
 
       SetCauchyStressArray <dim> (S_it, S_vect);
       VoigtHelper<dim>::transferBMatrixToSymVoigtBMatrix(*shapes_derivatives_filtered_it, B_tensor, nb_nodes_per_element);
       VoigtHelper<dim>::transferBMatrixToBL2(*shapes_derivatives_filtered_it, grad_u, B2_tensor, nb_nodes_per_element);
 
       B_tensor += B2_tensor;
 
       r_it.mul < true, false > (B_tensor, S_vect);
     }
 
     delete shapesd_filtered;
 
     /// compute @f$ k_e = \int_e \mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$
     Array<Real> * r_e = new Array<Real > (nb_element,
 					  bt_s_size, "r_e");
 
     model->getFEEngine().integrate(*bt_s, *r_e,
 			      bt_s_size,
 			      *it, ghost_type,
 			      elem_filter);
 
     delete bt_s;
 
     model->getFEEngine().assembleArray(*r_e, residual,
 				  model->getDOFSynchronizer().getLocalDOFEquationNumbers(),
 				  residual.getNbComponent(),
 				  *it, ghost_type, elem_filter, -1);
 
     delete r_e;
 
   }
   AKANTU_DEBUG_OUT();
 
 }
 
 /* -------------------------------------------------------------------------- */
 void Material::computeAllStressesFromTangentModuli(GhostType ghost_type) {
   AKANTU_DEBUG_IN();
 
   UInt spatial_dimension = model->getSpatialDimension();
 
   Mesh::type_iterator it = element_filter.firstType(spatial_dimension, ghost_type);
   Mesh::type_iterator last_type = element_filter.lastType(spatial_dimension, ghost_type);
   for(; it != last_type; ++it) {
     switch(spatial_dimension) {
     case 1: { computeAllStressesFromTangentModuli<1>(*it, ghost_type); break; }
     case 2: { computeAllStressesFromTangentModuli<2>(*it, ghost_type); break; }
     case 3: { computeAllStressesFromTangentModuli<3>(*it, ghost_type); break; }
     }
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template<UInt dim>
 void Material::computeAllStressesFromTangentModuli(const ElementType & type,
 						   GhostType ghost_type) {
   AKANTU_DEBUG_IN();
 
   const Array<Real> & shapes_derivatives = model->getFEEngine().getShapesDerivatives(type, ghost_type);
   Array<UInt> & elem_filter = element_filter(type, ghost_type);
   Array<Real> & gradu_vect = gradu(type, ghost_type);
 
   UInt nb_element                 = elem_filter.getSize();
   UInt nb_nodes_per_element       = Mesh::getNbNodesPerElement(type);
   UInt nb_quadrature_points       = model->getFEEngine().getNbQuadraturePoints(type, ghost_type);
 
 
   gradu_vect.resize(nb_quadrature_points * nb_element);
 
   Array<Real> & disp = model->getDisplacement();
 
   model->getFEEngine().gradientOnQuadraturePoints(disp, gradu_vect,
 					     dim, type, ghost_type, elem_filter);
 
   UInt tangent_moduli_size = getTangentStiffnessVoigtSize(dim);
 
   Array<Real> * tangent_moduli_tensors =
     new Array<Real>(nb_element*nb_quadrature_points, tangent_moduli_size * tangent_moduli_size,
 		    "tangent_moduli_tensors");
 
   tangent_moduli_tensors->clear();
   computeTangentModuli(type, *tangent_moduli_tensors, ghost_type);
 
   Array<Real> * shapesd_filtered =
     new Array<Real>(0, dim* nb_nodes_per_element, "filtered shapesd");
 
   FEEngine::filterElementalData(model->getFEEngine().getMesh(), shapes_derivatives, *shapesd_filtered,
 			   type, ghost_type, elem_filter);
 
   Array<Real> filtered_u(nb_element, nb_nodes_per_element * spatial_dimension);
 
   FEEngine::extractNodalToElementField(model->getFEEngine().getMesh(), disp, filtered_u,
 				  type, ghost_type, elem_filter);
 
   /// compute @f$\mathbf{D} \mathbf{B} \mathbf{u}@f$
   Array<Real>::matrix_iterator shapes_derivatives_filtered_it =
     shapesd_filtered->begin(dim, nb_nodes_per_element);
 
   Array<Real>::matrix_iterator D_it  = tangent_moduli_tensors->begin(tangent_moduli_size,
 								     tangent_moduli_size);
   Array<Real>::matrix_iterator sigma_it  = stress(type, ghost_type).begin(spatial_dimension,
 									  spatial_dimension);
   Array<Real>::vector_iterator u_it = filtered_u.begin(spatial_dimension * nb_nodes_per_element);
 
   Matrix<Real> B(tangent_moduli_size, spatial_dimension * nb_nodes_per_element);
   Vector<Real> Bu(tangent_moduli_size);
   Vector<Real> DBu(tangent_moduli_size);
 
   for (UInt e = 0; e < nb_element; ++e, ++u_it) {
     for (UInt q = 0; q < nb_quadrature_points; ++q, ++D_it, ++shapes_derivatives_filtered_it, ++sigma_it) {
       Vector<Real> & u = *u_it;
       Matrix<Real> & sigma = *sigma_it;
       Matrix<Real> & D = *D_it;
 
       VoigtHelper<dim>::transferBMatrixToSymVoigtBMatrix(*shapes_derivatives_filtered_it, B, nb_nodes_per_element);
 
       Bu.mul<false>(B, u);
       DBu.mul<false>(D, Bu);
 
       // Voigt notation to full symmetric tensor
       for (UInt i = 0; i < dim; ++i) sigma(i, i) = DBu(i);
       if(dim == 2) {
 	sigma(0,1) = sigma(1,0) = DBu(2);
       } else if(dim == 3) {
 	sigma(1,2) = sigma(2,1) = DBu(3);
 	sigma(0,2) = sigma(2,0) = DBu(4);
 	sigma(0,1) = sigma(1,0) = DBu(5);
       }
     }
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 
 /* -------------------------------------------------------------------------- */
 void Material::computePotentialEnergyByElements() {
   AKANTU_DEBUG_IN();
 
   Mesh::type_iterator it = element_filter.firstType(spatial_dimension);
   Mesh::type_iterator last_type = element_filter.lastType(spatial_dimension);
   for(; it != last_type; ++it) {
     computePotentialEnergy(*it);
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 void Material::computePotentialEnergy(ElementType el_type, GhostType ghost_type) {
   AKANTU_DEBUG_IN();
 
   if(!potential_energy.exists(el_type, ghost_type)) {
       UInt nb_element = element_filter(el_type, ghost_type).getSize();
       UInt nb_quadrature_points = model->getFEEngine().getNbQuadraturePoints(el_type, _not_ghost);
 
       potential_energy.alloc(nb_element * nb_quadrature_points, 1,
 			     el_type, ghost_type);
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 
 /* -------------------------------------------------------------------------- */
 Real Material::getPotentialEnergy() {
   AKANTU_DEBUG_IN();
   Real epot = 0.;
 
   computePotentialEnergyByElements();
 
   /// integrate the potential energy for each type of elements
   Mesh::type_iterator it = element_filter.firstType(spatial_dimension);
   Mesh::type_iterator last_type = element_filter.lastType(spatial_dimension);
   for(; it != last_type; ++it) {
 
     epot += model->getFEEngine().integrate(potential_energy(*it, _not_ghost), *it,
 				      _not_ghost, element_filter(*it, _not_ghost));
   }
 
   AKANTU_DEBUG_OUT();
   return epot;
 }
 
 /* -------------------------------------------------------------------------- */
 Real Material::getPotentialEnergy(ElementType & type, UInt index) {
   AKANTU_DEBUG_IN();
   Real epot = 0.;
 
   Vector<Real> epot_on_quad_points(model->getFEEngine().getNbQuadraturePoints(type));
 
   computePotentialEnergyByElement(type, index, epot_on_quad_points);
 
   epot = model->getFEEngine().integrate(epot_on_quad_points, type, element_filter(type)(index));
 
   AKANTU_DEBUG_OUT();
   return epot;
 }
 
 /* -------------------------------------------------------------------------- */
 Real Material::getEnergy(std::string type) {
   AKANTU_DEBUG_IN();
   if(type == "potential") return getPotentialEnergy();
   AKANTU_DEBUG_OUT();
   return 0.;
 }
 
 /* -------------------------------------------------------------------------- */
 Real Material::getEnergy(std::string energy_id, ElementType type, UInt index) {
   AKANTU_DEBUG_IN();
   if(energy_id == "potential") return getPotentialEnergy(type, index);
   AKANTU_DEBUG_OUT();
   return 0.;
 }
 
 /* -------------------------------------------------------------------------- */
 void Material::computeQuadraturePointsCoordinates(ElementTypeMapArray<Real> & quadrature_points_coordinates,
 						  const GhostType & ghost_type) const {
   AKANTU_DEBUG_IN();
 
   const Mesh & mesh = this->model->getMesh();
 
   Array<Real> nodes_coordinates(mesh.getNodes(), true);
   nodes_coordinates += this->model->getDisplacement();
 
   Mesh::type_iterator it = this->element_filter.firstType(spatial_dimension, ghost_type);
   Mesh::type_iterator last_type = this->element_filter.lastType(spatial_dimension, ghost_type);
   for(; it != last_type; ++it) {
     const Array<UInt> & elem_filter = this->element_filter(*it, ghost_type);
 
     UInt nb_element  = elem_filter.getSize();
     UInt nb_tot_quad = this->model->getFEEngine().getNbQuadraturePoints(*it, ghost_type) * nb_element;
 
     Array<Real> & quads = quadrature_points_coordinates(*it, ghost_type);
     quads.resize(nb_tot_quad);
 
     this->model->getFEEngine().interpolateOnQuadraturePoints(nodes_coordinates,
 							     quads, spatial_dimension,
 							     *it, ghost_type, elem_filter);
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 void Material::initElementalFieldInterpolation(const ElementTypeMapArray<Real> & interpolation_points_coordinates) {
   AKANTU_DEBUG_IN();
   const Mesh & mesh = model->getFEEngine().getMesh();
 
   ElementTypeMapArray<Real> quadrature_points_coordinates("quadrature_points_coordinates_for_interpolation", getID());
   mesh.initElementTypeMapArray(quadrature_points_coordinates, spatial_dimension, spatial_dimension);
 
   for (ghost_type_t::iterator gt = ghost_type_t::begin();
        gt != ghost_type_t::end(); ++gt) {
 
     GhostType ghost_type = *gt;
 
     computeQuadraturePointsCoordinates(quadrature_points_coordinates, ghost_type);
 
     Mesh::type_iterator it   = element_filter.firstType(spatial_dimension, ghost_type);
     Mesh::type_iterator last = element_filter.lastType(spatial_dimension, ghost_type);
     for (; it != last; ++it) {
       ElementType type = *it;
       UInt nb_element = mesh.getNbElement(type, ghost_type);
       if (nb_element == 0) continue;
 
       const Array<Real> & interp_points_coord = interpolation_points_coordinates(type, ghost_type);
       UInt nb_interpolation_points_per_elem = interp_points_coord.getSize() / nb_element;
 
       AKANTU_DEBUG_ASSERT(interp_points_coord.getSize() % nb_element == 0,
 			  "Number of interpolation points is wrong");
 
 #define AKANTU_INIT_INTERPOLATE_ELEMENTAL_FIELD(type)			\
       initElementalFieldInterpolation<type>(quadrature_points_coordinates(type, ghost_type), \
 					    interp_points_coord,	\
 					    nb_interpolation_points_per_elem, \
 					    ghost_type)			\
 
       AKANTU_BOOST_REGULAR_ELEMENT_SWITCH(AKANTU_INIT_INTERPOLATE_ELEMENTAL_FIELD);
 #undef AKANTU_INIT_INTERPOLATE_ELEMENTAL_FIELD
     }
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementType type>
 void Material::initElementalFieldInterpolation(const Array<Real> & quad_coordinates,
 					       const Array<Real> & interpolation_points_coordinates,
 					       const UInt nb_interpolation_points_per_elem,
 					       const GhostType ghost_type) {
   AKANTU_DEBUG_IN();
   Array<UInt> & elem_fil = element_filter(type, ghost_type);
   UInt nb_element = elem_fil.getSize();
   UInt nb_quad_per_element = model->getFEEngine().getNbQuadraturePoints(type, ghost_type);
 
   if(!interpolation_inverse_coordinates.exists(type, ghost_type))
     interpolation_inverse_coordinates.alloc(nb_element,
 					    nb_quad_per_element*nb_quad_per_element,
 					    type, ghost_type);
   else
     interpolation_inverse_coordinates(type, ghost_type).resize(nb_element);
 
   if(!interpolation_points_matrices.exists(type, ghost_type))
     interpolation_points_matrices.alloc(nb_element,
 					nb_interpolation_points_per_elem * nb_quad_per_element,
 					type, ghost_type);
   else
     interpolation_points_matrices(type, ghost_type).resize(nb_element);
 
   Array<Real> & interp_inv_coord = interpolation_inverse_coordinates(type, ghost_type);
   Array<Real> & interp_points_mat = interpolation_points_matrices(type, ghost_type);
 
   Matrix<Real> quad_coord_matrix(nb_quad_per_element, nb_quad_per_element);
 
   Array<Real>::const_matrix_iterator quad_coords_it =
     quad_coordinates.begin_reinterpret(spatial_dimension,
 				       nb_quad_per_element,
 				       nb_element);
 
   Array<Real>::const_matrix_iterator points_coords_begin =
     interpolation_points_coordinates.begin_reinterpret(spatial_dimension,
 						       nb_interpolation_points_per_elem,
 						       interpolation_points_coordinates.getSize() / nb_interpolation_points_per_elem);
 
   Array<Real>::matrix_iterator inv_quad_coord_it =
     interp_inv_coord.begin(nb_quad_per_element, nb_quad_per_element);
 
   Array<Real>::matrix_iterator inv_points_mat_it =
     interp_points_mat.begin(nb_interpolation_points_per_elem, nb_quad_per_element);
 
   /// loop over the elements of the current material and element type
   for (UInt el = 0; el < nb_element; ++el, ++inv_quad_coord_it,
 	 ++inv_points_mat_it, ++quad_coords_it) {
     /// matrix containing the quadrature points coordinates
     const Matrix<Real> & quad_coords = *quad_coords_it;
     /// matrix to store the matrix inversion result
     Matrix<Real> & inv_quad_coord_matrix = *inv_quad_coord_it;
 
     /// insert the quad coordinates in a matrix compatible with the interpolation
     buildElementalFieldInterpolationCoodinates<type>(quad_coords,
 						     quad_coord_matrix);
 
     /// invert the interpolation matrix
     inv_quad_coord_matrix.inverse(quad_coord_matrix);
 
 
     /// matrix containing the interpolation points coordinates
     const Matrix<Real> & points_coords = points_coords_begin[elem_fil(el)];
     /// matrix to store the interpolation points coordinates
     /// compatible with these functions
     Matrix<Real> & inv_points_coord_matrix = *inv_points_mat_it;
 
     /// insert the quad coordinates in a matrix compatible with the interpolation
     buildElementalFieldInterpolationCoodinates<type>(points_coords,
 						     inv_points_coord_matrix);
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 void Material::interpolateStress(ElementTypeMapArray<Real> & result,
 				 const GhostType ghost_type) {
   interpolateElementalField(stress, result, ghost_type);
 }
 
 /* -------------------------------------------------------------------------- */
 void Material::interpolateElementalField(const ElementTypeMapArray<Real> & field,
 					 ElementTypeMapArray<Real> & result,
 					 const GhostType ghost_type) {
   AKANTU_DEBUG_IN();
 
   Mesh::type_iterator it   = element_filter.firstType(spatial_dimension, ghost_type);
   Mesh::type_iterator last = element_filter.lastType(spatial_dimension, ghost_type);
   for (; it != last; ++it) {
     ElementType type = *it;
 
     Array<UInt> & elem_fil = element_filter(type, ghost_type);
     UInt nb_element = elem_fil.getSize();
     UInt nb_quad_per_element = model->getFEEngine().getNbQuadraturePoints(type, ghost_type);
 
     const Array<Real> & field_vec = field(type, ghost_type);
     Array<Real> & result_vec = result(type, ghost_type);
 
     Matrix<Real> coefficients(nb_quad_per_element, field_vec.getNbComponent());
 
     const Array<Real> & interp_inv_coord = interpolation_inverse_coordinates(type, ghost_type);
     const Array<Real> & interp_points_coord = interpolation_points_matrices(type, ghost_type);
 
     UInt nb_interpolation_points_per_elem
       = interp_points_coord.getNbComponent() / nb_quad_per_element;
 
     Array<Real>::const_matrix_iterator field_it
       = field_vec.begin_reinterpret(field_vec.getNbComponent(),
 				    nb_quad_per_element,
 				    nb_element);
 
     Array<Real>::const_matrix_iterator interpolation_points_coordinates_it =
       interp_points_coord.begin(nb_interpolation_points_per_elem, nb_quad_per_element);
 
     Array<Real>::matrix_iterator result_begin
       = result_vec.begin_reinterpret(field_vec.getNbComponent(),
 				     nb_interpolation_points_per_elem,
 				     result_vec.getSize() / nb_interpolation_points_per_elem);
 
     Array<Real>::const_matrix_iterator inv_quad_coord_it =
       interp_inv_coord.begin(nb_quad_per_element, nb_quad_per_element);
 
     /// loop over the elements of the current material and element type
     for (UInt el = 0; el < nb_element;
 	 ++el, ++field_it, ++inv_quad_coord_it, ++interpolation_points_coordinates_it) {
       /**
        * matrix containing the inversion of the quadrature points'
        * coordinates
        */
       const Matrix<Real> & inv_quad_coord_matrix = *inv_quad_coord_it;
 
       /**
        * multiply it by the field values over quadrature points to get
        * the interpolation coefficients
        */
       coefficients.mul<false, true>(inv_quad_coord_matrix, *field_it);
 
       /// matrix containing the points' coordinates
       const Matrix<Real> & coord = *interpolation_points_coordinates_it;
 
       /// multiply the coordinates matrix by the coefficients matrix and store the result
       Matrix<Real> res(result_begin[elem_fil(el)]);
       res.mul<true, true>(coefficients, coord);
     }
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 void Material::interpolateStressOnFacets(ElementTypeMapArray<Real> & result,
 					 const GhostType ghost_type) {
   interpolateElementalFieldOnFacets(stress, result, ghost_type);
 }
 
 /* -------------------------------------------------------------------------- */
 void Material::interpolateElementalFieldOnFacets(const ElementTypeMapArray<Real> & field,
 						 ElementTypeMapArray<Real> & result,
 						 const GhostType ghost_type) {
   AKANTU_DEBUG_IN();
 
   UInt sp2 = spatial_dimension * spatial_dimension;
 
   const Mesh & mesh = this->model->getMesh();
   const Mesh & mesh_facets = mesh.getMeshFacets();
 
   Mesh::type_iterator it   = element_filter.firstType(spatial_dimension, ghost_type);
   Mesh::type_iterator last = element_filter.lastType(spatial_dimension, ghost_type);
   for (; it != last; ++it) {
     ElementType type = *it;
 
     Array<UInt> & elem_fil = element_filter(type, ghost_type);
     UInt nb_element = elem_fil.getSize();
     UInt nb_quad_per_element = model->getFEEngine().getNbQuadraturePoints(type, ghost_type);
 
     const Array<Real> & field_vec = field(type, ghost_type);
 
     Matrix<Real> coefficients(nb_quad_per_element, field_vec.getNbComponent());
 
     const Array<Real> & interp_inv_coord = interpolation_inverse_coordinates(type, ghost_type);
     const Array<Real> & interp_points_coord = interpolation_points_matrices(type, ghost_type);
 
     UInt nb_interpolation_points_per_elem
       = interp_points_coord.getNbComponent() / nb_quad_per_element;
 
     Array<Real>::const_matrix_iterator field_it
       = field_vec.begin_reinterpret(field_vec.getNbComponent(),
 				    nb_quad_per_element,
 				    nb_element);
 
     Array<Real>::const_matrix_iterator interpolation_points_coordinates_it =
       interp_points_coord.begin(nb_interpolation_points_per_elem, nb_quad_per_element);
 
     Array<Real>::const_matrix_iterator inv_quad_coord_it =
       interp_inv_coord.begin(nb_quad_per_element, nb_quad_per_element);
 
     Matrix<Real> result_tmp(sp2, nb_interpolation_points_per_elem);
 
     const Array<Element> & facet_to_element =
       mesh_facets.getSubelementToElement(type, ghost_type);
     ElementType type_facet = Mesh::getFacetType(type);
     UInt nb_facet_per_elem = facet_to_element.getNbComponent();
     UInt nb_quad_per_facet = nb_interpolation_points_per_elem / nb_facet_per_elem;
     Element element_for_comparison(type, 0, ghost_type);
     const Array< std::vector<Element> > * element_to_facet = NULL;
     GhostType current_ghost_type = _casper;
     Array<Real> * result_vec = NULL;
 
     /// loop over the elements of the current material and element type
     for (UInt el = 0; el < nb_element;
 	 ++el, ++field_it, ++inv_quad_coord_it, ++interpolation_points_coordinates_it) {
       /**
        * matrix containing the inversion of the quadrature points'
        * coordinates
        */
       const Matrix<Real> & inv_quad_coord_matrix = *inv_quad_coord_it;
 
       /**
        * multiply it by the field values over quadrature points to get
        * the interpolation coefficients
        */
       coefficients.mul<false, true>(inv_quad_coord_matrix, *field_it);
 
       /// matrix containing the points' coordinates
       const Matrix<Real> & coord = *interpolation_points_coordinates_it;
 
       /// multiply the coordinates matrix by the coefficients matrix and store the result
       result_tmp.mul<true, true>(coefficients, coord);
 
       UInt global_el = elem_fil(el);
       element_for_comparison.element = global_el;
       for (UInt f = 0; f < nb_facet_per_elem; ++f) {
 	Element facet_elem = facet_to_element(global_el, f);
 	UInt global_facet = facet_elem.element;
 	if (facet_elem.ghost_type != current_ghost_type) {
 	  current_ghost_type = facet_elem.ghost_type;
 	  element_to_facet = &mesh_facets.getElementToSubelement(type_facet,
 								 current_ghost_type);
 	  result_vec = &result(type_facet, current_ghost_type);
 	}
 
 	bool is_second_element = (*element_to_facet)(global_facet)[0] != element_for_comparison;
 
 	for (UInt q = 0; q < nb_quad_per_facet; ++q) {
 	  Vector<Real> result_local(result_vec->storage()
 				    + (global_facet * nb_quad_per_facet + q) * result_vec->getNbComponent()
 				    + is_second_element * sp2,
 				    sp2);
 
 	  result_local = result_tmp(f * nb_quad_per_facet + q);
 	}
       }
     }
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 
 const Array<Real> & Material::getArray(const ID & vect_id, const ElementType & type, const GhostType & ghost_type) const {
   std::stringstream sstr;
   std::string ghost_id = "";
   if (ghost_type == _ghost) ghost_id = ":ghost";
   sstr << getID() << ":" << vect_id << ":" << type << ghost_id;
 
   ID fvect_id = sstr.str();
   try {
     return Memory::getArray<Real>(fvect_id);
   } catch(debug::Exception & e) {
     AKANTU_SILENT_EXCEPTION("The material " << name << "(" <<getID() << ") does not contain a vector " << vect_id << "(" << fvect_id << ") [" << e << "]");
   }
 }
 
 /* -------------------------------------------------------------------------- */
 Array<Real> & Material::getArray(const ID & vect_id, const ElementType & type, const GhostType & ghost_type) {
   std::stringstream sstr;
   std::string ghost_id = "";
   if (ghost_type == _ghost) ghost_id = ":ghost";
   sstr << getID() << ":" << vect_id << ":" << type << ghost_id;
 
   ID fvect_id = sstr.str();
   try {
     return Memory::getArray<Real>(fvect_id);
   } catch(debug::Exception & e) {
     AKANTU_SILENT_EXCEPTION("The material " << name << "(" << getID() << ") does not contain a vector " << vect_id << "(" << fvect_id << ") [" << e << "]");
   }
 }
 
 /* -------------------------------------------------------------------------- */
 const InternalField<Real> & Material::getInternal(const ID & int_id) const {
   std::map<ID, InternalField<Real> *>::const_iterator it = internal_vectors_real.find(getID() + ":" + int_id);
   if(it == internal_vectors_real.end()) {
     AKANTU_SILENT_EXCEPTION("The material " << name << "(" << getID()
 		     << ") does not contain an internal "
 		     << int_id << " (" << (getID() + ":" + int_id) << ")");
   }
   return *it->second;
 }
 
 /* -------------------------------------------------------------------------- */
 InternalField<Real> & Material::getInternal(const ID & int_id) {
   std::map<ID, InternalField<Real> *>::iterator it = internal_vectors_real.find(getID() + ":" + int_id);
   if(it == internal_vectors_real.end()) {
     AKANTU_SILENT_EXCEPTION("The material " << name << "(" << getID()
 			    << ") does not contain an internal "
 			    << int_id << " (" << (getID() + ":" + int_id) << ")");
   }
   return *it->second;
 }
 
 /* -------------------------------------------------------------------------- */
 void Material::addElements(const Array<Element> & elements_to_add) {
   AKANTU_DEBUG_IN();
   UInt mat_id = model->getInternalIndexFromID(getID());
   Array<Element>::const_iterator<Element> el_begin = elements_to_add.begin();
   Array<Element>::const_iterator<Element> el_end   = elements_to_add.end();
   for(;el_begin != el_end; ++el_begin) {
     const Element & element = *el_begin;
     Array<UInt> & mat_indexes = model->getMaterialByElement     (element.type, element.ghost_type);
     Array<UInt> & mat_loc_num = model->getMaterialLocalNumbering(element.type, element.ghost_type);
 
     UInt index = this->addElement(element.type, element.element, element.ghost_type);
     mat_indexes(element.element) = mat_id;
     mat_loc_num(element.element) = index;
   }
 
   this->resizeInternals();
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 void Material::removeElements(const Array<Element> & elements_to_remove) {
   AKANTU_DEBUG_IN();
 
   Array<Element>::const_iterator<Element> el_begin = elements_to_remove.begin();
   Array<Element>::const_iterator<Element> el_end   = elements_to_remove.end();
 
   if(el_begin==el_end)
     return;
 
   ElementTypeMapArray<UInt> material_local_new_numbering("remove mat filter elem", getID());
 
   Element element;
   for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) {
     GhostType ghost_type = *gt;
     element.ghost_type = ghost_type;
 
     ElementTypeMapArray<UInt>::type_iterator it  = element_filter.firstType(_all_dimensions, ghost_type, _ek_not_defined);
     ElementTypeMapArray<UInt>::type_iterator end = element_filter.lastType(_all_dimensions, ghost_type, _ek_not_defined);
 
     for(; it != end; ++it) {
       ElementType type = *it;
       element.type = type;
 
       Array<UInt> & elem_filter = this->element_filter(type, ghost_type);
       Array<UInt> & mat_loc_num = this->model->getMaterialLocalNumbering(type, ghost_type);
 
       if(!material_local_new_numbering.exists(type, ghost_type))
 	material_local_new_numbering.alloc(elem_filter.getSize(), 1, type, ghost_type);
       Array<UInt> & mat_renumbering = material_local_new_numbering(type, ghost_type);
 
       UInt nb_element = elem_filter.getSize();
       element.kind=(*el_begin).kind;
       Array<UInt> elem_filter_tmp;
 
       UInt new_id = 0;
       for (UInt el = 0; el < nb_element; ++el) {
 	element.element = elem_filter(el);
 
 	if(std::find(el_begin, el_end, element) == el_end) {
 	  elem_filter_tmp.push_back(element.element);
 
 	  mat_renumbering(el) = new_id;
 	  mat_loc_num(element.element) = new_id;
 	  ++new_id;
 	} else {
 	  mat_renumbering(el) = UInt(-1);
 	}
       }
 
       elem_filter.resize(elem_filter_tmp.getSize());
       elem_filter.copy(elem_filter_tmp);
     }
   }
 
   for (std::map<ID, InternalField<Real> *>::iterator it = internal_vectors_real.begin();
        it != internal_vectors_real.end();
        ++it) it->second->removeQuadraturePoints(material_local_new_numbering);
 
   for (std::map<ID, InternalField<UInt> *>::iterator it = internal_vectors_uint.begin();
        it != internal_vectors_uint.end();
        ++it) it->second->removeQuadraturePoints(material_local_new_numbering);
 
   for (std::map<ID, InternalField<bool> *>::iterator it = internal_vectors_bool.begin();
        it != internal_vectors_bool.end();
        ++it) it->second->removeQuadraturePoints(material_local_new_numbering);
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 void Material::resizeInternals() {
   AKANTU_DEBUG_IN();
   for (std::map<ID, InternalField<Real> *>::iterator it = internal_vectors_real.begin();
        it != internal_vectors_real.end();
        ++it) it->second->resize();
 
   for (std::map<ID, InternalField<UInt> *>::iterator it = internal_vectors_uint.begin();
        it != internal_vectors_uint.end();
        ++it) it->second->resize();
 
   for (std::map<ID, InternalField<bool> *>::iterator it = internal_vectors_bool.begin();
        it != internal_vectors_bool.end();
        ++it) it->second->resize();
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 void Material::onElementsAdded(__attribute__((unused)) const Array<Element> & element_list,
 			       __attribute__((unused)) const NewElementsEvent & event) {
   this->resizeInternals();
 }
 
 /* -------------------------------------------------------------------------- */
 void Material::onElementsRemoved(const Array<Element> & element_list,
 				 const ElementTypeMapArray<UInt> & new_numbering,
 				 __attribute__((unused)) const RemovedElementsEvent & event) {
   UInt my_num = model->getInternalIndexFromID(getID());
 
   ElementTypeMapArray<UInt> material_local_new_numbering("remove mat filter elem", getID());
 
   Array<Element>::const_iterator<Element> el_begin = element_list.begin();
   Array<Element>::const_iterator<Element> el_end   = element_list.end();
 
   for (ghost_type_t::iterator g = ghost_type_t::begin(); g != ghost_type_t::end(); ++g) {
     GhostType gt = *g;
 
     ElementTypeMapArray<UInt>::type_iterator it  = new_numbering.firstType(_all_dimensions, gt, _ek_not_defined);
     ElementTypeMapArray<UInt>::type_iterator end = new_numbering.lastType (_all_dimensions, gt, _ek_not_defined);
     for (; it != end; ++it) {
       ElementType type = *it;
       if(element_filter.exists(type, gt)){
 	Array<UInt> & elem_filter = element_filter(type, gt);
 
 	Array<UInt> & mat_indexes = this->model->getMaterialByElement     (*it, gt);
 	Array<UInt> & mat_loc_num = this->model->getMaterialLocalNumbering(*it, gt);
 	UInt nb_element = this->model->getMesh().getNbElement(type, gt);
 
 	// all materials will resize of the same size...
 	mat_indexes.resize(nb_element);
 	mat_loc_num.resize(nb_element);
 
 	if(!material_local_new_numbering.exists(type, gt))
 	  material_local_new_numbering.alloc(elem_filter.getSize(), 1, type, gt);
 
 	Array<UInt> & mat_renumbering = material_local_new_numbering(type, gt);
 	const Array<UInt> & renumbering = new_numbering(type, gt);
 	Array<UInt> elem_filter_tmp;
 	UInt ni = 0;
 	Element el;
 	el.type = type;
 	el.ghost_type = gt;
 	for (UInt i = 0; i < elem_filter.getSize(); ++i) {
 	  el.element = elem_filter(i);
 	  if(std::find(el_begin, el_end, el) == el_end) {
 	    UInt new_el = renumbering(el.element);
 	    AKANTU_DEBUG_ASSERT(new_el != UInt(-1), "A not removed element as been badly renumbered");
 	    elem_filter_tmp.push_back(new_el);
 	    mat_renumbering(i) = ni;
 
 	    mat_indexes(new_el) = my_num;
 	    mat_loc_num(new_el) = ni;
 	    ++ni;
 	  } else {
 	    mat_renumbering(i) = UInt(-1);
 	  }
 	}
 
 	elem_filter.resize(elem_filter_tmp.getSize());
 	elem_filter.copy(elem_filter);
       }
     }
   }
 
   for (std::map<ID, InternalField<Real> *>::iterator it = internal_vectors_real.begin();
        it != internal_vectors_real.end();
        ++it) it->second->removeQuadraturePoints(material_local_new_numbering);
 
   for (std::map<ID, InternalField<UInt> *>::iterator it = internal_vectors_uint.begin();
        it != internal_vectors_uint.end();
        ++it) it->second->removeQuadraturePoints(material_local_new_numbering);
 
   for (std::map<ID, InternalField<bool> *>::iterator it = internal_vectors_bool.begin();
        it != internal_vectors_bool.end();
        ++it) it->second->removeQuadraturePoints(material_local_new_numbering);
 }
 
 /* -------------------------------------------------------------------------- */
 void Material::onBeginningSolveStep(const AnalysisMethod & method) {
   this->savePreviousState();
 }
 
 /* -------------------------------------------------------------------------- */
 void Material::onEndSolveStep(const AnalysisMethod & method) {
   ElementTypeMapArray<UInt>::type_iterator it
     = this->element_filter.firstType(_all_dimensions, _not_ghost, _ek_not_defined);
   ElementTypeMapArray<UInt>::type_iterator end
     = element_filter.lastType(_all_dimensions, _not_ghost, _ek_not_defined);
 
   for(; it != end; ++it) {
     this->updateEnergies(*it, _not_ghost);
   }
 }
 /* -------------------------------------------------------------------------- */
 void Material::onDamageIteration() {
   this->savePreviousState();
 }
 
 /* -------------------------------------------------------------------------- */
 void Material::onDamageUpdate() {
   ElementTypeMapArray<UInt>::type_iterator it
     = this->element_filter.firstType(_all_dimensions, _not_ghost, _ek_not_defined);
   ElementTypeMapArray<UInt>::type_iterator end
     = element_filter.lastType(_all_dimensions, _not_ghost, _ek_not_defined);
 
   for(; it != end; ++it) {
 
     if(!this->potential_energy.exists(*it, _not_ghost)) {
       UInt nb_element = this->element_filter(*it, _not_ghost).getSize();
       UInt nb_quadrature_points = this->model->getFEEngine().getNbQuadraturePoints(*it, _not_ghost);
 
       this->potential_energy.alloc(nb_element * nb_quadrature_points, 1,
 				   *it, _not_ghost);
     }
     this->updateEnergiesAfterDamage(*it, _not_ghost);
   }
 }
 
 /* -------------------------------------------------------------------------- */
 void Material::onDump(){
   if(this->isFiniteDeformation())
     this->computeAllCauchyStresses(_not_ghost);
 }
 
 /* -------------------------------------------------------------------------- */
 void Material::printself(std::ostream & stream, int indent) const {
   std::string space;
   for(Int i = 0; i < indent; i++, space += AKANTU_INDENT);
 
   std::string type = getID().substr(getID().find_last_of(":") + 1);
 
   stream << space << "Material " << type << " [" << std::endl;
   Parsable::printself(stream, indent);
   stream << space << "]" << std::endl;
 }
 
 /* -------------------------------------------------------------------------- */
 inline ElementTypeMap<UInt> Material::getInternalDataPerElem(const ID & id,
-							     const ElementKind & element_kind,
-                                                             const ID & fe_engine_id) const {
+                                                      const ElementKind & element_kind,
+                                                      const ID & fe_engine_id) const {
 
   std::map<ID, InternalField<Real> *>::const_iterator internal_array =
     internal_vectors_real.find(this->getID()+":"+id);
 
   if (internal_array == internal_vectors_real.end() ||
       internal_array->second->getElementKind() != element_kind)
     AKANTU_EXCEPTION("Cannot find internal field " << id << " in material " << this->name);
 
   InternalField<Real> & internal = *internal_array->second;
   InternalField<Real>::type_iterator it = internal.firstType(internal.getSpatialDimension(),
 							     _not_ghost, element_kind);
   InternalField<Real>::type_iterator last_type = internal.lastType(internal.getSpatialDimension(),
 								   _not_ghost, element_kind);
 
   ElementTypeMap<UInt> res;
   for(; it != last_type; ++it) {
     UInt nb_quadrature_points = 0;
     nb_quadrature_points = model->getFEEngine(fe_engine_id).getNbQuadraturePoints(*it);
 
     res(*it) = internal.getNbComponent() * nb_quadrature_points;
   }
   return res;
 }
 
 /* -------------------------------------------------------------------------- */
 void Material::flattenInternal(const std::string & field_id,
 			       ElementTypeMapArray<Real> & internal_flat,
 			       const GhostType ghost_type,
 			       ElementKind element_kind) const {
   this->flattenInternalIntern(field_id, internal_flat,
 			      this->spatial_dimension,
 			      ghost_type, element_kind);
 }
 
 /* -------------------------------------------------------------------------- */
 void Material::flattenInternalIntern(const std::string & field_id,
 				     ElementTypeMapArray<Real> & internal_flat,
 				     UInt spatial_dimension,
 				     const GhostType ghost_type,
 				     ElementKind element_kind,
 				     const ElementTypeMapArray<UInt> * element_filter,
 				     const Mesh * mesh) const {
   typedef ElementTypeMapArray<UInt>::type_iterator iterator;
 
   if(element_filter == NULL) element_filter = &(this->element_filter);
   if(mesh == NULL)  mesh = &(this->model->mesh);
 
   iterator tit = element_filter->firstType(spatial_dimension,
 					  ghost_type,
 					  element_kind);
   iterator end = element_filter->lastType(spatial_dimension,
 					 ghost_type,
 					 element_kind);
 
   for (; tit != end; ++tit) {
     ElementType type = *tit;
 
     try {
       __attribute__((unused)) const Array<Real> & src_vect
 	= this->getArray(field_id, type, ghost_type);
     } catch(debug::Exception & e) {
       continue;
     }
 
     const Array<Real> & src_vect = this->getArray(field_id, type, ghost_type);
     const Array<UInt> & filter   = (*element_filter)(type, ghost_type);
 
     // total number of elements for a given type
     UInt nb_element = mesh->getNbElement(type,ghost_type);
     // number of filtered elements
     UInt nb_element_src = filter.getSize();
     // number of quadrature points per elem
     UInt nb_quad_per_elem = 0;
     // number of data per quadrature point
     UInt nb_data_per_quad = src_vect.getNbComponent();
 
     if (!internal_flat.exists(type,ghost_type)) {
       internal_flat.alloc(nb_element * nb_quad_per_elem,
 			  nb_data_per_quad,
 			  type,
 			  ghost_type);
     }
 
     if (nb_element_src == 0) continue;
     nb_quad_per_elem = (src_vect.getSize() / nb_element_src);
 
     // number of data per element
     UInt nb_data = nb_quad_per_elem * src_vect.getNbComponent();
 
     Array<Real> & dst_vect = internal_flat(type,ghost_type);
     dst_vect.resize(nb_element*nb_quad_per_elem);
 
     Array<UInt>::const_scalar_iterator it  = filter.begin();
     Array<UInt>::const_scalar_iterator end = filter.end();
     Array<Real>::const_vector_iterator it_src =
       src_vect.begin_reinterpret(nb_data, nb_element_src);
 
     Array<Real>::vector_iterator it_dst =
       dst_vect.begin_reinterpret(nb_data, nb_element);
 
     for (; it != end ; ++it,++it_src) {
       it_dst[*it] = *it_src;
     }
   }
 };
 /* -------------------------------------------------------------------------- */
 
 
 __END_AKANTU__
diff --git a/src/model/solid_mechanics/materials/material_embedded/material_reinforcement.cc b/src/model/solid_mechanics/materials/material_embedded/material_reinforcement.cc
index 7c430e4b4..750936fe1 100644
--- a/src/model/solid_mechanics/materials/material_embedded/material_reinforcement.cc
+++ b/src/model/solid_mechanics/materials/material_embedded/material_reinforcement.cc
@@ -1,731 +1,738 @@
 /**
  * @file   material_reinforcement.cc
  *
  * @author Lucas Frérot <lucas.frerot@epfl.ch>
  *
  * @date creation: Thu Mar 12 2015
  * @date last modification: Thu Mar 12 2015
  *
  * @brief  Reinforcement 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)
  *
  * 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_common.hh"
 #include "aka_voigthelper.hh"
 #include "material_reinforcement.hh"
 
 __BEGIN_AKANTU__
 
 /* -------------------------------------------------------------------------- */
 template<UInt dim>
 MaterialReinforcement<dim>::MaterialReinforcement(SolidMechanicsModel & model,
                                                   UInt spatial_dimension,
                                                   const Mesh & mesh,
                                                   FEEngine & fe_engine,
                                                   const ID & id) :
   Material(model, dim, mesh, fe_engine, id), // /!\ dim, not spatial_dimension !
   model(NULL),
   stress_embedded("stress_embedded", *this, 1, fe_engine, this->element_filter),
   gradu_embedded("gradu_embedded", *this, 1, fe_engine, this->element_filter),
   directing_cosines("directing_cosines", *this, 1, fe_engine, this->element_filter),
   pre_stress("pre_stress", *this, 1, fe_engine, this->element_filter),
   area(1.0),
   shape_derivatives() {
   AKANTU_DEBUG_IN();
   this->initialize(model);
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template<UInt dim>
 void MaterialReinforcement<dim>::initialize(SolidMechanicsModel & a_model) {
   this->model = dynamic_cast<EmbeddedInterfaceModel *>(&a_model);
   AKANTU_DEBUG_ASSERT(this->model != NULL, "MaterialReinforcement needs an EmbeddedInterfaceModel");
 
   this->registerParam("area", area, _pat_parsable | _pat_modifiable,
                       "Reinforcement cross-sectional area");
   this->registerParam("pre_stress", pre_stress, _pat_parsable | _pat_modifiable,
                       "Uniform pre-stress");
 
 
   // Reallocate the element filter
   this->element_filter.free();
   this->model->getInterfaceMesh().initElementTypeMapArray(this->element_filter,
                                                           1, 1, false, _ek_regular);
 }
 
 /* -------------------------------------------------------------------------- */
 
 template<UInt dim>
 MaterialReinforcement<dim>::~MaterialReinforcement() {
   AKANTU_DEBUG_IN();
 
   ElementTypeMap<ElementTypeMapArray<Real> *>::type_iterator it = shape_derivatives.firstType();
   ElementTypeMap<ElementTypeMapArray<Real> *>::type_iterator end = shape_derivatives.lastType();
 
   for (; it != end ; ++it) {
     delete shape_derivatives(*it, _not_ghost);
     delete shape_derivatives(*it, _ghost);
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 
 template<UInt dim>
 void MaterialReinforcement<dim>::initMaterial() {
   Material::initMaterial();
 
   stress_embedded.initialize(dim * dim);
   gradu_embedded.initialize(dim * dim);
   pre_stress.initialize(1);
 
   /// We initialise the stuff that is not going to change during the simulation
   this->allocBackgroundShapeDerivatives();
   this->initBackgroundShapeDerivatives();
   this->initDirectingCosines();
 }
 
 /* -------------------------------------------------------------------------- */
 
+/**
+ * Background shape derivatives need to be stored per background element
+ * types but also per embedded element type, which is why they are stored
+ * in an ElementTypeMap<ElementTypeMapArray<Real> *>. The outer ElementTypeMap
+ * refers to the embedded types, and the inner refers to the background types.
+ */
 template<UInt dim>
 void MaterialReinforcement<dim>::allocBackgroundShapeDerivatives() {
   AKANTU_DEBUG_IN();
 
   Mesh & interface_mesh = model->getInterfaceMesh();
   Mesh & mesh = model->getMesh();
 
   ghost_type_t::iterator int_ghost_it = ghost_type_t::begin();
 
   // Loop over interface ghosts
   for (; int_ghost_it != ghost_type_t::end() ; ++int_ghost_it) {
     Mesh::type_iterator interface_type_it = interface_mesh.firstType(1, *int_ghost_it);
     Mesh::type_iterator interface_type_end = interface_mesh.lastType(1, *int_ghost_it);
 
     // Loop over interface types
     for (; interface_type_it != interface_type_end ; ++interface_type_it) {
       Mesh::type_iterator background_type_it = mesh.firstType(dim, *int_ghost_it);
       Mesh::type_iterator background_type_end = mesh.lastType(dim, *int_ghost_it);
 
       // Loop over background types
       for (; background_type_it != background_type_end ; ++background_type_it) {
 	const ElementType & int_type = *interface_type_it;
 	const ElementType & back_type = *background_type_it;
 	const GhostType & int_ghost = *int_ghost_it;
 
 	std::string shaped_id = "embedded_shape_derivatives";
 
 	if (int_ghost == _ghost) shaped_id += ":ghost";
 
 	ElementTypeMapArray<Real> * shaped_etma = new ElementTypeMapArray<Real>(shaped_id, this->name);
 
 	UInt nb_points = Mesh::getNbNodesPerElement(back_type);
 	UInt nb_quad_points = model->getFEEngine("EmbeddedInterfaceFEEngine").getNbQuadraturePoints(int_type);
 	UInt nb_elements = element_filter(int_type, int_ghost).getSize();
 
         // Alloc the background ElementTypeMapArray
         shaped_etma->alloc(nb_elements * nb_quad_points,
             dim * nb_points,
             back_type);
 
         // Insert the background ElementTypeMapArray in the interface ElementTypeMap
         shape_derivatives(shaped_etma, int_type, int_ghost);
       }
     }
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 
 template<UInt dim>
 void MaterialReinforcement<dim>::initBackgroundShapeDerivatives() {
   AKANTU_DEBUG_IN();
 
   Mesh & mesh = model->getMesh();
 
   Mesh::type_iterator type_it = mesh.firstType(dim, _not_ghost);
   Mesh::type_iterator type_end = mesh.lastType(dim, _not_ghost);
 
   for (; type_it != type_end ; ++type_it) {
     computeBackgroundShapeDerivatives(*type_it, _not_ghost);
     //computeBackgroundShapeDerivatives(*type_it, _ghost);
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 
 template<UInt dim>
 void MaterialReinforcement<dim>::initDirectingCosines() {
   AKANTU_DEBUG_IN();
 
   Mesh & mesh = model->getInterfaceMesh();
 
   Mesh::type_iterator type_it = mesh.firstType(1, _not_ghost);
   Mesh::type_iterator type_end = mesh.lastType(1, _not_ghost);
 
   const UInt voigt_size = getTangentStiffnessVoigtSize(dim);
   directing_cosines.initialize(voigt_size * voigt_size);
 
   for (; type_it != type_end ; ++type_it) {
     computeDirectingCosines(*type_it, _not_ghost);
     //computeDirectingCosines(*type_it, _ghost);
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 
 template<UInt dim>
 void MaterialReinforcement<dim>::assembleStiffnessMatrix(GhostType ghost_type) {
   AKANTU_DEBUG_IN();
 
   Mesh & interface_mesh = model->getInterfaceMesh();
 
   Mesh::type_iterator type_it = interface_mesh.firstType(1, _not_ghost);
   Mesh::type_iterator type_end = interface_mesh.lastType(1, _not_ghost);
 
   for (; type_it != type_end ; ++type_it) {
     assembleStiffnessMatrix(*type_it, ghost_type);
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 
 template<UInt dim>
 void MaterialReinforcement<dim>::updateResidual(GhostType ghost_type) {
   AKANTU_DEBUG_IN();
 
   computeAllStresses(ghost_type);
   assembleResidual(ghost_type);
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 
 template<UInt dim>
 void MaterialReinforcement<dim>::assembleResidual(GhostType ghost_type) {
   AKANTU_DEBUG_IN();
 
   Mesh & interface_mesh = model->getInterfaceMesh();
 
   Mesh::type_iterator type_it = interface_mesh.firstType(1, _not_ghost);
   Mesh::type_iterator type_end = interface_mesh.lastType(1, _not_ghost);
 
   for (; type_it != type_end ; ++type_it) {
     assembleResidual(*type_it, ghost_type);
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template<UInt dim>
 void MaterialReinforcement<dim>::computeGradU(const ElementType & type, GhostType ghost_type) {
   AKANTU_DEBUG_IN();
 
   Array<UInt> & elem_filter = element_filter(type, ghost_type);
   UInt nb_element = elem_filter.getSize();
   UInt nb_quad_points = model->getFEEngine("EmbeddedInterfaceFEEngine").getNbQuadraturePoints(type);
 
   Array<Real> & gradu_vec = gradu_embedded(type, ghost_type);
 
   Mesh::type_iterator back_it = model->getMesh().firstType(dim, ghost_type);
   Mesh::type_iterator back_end = model->getMesh().lastType(dim, ghost_type);
 
   for (; back_it != back_end ; ++back_it) {
     UInt nodes_per_background_e = Mesh::getNbNodesPerElement(*back_it);
 
     Array<Real> & shapesd = shape_derivatives(type, ghost_type)->operator()(*back_it, ghost_type);
 
     Array<UInt> * background_filter = new Array<UInt>(nb_element, 1, "background_filter");
     filterInterfaceBackgroundElements(*background_filter, *back_it, type, ghost_type);
 
     Array<Real> * disp_per_element = new Array<Real>(0, dim * nodes_per_background_e, "disp_elem");
     FEEngine::extractNodalToElementField(model->getMesh(),
 	model->getDisplacement(),
 	*disp_per_element,
 	*back_it, ghost_type, *background_filter);
 
     Array<Real>::matrix_iterator disp_it = disp_per_element->begin(dim, nodes_per_background_e);
     Array<Real>::matrix_iterator disp_end = disp_per_element->end(dim, nodes_per_background_e);
 
     Array<Real>::matrix_iterator shapes_it = shapesd.begin(dim, nodes_per_background_e);
     Array<Real>::matrix_iterator grad_u_it = gradu_vec.begin(dim, dim);
 
     for (; disp_it != disp_end ; ++disp_it) {
       for (UInt i = 0; i < nb_quad_points; i++, ++shapes_it, ++grad_u_it) {
 	Matrix<Real> & B = *shapes_it;
 	Matrix<Real> & du = *grad_u_it;
 	Matrix<Real> & u = *disp_it;
 
 	du.mul<false, true>(u, B);
       }
     }
 
     delete background_filter;
     delete disp_per_element;
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template<UInt dim>
 void MaterialReinforcement<dim>::computeAllStresses(GhostType ghost_type) {
   AKANTU_DEBUG_IN();
 
   Mesh::type_iterator it = model->getInterfaceMesh().firstType();
   Mesh::type_iterator last_type = model->getInterfaceMesh().lastType();
 
   for(; it != last_type; ++it) {
     computeGradU(*it, ghost_type);
     computeStress(*it, ghost_type);
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template<UInt dim>
 void MaterialReinforcement<dim>::assembleResidual(const ElementType & type, GhostType ghost_type) {
   AKANTU_DEBUG_IN();
 
   Mesh & mesh = model->getMesh();
 
   Mesh::type_iterator type_it = mesh.firstType(dim, ghost_type);
   Mesh::type_iterator type_end = mesh.lastType(dim, ghost_type);
 
   for (; type_it != type_end ; ++type_it) {
-    assembleResidual(type, *type_it, ghost_type);
+    assembleResidualInterface(type, *type_it, ghost_type);
   }
 
 
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
+
+/**
+ * Computes and assemble the residual. Residual in reinforcement is computed as :
+ * \f[
+ * \vec{r} = A_s \int_S{\mathbf{B}^T\mathbf{C}^T \vec{\sigma_s}\,\mathrm{d}s}
+ * \f]
+ */
 template<UInt dim>
-void MaterialReinforcement<dim>::assembleResidual(const ElementType & interface_type,
-                                                  const ElementType & background_type,
-                                                  GhostType ghost_type) {
+void MaterialReinforcement<dim>::assembleResidualInterface(const ElementType & interface_type,
+                                                           const ElementType & background_type,
+                                                           GhostType ghost_type) {
   AKANTU_DEBUG_IN();
 
   UInt voigt_size = getTangentStiffnessVoigtSize(dim);
 
   Array<Real> & residual = const_cast<Array<Real> &>(model->getResidual());
 
   FEEngine & interface_engine = model->getFEEngine("EmbeddedInterfaceFEEngine");
   FEEngine & background_engine = model->getFEEngine();
 
   Array<UInt> & elem_filter = element_filter(interface_type, ghost_type);
 
   UInt nodes_per_background_e = Mesh::getNbNodesPerElement(background_type);
   UInt nb_quadrature_points = interface_engine.getNbQuadraturePoints(interface_type, ghost_type);
   UInt nb_element = elem_filter.getSize();
 
   UInt back_dof = dim * nodes_per_background_e;
 
   Array<Real> & shapesd = shape_derivatives(interface_type, ghost_type)->operator()(background_type, ghost_type);
 
   Array<Real> * integrant = new Array<Real>(nb_quadrature_points * nb_element,
 					    back_dof,
 					    "integrant");
 
   Array<Real>::vector_iterator integrant_it =
     integrant->begin(back_dof);
   Array<Real>::vector_iterator integrant_end =
     integrant->end(back_dof);
 
   Array<Real>::matrix_iterator B_it =
     shapesd.begin(dim, nodes_per_background_e);
   Array<Real>::matrix_iterator C_it =
     directing_cosines(interface_type, ghost_type).begin(voigt_size, voigt_size);
   Array<Real>::matrix_iterator sigma_it =
     stress_embedded(interface_type, ghost_type).begin(dim, dim);
 
   Vector<Real> sigma(voigt_size);
   Matrix<Real> Bvoigt(voigt_size, back_dof);
   Vector<Real> Ct_sigma(voigt_size);
 
   for (; integrant_it != integrant_end ; ++integrant_it,
 					 ++B_it,
 					 ++C_it,
 					 ++sigma_it) {
     VoigtHelper<dim>::transferBMatrixToSymVoigtBMatrix(*B_it, Bvoigt, nodes_per_background_e);
     Matrix<Real> & C = *C_it;
     Vector<Real> & BtCt_sigma = *integrant_it;
 
     stressTensorToVoigtVector(*sigma_it, sigma);
 
     Ct_sigma.mul<true>(C, sigma);
     BtCt_sigma.mul<true>(Bvoigt, Ct_sigma);
     BtCt_sigma *= area;
   }
 
 
   Array<Real> * residual_interface = new Array<Real>(nb_element, back_dof, "residual_interface");
   interface_engine.integrate(*integrant,
                              *residual_interface,
                              back_dof,
                              interface_type, ghost_type,
                              elem_filter);
 
   delete integrant;
 
   Array<UInt> * background_filter = new Array<UInt>(nb_element, 1, "background_filter");
 
   filterInterfaceBackgroundElements(*background_filter,
                                     background_type, interface_type,
                                     ghost_type);
   background_engine.assembleArray(*residual_interface, residual,
                                   model->getDOFSynchronizer().getLocalDOFEquationNumbers(),
                                   dim, background_type, ghost_type, *background_filter, -1.0);
   delete residual_interface;
   delete background_filter;
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template<UInt dim>
 void MaterialReinforcement<dim>::filterInterfaceBackgroundElements(Array<UInt> & filter,
 
                                                                    const ElementType & type,
                                                                    const ElementType & interface_type,
                                                                    GhostType ghost_type) {
   AKANTU_DEBUG_IN();
 
   filter.resize(0);
   filter.clear();
 
   Array<Element> & elements = model->getInterfaceAssociatedElements(interface_type, ghost_type);
   Array<UInt> & elem_filter = element_filter(interface_type, ghost_type);
 
   Array<UInt>::scalar_iterator
     filter_it = elem_filter.begin(),
     filter_end = elem_filter.end();
 
   for (; filter_it != filter_end ; ++filter_it) {
     Element & elem = elements(*filter_it);
     if (elem.type == type) filter.push_back(elem.element);
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 
 template<UInt dim>
 void MaterialReinforcement<dim>::computeDirectingCosines(const ElementType & type, GhostType ghost_type) {
   AKANTU_DEBUG_IN();
 
   Mesh & interface_mesh = this->model->getInterfaceMesh();
 
   const UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
   const UInt steel_dof = dim * nb_nodes_per_element;
   const UInt voigt_size = getTangentStiffnessVoigtSize(dim);
   const UInt nb_quad_points =
     model->getFEEngine("EmbeddedInterfaceFEEngine").getNbQuadraturePoints(type, ghost_type);
 
   Array<Real> node_coordinates(this->element_filter(type, ghost_type).getSize(), steel_dof);
 
   this->model->getFEEngine().template extractNodalToElementField<Real>(interface_mesh,
 								       interface_mesh.getNodes(),
 								       node_coordinates,
 								       type,
 								       ghost_type,
 								       this->element_filter(type, ghost_type));
 
   Array<Real>::matrix_iterator
     directing_cosines_it = directing_cosines(type, ghost_type).begin(voigt_size, voigt_size);
 
   Array<Real>::matrix_iterator node_coordinates_it = node_coordinates.begin(dim, nb_nodes_per_element);
   Array<Real>::matrix_iterator node_coordinates_end = node_coordinates.end(dim, nb_nodes_per_element);
 
   for (; node_coordinates_it != node_coordinates_end ; ++node_coordinates_it) {
     for (UInt i = 0 ; i < nb_quad_points ; i++, ++directing_cosines_it) {
       Matrix<Real> & nodes = *node_coordinates_it;
       Matrix<Real> & cosines = *directing_cosines_it;
 
       computeDirectingCosinesOnQuad(nodes, cosines);
     }
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 
 template<UInt dim>
 void MaterialReinforcement<dim>::assembleStiffnessMatrix(const ElementType & type, GhostType ghost_type) {
   AKANTU_DEBUG_IN();
 
   Mesh & mesh = model->getMesh();
 
   Mesh::type_iterator type_it = mesh.firstType(dim, ghost_type);
   Mesh::type_iterator type_end = mesh.lastType(dim, ghost_type);
 
   for (; type_it != type_end ; ++type_it) {
-    assembleStiffnessMatrix(type, *type_it, ghost_type);
+    assembleStiffnessMatrixInterface(type, *type_it, ghost_type);
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
+/**
+ * Computes the reinforcement stiffness matrix (Gomes & Awruch, 2001)
+ * \f[
+ * \mathbf{K}_e = \sum_{i=1}^R{A_i\int_{S_i}{\mathbf{B}^T
+ * \mathbf{C}_i^T \mathbf{D}_{s, i} \mathbf{C}_i \mathbf{B}\,\mathrm{d}s}}
+ * \f]
+ */
 template<UInt dim>
-void MaterialReinforcement<dim>::assembleStiffnessMatrix(const ElementType & interface_type,
-                                                         const ElementType & background_type,
-                                                         GhostType ghost_type) {
+void MaterialReinforcement<dim>::assembleStiffnessMatrixInterface(const ElementType & interface_type,
+                                                                  const ElementType & background_type,
+                                                                  GhostType ghost_type) {
   AKANTU_DEBUG_IN();
 
   UInt voigt_size = getTangentStiffnessVoigtSize(dim);
 
   SparseMatrix & K = const_cast<SparseMatrix &>(model->getStiffnessMatrix());
 
   FEEngine & background_engine = model->getFEEngine();
   FEEngine & interface_engine = model->getFEEngine("EmbeddedInterfaceFEEngine");
 
   Array<UInt> & elem_filter = element_filter(interface_type, ghost_type);
   Array<Real> & grad_u = gradu_embedded(interface_type, ghost_type);
 
   UInt nb_element = elem_filter.getSize();
   UInt nodes_per_background_e = Mesh::getNbNodesPerElement(background_type);
   UInt nb_quadrature_points = interface_engine.getNbQuadraturePoints(interface_type, ghost_type);
 
   UInt back_dof = dim * nodes_per_background_e;
 
   UInt integrant_size = back_dof;
 
   grad_u.resize(nb_quadrature_points * nb_element);
 
   Array<Real> * tangent_moduli = new Array<Real>(nb_element * nb_quadrature_points,
 						 1, "interface_tangent_moduli");
   tangent_moduli->clear();
   computeTangentModuli(interface_type, *tangent_moduli, ghost_type);
 
   Array<Real> & shapesd = shape_derivatives(interface_type, ghost_type)->operator()(background_type, ghost_type);
 
   Array<Real> * integrant = new Array<Real>(nb_element * nb_quadrature_points,
 					    integrant_size * integrant_size,
 					    "B^t*C^t*D*C*B");
   integrant->clear();
 
   /// Temporary matrices for integrant product
   Matrix<Real> Bvoigt(voigt_size, back_dof);
   Matrix<Real> DC(voigt_size, voigt_size);
   Matrix<Real> DCB(voigt_size, back_dof);
   Matrix<Real> CtDCB(voigt_size, back_dof);
 
   Array<Real>::scalar_iterator D_it = tangent_moduli->begin();
   Array<Real>::scalar_iterator D_end = tangent_moduli->end();
 
   Array<Real>::matrix_iterator C_it =
     directing_cosines(interface_type, ghost_type).begin(voigt_size, voigt_size);
   Array<Real>::matrix_iterator B_it =
     shapesd.begin(dim, nodes_per_background_e);
   Array<Real>::matrix_iterator integrant_it =
     integrant->begin(integrant_size, integrant_size);
 
   for (; D_it != D_end ; ++D_it, ++C_it, ++B_it, ++integrant_it) {
     Real & D = *D_it;
     Matrix<Real> & C = *C_it;
     Matrix<Real> & B = *B_it;
     Matrix<Real> & BtCtDCB = *integrant_it;
 
     VoigtHelper<dim>::transferBMatrixToSymVoigtBMatrix(B, Bvoigt, nodes_per_background_e);
 
     DC.clear();
     DC(0, 0) = D * area;
     DC *= C;
     DCB.mul<false, false>(DC, Bvoigt);
     CtDCB.mul<true, false>(C, DCB);
     BtCtDCB.mul<true, false>(Bvoigt, CtDCB);
   }
 
   delete tangent_moduli;
 
   Array<Real> * K_interface = new Array<Real>(nb_element, integrant_size * integrant_size, "K_interface");
   interface_engine.integrate(*integrant, *K_interface,
                              integrant_size * integrant_size,
                              interface_type, ghost_type,
                              elem_filter);
 
   delete integrant;
 
   Array<UInt> * background_filter = new Array<UInt>(nb_element, 1, "background_filter");
 
   filterInterfaceBackgroundElements(*background_filter,
                                     background_type, interface_type,
                                     ghost_type);
 
   background_engine.assembleMatrix(*K_interface, K, dim, background_type, ghost_type, *background_filter);
 
   delete K_interface;
   delete background_filter;
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 
 /// In this function, type and ghost_type refer to background elements
 template<UInt dim>
 void MaterialReinforcement<dim>::computeBackgroundShapeDerivatives(const ElementType & type, GhostType ghost_type) {
   AKANTU_DEBUG_IN();
 
   Mesh & interface_mesh = model->getInterfaceMesh();
 
   FEEngine & engine = model->getFEEngine();
   FEEngine & interface_engine = model->getFEEngine("EmbeddedInterfaceFEEngine");
 
   Mesh::type_iterator interface_type = interface_mesh.firstType();
   Mesh::type_iterator interface_last = interface_mesh.lastType();
 
   for (; interface_type != interface_last ; ++interface_type) {
     Array<UInt> & filter = element_filter(*interface_type, ghost_type);
 
     const UInt nb_elements = filter.getSize();
     const UInt nb_nodes = Mesh::getNbNodesPerElement(type);
     const UInt nb_quad_per_element = interface_engine.getNbQuadraturePoints(*interface_type);
 
     Array<Real> quad_pos(nb_quad_per_element * nb_elements, dim, "interface_quad_points");
     quad_pos.resize(nb_quad_per_element * nb_elements);
     interface_engine.interpolateOnQuadraturePoints(interface_mesh.getNodes(),
 						   quad_pos, dim, *interface_type,
 						   ghost_type, filter);
 
 
     Array<Real> & background_shapesd = shape_derivatives(*interface_type, ghost_type)->operator()(type, ghost_type);
     background_shapesd.clear();
 
     Array<UInt> * background_elements = new Array<UInt>(nb_elements, 1, "computeBackgroundShapeDerivatives:background_filter");
     filterInterfaceBackgroundElements(*background_elements, type, *interface_type, ghost_type);
 
     Array<UInt>::scalar_iterator
       back_it = background_elements->begin(),
       back_end = background_elements->end();
 
     Array<Real>::matrix_iterator shapesd_it = background_shapesd.begin(dim, nb_nodes);
 
     Array<Real>::vector_iterator quad_pos_it = quad_pos.begin(dim);
 
 
     for (; back_it != back_end ; ++back_it) {
       for (UInt i = 0 ; i < nb_quad_per_element ; i++, ++shapesd_it, ++quad_pos_it)
-	engine.computeShapeDerivatives(*quad_pos_it, *back_it, type, *shapesd_it, ghost_type);
+        engine.computeShapeDerivatives(*quad_pos_it, *back_it, type, *shapesd_it, ghost_type);
     }
 
     delete background_elements;
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 template<UInt dim>
 Real MaterialReinforcement<dim>::getEnergy(std::string id) {
   AKANTU_DEBUG_IN();
   if (id == "potential") {
     Real epot = 0.;
 
     computePotentialEnergyByElements();
 
     Mesh::type_iterator
       it = element_filter.firstType(spatial_dimension),
       end = element_filter.lastType(spatial_dimension);
 
     for (; it != end ; ++it) {
       FEEngine & interface_engine = model->getFEEngine("EmbeddedInterfaceFEEngine");
       epot += interface_engine.integrate(potential_energy(*it, _not_ghost),
                                          *it, _not_ghost,
                                          element_filter(*it, _not_ghost));
       epot *= area;
     }
 
     return epot;
   }
 
   AKANTU_DEBUG_OUT();
   return 0;
 }
 
 /* -------------------------------------------------------------------------- */
 template<UInt dim>
 ElementTypeMap<UInt> MaterialReinforcement<dim>::getInternalDataPerElem(const ID & field_name,
-									const ElementKind & kind,
-									const ID & fe_engine_id) {
-  if (field_name == "stress_embedded")
-    return Material::getInternalDataPerElem(field_name, kind, "EmbeddedInterfaceFEEngine");
-  else {
-    return Material::getInternalDataPerElem(field_name, kind, fe_engine_id);
-  }
+                                                                        const ElementKind & kind,
+                                                                        const ID & fe_engine_id) const {
+  return Material::getInternalDataPerElem(field_name, kind, "EmbeddedInterfaceFEEngine");
 }
 
 
 /* -------------------------------------------------------------------------- */
 // Author is Guillaume Anciaux, see material.cc
 template<UInt dim>
 void MaterialReinforcement<dim>::flattenInternal(const std::string & field_id,
-						 ElementTypeMapArray<Real> & internal_flat,
-						 const GhostType ghost_type,
-						 ElementKind element_kind) {
+                                                 ElementTypeMapArray<Real> & internal_flat,
+                                                 const GhostType ghost_type,
+                                                 ElementKind element_kind) const {
   AKANTU_DEBUG_IN();
 
-  // si c'est pour tous les champ qu'il faut spatial_dimension = 1 vire juste le if statement
-  if(field_id == "stress_embedded") {
-    Material::flattenInternalIntern(field_id,
-				    internal_flat,
-				    1,
-				    ghost_type,
-				    _ek_not_defined,
-				    &(this->element_filter),
-				    &(this->model->getInterfaceMesh()));
-
-  } else {
-    Material::flattenInternal(field_id,
-			      internal_flat,
-			      ghost_type,
-			      _ek_not_defined);
-  }
+  Material::flattenInternalIntern(field_id,
+                                  internal_flat,
+                                  1,
+                                  ghost_type,
+                                  _ek_not_defined,
+                                  &(this->element_filter),
+                                  &(this->model->getInterfaceMesh()));
 
   AKANTU_DEBUG_OUT();
 }
 
 
 /* -------------------------------------------------------------------------- */
 
 INSTANTIATE_MATERIAL(MaterialReinforcement);
 
 /* -------------------------------------------------------------------------- */
 
 __END_AKANTU__
diff --git a/src/model/solid_mechanics/materials/material_embedded/material_reinforcement.hh b/src/model/solid_mechanics/materials/material_embedded/material_reinforcement.hh
index bc6f43aa6..e6b5ffeb6 100644
--- a/src/model/solid_mechanics/materials/material_embedded/material_reinforcement.hh
+++ b/src/model/solid_mechanics/materials/material_embedded/material_reinforcement.hh
@@ -1,232 +1,197 @@
 /**
  * @file   material_reinforcement.hh
  *
  * @author Lucas Frérot <lucas.frerot@epfl.ch>
  *
  * @date creation: Thu Mar 12 2015
  * @date last modification: Thu Mar 12 2015
  *
  * @brief  Reinforcement 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)
  *
  * 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/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 
 #ifndef __AKANTU_MATERIAL_REINFORCEMENT_HH__
 #define __AKANTU_MATERIAL_REINFORCEMENT_HH__
 
 #include "aka_common.hh"
 
 #include "material.hh"
 #include "embedded_interface_model.hh"
 #include "embedded_internal_field.hh"
 
 /* -------------------------------------------------------------------------- */
 
 __BEGIN_AKANTU__
 
 /**
  * @brief Material used to represent embedded reinforcements
  *
  * This class is used for computing the reinforcement stiffness matrix
  * along with the reinforcement residual. Room is made for constitutive law,
  * but actual use of contitutive laws is made in MaterialReinforcementTemplate.
  *
  * Be careful with the dimensions in this class :
  *  -  this->spatial_dimension is always 1
  *  -  the template parameter dim is the dimension of the problem
  */
 template<UInt dim>
 class MaterialReinforcement : virtual public Material {
 
   /* ------------------------------------------------------------------------ */
   /* Constructors/Destructors                                                 */
   /* ------------------------------------------------------------------------ */
 public:
   /// Constructor
   MaterialReinforcement(SolidMechanicsModel & model,
                         UInt spatial_dimension,
                         const Mesh & mesh,
                         FEEngine & fe_engine,
                         const ID & id = "");
 
   /// Destructor
   virtual ~MaterialReinforcement();
 
 protected:
   void initialize(SolidMechanicsModel & a_model);
 
   /* ------------------------------------------------------------------------ */
   /* Methods                                                                  */
   /* ------------------------------------------------------------------------ */
 public:
   /// Init the material
   virtual void initMaterial();
 
   /// Init the background shape derivatives
   void initBackgroundShapeDerivatives();
 
   /// Init the cosine matrices
   void initDirectingCosines();
 
   /// Assemble stiffness matrix
   virtual void assembleStiffnessMatrix(GhostType ghost_type);
 
   /// Update the residual
   virtual void updateResidual(GhostType ghost_type = _not_ghost);
 
   /// Assembled the residual
   virtual void assembleResidual(GhostType ghost_type);
 
   /// Compute all the stresses !
   virtual void computeAllStresses(GhostType ghost_type);
 
   /// Compute the stiffness parameter for elements of a type
   virtual void computeTangentModuli(const ElementType & type,
                                     Array<Real> & tangent,
                                     GhostType ghost_type) = 0;
 
   /// Compute energy
   virtual Real getEnergy(std::string id);
 
-  ElementTypeMap<UInt> getInternalDataPerElem(const ID & field_name,
-					      const ElementKind & kind,
-					      const ID & fe_engine_id);
+  virtual ElementTypeMap<UInt> getInternalDataPerElem(const ID & field_name,
+                                                      const ElementKind & kind,
+                                                      const ID & fe_engine_id) const;
 
   /// Reimplementation of Material's function to accomodate for interface mesh
   virtual void flattenInternal(const std::string & field_id,
-			       ElementTypeMapArray<Real> & internal_flat,
+                               ElementTypeMapArray<Real> & internal_flat,
                                const GhostType ghost_type = _not_ghost,
-                               ElementKind element_kind = _ek_not_defined);
+                               ElementKind element_kind = _ek_not_defined) const;
 
   /* ------------------------------------------------------------------------ */
   /* Protected methods                                                        */
   /* ------------------------------------------------------------------------ */
 protected:
-  /**
-   * @brief Allocate the background shape derivatives
-   *
-   * Background shape derivatives need to be stored per background element
-   * types but also per embedded element type, which is why they are stored
-   * in an ElementTypeMap<ElementTypeMapArray<Real> *>. The outer ElementTypeMap
-   * refers to the embedded types, and the inner refers to the background types.
-   */
+  /// Allocate the background shape derivatives
   void allocBackgroundShapeDerivatives();
 
   /// Compute the directing cosines matrix for one element type
   void computeDirectingCosines(const ElementType & type, GhostType ghost_type);
 
-  /**
-   * @brief Compute the directing cosines matrix on quadrature points.
-   *
-   * The structure of the directing cosines matrix is :
-   * \f{eqnarray*}{
-   *  C_{1,\cdot} & = & (l^2, m^2, n^2, lm, mn, ln) \\
-   *  C_{i,j} & = & 0
-   * \f}
-   *
-   * with :
-   * \f[
-   * (l, m, n) = \frac{1}{\|\frac{\mathrm{d}\vec{r}(s)}{\mathrm{d}s}\|} \cdot \frac{\mathrm{d}\vec{r}(s)}{\mathrm{d}s}
-   * \f]
-   */
+  /// Compute the directing cosines matrix on quadrature points.
   inline void computeDirectingCosinesOnQuad(const Matrix<Real> & nodes,
                                             Matrix<Real> & cosines);
 
   /// Assemble the stiffness matrix for an element type (typically _segment_2)
   void assembleStiffnessMatrix(const ElementType & type, GhostType ghost_type);
 
-  /**
-   * @brief Assemble the stiffness matrix for background & interface types
-   *
-   * Computes the reinforcement stiffness matrix (Gomes & Awruch, 2001)
-   * \f[
-   * \mathbf{K}_e = \sum_{i=1}^R{A_i\int_{S_i}{\mathbf{B}^T
-   * \mathbf{C}_i^T \mathbf{D}_{s, i} \mathbf{C}_i \mathbf{B}\,\mathrm{d}s}}
-   * \f]
-   */
-  void assembleStiffnessMatrix(const ElementType & interface_type,
-                               const ElementType & background_type,
-                               GhostType ghost_type);
+  /// Assemble the stiffness matrix for background & interface types
+  void assembleStiffnessMatrixInterface(const ElementType & interface_type,
+                                        const ElementType & background_type,
+                                        GhostType ghost_type);
 
   /// Compute the background shape derivatives for a type
   void computeBackgroundShapeDerivatives(const ElementType & type, GhostType ghost_type);
 
   /// Filter elements crossed by interface of a type
   void filterInterfaceBackgroundElements(Array<UInt> & filter,
                                          const ElementType & type,
                                          const ElementType & interface_type,
                                          GhostType ghost_type);
 
   /// Assemble the residual of one type of element (typically _segment_2)
   void assembleResidual(const ElementType & type, GhostType ghost_type);
 
-  /**
-   * @brief Assemble the residual for a pair of elements
-   *
-   * Computes and assemble the residual. Residual in reinforcement is computed as :
-   * \f[
-   * \vec{r} = A_s \int_S{\mathbf{B}^T\mathbf{C}^T \vec{\sigma_s}\,\mathrm{d}s}
-   * \f]
-   */
-  void assembleResidual(const ElementType & interface_type,
-                        const ElementType & background_type,
-                        GhostType ghost_type);
+  /// Assemble the residual for a pair of elements
+  void assembleResidualInterface(const ElementType & interface_type,
+                                 const ElementType & background_type,
+                                 GhostType ghost_type);
 
   // TODO figure out why voigt size is 4 in 2D
   inline void stressTensorToVoigtVector(const Matrix<Real> & tensor, Vector<Real> & vector);
   inline void strainTensorToVoigtVector(const Matrix<Real> & tensor, Vector<Real> & vector);
 
   /// Compute gradu on the interface quadrature points
   virtual void computeGradU(const ElementType & type, GhostType ghost_type);
 
   /* ------------------------------------------------------------------------ */
   /* Class Members                                                            */
   /* ------------------------------------------------------------------------ */
 protected:
   /// Embedded model
   EmbeddedInterfaceModel * model;
 
   /// Stress in the reinforcement
   InternalField<Real> stress_embedded;
 
   /// Gradu of concrete on reinforcement
   InternalField<Real> gradu_embedded;
 
   /// C matrix on quad
   InternalField<Real> directing_cosines;
 
   /// Prestress on quad
   InternalField<Real> pre_stress;
 
   /// Cross-sectional area
   Real area;
 
   /// Background mesh shape derivatives
   ElementTypeMap< ElementTypeMapArray<Real> * > shape_derivatives;
 
 };
 
 #include "material_reinforcement_inline_impl.cc"
 
 __END_AKANTU__
 
 #endif // __AKANTU_MATERIAL_REINFORCEMENT_HH__
diff --git a/src/model/solid_mechanics/materials/material_embedded/material_reinforcement_inline_impl.cc b/src/model/solid_mechanics/materials/material_embedded/material_reinforcement_inline_impl.cc
index 8bde1600e..d6b2e756e 100644
--- a/src/model/solid_mechanics/materials/material_embedded/material_reinforcement_inline_impl.cc
+++ b/src/model/solid_mechanics/materials/material_embedded/material_reinforcement_inline_impl.cc
@@ -1,99 +1,111 @@
 /**
  * @file   material_reinforcement_inline_impl.cc
  *
  * @author Lucas Frérot <lucas.frerot@epfl.ch>
  *
  * @date creation: Thu Mar 12 2015
  * @date last modification: Thu Mar 12 2015
  *
  * @brief  Reinforcement 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)
  *
  * 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/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 
+/**
+ * The structure of the directing cosines matrix is :
+ * \f{eqnarray*}{
+ *  C_{1,\cdot} & = & (l^2, m^2, n^2, lm, mn, ln) \\
+ *  C_{i,j} & = & 0
+ * \f}
+ *
+ * with :
+ * \f[
+ * (l, m, n) = \frac{1}{\|\frac{\mathrm{d}\vec{r}(s)}{\mathrm{d}s}\|} \cdot \frac{\mathrm{d}\vec{r}(s)}{\mathrm{d}s}
+ * \f]
+ */
 template<UInt dim>
 inline void MaterialReinforcement<dim>::computeDirectingCosinesOnQuad(const Matrix<Real> & nodes,
                                                                  Matrix<Real> & cosines) {
   AKANTU_DEBUG_IN();
 
   AKANTU_DEBUG_ASSERT(nodes.cols() == 2, "Higher order reinforcement elements not implemented");
 
   const Vector<Real> & a = nodes(0), b = nodes(1);
 
   cosines.clear();
 
   Real sq_length = 0.;
   for (UInt i = 0 ; i < dim ; i++) {
     sq_length += Math::pow<2, Real>(b(i) - a(i));
   }
 
   // Fill the first row of cosine matrix
   for (UInt i = 0 ; i < dim ; i++) {
     cosines(0, i) = Math::pow<2, Real>(b(i) - a(i)) / sq_length;
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 
 template<UInt dim>
 inline void MaterialReinforcement<dim>::stressTensorToVoigtVector(const Matrix<Real> & tensor,
                                                                   Vector<Real> & vector) {
   AKANTU_DEBUG_IN();
   
   for (UInt i = 0; i < dim; i++) {
     vector(i) = tensor(i, i);
   }
 
   if (dim == 2) {
     vector(2) = tensor(0, 1);
   } else if (dim == 3) {
     vector(3) = tensor(0, 1);
     vector(4) = tensor(0, 2);
     vector(5) = tensor(1, 2);
   }
   
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 
 template<UInt dim>
 inline void MaterialReinforcement<dim>::strainTensorToVoigtVector(const Matrix<Real> & tensor,
                                                                   Vector<Real> & vector) {
   AKANTU_DEBUG_IN();
   
   for (UInt i = 0; i < dim; i++) {
     vector(i) = tensor(i, i);
   }
 
   if (dim == 2) {
     vector(2) = 2 * tensor(0, 1);
   } else if (dim == 3) {
     vector(3) = 2 * tensor(0, 1);
     vector(4) = 2 * tensor(0, 2);
     vector(5) = 2 * tensor(1, 2);
   }
   
   AKANTU_DEBUG_OUT();
 }
diff --git a/test/test_model/test_solid_mechanics_model/test_embedded_interface/CMakeLists.txt b/test/test_model/test_solid_mechanics_model/test_embedded_interface/CMakeLists.txt
index 31ac629ff..3081db0b6 100644
--- a/test/test_model/test_solid_mechanics_model/test_embedded_interface/CMakeLists.txt
+++ b/test/test_model/test_solid_mechanics_model/test_embedded_interface/CMakeLists.txt
@@ -1,49 +1,50 @@
 #===============================================================================
 # @file   CMakeLists.txt
 #
 # @author Lucas Frérot <lucas.frerot@epfl.ch>
 #
 # @date creation: Wed Mar 25 2015
 # @date last modification: Wed Mar 25 2015
 #
 # @brief  configuration for embedded interface tests
 #
 # @section LICENSE
 #
 # Copyright (©) 2015 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/>.
 #
 # @section DESCRIPTION
 #
 #===============================================================================
 
 register_test(test_embedded_element_matrix
   SOURCES test_embedded_element_matrix.cc
   FILES_TO_COPY triangle.msh embedded_element.dat
   PACKAGE embedded implicit
   )
 
 register_test(test_embedded_interface_model
   SOURCES test_embedded_interface_model.cc
   FILES_TO_COPY embedded_mesh.msh material.dat matrix
+  DIRECTORIES_TO_CREATE paraview
   PACKAGE embedded implicit
   )
 
 register_test(test_embedded_interface_model_prestress
   SOURCES test_embedded_interface_model_prestress.cc
   FILES_TO_COPY embedded_mesh_prestress.msh embedded_mesh_prestress_reinforcement.msh prestress.dat
   PACKAGE embedded implicit
   )
diff --git a/test/test_model/test_solid_mechanics_model/test_embedded_interface/test_embedded_interface_model.cc b/test/test_model/test_solid_mechanics_model/test_embedded_interface/test_embedded_interface_model.cc
index 47e1b4cb2..befd92f17 100644
--- a/test/test_model/test_solid_mechanics_model/test_embedded_interface/test_embedded_interface_model.cc
+++ b/test/test_model/test_solid_mechanics_model/test_embedded_interface/test_embedded_interface_model.cc
@@ -1,98 +1,105 @@
 /**
  * @file  test_embedded_interface_model.cc
  *
  * @author Lucas Frerot <lucas.frerot@epfl.ch>
  *
  * @date   lun. 09 29 16:03:10 2014
  *
  * @brief Embedded model test based on potential energy
  *
  * @section LICENSE
  *
  * Copyright (©) 2010-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 <iostream>
 
 #include "aka_common.hh"
 #include "embedded_interface_model.hh"
 
 using namespace akantu;
 
 int main (int argc, char * argv[]) {
   debug::setDebugLevel(dblWarning);
   initialize("material.dat", argc, argv);
 
   UInt dim = 2;
   Math::setTolerance(1e-7);
   
   // Mesh here is a 1x1 patch
   Mesh mesh(dim);
   mesh.read("embedded_mesh.msh");
 
   Array<Real> nodes_vec(2, dim, "reinforcement_nodes");
   nodes_vec.storage()[0] = 0; nodes_vec.storage()[1] = 0.5;
   nodes_vec.storage()[2] = 1; nodes_vec.storage()[3] = 0.5;
 
   Array<UInt> conn_vec(1, 2, "reinforcement_connectivity");
   conn_vec.storage()[0] = 0; conn_vec.storage()[1] = 1;
 
   Array<std::string> names_vec(1, 1, "reinforcement", "reinforcement_names");
 
   Mesh reinforcement_mesh(dim, "reinforcement_mesh");
   reinforcement_mesh.getNodes().copy(nodes_vec);
   reinforcement_mesh.addConnectivityType(_segment_2);
   reinforcement_mesh.getConnectivity(_segment_2).copy(conn_vec);
   reinforcement_mesh.registerData<std::string>("physical_names").alloc(1, 1, _segment_2);
   reinforcement_mesh.getData<std::string>("physical_names")(_segment_2).copy(names_vec);
 
   EmbeddedInterfaceModel model(mesh, reinforcement_mesh, dim);
   model.initFull(EmbeddedInterfaceModelOptions(_static));
 
   Array<Real> & nodes  = mesh.getNodes();
   Array<Real> & forces = model.getForce();
   Array<bool> & bound  = model.getBlockedDOFs();
 
   forces(2, 0) = -250;
   forces(5, 0) = -500;
   forces(8, 0) = -250;
 
   for (UInt i = 0 ; i < mesh.getNbNodes() ; i++) {
     if (Math::are_float_equal(nodes(i, 0), 0.))
       bound(i, 0) = true;
     if (Math::are_float_equal(nodes(i, 1), 0.))
       bound(i, 1) = true;
   }
 
+  model.addDumpFieldVector("displacement");
+  model.addDumpFieldTensor("stress");
+
+  model.setBaseNameToDumper("reinforcement", "reinforcement");
+  model.addDumpFieldTensorToDumper("reinforcement", "stress_embedded");
+
   // Assemble the global stiffness matrix
   model.assembleStiffnessMatrix();
   model.updateResidual();
 
   model.getStiffnessMatrix().saveMatrix("matrix_test");
 
   model.solveStatic<_scm_newton_raphson_tangent_not_computed, _scc_residual>(1e-7, 1);
   model.updateResidual();
+  model.dump();
 
   Real pot_energy = model.getEnergy("potential");
 
   if (std::abs(pot_energy - 7.37343e-06) > 1e-5)
     return EXIT_FAILURE;
 
   finalize();
   return 0;
 }