diff --git a/python/swig/akantu.i b/python/swig/akantu.i
index abfd3456c..a18640324 100644
--- a/python/swig/akantu.i
+++ b/python/swig/akantu.i
@@ -1,84 +1,85 @@
 /**
  * @file   akantu.i
  *
  * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
  * @author Fabian Barras <fabian.barras@epfl.ch>
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Fri Dec 12 2014
  * @date last modification: Mon Nov 23 2015
  *
  * @brief  Main swig file for akantu' python interface
  *
  * @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/>.
  *
  */
 
 %module akantu
 
 %exception {
   try {
     $action
   } catch (akantu::debug::Exception e) {
     PyErr_SetString(PyExc_IndexError,e.what());
     return NULL;
   }
  }
 
 #define __attribute__(x)
 
 %ignore akantu::operator <<;
 
 %include "aka_common.i"
 %include "aka_csr.i"
 %include "aka_array.i"
 
 %define print_self(MY_CLASS)
   %extend akantu::MY_CLASS {
     std::string __str__() {
       std::stringstream sstr;
       sstr << *($self);
       return sstr.str();
     }
  }
 %enddef
 
 %include "mesh.i"
 %include "mesh_utils.i"
+%include "fe_engine.i"
 %include "model.i"
 %include "communicator.i"
 
 %include "solid_mechanics_model.i"
 #if defined(AKANTU_COHESIVE_ELEMENT)
 %include "solid_mechanics_model_cohesive.i"
 #endif
 
 #if defined(AKANTU_HEAT_TRANSFER)
 %include "heat_transfer_model.i"
 #endif
 
 
 #if defined(AKANTU_STRUCTURAL_MECHANICS)
 %include "load_functions.i"
 %include "structural_mechanics_model.i"
 #endif
 
 %pythoncode %{
   __initialize()
 %}
diff --git a/python/swig/fe_engine.i b/python/swig/fe_engine.i
new file mode 100644
index 000000000..51d7b6b91
--- /dev/null
+++ b/python/swig/fe_engine.i
@@ -0,0 +1,53 @@
+/**
+ * @file   fe_engine.i
+ *
+ * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
+ * @author Aurelia Isabel Cuba Ramos <aurelia.cubaramos@epfl.ch>
+ * @author Nicolas Richart <nicolas.richart@epfl.ch>
+ *
+ * @date creation: Fri Dec 12 2014
+ * @date last modification: Wed Nov 11 2015
+ *
+ * @brief  FEEngine wrapper
+ *
+ * @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/>.
+ *
+ */
+
+%{
+  #include "mesh.hh"
+  %}
+
+namespace akantu {
+  %ignore FEEngine::getIGFEMElementTypes;
+  %ignore FEEngine::interpolateOnIntegrationPoints(const Array<Real> &,ElementTypeMapArray<Real> &,const ElementTypeMapArray<UInt> *) const;
+  %ignore FEEngine::interpolateOnIntegrationPoints(const Array<Real> &,ElementTypeMapArray<Real> &) const;
+  %ignore FEEngine::interpolateOnIntegrationPoints(const Array<Real> &,Array<Real> &,UInt,const ElementType&,const GhostType &,const Array< UInt > &) const;
+  %ignore FEEngine::interpolateOnIntegrationPoints(const Array<Real> &,Array<Real> &,UInt,const ElementType&,const GhostType &) const;
+  %ignore FEEngine::onNodesAdded;
+  %ignore FEEngine::onNodesRemoved;
+  %ignore FEEngine::onElementsAdded;
+  %ignore FEEngine::onElementsChanged;
+  %ignore FEEngine::onElementsRemoved;
+  %ignore FEEngine::elementTypes;
+  %ignore
+}
+
+%include "sparse_matrix.i"
+%include "fe_engine.hh"
diff --git a/python/swig/mesh.i b/python/swig/mesh.i
index 08016e7f9..0d7db0204 100644
--- a/python/swig/mesh.i
+++ b/python/swig/mesh.i
@@ -1,166 +1,166 @@
 /**
  * @file   mesh.i
  *
  * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
  * @author Fabian Barras <fabian.barras@epfl.ch>
  * @author Aurelia Isabel Cuba Ramos <aurelia.cubaramos@epfl.ch>
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Fri Dec 12 2014
  * @date last modification: Wed Jan 13 2016
  *
  * @brief  mesh wrapper
  *
  * @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/>.
  *
  */
 
 %{
 #include "mesh.hh"
 #include "node_group.hh"
 #include "solid_mechanics_model.hh"
 #include "python_functor.hh"
 #include "mesh_utils.hh"
   
 using akantu::IntegrationPoint;
 using akantu::Vector;
 using akantu::ElementTypeMapArray;
 using akantu::MatrixProxy;
 using akantu::Matrix;
 using akantu::UInt;
 using akantu::Real;
 using akantu::Array;
 using akantu::SolidMechanicsModel;
 %}
 
 
 namespace akantu {
   %ignore NewNodesEvent;
   %ignore RemovedNodesEvent;
   %ignore NewElementsEvent;
   %ignore RemovedElementsEvent;
   %ignore MeshEventHandler;
   %ignore MeshEvent< UInt >;
   %ignore MeshEvent< Element >;
   %ignore Mesh::extractNodalCoordinatesFromPBCElement;
   %ignore Mesh::getGroupDumer;
   %ignore Mesh::getFacetLocalConnectivity;
   %ignore Mesh::getAllFacetTypes;
   %ignore Mesh::getCommunicator;
   %ignore GroupManager::getElementGroups;
   %ignore Dumpable::addDumpFieldExternalReal;
 }
 
 print_self(Mesh)
 
 // Swig considers enums to be ints, and it creates a conflict with two versions of getNbElement()
 %rename(getNbElementByDimension)
 akantu::Mesh::getNbElement(const UInt spatial_dimension = _all_dimensions,
                            const GhostType& ghost_type = _not_ghost,
                            const ElementKind& kind = _ek_not_defined) const;
 
 %extend akantu::Mesh {
 
   PyObject * getElementGroups(){
     return akantu::PythonFunctor::convertToPython($self->getElementGroups());
   }
 
   void resizeMesh(UInt nb_nodes, UInt nb_element, const ElementType & type) {
     Array<Real> & nodes = const_cast<Array<Real> &>($self->getNodes());
     nodes.resize(nb_nodes);
 
     $self->addConnectivityType(type);
     Array<UInt> & connectivity = const_cast<Array<UInt> &>($self->getConnectivity(type));
     connectivity.resize(nb_element);
   }
 }
 
 %extend akantu::GroupManager {
   void createGroupsFromStringMeshData(const std::string & dataset_name) {
     if (dataset_name == "physical_names"){
       AKANTU_EXCEPTION("Deprecated behavior: no need to call 'createGroupsFromStringMeshData' for physical names");
     }
     $self->createGroupsFromMeshData<std::string>(dataset_name);
   }
 
   void createGroupsFromUIntMeshData(const std::string & dataset_name) {
     $self->createGroupsFromMeshData<akantu::UInt>(dataset_name);
   }
 }
 
 %extend akantu::NodeGroup {
-  akantu::Array<akantu::Real> & getGroupedNodes(akantu::Array<akantu::Real, true> & surface_array, Mesh & mesh) {
+    akantu::Array<akantu::Real> & getGroupedNodes(akantu::Array<akantu::Real, true> & surface_array, Mesh & mesh) {
     akantu::Array<akantu::UInt> group_node = $self->getNodes();
     akantu::Array<akantu::Real> & full_array = mesh.getNodes();
     surface_array.resize(group_node.size());
 
     for (UInt i = 0; i < group_node.size(); ++i) {
       for (UInt cmp = 0; cmp < full_array.getNbComponent(); ++cmp) {
 
         surface_array(i,cmp) = full_array(group_node(i),cmp);
       }
     }
 
     akantu::Array<akantu::Real> & res(surface_array);
     return res;
   }
 
   akantu::Array<akantu::Real> & getGroupedArray(akantu::Array<akantu::Real, true> & surface_array, akantu::SolidMechanicsModel & model, int type) {
     akantu::Array<akantu::Real> * full_array;
 
     switch (type) {
 
     case 0 : full_array = new akantu::Array<akantu::Real>(model.getDisplacement());
       break;
     case 1 : full_array = new akantu::Array<akantu::Real>(model.getVelocity());
       break;
     case 2 : full_array = new akantu::Array<akantu::Real>(model.getForce());
       break;
     }
     akantu::Array<akantu::UInt> group_node = $self->getNodes();
     surface_array.resize(group_node.size());
 
     for (UInt i = 0; i < group_node.size(); ++i) {
       for (UInt cmp = 0; cmp < full_array->getNbComponent(); ++cmp) {
 
         surface_array(i,cmp) = (*full_array)(group_node(i),cmp);
       }
     }
 
     akantu::Array<akantu::Real> & res(surface_array);
     return res;
   }
 }
 
 %include "group_manager.hh"
 %include "node_group.hh"
 %include "dumper_iohelper.hh"
 %include "dumpable_iohelper.hh"
 %include "element_group.hh"
 %include "mesh.hh"
 %include "mesh_utils.hh"
 
 namespace akantu{
 %extend Dumpable {
     void addDumpFieldExternalReal(const std::string & field_id,
                                   const Array<Real> & field){
       $self->addDumpFieldExternal<Real>(field_id,field);
     }
   }
  }
diff --git a/python/swig/model.i b/python/swig/model.i
index 8f0353180..4cb6acdba 100644
--- a/python/swig/model.i
+++ b/python/swig/model.i
@@ -1,132 +1,116 @@
 /**
  * @file   model.i
  *
  * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
  * @author Aurelia Isabel Cuba Ramos <aurelia.cubaramos@epfl.ch>
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Fri Dec 12 2014
  * @date last modification: Wed Nov 11 2015
  *
  * @brief  model wrapper
  *
  * @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/>.
  *
  */
 
 %{
   #include "boundary_condition_python_functor.hh"
   #include "model_solver.hh"
   #include "non_linear_solver.hh"
   %}
 
 
 namespace akantu {
   %ignore Model::createSynchronizerRegistry;
   %ignore Model::getSynchronizerRegistry;
   %ignore Model::createParallelSynch;
   %ignore Model::getDOFSynchronizer;
   %ignore Model::registerFEEngineObject;
   %ignore Model::unregisterFEEngineObject;
-  %ignore Model::getFEEngineBoundary;
+  //  %ignore Model::getFEEngineBoundary;
   //  %ignore Model::getFEEngine;
   %ignore Model::getFEEngineClass;
   %ignore Model::getFEEngineClassBoundary;
   %ignore Model::setParser;
   %ignore Model::updateDataForNonLocalCriterion;
   %ignore IntegrationPoint::operator=;
-
-  %ignore FEEngine::getNbIntegrationPoints;
-  %ignore FEEngine::getShapes;
-  %ignore FEEngine::getShapesDerivatives;
-  %ignore FEEngine::getIntegrationPoints;
-  %ignore FEEngine::getIGFEMElementTypes;
-  %ignore FEEngine::interpolateOnIntegrationPoints(const Array<Real> &,ElementTypeMapArray<Real> &,const ElementTypeMapArray<UInt> *) const;
-  %ignore FEEngine::interpolateOnIntegrationPoints(const Array<Real> &,ElementTypeMapArray<Real> &) const;
-  %ignore FEEngine::interpolateOnIntegrationPoints(const Array<Real> &,Array<Real> &,UInt,const ElementType&,const GhostType &,const Array< UInt > &) const;
-  %ignore FEEngine::interpolateOnIntegrationPoints(const Array<Real> &,Array<Real> &,UInt,const ElementType&,const GhostType &) const;
-  %ignore FEEngine::onNodesAdded;
-  %ignore FEEngine::onNodesRemoved;
-  %ignore FEEngine::onElementsAdded;
-  %ignore FEEngine::onElementsChanged;
-  %ignore FEEngine::onElementsRemoved;
-  %ignore FEEngine::elementTypes;
 }
 
 %include "sparse_matrix.i"
 %include "fe_engine.hh"
 
 %rename(FreeBoundaryDirichlet) akantu::BC::Dirichlet::FreeBoundary;
 %rename(FreeBoundaryNeumann) akantu::BC::Neumann::FreeBoundary;
 %rename(PythonBoundary) akantu::BC::Dirichlet::PythonFunctor;
 
 %include "boundary_condition_functor.hh"
 %include "boundary_condition.hh"
 %include "boundary_condition_python_functor.hh"
 %include "communication_buffer.hh"
 %include "data_accessor.hh"
 //%include "synchronizer.hh"
 //%include "synchronizer_registry.hh"
 %include "model.hh"
 %include "non_linear_solver.hh"
 %include "model_options.hh"
 
 %extend akantu::Model {
   void initFullImpl(
        const akantu::ModelOptions & options = akantu::ModelOptions()){
      $self->initFull(options);
   };
 
   %insert("python") %{
     def initFull(self, *args, **kwargs):
         if len(args) == 0:
             import importlib as __aka_importlib
             _module =  __aka_importlib.import_module(self.__module__)
             _cls = getattr(_module, "{0}Options".format(self.__class__.__name__))
             options = _cls()
             if len(kwargs) > 0:
                 for key, val in kwargs.items():
                     if key[0] == '_':
                         key = key[1:]
                     setattr(options, key, val)
         else:
             options = args[0]
 
         self.initFullImpl(options)
   %}
 
 
   void solveStep(const akantu::ID & id=""){
     $self->solveStep(id);
   }
 
   akantu::NonLinearSolver & getNonLinearSolver(const akantu::ID & id=""){
    return $self->getNonLinearSolver(id);
   }
  }
 
 %extend akantu::NonLinearSolver {
   void set(const std::string & id, akantu::Real val){
     if (id == "max_iterations")
       $self->set(id, int(val));
     else if (id == "convergence_type")
       $self->set(id, akantu::SolveConvergenceCriteria(UInt(val)));
     else $self->set(id, val);
   }
  }
diff --git a/src/fe_engine/fe_engine.hh b/src/fe_engine/fe_engine.hh
index adcd28bb9..a4056daed 100644
--- a/src/fe_engine/fe_engine.hh
+++ b/src/fe_engine/fe_engine.hh
@@ -1,385 +1,354 @@
 /**
  * @file   fe_engine.hh
  *
  * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Fri Jun 18 2010
  * @date last modification: Tue Feb 20 2018
  *
  * @brief  FEM class
  *
  * @section LICENSE
  *
  * Copyright (©)  2010-2018 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_memory.hh"
 #include "element_type_map.hh"
 #include "mesh_events.hh"
 /* -------------------------------------------------------------------------- */
 
 #ifndef __AKANTU_FE_ENGINE_HH__
 #define __AKANTU_FE_ENGINE_HH__
 
 namespace akantu {
 class Mesh;
 class Integrator;
 class ShapeFunctions;
 class DOFManager;
 class Element;
 } // namespace akantu
 
 /* -------------------------------------------------------------------------- */
 namespace akantu {
 /* -------------------------------------------------------------------------- */
 
 /**
  * The  generic  FEEngine class  derived  in  a  FEEngineTemplate class
  * containing  the
  * shape functions and the integration method
  */
 class FEEngine : protected Memory, public MeshEventHandler {
   /* ------------------------------------------------------------------------ */
   /* Constructors/Destructors                                                 */
   /* ------------------------------------------------------------------------ */
 public:
   FEEngine(Mesh & mesh, UInt spatial_dimension = _all_dimensions,
            const ID & id = "fem", MemoryID memory_id = 0);
 
   ~FEEngine() override;
 
   /* ------------------------------------------------------------------------ */
   /* Methods                                                                  */
   /* ------------------------------------------------------------------------ */
 public:
   /// pre-compute all the shape functions, their derivatives and the jacobians
   virtual void
   initShapeFunctions(const GhostType & ghost_type = _not_ghost) = 0;
 
   /// extract the nodal values and store them per element
   template <typename T>
   static void extractNodalToElementField(
       const Mesh & mesh, const Array<T> & nodal_f, Array<T> & elemental_f,
       const ElementType & type, const GhostType & ghost_type = _not_ghost,
       const Array<UInt> & filter_elements = empty_filter);
 
   /// filter a field
   template <typename T>
   static void
   filterElementalData(const Mesh & mesh, const Array<T> & quad_f,
                       Array<T> & filtered_f, const ElementType & type,
                       const GhostType & ghost_type = _not_ghost,
                       const Array<UInt> & filter_elements = empty_filter);
 
   /* ------------------------------------------------------------------------ */
   /* Integration method bridges                                               */
   /* ------------------------------------------------------------------------ */
   /// integrate f for all elements of type "type"
   virtual void
   integrate(const Array<Real> & f, Array<Real> & intf,
             UInt nb_degree_of_freedom, const ElementType & type,
             const GhostType & ghost_type = _not_ghost,
             const Array<UInt> & filter_elements = empty_filter) const = 0;
 
   /// integrate a scalar value f on all elements of type "type"
   virtual Real
   integrate(const Array<Real> & f, const ElementType & type,
             const GhostType & ghost_type = _not_ghost,
             const Array<UInt> & filter_elements = empty_filter) const = 0;
 
   /// integrate f for all integration points of type "type" but don't sum over
   /// all integration points
   virtual void integrateOnIntegrationPoints(
       const Array<Real> & f, Array<Real> & intf, UInt nb_degree_of_freedom,
       const ElementType & type, const GhostType & ghost_type = _not_ghost,
       const Array<UInt> & filter_elements = empty_filter) const = 0;
 
   /// integrate one element scalar value on all elements of type "type"
   virtual Real integrate(const Vector<Real> & f, const ElementType & type,
                          UInt index,
                          const GhostType & ghost_type = _not_ghost) const = 0;
 
 /* ------------------------------------------------------------------------ */
 /* compatibility with old FEEngine fashion */
 /* ------------------------------------------------------------------------ */
-#ifndef SWIG
   /// get the number of integration points
   virtual UInt
   getNbIntegrationPoints(const ElementType & type,
                          const GhostType & ghost_type = _not_ghost) const = 0;
+
+#ifndef SWIG
   /// get the precomputed shapes
   const virtual Array<Real> &
   getShapes(const ElementType & type, const GhostType & ghost_type = _not_ghost,
             UInt id = 0) const = 0;
 
   /// get the derivatives of shapes
   const virtual Array<Real> &
   getShapesDerivatives(const ElementType & type,
                        const GhostType & ghost_type = _not_ghost,
                        UInt id = 0) const = 0;
 
   /// get integration points
   const virtual Matrix<Real> &
   getIntegrationPoints(const ElementType & type,
                        const GhostType & ghost_type = _not_ghost) const = 0;
 #endif
   /* ------------------------------------------------------------------------ */
   /* Shape method bridges                                                     */
   /* ------------------------------------------------------------------------ */
   /// Compute the gradient nablauq on the integration points of an element type
   /// from nodal values u
   virtual void gradientOnIntegrationPoints(
       const Array<Real> & u, Array<Real> & nablauq,
       const UInt nb_degree_of_freedom, const ElementType & type,
       const GhostType & ghost_type = _not_ghost,
       const Array<UInt> & filter_elements = empty_filter) const = 0;
 
   /// Interpolate a nodal field u at the integration points of an element type
   /// -> uq
   virtual void interpolateOnIntegrationPoints(
       const Array<Real> & u, Array<Real> & uq, UInt nb_degree_of_freedom,
       const ElementType & type, const GhostType & ghost_type = _not_ghost,
       const Array<UInt> & filter_elements = empty_filter) const = 0;
 
   /// Interpolate a nodal field u at the integration points of many element
   /// types -> uq
   virtual void interpolateOnIntegrationPoints(
       const Array<Real> & u, ElementTypeMapArray<Real> & uq,
       const ElementTypeMapArray<UInt> * filter_elements = nullptr) const = 0;
 
   /// pre multiplies a tensor by the shapes derivaties
   virtual void
   computeBtD(const Array<Real> & Ds, Array<Real> & BtDs,
              const ElementType & type, const GhostType & ghost_type,
              const Array<UInt> & filter_elements = empty_filter) const = 0;
 
   /// left and right  multiplies a tensor by the shapes derivaties
   virtual void
   computeBtDB(const Array<Real> & Ds, Array<Real> & BtDBs, UInt order_d,
               const ElementType & type, const GhostType & ghost_type,
               const Array<UInt> & filter_elements = empty_filter) const = 0;
 
   /// left multiples a vector by the shape functions
   virtual void computeNtb(const Array<Real> & bs, Array<Real> & Ntbs,
                           const ElementType & type,
                           const GhostType & ghost_type,
                           const Array<UInt> & filter_elements) const = 0;
 
   /// Compute the interpolation point position in the global coordinates for
   /// many element types
   virtual void computeIntegrationPointsCoordinates(
       ElementTypeMapArray<Real> & integration_points_coordinates,
       const ElementTypeMapArray<UInt> * filter_elements = nullptr) const = 0;
 
   /// Compute the interpolation point position in the global coordinates for an
   /// element type
   virtual void computeIntegrationPointsCoordinates(
       Array<Real> & integration_points_coordinates, const ElementType & type,
       const GhostType & ghost_type = _not_ghost,
       const Array<UInt> & filter_elements = empty_filter) const = 0;
 
   /// Build pre-computed matrices for interpolation of field form integration
   /// points at other given positions (interpolation_points)
   virtual void initElementalFieldInterpolationFromIntegrationPoints(
       const ElementTypeMapArray<Real> & interpolation_points_coordinates,
       ElementTypeMapArray<Real> & interpolation_points_coordinates_matrices,
       ElementTypeMapArray<Real> & integration_points_coordinates_inv_matrices,
       const ElementTypeMapArray<UInt> * element_filter) const = 0;
 
   /// interpolate field at given position (interpolation_points) from given
   /// values of this field at integration points (field)
   virtual void interpolateElementalFieldFromIntegrationPoints(
       const ElementTypeMapArray<Real> & field,
       const ElementTypeMapArray<Real> & interpolation_points_coordinates,
       ElementTypeMapArray<Real> & result, const GhostType ghost_type,
       const ElementTypeMapArray<UInt> * element_filter) const = 0;
 
   /// Interpolate field at given position from given values of this field at
   /// integration points (field)
   /// using matrices precomputed with
   /// initElementalFieldInterplationFromIntegrationPoints
   virtual void interpolateElementalFieldFromIntegrationPoints(
       const ElementTypeMapArray<Real> & field,
       const ElementTypeMapArray<Real> &
           interpolation_points_coordinates_matrices,
       const ElementTypeMapArray<Real> &
           integration_points_coordinates_inv_matrices,
       ElementTypeMapArray<Real> & result, const GhostType ghost_type,
       const ElementTypeMapArray<UInt> * element_filter) const = 0;
 
   /// interpolate on a phyiscal point inside an element
   virtual void interpolate(const Vector<Real> & real_coords,
                            const Matrix<Real> & nodal_values,
                            Vector<Real> & interpolated,
                            const Element & element) const = 0;
 
   /// compute the shape on a provided point
   virtual void
   computeShapes(const Vector<Real> & real_coords, UInt elem,
                 const ElementType & type, Vector<Real> & shapes,
                 const GhostType & ghost_type = _not_ghost) const = 0;
 
   /// compute the shape derivatives on a provided point
   virtual void
   computeShapeDerivatives(const Vector<Real> & real__coords, UInt element,
                           const ElementType & type,
                           Matrix<Real> & shape_derivatives,
                           const GhostType & ghost_type = _not_ghost) const = 0;
 
   /* ------------------------------------------------------------------------ */
   /* Other methods                                                            */
   /* ------------------------------------------------------------------------ */
 
   /// pre-compute normals on integration points
   virtual void computeNormalsOnIntegrationPoints(
       const GhostType & ghost_type = _not_ghost) = 0;
 
   /// pre-compute normals on integration points
   virtual void computeNormalsOnIntegrationPoints(
       __attribute__((unused)) const Array<Real> & field,
       __attribute__((unused)) const GhostType & ghost_type = _not_ghost) {
     AKANTU_TO_IMPLEMENT();
   }
 
   /// pre-compute normals on integration points
   virtual void computeNormalsOnIntegrationPoints(
       __attribute__((unused)) const Array<Real> & field,
       __attribute__((unused)) Array<Real> & normal,
       __attribute__((unused)) const ElementType & type,
       __attribute__((unused)) const GhostType & ghost_type = _not_ghost) const {
     AKANTU_TO_IMPLEMENT();
   }
 
-  // #ifdef AKANTU_STRUCTURAL_MECHANICS
-
-  //   /// assemble a field as a matrix for structural elements (ex. rho to mass
-  //   /// matrix)
-  //   virtual void assembleFieldMatrix(
-  //       __attribute__((unused)) const Array<Real> & field_1,
-  //       __attribute__((unused)) UInt nb_degree_of_freedom,
-  //       __attribute__((unused)) SparseMatrix & M,
-  //       __attribute__((unused)) Array<Real> * n,
-  //       __attribute__((unused)) ElementTypeMapArray<Real> & rotation_mat,
-  //       __attribute__((unused)) ElementType type,
-  //       __attribute__((unused)) const GhostType & ghost_type) const {
-
-  //     AKANTU_TO_IMPLEMENT();
-  //   }
-  //   /// compute shapes function in a matrix for structural elements
-  //   virtual void computeShapesMatrix(
-  //       __attribute__((unused)) const ElementType & type,
-  //       __attribute__((unused)) UInt nb_degree_of_freedom,
-  //       __attribute__((unused)) UInt nb_nodes_per_element,
-  //       __attribute__((unused)) Array<Real> * n, __attribute__((unused)) UInt
-  //       id,
-  //       __attribute__((unused)) UInt degree_to_interpolate,
-  //       __attribute__((unused)) UInt degree_interpolated,
-  //       __attribute__((unused)) const bool sign,
-  //       __attribute__((unused)) const GhostType & ghost_type) const {
-
-  //     AKANTU_TO_IMPLEMENT();
-  //   }
-
-  // #endif
-
   /// function to print the containt of the class
   virtual void printself(std::ostream & stream, int indent = 0) const;
 
 private:
   /// initialise the class
   void init();
 
   /* ------------------------------------------------------------------------ */
   /* Accessors                                                                */
   /* ------------------------------------------------------------------------ */
 public:
   using ElementTypesIteratorHelper =
       ElementTypeMapArray<Real, ElementType>::ElementTypesIteratorHelper;
 
   ElementTypesIteratorHelper elementTypes(UInt dim = _all_dimensions,
                                           GhostType ghost_type = _not_ghost,
                                           ElementKind kind = _ek_regular) const;
 
   /// get the dimension of the element handeled by this fe_engine object
   AKANTU_GET_MACRO(ElementDimension, element_dimension, UInt);
 
   /// get the mesh contained in the fem object
   AKANTU_GET_MACRO(Mesh, mesh, const Mesh &);
   /// get the mesh contained in the fem object
   AKANTU_GET_MACRO_NOT_CONST(Mesh, mesh, Mesh &);
 
   /// get the in-radius of an element
   static inline Real getElementInradius(const Matrix<Real> & coord,
                                         const ElementType & type);
 
   /// get the normals on integration points
   AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(NormalsOnIntegrationPoints,
                                          normals_on_integration_points, Real);
 
   /// get cohesive element type for a given facet type
   static inline ElementType
   getCohesiveElementType(const ElementType & type_facet);
 
   /// get igfem element type for a given regular type
   static inline Vector<ElementType>
   getIGFEMElementTypes(const ElementType & type);
 
   /// get the interpolation element associated to an element type
   static inline InterpolationType
   getInterpolationType(const ElementType & el_type);
 
   /// get the shape function class (probably useless: see getShapeFunction in
   /// fe_engine_template.hh)
   virtual const ShapeFunctions & getShapeFunctionsInterface() const = 0;
   /// get the integrator class (probably useless: see getIntegrator in
   /// fe_engine_template.hh)
   virtual const Integrator & getIntegratorInterface() const = 0;
 
   /* ------------------------------------------------------------------------ */
   /* Class Members                                                            */
   /* ------------------------------------------------------------------------ */
 protected:
   /// spatial dimension of the problem
   UInt element_dimension;
 
   /// the mesh on which all computation are made
   Mesh & mesh;
 
   /// normals at integration points
   ElementTypeMapArray<Real> normals_on_integration_points;
 };
 
 /* -------------------------------------------------------------------------- */
 /* inline functions                                                           */
 /* -------------------------------------------------------------------------- */
 
 /// standard output stream operator
 inline std::ostream & operator<<(std::ostream & stream,
                                  const FEEngine & _this) {
   _this.printself(stream);
   return stream;
 }
 
 } // namespace akantu
 
 #include "fe_engine_inline_impl.cc"
 #include "fe_engine_template.hh"
 
 #endif /* __AKANTU_FE_ENGINE_HH__ */