diff --git a/src/fe_engine/integrator_gauss_inline_impl.cc b/src/fe_engine/integrator_gauss_inline_impl.cc
index e1d54b647..ca419deda 100644
--- a/src/fe_engine/integrator_gauss_inline_impl.cc
+++ b/src/fe_engine/integrator_gauss_inline_impl.cc
@@ -1,712 +1,713 @@
 /**
  * @file   integrator_gauss_inline_impl.cc
  *
  * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Tue Feb 15 2011
  * @date last modification: Thu Nov 19 2015
  *
  * @brief  inline function of gauss integrator
  *
  * @section LICENSE
  *
  * Copyright (©)  2010-2012, 2014,  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/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #include "fe_engine.hh"
 #include "mesh_iterators.hh"
 #if defined(AKANTU_DEBUG_TOOLS)
 #include "aka_debug_tools.hh"
 #endif
 /* -------------------------------------------------------------------------- */
 
 namespace akantu {
 namespace debug {
 struct IntegratorGaussException : public Exception {};
 }
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind, class IntegrationOrderFunctor>
 template <ElementType type>
 inline void IntegratorGauss<kind, IntegrationOrderFunctor>::integrateOnElement(
     const Array<Real> & f, Real * intf, UInt nb_degree_of_freedom,
     const UInt elem, const GhostType & ghost_type) const {
   Array<Real> & jac_loc = jacobians(type, ghost_type);
 
   UInt nb_quadrature_points = ElementClass<type>::getNbQuadraturePoints();
   AKANTU_DEBUG_ASSERT(f.getNbComponent() == nb_degree_of_freedom,
-                      "The vector f do not have the good number of component.");
+		      "The vector f do not have the good number of component.");
 
   Real * f_val = f.storage() + elem * f.getNbComponent();
   Real * jac_val = jac_loc.storage() + elem * nb_quadrature_points;
 
   integrate(f_val, jac_val, intf, nb_degree_of_freedom, nb_quadrature_points);
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind, class IntegrationOrderFunctor>
 template <ElementType type>
 inline Real IntegratorGauss<kind, IntegrationOrderFunctor>::integrate(
     const Vector<Real> & in_f, UInt index, const GhostType & ghost_type) const {
   const Array<Real> & jac_loc = jacobians(type, ghost_type);
 
   UInt nb_quadrature_points =
       GaussIntegrationElement<type>::getNbQuadraturePoints();
   AKANTU_DEBUG_ASSERT(in_f.size() == nb_quadrature_points,
-                      "The vector f do not have nb_quadrature_points entries.");
+		      "The vector f do not have nb_quadrature_points entries.");
 
   Real * jac_val = jac_loc.storage() + index * nb_quadrature_points;
   Real intf;
 
   integrate(in_f.storage(), jac_val, &intf, 1, nb_quadrature_points);
 
   return intf;
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind, class IntegrationOrderFunctor>
 inline void IntegratorGauss<kind, IntegrationOrderFunctor>::integrate(
     Real * f, Real * jac, Real * inte, UInt nb_degree_of_freedom,
     UInt nb_quadrature_points) const {
   memset(inte, 0, nb_degree_of_freedom * sizeof(Real));
 
   Real * cjac = jac;
   for (UInt q = 0; q < nb_quadrature_points; ++q) {
     for (UInt dof = 0; dof < nb_degree_of_freedom; ++dof) {
       inte[dof] += *f * *cjac;
       ++f;
     }
     ++cjac;
   }
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind, class IntegrationOrderFunctor>
 template <ElementType type>
 inline const Matrix<Real> &
 IntegratorGauss<kind, IntegrationOrderFunctor>::getIntegrationPoints(
     const GhostType & ghost_type) const {
   AKANTU_DEBUG_ASSERT(
       quadrature_points.exists(type, ghost_type),
       "Quadrature points for type "
-          << quadrature_points.printType(type, ghost_type)
-          << " have not been initialized."
-          << " Did you use 'computeQuadraturePoints' function ?");
+	  << quadrature_points.printType(type, ghost_type)
+	  << " have not been initialized."
+	  << " Did you use 'computeQuadraturePoints' function ?");
   return quadrature_points(type, ghost_type);
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind, class IntegrationOrderFunctor>
 template <ElementType type>
 inline UInt
 IntegratorGauss<kind, IntegrationOrderFunctor>::getNbIntegrationPoints(
     const GhostType & ghost_type) const {
   AKANTU_DEBUG_ASSERT(
       quadrature_points.exists(type, ghost_type),
       "Quadrature points for type "
-          << quadrature_points.printType(type, ghost_type)
-          << " have not been initialized."
-          << " Did you use 'computeQuadraturePoints' function ?");
+	  << quadrature_points.printType(type, ghost_type)
+	  << " have not been initialized."
+	  << " Did you use 'computeQuadraturePoints' function ?");
   return quadrature_points(type, ghost_type).cols();
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind, class IntegrationOrderFunctor>
 template <ElementType type, UInt polynomial_degree>
 inline Matrix<Real>
 IntegratorGauss<kind, IntegrationOrderFunctor>::getIntegrationPoints() const {
   return GaussIntegrationElement<type,
-                                 polynomial_degree>::getQuadraturePoints();
+				 polynomial_degree>::getQuadraturePoints();
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind, class IntegrationOrderFunctor>
 template <ElementType type, UInt polynomial_degree>
 inline Vector<Real>
 IntegratorGauss<kind, IntegrationOrderFunctor>::getIntegrationWeights() const {
   return GaussIntegrationElement<type, polynomial_degree>::getWeights();
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind, class IntegrationOrderFunctor>
 template <ElementType type>
 inline void
 IntegratorGauss<kind, IntegrationOrderFunctor>::computeQuadraturePoints(
     const GhostType & ghost_type) {
   Matrix<Real> & quads = quadrature_points(type, ghost_type);
   const UInt polynomial_degree =
       IntegrationOrderFunctor::template getOrder<type>();
   quads =
       GaussIntegrationElement<type, polynomial_degree>::getQuadraturePoints();
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind, class IntegrationOrderFunctor>
 template <ElementType type>
 inline void IntegratorGauss<kind, IntegrationOrderFunctor>::
     computeJacobianOnQuadPointsByElement(const Matrix<Real> & node_coords,
-                                         const Matrix<Real> & quad,
-                                         Vector<Real> & jacobians) const {
+					 const Matrix<Real> & quad,
+					 Vector<Real> & jacobians) const {
   // jacobian
   ElementClass<type>::computeJacobian(quad, node_coords, jacobians);
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind, class IntegrationOrderFunctor>
 IntegratorGauss<kind, IntegrationOrderFunctor>::IntegratorGauss(
     const Mesh & mesh, const ID & id, const MemoryID & memory_id)
     : Integrator(mesh, id, memory_id) {
   AKANTU_DEBUG_IN();
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind, class IntegrationOrderFunctor>
 template <ElementType type>
 void IntegratorGauss<kind, IntegrationOrderFunctor>::checkJacobians(
     const GhostType & ghost_type) const {
   AKANTU_DEBUG_IN();
 
   UInt nb_quadrature_points = this->quadrature_points(type, ghost_type).cols();
 
   UInt nb_element = mesh.getConnectivity(type, ghost_type).size();
 
   Real * jacobians_val = jacobians(type, ghost_type).storage();
 
   for (UInt i = 0; i < nb_element * nb_quadrature_points;
        ++i, ++jacobians_val) {
     if (*jacobians_val < 0)
       AKANTU_CUSTOM_EXCEPTION_INFO(debug::IntegratorGaussException{},
-          "Negative jacobian computed,"
-          << " possible problem in the element node ordering (Quadrature Point "
-          << i % nb_quadrature_points << ":" << i / nb_quadrature_points << ":"
-          << type << ":" << ghost_type << ")");
+	  "Negative jacobian computed,"
+	  << " possible problem in the element node ordering (Quadrature Point "
+	  << i % nb_quadrature_points << ":" << i / nb_quadrature_points << ":"
+	  << type << ":" << ghost_type << ")");
   }
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind, class IntegrationOrderFunctor>
 template <ElementType type>
 void IntegratorGauss<kind, IntegrationOrderFunctor>::
     computeJacobiansOnIntegrationPoints(
-        const Array<Real> & nodes, const Matrix<Real> & quad_points,
-        Array<Real> & jacobians, const GhostType & ghost_type,
-        const Array<UInt> & filter_elements) const {
+	const Array<Real> & nodes, const Matrix<Real> & quad_points,
+	Array<Real> & jacobians, const GhostType & ghost_type,
+	const Array<UInt> & filter_elements) const {
   AKANTU_DEBUG_IN();
 
   UInt spatial_dimension = mesh.getSpatialDimension();
   UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
   UInt nb_quadrature_points = quad_points.cols();
 
   UInt nb_element = mesh.getNbElement(type, ghost_type);
   jacobians.resize(nb_element * nb_quadrature_points);
 
   auto jacobians_it =
       jacobians.begin_reinterpret(nb_quadrature_points, nb_element);
   auto jacobians_begin = jacobians_it;
 
   Array<Real> x_el(0, spatial_dimension * nb_nodes_per_element);
   FEEngine::extractNodalToElementField(mesh, nodes, x_el, type, ghost_type,
-                                       filter_elements);
+				       filter_elements);
 
   auto x_it = x_el.begin(spatial_dimension, nb_nodes_per_element);
 
   nb_element = x_el.size();
 
   //  Matrix<Real> local_coord(spatial_dimension, nb_nodes_per_element);
   for (UInt elem = 0; elem < nb_element; ++elem, ++x_it) {
     const Matrix<Real> & x = *x_it;
     Vector<Real> & J = *jacobians_it;
     computeJacobianOnQuadPointsByElement<type>(x, quad_points, J);
 
     if (filter_elements == empty_filter) {
       ++jacobians_it;
     } else {
       jacobians_it = jacobians_begin + filter_elements(elem);
     }
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 #if defined(AKANTU_STRUCTURAL_MECHANICS)
 template <>
 template <ElementType type>
 void IntegratorGauss<_ek_structural, DefaultIntegrationOrderFunctor>::
     computeJacobiansOnIntegrationPoints(
-        const Array<Real> & nodes, const Matrix<Real> & quad_points,
-        Array<Real> & jacobians, const GhostType & ghost_type,
-        const Array<UInt> & filter_elements) const {
+	const Array<Real> & nodes, const Matrix<Real> & quad_points,
+	Array<Real> & jacobians, const GhostType & ghost_type,
+	const Array<UInt> & filter_elements) const {
   AKANTU_DEBUG_IN();
 
   const UInt spatial_dimension = mesh.getSpatialDimension();
   const UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
   const UInt nb_quadrature_points = quad_points.cols();
   const UInt nb_dofs = ElementClass<type>::getNbDegreeOfFreedom();
 
   UInt nb_element = mesh.getNbElement(type, ghost_type);
   jacobians.resize(nb_element * nb_quadrature_points);
 
   auto jacobians_it =
       jacobians.begin_reinterpret(nb_quadrature_points, nb_element);
   auto jacobians_begin = jacobians_it;
 
   Array<Real> x_el(0, spatial_dimension * nb_nodes_per_element);
   FEEngine::extractNodalToElementField(mesh, nodes, x_el, type, ghost_type,
-                                       filter_elements);
+				       filter_elements);
 
   auto x_it = x_el.begin(spatial_dimension, nb_nodes_per_element);
 
   nb_element = x_el.size();
 
   const bool has_extra_normal =
       mesh.hasData<Real>("extra_normal", type, ghost_type);
   Array<Real>::vector_iterator extra_normal, extra_normal_begin;
-  if (has_extra_normal)
+  if (has_extra_normal) {
     extra_normal = mesh.getData<Real>("extra_normal", type, ghost_type)
-                       .begin(spatial_dimension);
-  extra_normal_begin = extra_normal;
+		       .begin(spatial_dimension);
+    extra_normal_begin = extra_normal;
+  }
 
   //  Matrix<Real> local_coord(spatial_dimension, nb_nodes_per_element);
   for (UInt elem = 0; elem < nb_element; ++elem, ++x_it) {
     const Matrix<Real> & X = *x_it;
     Vector<Real> & J = *jacobians_it;
     Matrix<Real> R(nb_dofs, nb_dofs);
 
     if (has_extra_normal)
       ElementClass<type>::computeRotationMatrix(R, X, *extra_normal);
     else
       ElementClass<type>::computeRotationMatrix(R, X, Vector<Real>(X.rows()));
     // Extracting relevant lines
     auto x = (R.block(0, 0, spatial_dimension, spatial_dimension) * X)
-                 .block(0, 0, ElementClass<type>::getNaturalSpaceDimension(),
-                        ElementClass<type>::getNbNodesPerElement());
+		 .block(0, 0, ElementClass<type>::getNaturalSpaceDimension(),
+			ElementClass<type>::getNbNodesPerElement());
 
     computeJacobianOnQuadPointsByElement<type>(x, quad_points, J);
 
     if (filter_elements == empty_filter) {
       ++jacobians_it;
       ++extra_normal;
     } else {
       jacobians_it = jacobians_begin + filter_elements(elem);
       extra_normal = extra_normal_begin + filter_elements(elem);
     }
   }
 
   AKANTU_DEBUG_OUT();
 }
 #endif
-  
+
 /* -------------------------------------------------------------------------- */
 #if defined(AKANTU_COHESIVE_ELEMENT)
 template <>
 template <ElementType type>
 void IntegratorGauss<_ek_cohesive, DefaultIntegrationOrderFunctor>::
     computeJacobiansOnIntegrationPoints(
-        const Array<Real> & nodes, const Matrix<Real> & quad_points,
-        Array<Real> & jacobians, const GhostType & ghost_type,
-        const Array<UInt> & filter_elements) const {
+	const Array<Real> & nodes, const Matrix<Real> & quad_points,
+	Array<Real> & jacobians, const GhostType & ghost_type,
+	const Array<UInt> & filter_elements) const {
   AKANTU_DEBUG_IN();
 
   UInt spatial_dimension = mesh.getSpatialDimension();
   UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
   UInt nb_quadrature_points = quad_points.cols();
 
   UInt nb_element = mesh.getNbElement(type, ghost_type);
 
   jacobians.resize(nb_element * nb_quadrature_points);
 
   auto jacobians_it =
       jacobians.begin_reinterpret(nb_quadrature_points, nb_element);
   auto jacobians_begin = jacobians_it;
 
   Vector<Real> weights = GaussIntegrationElement<type>::getWeights();
 
   Array<Real> x_el(0, spatial_dimension * nb_nodes_per_element);
   FEEngine::extractNodalToElementField(mesh, nodes, x_el, type, ghost_type,
-                                       filter_elements);
+				       filter_elements);
 
   auto x_it = x_el.begin(spatial_dimension, nb_nodes_per_element);
 
   UInt nb_nodes_per_subelement = nb_nodes_per_element / 2;
   Matrix<Real> x(spatial_dimension, nb_nodes_per_subelement);
 
   nb_element = x_el.size();
 
   auto compute = [&](const auto & el) {
     Vector<Real> J(jacobians_begin[el]);
 
     for (UInt s = 0; s < spatial_dimension; ++s)
       for (UInt n = 0; n < nb_nodes_per_subelement; ++n)
-        x(s, n) =
-            ((*x_it)(s, n) + (*x_it)(s, n + nb_nodes_per_subelement)) * .5;
+	x(s, n) =
+	    ((*x_it)(s, n) + (*x_it)(s, n + nb_nodes_per_subelement)) * .5;
 
     if (type == _cohesive_1d_2)
       J(0) = 1;
     else
       this->computeJacobianOnQuadPointsByElement<type>(x, quad_points, J);
   };
 
   for_each_elements(nb_element, filter_elements, compute);
 
 
   AKANTU_DEBUG_OUT();
 }
 #endif
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind, class IntegrationOrderFunctor>
 template <ElementType type>
 void IntegratorGauss<kind, IntegrationOrderFunctor>::
     precomputeJacobiansOnQuadraturePoints(const Array<Real> & nodes,
-                                          const GhostType & ghost_type) {
+					  const GhostType & ghost_type) {
   AKANTU_DEBUG_IN();
 
   Array<Real> & jacobians_tmp = jacobians.alloc(0, 1, type, ghost_type);
 
   this->computeJacobiansOnIntegrationPoints<type>(
       nodes, quadrature_points(type, ghost_type), jacobians_tmp, ghost_type);
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind, class IntegrationOrderFunctor>
 template <ElementType type, UInt polynomial_degree>
 void IntegratorGauss<kind, IntegrationOrderFunctor>::multiplyJacobiansByWeights(
     Array<Real> & jacobians) const {
   AKANTU_DEBUG_IN();
 
   UInt nb_quadrature_points =
       GaussIntegrationElement<type, polynomial_degree>::getNbQuadraturePoints();
   UInt nb_element = jacobians.size() / nb_quadrature_points;
 
   Vector<Real> weights =
       GaussIntegrationElement<type, polynomial_degree>::getWeights();
 
   Array<Real>::vector_iterator jacobians_it =
       jacobians.begin_reinterpret(nb_quadrature_points, nb_element);
   Array<Real>::vector_iterator jacobians_end =
       jacobians.end_reinterpret(nb_quadrature_points, nb_element);
 
   for (; jacobians_it != jacobians_end; ++jacobians_it) {
     Vector<Real> & J = *jacobians_it;
     J *= weights;
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind, class IntegrationOrderFunctor>
 void IntegratorGauss<kind, IntegrationOrderFunctor>::integrate(
     const Array<Real> & in_f, Array<Real> & intf, UInt nb_degree_of_freedom,
     const Array<Real> & jacobians, UInt nb_element) const {
   AKANTU_DEBUG_IN();
 
   intf.resize(nb_element);
   if (nb_element == 0)
     return;
 
   UInt nb_points = jacobians.size() / nb_element;
 
   Array<Real>::const_matrix_iterator J_it;
   Array<Real>::matrix_iterator inte_it;
   Array<Real>::const_matrix_iterator f_it;
 
   f_it = in_f.begin_reinterpret(nb_degree_of_freedom, nb_points, nb_element);
   inte_it = intf.begin_reinterpret(nb_degree_of_freedom, 1, nb_element);
   J_it = jacobians.begin_reinterpret(nb_points, 1, nb_element);
 
   for (UInt el = 0; el < nb_element; ++el, ++J_it, ++f_it, ++inte_it) {
     const Matrix<Real> & f = *f_it;
     const Matrix<Real> & J = *J_it;
     Matrix<Real> & inte_f = *inte_it;
 
     inte_f.mul<false, false>(f, J);
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind, class IntegrationOrderFunctor>
 template <ElementType type>
 void IntegratorGauss<kind, IntegrationOrderFunctor>::integrate(
     const Array<Real> & in_f, Array<Real> & intf, UInt nb_degree_of_freedom,
     const GhostType & ghost_type, const Array<UInt> & filter_elements) const {
   AKANTU_DEBUG_IN();
 
   AKANTU_DEBUG_ASSERT(jacobians.exists(type, ghost_type),
-                      "No jacobians for the type "
-                          << jacobians.printType(type, ghost_type));
+		      "No jacobians for the type "
+			  << jacobians.printType(type, ghost_type));
 
   const Array<Real> & jac_loc = jacobians(type, ghost_type);
   if (filter_elements != empty_filter) {
     UInt nb_element = filter_elements.size();
     Array<Real> * filtered_J = new Array<Real>(0, jac_loc.getNbComponent());
     FEEngine::filterElementalData(mesh, jac_loc, *filtered_J, type, ghost_type,
-                                  filter_elements);
+				  filter_elements);
     this->integrate(in_f, intf, nb_degree_of_freedom, *filtered_J, nb_element);
     delete filtered_J;
   } else {
     UInt nb_element = mesh.getNbElement(type, ghost_type);
     this->integrate(in_f, intf, nb_degree_of_freedom, jac_loc, nb_element);
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind, class IntegrationOrderFunctor>
 template <ElementType type, UInt polynomial_degree>
 void IntegratorGauss<kind, IntegrationOrderFunctor>::integrate(
     const Array<Real> & in_f, Array<Real> & intf, UInt nb_degree_of_freedom,
     const GhostType & ghost_type) const {
   AKANTU_DEBUG_IN();
 
   Matrix<Real> quads = this->getIntegrationPoints<type, polynomial_degree>();
 
   Array<Real> jacobians;
   this->computeJacobiansOnIntegrationPoints<type>(mesh.getNodes(), quads,
-                                                  jacobians, ghost_type);
+						  jacobians, ghost_type);
   this->multiplyJacobiansByWeights<type, polynomial_degree>(jacobians);
 
   this->integrate(in_f, intf, nb_degree_of_freedom, jacobians,
-                  mesh.getNbElement(type, ghost_type));
+		  mesh.getNbElement(type, ghost_type));
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind, class IntegrationOrderFunctor>
 template <ElementType type, UInt polynomial_degree>
 Real IntegratorGauss<kind, IntegrationOrderFunctor>::integrate(
     const Array<Real> & in_f, const GhostType & ghost_type) const {
   AKANTU_DEBUG_IN();
 
   Array<Real> intfv(0, 1);
   integrate<type, polynomial_degree>(in_f, intfv, 1, ghost_type);
 
   Real res = Math::reduce(intfv);
 
   AKANTU_DEBUG_OUT();
   return res;
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind, class IntegrationOrderFunctor>
 template <ElementType type>
 Real IntegratorGauss<kind, IntegrationOrderFunctor>::integrate(
     const Array<Real> & in_f, const GhostType & ghost_type,
     const Array<UInt> & filter_elements) const {
   AKANTU_DEBUG_IN();
 
   AKANTU_DEBUG_ASSERT(jacobians.exists(type, ghost_type),
-                      "No jacobians for the type "
-                          << jacobians.printType(type, ghost_type));
+		      "No jacobians for the type "
+			  << jacobians.printType(type, ghost_type));
 
   Array<Real> intfv(0, 1);
   integrate<type>(in_f, intfv, 1, ghost_type, filter_elements);
 
   Real res = Math::reduce(intfv);
 
   AKANTU_DEBUG_OUT();
   return res;
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind, class IntegrationOrderFunctor>
 void IntegratorGauss<kind, IntegrationOrderFunctor>::
     integrateOnIntegrationPoints(const Array<Real> & in_f, Array<Real> & intf,
-                                 UInt nb_degree_of_freedom,
-                                 const Array<Real> & jacobians,
-                                 UInt nb_element) const {
+				 UInt nb_degree_of_freedom,
+				 const Array<Real> & jacobians,
+				 UInt nb_element) const {
   AKANTU_DEBUG_IN();
 
   UInt nb_points = jacobians.size() / nb_element;
 
   Array<Real>::const_scalar_iterator J_it;
   Array<Real>::vector_iterator inte_it;
   Array<Real>::const_vector_iterator f_it;
 
   intf.resize(nb_element * nb_points);
 
   J_it = jacobians.begin();
   f_it = in_f.begin(nb_degree_of_freedom);
   inte_it = intf.begin(nb_degree_of_freedom);
 
   for (UInt el = 0; el < nb_element; ++el, ++J_it, ++f_it, ++inte_it) {
     const Real & J = *J_it;
     const Vector<Real> & f = *f_it;
     Vector<Real> & inte_f = *inte_it;
 
     inte_f = f;
     inte_f *= J;
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind, class IntegrationOrderFunctor>
 template <ElementType type>
 void IntegratorGauss<kind, IntegrationOrderFunctor>::
     integrateOnIntegrationPoints(const Array<Real> & in_f, Array<Real> & intf,
-                                 UInt nb_degree_of_freedom,
-                                 const GhostType & ghost_type,
-                                 const Array<UInt> & filter_elements) const {
+				 UInt nb_degree_of_freedom,
+				 const GhostType & ghost_type,
+				 const Array<UInt> & filter_elements) const {
   AKANTU_DEBUG_IN();
 
   AKANTU_DEBUG_ASSERT(jacobians.exists(type, ghost_type),
-                      "No jacobians for the type "
-                          << jacobians.printType(type, ghost_type));
+		      "No jacobians for the type "
+			  << jacobians.printType(type, ghost_type));
 
   const Array<Real> & jac_loc = this->jacobians(type, ghost_type);
 
   if (filter_elements != empty_filter) {
 
     UInt nb_element = filter_elements.size();
     Array<Real> * filtered_J = new Array<Real>(0, jac_loc.getNbComponent());
     FEEngine::filterElementalData(mesh, jac_loc, *filtered_J, type, ghost_type,
-                                  filter_elements);
+				  filter_elements);
 
     this->integrateOnIntegrationPoints(in_f, intf, nb_degree_of_freedom,
-                                       *filtered_J, nb_element);
+				       *filtered_J, nb_element);
   } else {
     UInt nb_element = mesh.getNbElement(type, ghost_type);
     this->integrateOnIntegrationPoints(in_f, intf, nb_degree_of_freedom,
-                                       jac_loc, nb_element);
+				       jac_loc, nb_element);
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind, class IntegrationOrderFunctor>
 void IntegratorGauss<kind, IntegrationOrderFunctor>::onElementsAdded(
     const Array<Element> & new_elements) {
 
   const auto & nodes = mesh.getNodes();
 
   for (auto elements_range : MeshElementsByTypes(new_elements)) {
     auto type = elements_range.getType();
     auto ghost_type = elements_range.getGhostType();
 
     if (mesh.getKind(type) != kind)
       continue;
 
     auto & elements = elements_range.getElements();
 
     if (not jacobians.exists(type, ghost_type))
       jacobians.alloc(0, 1, type, ghost_type);
 
     this->computeJacobiansOnIntegrationPoints(
-        nodes, quadrature_points(type, ghost_type), jacobians(type, ghost_type),
-        type, ghost_type, elements);
+	nodes, quadrature_points(type, ghost_type), jacobians(type, ghost_type),
+	type, ghost_type, elements);
   }
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind, class IntegrationOrderFunctor>
 template <ElementType type>
 inline void IntegratorGauss<kind, IntegrationOrderFunctor>::initIntegrator(
     const Array<Real> & nodes, const GhostType & ghost_type) {
   computeQuadraturePoints<type>(ghost_type);
   precomputeJacobiansOnQuadraturePoints<type>(nodes, ghost_type);
   checkJacobians<type>(ghost_type);
   constexpr UInt polynomial_degree =
       IntegrationOrderFunctor::template getOrder<type>();
   multiplyJacobiansByWeights<type, polynomial_degree>(
       this->jacobians(type, ghost_type));
 }
 
 namespace integrator {
   namespace details {
     template <ElementKind kind> struct GaussIntegratorInitHelper {};
 
 #define INIT_INTEGRATOR(type)                                                  \
   _int.template initIntegrator<type>(nodes, ghost_type)
 
 #define AKANTU_GAUSS_INTERGRATOR_INIT_HELPER(kind)                             \
   template <> struct GaussIntegratorInitHelper<kind> {                         \
     template <ElementKind k, class IOF>                                        \
     static void call(IntegratorGauss<k, IOF> & _int,                           \
-                     const Array<Real> & nodes, const ElementType & type,      \
-                     const GhostType & ghost_type) {                           \
+		     const Array<Real> & nodes, const ElementType & type,      \
+		     const GhostType & ghost_type) {                           \
       AKANTU_BOOST_KIND_ELEMENT_SWITCH(INIT_INTEGRATOR, kind);                 \
     }                                                                          \
   };
 
     AKANTU_BOOST_ALL_KIND(AKANTU_GAUSS_INTERGRATOR_INIT_HELPER)
 
 #undef AKANTU_GAUSS_INTERGRATOR_INIT_HELPER
 #undef INIT_INTEGRATOR
   } // namespace details
 } // namespace integrator
 
 template <ElementKind kind, class IntegrationOrderFunctor>
 inline void IntegratorGauss<kind, IntegrationOrderFunctor>::initIntegrator(
     const Array<Real> & nodes, const ElementType & type,
     const GhostType & ghost_type) {
   integrator::details::GaussIntegratorInitHelper<kind>::call(*this, nodes, type,
-                                                             ghost_type);
+							     ghost_type);
 }
 
 namespace integrator {
   namespace details {
     template <ElementKind kind> struct GaussIntegratorComputeJacobiansHelper {};
 
 #define AKANTU_COMPUTE_JACOBIANS(type)                                         \
   _int.template computeJacobiansOnIntegrationPoints<type>(                     \
       nodes, quad_points, jacobians, ghost_type, filter_elements);
 
 #define AKANTU_GAUSS_INTERGRATOR_COMPUTE_JACOBIANS(kind)                       \
   template <> struct GaussIntegratorComputeJacobiansHelper<kind> {             \
     template <ElementKind k, class IOF>                                        \
     static void                                                                \
     call(const IntegratorGauss<k, IOF> & _int, const Array<Real> & nodes,      \
-         const Matrix<Real> & quad_points, Array<Real> & jacobians,            \
-         const ElementType & type, const GhostType & ghost_type,               \
-         const Array<UInt> & filter_elements) {                                \
+	 const Matrix<Real> & quad_points, Array<Real> & jacobians,            \
+	 const ElementType & type, const GhostType & ghost_type,               \
+	 const Array<UInt> & filter_elements) {                                \
       AKANTU_BOOST_KIND_ELEMENT_SWITCH(AKANTU_COMPUTE_JACOBIANS, kind);        \
     }                                                                          \
   };
 
     AKANTU_BOOST_ALL_KIND(AKANTU_GAUSS_INTERGRATOR_COMPUTE_JACOBIANS)
 
 #undef AKANTU_GAUSS_INTERGRATOR_COMPUTE_JACOBIANS
 #undef AKANTU_COMPUTE_JACOBIANS
   } // namespace details
 } // namespace integrator
 
 template <ElementKind kind, class IntegrationOrderFunctor>
 void IntegratorGauss<kind, IntegrationOrderFunctor>::
     computeJacobiansOnIntegrationPoints(
-        const Array<Real> & nodes, const Matrix<Real> & quad_points,
-        Array<Real> & jacobians, const ElementType & type,
-        const GhostType & ghost_type,
-        const Array<UInt> & filter_elements) const {
+	const Array<Real> & nodes, const Matrix<Real> & quad_points,
+	Array<Real> & jacobians, const ElementType & type,
+	const GhostType & ghost_type,
+	const Array<UInt> & filter_elements) const {
   integrator::details::GaussIntegratorComputeJacobiansHelper<kind>::call(
       *this, nodes, quad_points, jacobians, type, ghost_type, filter_elements);
 }
 
 } // namespace akantu