diff --git a/examples/python/custom-material/bi-material.py b/examples/python/custom-material/bi-material.py
index 37bf307b4..4d03347a4 100644
--- a/examples/python/custom-material/bi-material.py
+++ b/examples/python/custom-material/bi-material.py
@@ -1,184 +1,184 @@
 from __future__ import print_function
 # ------------------------------------------------------------- #
 import akantu as aka
 import subprocess
 import numpy as np
 import time
 # ------------------------------------------------------------- #
 
 
 class LocalElastic:
 
     # declares all the internals
     def initMaterial(self, internals, params):
         self.E = params['E']
         self.nu = params['nu']
         self.rho = params['rho']
         # print(self.__dict__)
         # First Lame coefficient
         self.lame_lambda = self.nu * self.E / (
             (1. + self.nu) * (1. - 2. * self.nu))
         # Second Lame coefficient (shear modulus)
         self.lame_mu = self.E / (2. * (1. + self.nu))
 
         all_factor = internals['factor']
         all_quad_coords = internals['quad_coordinates']
 
         for elem_type in all_factor.keys():
             factor = all_factor[elem_type]
             quad_coords = all_quad_coords[elem_type]
 
             factor[:] = 1.
             factor[quad_coords[:, 1] < 0.5] = .5
 
     # declares all the internals
     @staticmethod
     def registerInternals():
         return ['potential', 'factor']
 
     # declares all the internals
     @staticmethod
     def registerInternalSizes():
         return [1, 1]
 
     # declares all the parameters that could be parsed
     @staticmethod
     def registerParam():
         return ['E', 'nu']
 
     # declares all the parameters that are needed
     def getPushWaveSpeed(self, params):
         return np.sqrt((self.lame_lambda + 2 * self.lame_mu) / self.rho)
 
     # compute small deformation tensor
     @staticmethod
     def computeEpsilon(grad_u):
         return 0.5 * (grad_u + np.einsum('aij->aji', grad_u))
 
     # constitutive law
     def computeStress(self, grad_u, sigma, internals, params):
         n_quads = grad_u.shape[0]
         grad_u = grad_u.reshape((n_quads, 2, 2))
         factor = internals['factor'].reshape(n_quads)
         epsilon = self.computeEpsilon(grad_u)
         sigma = sigma.reshape((n_quads, 2, 2))
         trace = np.einsum('aii->a', grad_u)
 
         sigma[:, :, :] = (
             np.einsum('a,ij->aij', trace,
                       self.lame_lambda * np.eye(2))
             + 2. * self.lame_mu * epsilon)
 
         # print(sigma.reshape((n_quads, 4)))
         # print(grad_u.reshape((n_quads, 4)))
         sigma[:, :, :] = np.einsum('aij, a->aij', sigma, factor)
 
     # constitutive law tangent modulii
     def computeTangentModuli(self, grad_u, tangent, internals, params):
         n_quads = tangent.shape[0]
         tangent = tangent.reshape(n_quads, 3, 3)
         factor = internals['factor'].reshape(n_quads)
 
         Miiii = self.lame_lambda + 2 * self.lame_mu
         Miijj = self.lame_lambda
         Mijij = self.lame_mu
 
         tangent[:, 0, 0] = Miiii
         tangent[:, 1, 1] = Miiii
         tangent[:, 0, 1] = Miijj
         tangent[:, 1, 0] = Miijj
         tangent[:, 2, 2] = Mijij
         tangent[:, :, :] = np.einsum('aij, a->aij', tangent, factor)
 
     # computes the energy density
     def getEnergyDensity(self, energy_type, energy_density,
                          grad_u, stress, internals, params):
 
         nquads = stress.shape[0]
         stress = stress.reshape(nquads, 2, 2)
         grad_u = grad_u.reshape((nquads, 2, 2))
 
         if energy_type != 'potential':
             raise RuntimeError('not known energy')
 
         epsilon = self.computeEpsilon(grad_u)
 
         energy_density[:, 0] = (
             0.5 * np.einsum('aij,aij->a', stress, epsilon))
 
 
 # applies manually the boundary conditions
 def applyBC(model):
 
     nbNodes = model.getMesh().getNbNodes()
     position = model.getMesh().getNodes()
     displacement = model.getDisplacement()
     blocked_dofs = model.getBlockedDOFs()
 
     width = 1.
     height = 1.
     epsilon = 1e-8
     for node in range(0, nbNodes):
         if((np.abs(position[node, 0]) < epsilon) or  # left side
            (np.abs(position[node, 0] - width) < epsilon)):  # right side
             blocked_dofs[node, 0] = True
             displacement[node, 0] = 0 * position[node, 0] + 0.
 
         if(np.abs(position[node, 1]) < epsilon):  # lower side
             blocked_dofs[node, 1] = True
             displacement[node, 1] = - 1.
         if(np.abs(position[node, 1] - height) < epsilon):  # upper side
             blocked_dofs[node, 1] = True
             displacement[node, 1] = 1.
 
 # main parameters
 spatial_dimension = 2
 mesh_file = 'square.msh'
 
 # call gmsh to generate the mesh
-ret = subprocess.call(['gmsh', '-2', 'square.geo', '-optimize', 'square.msh'])
+ret = subprocess.call('gmsh -format msh2 -2 square.geo -optimize square.msh', shell=True)
 if ret != 0:
     raise Exception(
         'execution of GMSH failed: do you have it installed ?')
 
 time.sleep(1)
 
 # read mesh
 mesh = aka.Mesh(spatial_dimension)
 mesh.read(mesh_file)
 
 # create the custom material
 mat = LocalElastic()
 aka.registerNewPythonMaterial(mat, "local_elastic")
 
 # parse input file
 aka.parseInput('material.dat')
 
 # init the SolidMechanicsModel
 model = aka.SolidMechanicsModel(mesh)
 model.initFull(_analysis_method=aka._static)
 
 # configure the solver
 solver = model.getNonLinearSolver()
 solver.set("max_iterations", 2)
 solver.set("threshold", 1e-3)
 solver.set("convergence_type", aka._scc_solution)
 
 # prepare the dumper
 model.setBaseName("bimaterial")
 model.addDumpFieldVector("displacement")
 model.addDumpFieldVector("internal_force")
 model.addDumpFieldVector("external_force")
 model.addDumpField("strain")
 model.addDumpField("stress")
 model.addDumpField("factor")
 model.addDumpField("blocked_dofs")
 
 # Boundary conditions
 applyBC(model)
 
 # solve the problem
 model.solveStep()
 
 # dump paraview files
 model.dump()
diff --git a/examples/python/dynamics/dynamics.py b/examples/python/dynamics/dynamics.py
index c8be57fec..5b7e3ef84 100644
--- a/examples/python/dynamics/dynamics.py
+++ b/examples/python/dynamics/dynamics.py
@@ -1,129 +1,131 @@
 #!/usr/bin/env python3
 
 from __future__ import print_function
 ################################################################
 import os
 import subprocess
 import numpy as np
 import akantu
 ################################################################
 
 
-class FixedValue:
+class MyFixedValue(akantu.FixedValue):
 
     def __init__(self, value, axis):
+        super().__init__(value, axis)
         self.value = value
-        self.axis = axis
+        self.axis = int(axis)
 
-    def operator(self, node, flags, disp, coord):
+    def __call__(self, node, flags, disp, coord):
         # sets the displacement to the desired value in the desired axis
         disp[self.axis] = self.value
         # sets the blocked dofs vector to true in the desired axis
         flags[self.axis] = True
 
 ################################################################
 
 
 def main():
 
     spatial_dimension = 2
 
     akantu.parseInput('material.dat')
 
     mesh_file = 'bar.msh'
     max_steps = 250
     time_step = 1e-3
 
     # if mesh was not created the calls gmsh to generate it
     if not os.path.isfile(mesh_file):
-        ret = subprocess.call('gmsh -2 bar.geo bar.msh', shell=True)
+        ret = subprocess.call('gmsh -format msh2 -2 bar.geo bar.msh', shell=True)
         if ret != 0:
             raise Exception(
                 'execution of GMSH failed: do you have it installed ?')
 
     ################################################################
     # Initialization
     ################################################################
     mesh = akantu.Mesh(spatial_dimension)
     mesh.read(mesh_file)
 
     model = akantu.SolidMechanicsModel(mesh)
 
     model.initFull(_analysis_method=akantu._explicit_lumped_mass)
     # model.initFull(_analysis_method=akantu._implicit_dynamic)
 
     model.setBaseName("waves")
     model.addDumpFieldVector("displacement")
     model.addDumpFieldVector("acceleration")
     model.addDumpFieldVector("velocity")
     model.addDumpFieldVector("internal_force")
     model.addDumpFieldVector("external_force")
     model.addDumpField("strain")
     model.addDumpField("stress")
     model.addDumpField("blocked_dofs")
 
     ################################################################
     # boundary conditions
     ################################################################
 
-    model.applyDirichletBC(FixedValue(0, akantu._x), "XBlocked")
-    model.applyDirichletBC(FixedValue(0, akantu._y), "YBlocked")
+    model.applyBC(MyFixedValue(0, akantu._x), "XBlocked")
+    model.applyBC(MyFixedValue(0, akantu._y), "YBlocked")
 
     ################################################################
     # initial conditions
     ################################################################
 
     displacement = model.getDisplacement()
     nb_nodes = mesh.getNbNodes()
     position = mesh.getNodes()
 
     pulse_width = 1
     A = 0.01
     for i in range(0, nb_nodes):
         # Sinus * Gaussian
         x = position[i, 0] - 5.
         L = pulse_width
         k = 0.1 * 2 * np.pi * 3 / L
         displacement[i, 0] = A * \
             np.sin(k * x) * np.exp(-(k * x) * (k * x) / (L * L))
+        displacement[i, 1] = 0
 
     ################################################################
     # timestep value computation
     ################################################################
     time_factor = 0.8
     stable_time_step = model.getStableTimeStep() * time_factor
 
     print("Stable Time Step = {0}".format(stable_time_step))
     print("Required Time Step = {0}".format(time_step))
 
     time_step = stable_time_step * time_factor
 
     model.setTimeStep(time_step)
 
     ################################################################
     # loop for evolution of motion dynamics
     ################################################################
     model.assembleInternalForces()
 
     print("step,step * time_step,epot,ekin,epot + ekin")
     for step in range(0, max_steps + 1):
 
         model.solveStep()
 
         if step % 10 == 0:
             model.dump()
 
         epot = model.getEnergy('potential')
         ekin = model.getEnergy('kinetic')
 
         # output energy calculation to screen
         print("{0},{1},{2},{3},{4}".format(step, step * time_step,
                                            epot, ekin,
                                            (epot + ekin)))
 
     return
 
 
 ################################################################
 if __name__ == "__main__":
     main()
diff --git a/python/pybind11/py_aka_array.cc b/python/pybind11/py_aka_array.cc
index 3c5ee7a13..a5b17996e 100644
--- a/python/pybind11/py_aka_array.cc
+++ b/python/pybind11/py_aka_array.cc
@@ -1,215 +1,203 @@
 /* -------------------------------------------------------------------------- */
 #include "aka_array.hh"
 /* -------------------------------------------------------------------------- */
 #include <pybind11/numpy.h>
 #include <pybind11/pybind11.h>
 /* -------------------------------------------------------------------------- */
 
 namespace py = pybind11;
 namespace _aka = akantu;
 
 namespace akantu {
+
 template <typename T> class ArrayProxy : public Array<T> {
 protected:
   // deallocate the memory
   void deallocate() override final {}
 
   // allocate the memory
   void allocate(__attribute__((unused)) UInt size,
                 __attribute__((unused)) UInt nb_component) override final {}
 
   // allocate and initialize the memory
   void allocate(__attribute__((unused)) UInt size,
                 __attribute__((unused)) UInt nb_component,
                 __attribute__((unused)) const T & value) override final {}
 
 public:
   ArrayProxy(T * data, UInt size, UInt nb_component) {
     this->values = data;
     this->size_ = size;
     this->nb_component = nb_component;
   }
 
   ArrayProxy(const Array<T> & src) {
     this->values = src.storage();
     this->size_ = src.size();
     this->nb_component = src.getNbComponent();
   }
 
+  ~ArrayProxy() { this->values = nullptr; }
+
   void resize(__attribute__((unused)) UInt size,
               __attribute__((unused)) const T & val) override final {
     AKANTU_EXCEPTION("cannot resize a temporary array");
   }
 
   void resize(__attribute__((unused)) UInt new_size) override final {
     AKANTU_EXCEPTION("cannot resize a temporary array");
   }
 
   void reserve(__attribute__((unused)) UInt new_size) override final {
     AKANTU_EXCEPTION("cannot resize a temporary array");
   }
 };
 
 } // namespace akantu
 
 namespace pybind11 {
 namespace detail {
+
+  template <typename U>
+  using array_type = array_t<U, array::c_style | array::forcecast>;
+
   /* ------------------------------------------------------------------------ */
   template <typename T>
   py::handle aka_array_cast(const _aka::Array<T> & src,
                             py::handle base = handle(), bool writeable = true) {
-    constexpr ssize_t elem_size = sizeof(T);
     array a;
-    a = array({src.size(), src.getNbComponent()},
-              {elem_size * src.getNbComponent(), elem_size}, src.storage(),
-              base);
+    a = array_type<T>({src.size(), src.getNbComponent()}, src.storage(), base);
 
     if (not writeable)
       array_proxy(a.ptr())->flags &= ~detail::npy_api::NPY_ARRAY_WRITEABLE_;
 
     return a.release();
   }
 
   template <typename T> struct type_caster<_aka::Array<T>> {
   protected:
-    using py_array = array_t<T, array::forcecast | array::c_style>;
     using type = _aka::Array<T>;
+    type value;
 
   public:
-    PYBIND11_TYPE_CASTER(type, _("Array<T>"));
-
-    // static constexpr auto name =
-    //     _("numpy.ndarray[") + npy_format_descriptor<T>::name +
-    //     _(", flags.writeable") + _(", flags.c_contiguous") + _("]");
+    static PYBIND11_DESCR name() { return type_descr(_("Array<T>")); };
 
     /**
      * Conversion part 1 (Python->C++)
      */
     bool load(handle src, bool convert) {
-      bool need_copy = not isinstance<py_array>(src);
+      bool need_copy = not isinstance<array_type<T>>(src);
 
-      auto && fits = [&](auto && a) {
-        auto && dims = copy_or_ref.ndim();
+      auto && fits = [&](auto && aref) {
+        auto && dims = aref.ndim();
         if (dims < 1 || dims > 2)
           return false;
 
         return true;
       };
 
       if (not need_copy) {
         // We don't need a converting copy, but we also need to check whether
         // the strides are compatible with the Ref's stride requirements
-        py_array aref = reinterpret_borrow<py_array>(src);
+        auto aref = py::cast<array_type<T>>(src);
 
         if (not fits(aref)) {
           return false;
         }
-
         copy_or_ref = std::move(aref);
       } else {
         if (not convert) {
           return false;
         }
 
-        auto copy = py_array::ensure(src);
+        auto copy = array_type<T>::ensure(src);
         if (not copy) {
           return false;
         }
 
         if (not fits(copy)) {
           return false;
         }
-        copy_or_ref = std::move(py_array::ensure(src));
+        copy_or_ref = std::move(array_type<T>::ensure(src));
         loader_life_support::add_patient(copy_or_ref);
       }
 
-      array_proxy.reset(new _aka::ArrayProxy<T>> (copy_or_ref.mutable_data(),
-                                                  copy_or_ref.shape(0),
-                                                  copy_or_ref.shape(1)));
+      array_proxy = std::make_unique<_aka::ArrayProxy<T>>(
+          copy_or_ref.mutable_data(), copy_or_ref.shape(0),
+          copy_or_ref.shape(1));
       return true;
     }
 
-    // operator _aka::Array<T> *() { return array_proxy.get(); }
-    // operator _aka::Array<T> &() { return *array_proxy; }
-    // template <typename _T>
-    // using cast_op_type = pybind11::detail::cast_op_type<_T>;
+    operator _aka::Array<T> *() { return array_proxy.get(); }
+    operator _aka::Array<T> &() { return *array_proxy; }
+    template <typename _T>
+    using cast_op_type = pybind11::detail::cast_op_type<_T>;
 
     /**
      * Conversion part 2 (C++ -> Python)
      */
     static handle cast(const _aka::Array<T> & src, return_value_policy policy,
                        handle parent) {
       switch (policy) {
       case return_value_policy::copy:
         return aka_array_cast<T>(src);
       case return_value_policy::reference_internal:
         return aka_array_cast<T>(src, parent);
       case return_value_policy::reference:
       case return_value_policy::automatic:
       case return_value_policy::automatic_reference:
         return aka_array_cast<T>(src, none());
       default:
         pybind11_fail("Invalid return_value_policy for ArrayProxy type");
       }
     }
 
   protected:
-    std::unique_ptr<_aka::Array<T>> array_proxy;
-    py_array copy_or_ref;
+    std::unique_ptr<_aka::ArrayProxy<T>> array_proxy;
+    array_type<T> copy_or_ref;
   };
 
   /**
    * Type caster for Vector classes
    */
   template <template <typename> class V, typename T> struct type_caster<V<T>> {
     using type = V<T>;
-    template <typename U>
-    using array_type = array_t<U, array::c_style | array::forcecast>;
 
   public:
     PYBIND11_TYPE_CASTER(type, _("Vector<T>"));
 
     /**
      * Conversion part 1 (Python->C++): convert a PyObject into a Vector
      * instance or return false upon failure. The second argument
      * indicates whether implicit conversions should be applied.
      */
     bool load(handle src, bool convert) {
       if (!convert && !array_type<typename type::value_type>::check_(src))
         return false;
 
       auto buf = array_type<typename type::value_type>::ensure(src);
-      value.move(_aka::VectorProxy<T>(buf));
+      value = std::move(_aka::VectorProxy<T>(buf));
 
       return true;
     }
 
     /**
      * Conversion part 2 (C++ -> Python): convert a grid instance into
      * a Python object. The second and third arguments are used to
      * indicate the return value policy and parent object (for
      * ``return_value_policy::reference_internal``) and are generally
      * ignored by implicit casters.
      *
      * TODO: do not ignore policy (see pybind11/eigen.h)
      */
     static handle cast(const type & src, return_value_policy /* policy */,
                        handle /*parent*/) {
       // none() passed as parent to get a correct nocopy
       auto a = array_type<typename type::value_type>(src.size(), src.storage(),
                                                      none());
       return a.release();
     }
   };
 
 } // namespace detail
 } // namespace pybind11
-
-namespace {
-py::module & register_arrays(py::module & mod) {
-
-  // py::class_<_aka::Vector<bool>>(mod, "VectorBool");
-  // py::class_<_aka::Vector<double>>(mod, "VectorReal");
-  return mod;
-}
-} // namespace
diff --git a/python/pybind11/py_aka_boundary_conditions.cc b/python/pybind11/py_aka_boundary_conditions.cc
index c40670a23..fda58a839 100644
--- a/python/pybind11/py_aka_boundary_conditions.cc
+++ b/python/pybind11/py_aka_boundary_conditions.cc
@@ -1,66 +1,102 @@
 /* -------------------------------------------------------------------------- */
 #include "aka_common.hh"
 /* -------------------------------------------------------------------------- */
 #include <pybind11/operators.h>
 #include <pybind11/pybind11.h>
 #include <pybind11/stl.h>
 /* -------------------------------------------------------------------------- */
 #include "boundary_condition_functor.hh"
 /* -------------------------------------------------------------------------- */
 namespace py = pybind11;
 namespace _aka = akantu;
 
 namespace {
 
 /* -------------------------------------------------------------------------- */
 
-class PyDirichletFunctor : public _aka::BC::DirichletFunctor {
+template <typename daughter = _aka::BC::DirichletFunctor>
+class PyDirichletFunctor : public daughter {
 public:
   /* Inherit the constructors */
-  using _aka::BC::DirichletFunctor::DirichletFunctor;
+  using daughter::daughter;
 
   /* Trampoline (need one for each virtual function) */
   void operator()(_aka::UInt node, _aka::Vector<bool> & flags,
                   _aka::Vector<_aka::Real> & primal,
                   const _aka::Vector<_aka::Real> & coord) const override {
 
-    PYBIND11_OVERLOAD_NAME(void, _aka::BC::DirichletFunctor,
-                           "__call__", operator(), node, flags, primal, coord);
+    PYBIND11_OVERLOAD_NAME(void, daughter, "__call__", operator(), node, flags,
+                           primal, coord);
   }
 };
 /* -------------------------------------------------------------------------- */
 
-class PyNeumannFunctor : public _aka::BC::NeumannFunctor {
+template <typename daughter = _aka::BC::NeumannFunctor>
+class PyNeumannFunctor : public daughter {
 public:
   /* Inherit the constructors */
-  using _aka::BC::NeumannFunctor::NeumannFunctor;
+  using daughter::daughter;
 
   /* Trampoline (need one for each virtual function) */
   void operator()(const _aka::IntegrationPoint & quad_point,
                   _aka::Vector<_aka::Real> & dual,
                   const _aka::Vector<_aka::Real> & coord,
                   const _aka::Vector<_aka::Real> & normals) const override {
 
-    PYBIND11_OVERLOAD_PURE_NAME(void, _aka::BC::NeumannFunctor,
-                           "__call__", operator(), quad_point, dual, coord,
-                           normals);
+    PYBIND11_OVERLOAD_PURE_NAME(void, daughter, "__call__", operator(),
+                                quad_point, dual, coord, normals);
   }
 };
+
+/* -------------------------------------------------------------------------- */
+
+template <typename Functor, typename Constructor>
+decltype(auto) declareDirichletFunctor(py::module mod, const char * name,
+                                       Constructor && cons) {
+  py::class_<Functor, PyDirichletFunctor<Functor>, _aka::BC::DirichletFunctor>(
+      mod, name)
+      .def(cons);
+}
+template <typename Functor, typename Constructor>
+decltype(auto) declareNeumannFunctor(py::module mod, const char * name,
+                                     Constructor && cons) {
+  py::class_<Functor, PyNeumannFunctor<Functor>, _aka::BC::NeumannFunctor>(mod,
+                                                                           name)
+      .def(cons);
+}
+
 /* -------------------------------------------------------------------------- */
 
 py::module & register_boundary_conditions(py::module & mod) {
 
-  py::class_<_aka::BC::Functor>(mod, "BCFunctor");
-  py::class_<_aka::BC::DirichletFunctor, PyDirichletFunctor, _aka::BC::Functor>(
-      mod, "DirichletFunctor")
+  py::class_<_aka::BC::Functor>(mod, "_aka::BCFunctor");
+  py::class_<_aka::BC::DirichletFunctor, PyDirichletFunctor<>,
+             _aka::BC::Functor>(mod, "DirichletFunctor")
       .def(py::init())
       .def(py::init<_aka::SpatialDirection>());
 
-  py::class_<_aka::BC::NeumannFunctor, PyNeumannFunctor, _aka::BC::Functor>(
+  py::class_<_aka::BC::NeumannFunctor, PyNeumannFunctor<>, _aka::BC::Functor>(
       mod, "NeumannFunctor")
       .def(py::init());
 
-  return mod;
-} // namespace
+  declareDirichletFunctor<_aka::BC::FixedValue>(
+      mod, "FixedValue", py::init<_aka::Real, _aka::BC::Axis>());
 
+  declareDirichletFunctor<_aka::BC::IncrementValue>(
+      mod, "IncrementValue", py::init<_aka::Real, _aka::BC::Axis>());
+
+  // declareDirichletFunctor<_aka::BC::Increment>(
+  //     mod, "Increment", py::init<_aka::Vector<_aka::Real> &>());
+
+  // declareNeumannFunctor<_aka::BC::FromHigherDim>(
+  //     mod, "FromHigherDim", py::init<_aka::Matrix<_aka::Real> &>());
+
+  // declareNeumannFunctor<_aka::BC::FromSameDim>(
+  //     mod, "FromSameDim", py::init<_aka::Vector<_aka::Real> &>());
+
+  declareNeumannFunctor<_aka::BC::FreeBoundary>(mod, "FreeBoundary",
+                                                py::init());
+
+  return mod;
+}
 } // namespace
diff --git a/python/pybind11/py_aka_common.cc b/python/pybind11/py_aka_common.cc
index 0729b7505..a15c7463c 100644
--- a/python/pybind11/py_aka_common.cc
+++ b/python/pybind11/py_aka_common.cc
@@ -1,118 +1,121 @@
 /* -------------------------------------------------------------------------- */
 #include "aka_common.hh"
 /* -------------------------------------------------------------------------- */
 #include <pybind11/operators.h>
 #include <pybind11/pybind11.h>
 #include <pybind11/stl.h>
 /* -------------------------------------------------------------------------- */
 #include "integration_point.hh"
 #include "mesh.hh"
 #include "parser.hh"
 /* -------------------------------------------------------------------------- */
-namespace py = pybind11;
-namespace _aka = akantu;
 
 namespace {
+namespace py = pybind11;
+namespace _aka = akantu;
 
 /* -------------------------------------------------------------------------- */
 
 py::module & register_enums(py::module & mod) {
   py::enum_<_aka::SpatialDirection>(mod, "SpatialDirection")
       .value("_x", _aka::_x)
       .value("_y", _aka::_y)
       .value("_z", _aka::_z)
       .export_values();
 
   py::enum_<_aka::AnalysisMethod>(mod, "AnalysisMethod")
       .value("_static", _aka::_static)
       .value("_implicit_dynamic", _aka::_implicit_dynamic)
       .value("_explicit_lumped_mass", _aka::_explicit_lumped_mass)
       .value("_explicit_lumped_capacity", _aka::_explicit_lumped_capacity)
       .value("_explicit_consistent_mass", _aka::_explicit_consistent_mass)
       .export_values();
 
   py::enum_<_aka::NonLinearSolverType>(mod, "NonLinearSolverType")
       .value("_nls_linear", _aka::_nls_linear)
       .value("_nls_newton_raphson", _aka::_nls_newton_raphson)
       .value("_nls_newton_raphson_modified", _aka::_nls_newton_raphson_modified)
       .value("_nls_lumped", _aka::_nls_lumped)
       .value("_nls_auto", _aka::_nls_auto)
       .export_values();
 
   py::enum_<_aka::TimeStepSolverType>(mod, "TimeStepSolverType")
       .value("_tsst_static", _aka::_tsst_static)
       .value("_tsst_dynamic", _aka::_tsst_dynamic)
       .value("_tsst_dynamic_lumped", _aka::_tsst_dynamic_lumped)
       .value("_tsst_not_defined", _aka::_tsst_not_defined)
       .export_values();
 
   py::enum_<_aka::IntegrationSchemeType>(mod, "IntegrationSchemeType")
       .value("_ist_pseudo_time", _aka::_ist_pseudo_time)
       .value("_ist_forward_euler", _aka::_ist_forward_euler)
       .value("_ist_trapezoidal_rule_1", _aka::_ist_trapezoidal_rule_1)
       .value("_ist_backward_euler", _aka::_ist_backward_euler)
       .value("_ist_central_difference", _aka::_ist_central_difference)
       .value("_ist_fox_goodwin", _aka::_ist_fox_goodwin)
       .value("_ist_trapezoidal_rule_2", _aka::_ist_trapezoidal_rule_2)
       .value("_ist_linear_acceleration", _aka::_ist_linear_acceleration)
       .value("_ist_newmark_beta", _aka::_ist_newmark_beta)
       .value("_ist_generalized_trapezoidal", _aka::_ist_generalized_trapezoidal)
       .export_values();
 
   py::enum_<_aka::SolveConvergenceCriteria>(mod, "SolveConvergenceCriteria")
       .value("_scc_residual", _aka::_scc_residual)
       .value("_scc_solution", _aka::_scc_solution)
       .value("_scc_residual_mass_wgh", _aka::_scc_residual_mass_wgh)
       .export_values();
 
   py::enum_<_aka::CohesiveMethod>(mod, "CohesiveMethod")
       .value("_intrinsic", _aka::_intrinsic)
       .value("_extrinsic", _aka::_extrinsic)
       .export_values();
 
   py::enum_<_aka::GhostType>(mod, "GhostType")
       .value("_not_ghost", _aka::_not_ghost)
       .value("_ghost", _aka::_ghost)
       .value("_casper", _aka::_casper)
       .export_values();
 
   py::enum_<_aka::MeshIOType>(mod, "MeshIOType")
       .value("_miot_auto", _aka::_miot_auto)
       .value("_miot_gmsh", _aka::_miot_gmsh)
       .value("_miot_gmsh_struct", _aka::_miot_gmsh_struct)
       .value("_miot_diana", _aka::_miot_diana)
       .value("_miot_abaqus", _aka::_miot_abaqus)
       .export_values();
 
   py::enum_<_aka::ModelType>(mod, "ModelType")
 
       .value("_model", _aka::ModelType::_model)
       .value("_solid_mechanics_model", _aka::ModelType::_solid_mechanics_model)
       .value("_solid_mechanics_model_cohesive",
              _aka::ModelType::_solid_mechanics_model_cohesive)
       .value("_heat_transfer_model", _aka::ModelType::_heat_transfer_model)
       .value("_structural_mechanics_model",
              _aka::ModelType::_structural_mechanics_model)
       .value("_embedded_model", _aka::ModelType::_embedded_model)
       .export_values();
 
   mod.def("parseInput",
           [](const std::string & input_file) {
             _aka::getStaticParser().parse(input_file);
           },
           "Parse an Akantu input file");
 
   py::class_<_aka::Mesh>(mod, "Mesh")
       .def(py::init<_aka::UInt, const _aka::ID &, const _aka::MemoryID &>(),
            py::arg("spatial_dimension"), py::arg("id") = "mesh",
            py::arg("memory_id") = 0)
       .def("read", &_aka::Mesh::read, py::arg("filename"),
-           py::arg("mesh_io_type") = _aka::_miot_auto,
-           "read the mesh from a file");
+           py::arg("mesh_io_type") = _aka::_miot_auto, "read the mesh from a file")
+      .def("getNodes",
+           [](_aka::Mesh & self) -> _aka::Array<_aka::Real> { return self.getNodes(); },
+           py::return_value_policy::reference)
+      .def("getNbNodes", &_aka::Mesh::getNbNodes);
 
   py::class_<_aka::IntegrationPoint>(mod, "IntegrationPoint");
 
   return mod;
 } // namespace
 
 } // namespace
diff --git a/python/pybind11/py_aka_solid_mechanics_model.cc b/python/pybind11/py_aka_solid_mechanics_model.cc
index 289ed4f6a..6da0b2e88 100644
--- a/python/pybind11/py_aka_solid_mechanics_model.cc
+++ b/python/pybind11/py_aka_solid_mechanics_model.cc
@@ -1,82 +1,122 @@
 /* -------------------------------------------------------------------------- */
 #include "aka_common.hh"
 /* -------------------------------------------------------------------------- */
 #include <pybind11/operators.h>
 #include <pybind11/pybind11.h>
 #include <pybind11/stl.h>
 /* -------------------------------------------------------------------------- */
 #include "non_linear_solver.hh"
 #include "solid_mechanics_model.hh"
 /* -------------------------------------------------------------------------- */
 namespace py = pybind11;
 namespace _aka = akantu;
 
 namespace {
 
 /* -------------------------------------------------------------------------- */
+#define def_deprecated(func_name, mesg)                                        \
+  def(func_name, [](py::args, py::kwargs) { AKANTU_ERROR(mesg); })
+
+#define def_function_nocopy(func_name)                                         \
+  def(#func_name,                                                              \
+      [](_aka::SolidMechanicsModel & self) -> decltype(auto) {                 \
+        return self.func_name();                                               \
+      },                                                                       \
+      py::return_value_policy::reference)
+
+#define def_function(func_name)                                                \
+  def(#func_name, [](_aka::SolidMechanicsModel & self) -> decltype(auto) {     \
+    return self.func_name();                                                   \
+  })
+/* -------------------------------------------------------------------------- */
 
 py::module & register_solid_mechanics_models(py::module & mod) {
 
-  py::class_<_aka::SolidMechanicsModelOptions>(mod,
-                                               "SolidMechanicsModelOptions")
+  py::class_<_aka::SolidMechanicsModelOptions>(
+      mod, "_aka::SolidMechanicsModelOptions")
       .def(py::init<_aka::AnalysisMethod>(),
            py::arg("analysis_method") = _aka::_explicit_lumped_mass);
 
   py::class_<_aka::Model>(mod, "Model")
       .def("setBaseName", &_aka::Model::setBaseName)
       .def("addDumpFieldVector", &_aka::Model::addDumpFieldVector)
       .def("addDumpField", &_aka::Model::addDumpField)
       .def("dump", &_aka::Model::dump);
 
   py::class_<_aka::NonLinearSolver>(mod, "NonLinearSolver")
       .def("set",
            [](_aka::NonLinearSolver & self, const std::string & id,
               const _aka::Real & val) {
              if (id == "max_iterations")
                self.set(id, int(val));
              else if (id == "convergence_type")
                self.set(id, akantu::SolveConvergenceCriteria(_aka::UInt(val)));
              else
                self.set(id, val);
            })
       .def("set", &_aka::NonLinearSolver::set<_aka::SolveConvergenceCriteria>);
 
   py::class_<_aka::ModelSolver>(mod, "ModelSolver")
       .def("getNonLinearSolver",
            (_aka::NonLinearSolver & (_aka::ModelSolver::*)(const _aka::ID &)) &
                _aka::ModelSolver::getNonLinearSolver,
            py::arg("solver_id") = "", py::return_value_policy::reference)
       .def("solveStep", &_aka::ModelSolver::solveStep,
            py::arg("solver_id") = "");
 
   py::class_<_aka::SolidMechanicsModel, _aka::Model, _aka::ModelSolver>(
-      mod, "SolidMechanicsModel")
+      mod, "_aka::SolidMechanicsModel")
       .def(py::init<_aka::Mesh &, _aka::UInt, const _aka::ID &,
                     const _aka::MemoryID &, const _aka::ModelType>(),
            py::arg("mesh"),
            py::arg("spatial_dimension") = _aka::_all_dimensions,
            py::arg("id") = "solid_mechanics_model", py::arg("memory_id") = 0,
            py::arg("model_type") = _aka::ModelType::_solid_mechanics_model)
       .def("initFull",
            [](_aka::SolidMechanicsModel & self,
               const _aka::SolidMechanicsModelOptions & options) {
              self.initFull(options);
            },
-           py::arg("options") = _aka::SolidMechanicsModelOptions())
+           py::arg("_analysis_method") = _aka::SolidMechanicsModelOptions())
+      .def("initFull",
+           [](_aka::SolidMechanicsModel & self,
+              const _aka::AnalysisMethod & _analysis_method) {
+             self.initFull(_aka::SolidMechanicsModelOptions(_analysis_method));
+           },
+           py::arg("_analysis_method"))
+      .def_deprecated("applyDirichletBC", "Deprecated: use applyBC")
       .def("applyBC",
            [](_aka::SolidMechanicsModel & self,
               _aka::BC::DirichletFunctor & func,
               const std::string & element_group) {
              self.applyBC(func, element_group);
            })
       .def("applyBC",
            [](_aka::SolidMechanicsModel & self, _aka::BC::NeumannFunctor & func,
               const std::string & element_group) {
              self.applyBC(func, element_group);
            })
-      .def("getExternalForce", &_aka::SolidMechanicsModel::getExternalForce);
+      .def("setTimeStep", &_aka::SolidMechanicsModel::setTimeStep,
+           py::arg("time_step"), py::arg("solver_id") = "")
+      .def("getEnergy",
+           py::overload_cast<const std::string &>(
+               &_aka::SolidMechanicsModel::getEnergy),
+           py::arg("energy_id"))
+      .def_function(assembleStiffnessMatrix)
+      .def_function(assembleInternalForces)
+      .def_function(getStableTimeStep)
+      .def_function_nocopy(getExternalForce)
+      .def_function_nocopy(getDisplacement)
+      .def_function_nocopy(getPreviousDisplacement)
+      .def_function_nocopy(getIncrement)
+      .def_function_nocopy(getMass)
+      .def_function_nocopy(getVelocity)
+      .def_function_nocopy(getAcceleration)
+      .def_function_nocopy(getInternalForce)
+      .def_function_nocopy(getBlockedDOFs)
+      .def_function_nocopy(getIncrementFlag);
 
   return mod;
 } // namespace
 
 } // namespace
diff --git a/python/pybind11/py_akantu.cc b/python/pybind11/py_akantu.cc
index 188a60d3a..9b4ba99c5 100644
--- a/python/pybind11/py_akantu.cc
+++ b/python/pybind11/py_akantu.cc
@@ -1,19 +1,18 @@
 #include <pybind11/pybind11.h>
 #include <pybind11/numpy.h>
 
 namespace py = pybind11;
 
 #include "py_aka_common.cc"
 #include "py_aka_array.cc"
 #include "py_aka_boundary_conditions.cc"
 #include "py_aka_solid_mechanics_model.cc"
 
 PYBIND11_MODULE(py11_akantu, mod) {
   mod.doc() = "Akantu python interface";
 
   register_enums(mod);
-  register_arrays(mod);
   register_boundary_conditions(mod);
   register_solid_mechanics_models(mod);
 
 } // Module akantu
diff --git a/src/io/parser/parameter_registry_tmpl.hh b/src/io/parser/parameter_registry_tmpl.hh
index 21ca58755..b183f509d 100644
--- a/src/io/parser/parameter_registry_tmpl.hh
+++ b/src/io/parser/parameter_registry_tmpl.hh
@@ -1,412 +1,410 @@
 /**
  * @file   parameter_registry_tmpl.hh
  *
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Wed May 04 2016
  * @date last modification: Tue Jan 30 2018
  *
  * @brief  implementation of the templated part of ParameterRegistry class and
  * the derivated ones
  *
  * @section LICENSE
  *
  * Copyright (©) 2016-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_error.hh"
 #include "aka_iterators.hh"
 #include "parameter_registry.hh"
 #include "parser.hh"
 /* -------------------------------------------------------------------------- */
 #include <algorithm>
 #include <string>
 #include <vector>
 /* -------------------------------------------------------------------------- */
 
 #ifndef __AKANTU_PARAMETER_REGISTRY_TMPL_HH__
 #define __AKANTU_PARAMETER_REGISTRY_TMPL_HH__
 
 namespace akantu {
 
 namespace debug {
   class ParameterException : public Exception {
   public:
     ParameterException(const std::string & name, const std::string & message)
         : Exception(message), name(name) {}
     const std::string & name;
   };
 
   class ParameterUnexistingException : public ParameterException {
   public:
     ParameterUnexistingException(const std::string & name,
                                  const ParameterRegistry & registery)
-        : ParameterException(
-              name, "Parameter " + name + " does not exists in this scope") {
+        : ParameterException(name, "Parameter " + name +
+                                       " does not exists in this scope") {
       auto && params = registery.listParameters();
       this->_info =
           std::accumulate(params.begin(), params.end(),
                           this->_info + "\n Possible parameters are: ",
                           [](auto && str, auto && param) {
                             static auto first = true;
                             auto && ret = str + (first ? " " : ", ") + param;
                             first = false;
                             return ret;
                           });
     }
   };
 
   class ParameterAccessRightException : public ParameterException {
   public:
     ParameterAccessRightException(const std::string & name,
                                   const std::string & perm)
         : ParameterException(name, "Parameter " + name + " is not " + perm) {}
   };
 
   class ParameterWrongTypeException : public ParameterException {
   public:
     ParameterWrongTypeException(const std::string & name,
                                 const std::type_info & wrong_type,
                                 const std::type_info & type)
-        : ParameterException(name,
-                             "Parameter " + name +
-                                 " type error, cannot convert " +
-                                 debug::demangle(type.name()) + " to " +
-                                 debug::demangle(wrong_type.name())) {}
+        : ParameterException(name, "Parameter " + name +
+                                       " type error, cannot convert " +
+                                       debug::demangle(type.name()) + " to " +
+                                       debug::demangle(wrong_type.name())) {}
   };
 } // namespace debug
 /* -------------------------------------------------------------------------- */
 template <typename T>
 const ParameterTyped<T> & Parameter::getParameterTyped() const {
   try {
     const auto & tmp = dynamic_cast<const ParameterTyped<T> &>(*this);
     return tmp;
   } catch (std::bad_cast &) {
     AKANTU_CUSTOM_EXCEPTION(
         debug::ParameterWrongTypeException(name, typeid(T), this->type()));
   }
 }
 
 /* -------------------------------------------------------------------------- */
 template <typename T> ParameterTyped<T> & Parameter::getParameterTyped() {
   try {
     auto & tmp = dynamic_cast<ParameterTyped<T> &>(*this);
     return tmp;
   } catch (std::bad_cast &) {
     AKANTU_CUSTOM_EXCEPTION(
         debug::ParameterWrongTypeException(name, typeid(T), this->type()));
   }
 }
 
 /* ------------------------------------------------------------------------ */
 template <typename T, typename V> void Parameter::set(const V & value) {
   if (!(isWritable()))
     AKANTU_CUSTOM_EXCEPTION(
         debug::ParameterAccessRightException(name, "writable"));
   ParameterTyped<T> & typed_param = getParameterTyped<T>();
   typed_param.setTyped(value);
 }
 
 /* ------------------------------------------------------------------------ */
 inline void Parameter::setAuto(__attribute__((unused))
                                const ParserParameter & value) {
   if (!(isParsable()))
     AKANTU_CUSTOM_EXCEPTION(
         debug::ParameterAccessRightException(name, "parsable"));
 }
 
 /* -------------------------------------------------------------------------- */
 template <typename T> const T & Parameter::get() const {
   if (!(isReadable()))
     AKANTU_CUSTOM_EXCEPTION(
         debug::ParameterAccessRightException(name, "readable"));
   const ParameterTyped<T> & typed_param = getParameterTyped<T>();
   return typed_param.getTyped();
 }
 
 /* -------------------------------------------------------------------------- */
 template <typename T> T & Parameter::get() {
   ParameterTyped<T> & typed_param = getParameterTyped<T>();
   if (!(isReadable()) || !(this->isWritable()))
     AKANTU_CUSTOM_EXCEPTION(
         debug::ParameterAccessRightException(name, "accessible"));
   return typed_param.getTyped();
 }
 
 /* -------------------------------------------------------------------------- */
 template <typename T> inline Parameter::operator T() const {
   return this->get<T>();
 }
 
 /* -------------------------------------------------------------------------- */
 template <typename T>
 ParameterTyped<T>::ParameterTyped(std::string name, std::string description,
                                   ParameterAccessType param_type, T & param)
     : Parameter(name, description, param_type), param(param) {}
 
 /* -------------------------------------------------------------------------- */
 template <typename T>
 template <typename V>
 void ParameterTyped<T>::setTyped(const V & value) {
   param = value;
 }
 
 /* -------------------------------------------------------------------------- */
 template <typename T>
 inline void ParameterTyped<T>::setAuto(const ParserParameter & value) {
   Parameter::setAuto(value);
   param = static_cast<T>(value);
 }
 
 /* -------------------------------------------------------------------------- */
 template <>
 inline void
 ParameterTyped<std::string>::setAuto(const ParserParameter & value) {
   Parameter::setAuto(value);
   param = value.getValue();
 }
 
 /* -------------------------------------------------------------------------- */
 template <>
 inline void
 ParameterTyped<Vector<Real>>::setAuto(const ParserParameter & in_param) {
   Parameter::setAuto(in_param);
   Vector<Real> tmp = in_param;
   if (param.size() == 0) {
     param = tmp;
   } else {
     for (UInt i = 0; i < param.size(); ++i) {
       param(i) = tmp(i);
     }
   }
 }
 
 /* -------------------------------------------------------------------------- */
 template <>
 inline void
 ParameterTyped<Matrix<Real>>::setAuto(const ParserParameter & in_param) {
   Parameter::setAuto(in_param);
   Matrix<Real> tmp = in_param;
   if (param.size() == 0) {
     param = tmp;
   } else {
     for (UInt i = 0; i < param.rows(); ++i) {
       for (UInt j = 0; j < param.cols(); ++j) {
         param(i, j) = tmp(i, j);
       }
     }
   }
 }
 
 /* -------------------------------------------------------------------------- */
 template <typename T> const T & ParameterTyped<T>::getTyped() const {
   return param;
 }
 /* -------------------------------------------------------------------------- */
 template <typename T> T & ParameterTyped<T>::getTyped() { return param; }
 
 /* -------------------------------------------------------------------------- */
 template <typename T>
 inline void ParameterTyped<T>::printself(std::ostream & stream) const {
   Parameter::printself(stream);
   stream << param << "\n";
 }
 
 /* -------------------------------------------------------------------------- */
 template <typename T> class ParameterTyped<std::vector<T>> : public Parameter {
 public:
   ParameterTyped(std::string name, std::string description,
                  ParameterAccessType param_type, std::vector<T> & param)
       : Parameter(name, description, param_type), param(param) {}
 
   /* ------------------------------------------------------------------------ */
   template <typename V> void setTyped(const V & value) { param = value; }
   void setAuto(const ParserParameter & value) override {
     Parameter::setAuto(value);
     param.clear();
     const std::vector<T> & tmp = value;
     for (auto && z : tmp) {
       param.emplace_back(z);
     }
   }
 
   std::vector<T> & getTyped() { return param; }
   const std::vector<T> & getTyped() const { return param; }
 
   void printself(std::ostream & stream) const override {
     Parameter::printself(stream);
     stream << "[ ";
     for (auto && v : param)
       stream << v << " ";
     stream << "]\n";
   }
 
   inline const std::type_info & type() const override {
     return typeid(std::vector<T>);
   }
 
 private:
   /// Value of parameter
   std::vector<T> & param;
 };
 
 /* ------o--------------------------------------------------------------------
  */
 template <>
 inline void ParameterTyped<bool>::printself(std::ostream & stream) const {
   Parameter::printself(stream);
   stream << std::boolalpha << param << "\n";
 }
 
 /* -------------------------------------------------------------------------- */
 template <typename T>
 void ParameterRegistry::registerParam(std::string name, T & variable,
                                       ParameterAccessType type,
                                       const std::string & description) {
   auto it = params.find(name);
   if (it != params.end())
     AKANTU_CUSTOM_EXCEPTION(debug::ParameterException(
         name, "Parameter named " + name + " already registered."));
   auto * param = new ParameterTyped<T>(name, description, type, variable);
   params[name] = param;
 }
 
 /* -------------------------------------------------------------------------- */
 template <typename T>
 void ParameterRegistry::registerParam(std::string name, T & variable,
                                       const T & default_value,
                                       ParameterAccessType type,
                                       const std::string & description) {
   variable = default_value;
   registerParam(name, variable, type, description);
 }
 
 /* -------------------------------------------------------------------------- */
 template <typename T, typename V>
 void ParameterRegistry::setMixed(const std::string & name, const V & value) {
   auto it = params.find(name);
   if (it == params.end()) {
     if (consisder_sub) {
       for (auto it = sub_registries.begin(); it != sub_registries.end(); ++it) {
         it->second->setMixed<T>(name, value);
       }
     } else {
       AKANTU_CUSTOM_EXCEPTION(debug::ParameterUnexistingException(name, *this));
     }
   } else {
     Parameter & param = *(it->second);
     param.set<T>(value);
   }
 }
 
 /* -------------------------------------------------------------------------- */
 template <typename T>
 void ParameterRegistry::set(const std::string & name, const T & value) {
   this->template setMixed<T>(name, value);
 }
 
 /* -------------------------------------------------------------------------- */
 template <typename T> T & ParameterRegistry::get_(const std::string & name) {
   auto it = params.find(name);
   if (it == params.end()) {
     if (consisder_sub) {
       for (auto it = sub_registries.begin(); it != sub_registries.end(); ++it) {
         try {
           return it->second->get_<T>(name);
         } catch (...) {
         }
       }
     }
 
     // nothing was found not even in sub registries
     AKANTU_CUSTOM_EXCEPTION(debug::ParameterUnexistingException(name, *this));
   }
 
   Parameter & param = *(it->second);
   return param.get<T>();
 }
 
 /* -------------------------------------------------------------------------- */
 const Parameter & ParameterRegistry::get(const std::string & name) const {
   auto it = params.find(name);
   if (it == params.end()) {
     if (consisder_sub) {
       for (auto it = sub_registries.begin(); it != sub_registries.end(); ++it) {
         try {
           return it->second->get(name);
         } catch (...) {
         }
       }
     }
 
     // nothing was found not even in sub registries
     AKANTU_CUSTOM_EXCEPTION(debug::ParameterUnexistingException(name, *this));
   }
 
   Parameter & param = *(it->second);
   return param;
 }
 
 /* -------------------------------------------------------------------------- */
 Parameter & ParameterRegistry::get(const std::string & name) {
   auto it = params.find(name);
   if (it == params.end()) {
     if (consisder_sub) {
       for (auto it = sub_registries.begin(); it != sub_registries.end(); ++it) {
         try {
           return it->second->get(name);
         } catch (...) {
         }
       }
     }
 
     // nothing was found not even in sub registries
     AKANTU_CUSTOM_EXCEPTION(debug::ParameterUnexistingException(name, *this));
   }
 
   Parameter & param = *(it->second);
   return param;
 }
 
-
 /* -------------------------------------------------------------------------- */
 namespace {
   namespace details {
     template <class T, class R, class Enable = void> struct CastHelper {
       static R convert(const T &) { throw std::bad_cast(); }
     };
 
     template <class T, class R>
     struct CastHelper<T, R,
                       std::enable_if_t<std::is_convertible<T, R>::value>> {
       static R convert(const T & val) { return val; }
     };
   } // namespace details
 } // namespace
 
 template <typename T> inline ParameterTyped<T>::operator Real() const {
   if (not isReadable())
     AKANTU_CUSTOM_EXCEPTION(
         debug::ParameterAccessRightException(name, "accessible"));
   return details::CastHelper<T, Real>::convert(param);
 }
 
 } // namespace akantu
 
 #endif /* __AKANTU_PARAMETER_REGISTRY_TMPL_HH__ */
diff --git a/src/model/boundary_condition_functor.hh b/src/model/boundary_condition_functor.hh
index c7a9727d6..7175c7c34 100644
--- a/src/model/boundary_condition_functor.hh
+++ b/src/model/boundary_condition_functor.hh
@@ -1,207 +1,206 @@
 /**
  * @file   boundary_condition_functor.hh
  *
  * @author Dana Christen <dana.christen@gmail.com>
  * @author David Simon Kammer <david.kammer@epfl.ch>
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Fri May 03 2013
  * @date last modification: Tue Feb 20 2018
  *
  * @brief  Definitions of the functors to apply boundary conditions
  *
  * @section LICENSE
  *
  * Copyright (©) 2014-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_common.hh"
 #include "fe_engine.hh"
 #include "integration_point.hh"
 /* -------------------------------------------------------------------------- */
 
 #ifndef __AKANTU_BOUNDARY_CONDITION_FUNCTOR_HH__
 #define __AKANTU_BOUNDARY_CONDITION_FUNCTOR_HH__
 
 /* -------------------------------------------------------------------------- */
 
 namespace akantu {
 
 /* -------------------------------------------------------------------------- */
 namespace BC {
   using Axis = ::akantu::SpatialDirection;
 
   /* ---------------------------------------------------------------------- */
   struct Functor {
     enum Type { _dirichlet, _neumann };
   };
 
   /* ---------------------------------------------------------------------- */
   class DirichletFunctor : public Functor {
   public:
-
     DirichletFunctor() = default;
     explicit DirichletFunctor(Axis ax) : axis(ax) {}
 
     virtual void operator()(__attribute__((unused)) UInt node,
                             __attribute__((unused)) Vector<bool> & flags,
                             __attribute__((unused)) Vector<Real> & primal,
                             __attribute__((unused))
                             const Vector<Real> & coord) const {
       AKANTU_TO_IMPLEMENT();
     }
 
   public:
     static const Type type = _dirichlet;
 
   protected:
     Axis axis{_x};
   };
 
   /* ---------------------------------------------------------------------- */
   class FlagOnly : public DirichletFunctor {
   public:
     explicit FlagOnly(Axis ax = _x) : DirichletFunctor(ax) {}
 
   public:
     inline void operator()(UInt node, Vector<bool> & flags,
                            Vector<Real> & primal,
                            const Vector<Real> & coord) const;
   };
 
   /* ---------------------------------------------------------------------- */
-  class FreeBoundary : public DirichletFunctor {
-  public:
-    explicit FreeBoundary(Axis ax = _x) : DirichletFunctor(ax) {}
+  // class FreeBoundary : public DirichletFunctor {
+  // public:
+  //   explicit FreeBoundary(Axis ax = _x) : DirichletFunctor(ax) {}
 
-  public:
-    inline void operator()(UInt node, Vector<bool> & flags,
-                           Vector<Real> & primal,
-                           const Vector<Real> & coord) const;
-  };
+  // public:
+  //   inline void operator()(UInt node, Vector<bool> & flags,
+  //                          Vector<Real> & primal,
+  //                          const Vector<Real> & coord) const;
+  // };
 
   /* ---------------------------------------------------------------------- */
   class FixedValue : public DirichletFunctor {
   public:
     FixedValue(Real val, Axis ax = _x) : DirichletFunctor(ax), value(val) {}
 
   public:
     inline void operator()(UInt node, Vector<bool> & flags,
                            Vector<Real> & primal,
                            const Vector<Real> & coord) const;
 
   protected:
     Real value;
   };
 
   /* ---------------------------------------------------------------------- */
   class IncrementValue : public DirichletFunctor {
   public:
     IncrementValue(Real val, Axis ax = _x) : DirichletFunctor(ax), value(val) {}
 
   public:
     inline void operator()(UInt node, Vector<bool> & flags,
                            Vector<Real> & primal,
                            const Vector<Real> & coord) const;
 
     inline void setIncrement(Real val) { this->value = val; }
 
   protected:
     Real value;
   };
 
   /* ---------------------------------------------------------------------- */
   class Increment : public DirichletFunctor {
   public:
     explicit Increment(const Vector<Real> & val)
         : DirichletFunctor(_x), value(val) {}
 
   public:
     inline void operator()(UInt node, Vector<bool> & flags,
                            Vector<Real> & primal,
                            const Vector<Real> & coord) const;
 
     inline void setIncrement(const Vector<Real> & val) { this->value = val; }
 
   protected:
     Vector<Real> value;
   };
 
   /* ------------------------------------------------------------------------ */
   /* Neumann                                                                  */
   /* ------------------------------------------------------------------------ */
   class NeumannFunctor : public Functor {
 
   protected:
     NeumannFunctor() = default;
 
   public:
     virtual void operator()(const IntegrationPoint & quad_point,
                             Vector<Real> & dual, const Vector<Real> & coord,
                             const Vector<Real> & normals) const = 0;
 
     virtual ~NeumannFunctor() = default;
 
   public:
     static const Type type = _neumann;
   };
 
   /* ---------------------------------------------------------------------- */
   class FromHigherDim : public NeumannFunctor {
   public:
     explicit FromHigherDim(const Matrix<Real> & mat) : bc_data(mat) {}
     ~FromHigherDim() override = default;
 
   public:
     inline void operator()(const IntegrationPoint & quad_point,
                            Vector<Real> & dual, const Vector<Real> & coord,
                            const Vector<Real> & normals) const override;
 
   protected:
     Matrix<Real> bc_data;
   };
 
   /* ---------------------------------------------------------------------- */
   class FromSameDim : public NeumannFunctor {
   public:
     explicit FromSameDim(const Vector<Real> & vec) : bc_data(vec) {}
     ~FromSameDim() override = default;
 
   public:
     inline void operator()(const IntegrationPoint & quad_point,
                            Vector<Real> & dual, const Vector<Real> & coord,
                            const Vector<Real> & normals) const override;
 
   protected:
     Vector<Real> bc_data;
   };
 
   /* ---------------------------------------------------------------------- */
-  // class FreeBoundary : public NeumannFunctor {
-  // public:
-  //   inline void operator()(const IntegrationPoint & quad_point,
-  //                          Vector<Real> & dual, const Vector<Real> & coord,
-  //                          const Vector<Real> & normals) const override;
-  // };
+  class FreeBoundary : public NeumannFunctor {
+  public:
+    inline void operator()(const IntegrationPoint & quad_point,
+                           Vector<Real> & dual, const Vector<Real> & coord,
+                           const Vector<Real> & normals) const override;
+  };
 } // namespace BC
 } // namespace akantu
 
 #include "boundary_condition_functor_inline_impl.cc"
 
 #endif /* __AKANTU_BOUNDARY_CONDITION_FUNCTOR_HH__ */
diff --git a/src/model/boundary_condition_functor_inline_impl.cc b/src/model/boundary_condition_functor_inline_impl.cc
index 0d960be68..8cc42b5c0 100644
--- a/src/model/boundary_condition_functor_inline_impl.cc
+++ b/src/model/boundary_condition_functor_inline_impl.cc
@@ -1,150 +1,150 @@
 /**
  * @file   boundary_condition_functor_inline_impl.cc
  *
  * @author Dana Christen <dana.christen@gmail.com>
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Fri May 03 2013
  * @date last modification: Mon Feb 19 2018
  *
  * @brief  implementation of the BC::Functors
  *
  * @section LICENSE
  *
  * Copyright (©) 2014-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 "boundary_condition_functor.hh"
 /* -------------------------------------------------------------------------- */
 
 #ifndef __AKANTU_BOUNDARY_CONDITION_FUNCTOR_INLINE_IMPL_CC__
 #define __AKANTU_BOUNDARY_CONDITION_FUNCTOR_INLINE_IMPL_CC__
 
 /* -------------------------------------------------------------------------- */
 #define DIRICHLET_SANITY_CHECK                                                 \
   AKANTU_DEBUG_ASSERT(                                                         \
       coord.size() <= flags.size(),                                            \
       "The coordinates and flags vectors given to the boundary"                \
           << " condition functor have different sizes!");                      \
   AKANTU_DEBUG_ASSERT(                                                         \
       primal.size() <= coord.size(),                                           \
       "The primal vector and coordinates vector given"                         \
           << " to the boundary condition functor have different sizes!");
 
 #define NEUMANN_SANITY_CHECK                                                   \
   AKANTU_DEBUG_ASSERT(                                                         \
       coord.size() <= normals.size(),                                          \
       "The coordinates and normals vectors given to the"                       \
           << " boundary condition functor have different sizes!");             \
   AKANTU_DEBUG_ASSERT(                                                         \
       dual.size() <= coord.size(),                                             \
       "The dual vector and coordinates vector given to"                        \
           << " the boundary condition functor have different sizes!");
 
 namespace akantu {
 namespace BC {
   /* ---------------------------------------------------------------------- */
   inline void FlagOnly::
   operator()(__attribute__((unused)) UInt node, Vector<bool> & flags,
              __attribute__((unused)) Vector<Real> & primal,
              __attribute__((unused)) const Vector<Real> & coord) const {
 
     DIRICHLET_SANITY_CHECK;
 
     flags(this->axis) = true;
   }
 
   /* ---------------------------------------------------------------------- */
-  inline void FreeBoundary::
-  operator()(__attribute__((unused)) UInt node, Vector<bool> & flags,
-             __attribute__((unused)) Vector<Real> & primal,
-             __attribute__((unused)) const Vector<Real> & coord) const {
+  // inline void FreeBoundary::
+  // operator()(__attribute__((unused)) UInt node, Vector<bool> & flags,
+  //            __attribute__((unused)) Vector<Real> & primal,
+  //            __attribute__((unused)) const Vector<Real> & coord) const {
 
-    DIRICHLET_SANITY_CHECK;
+  //   DIRICHLET_SANITY_CHECK;
 
-    flags(this->axis) = false;
-  }
+  //   flags(this->axis) = false;
+  // }
 
   /* ---------------------------------------------------------------------- */
   inline void FixedValue::operator()(__attribute__((unused)) UInt node,
                                      Vector<bool> & flags,
                                      Vector<Real> & primal,
                                      __attribute__((unused))
                                      const Vector<Real> & coord) const {
     DIRICHLET_SANITY_CHECK;
     flags(this->axis) = true;
     primal(this->axis) = value;
   }
 
   /* ---------------------------------------------------------------------- */
   inline void IncrementValue::operator()(__attribute__((unused)) UInt node,
                                          Vector<bool> & flags,
                                          Vector<Real> & primal,
                                          __attribute__((unused))
                                          const Vector<Real> & coord) const {
     DIRICHLET_SANITY_CHECK;
     flags(this->axis) = true;
     primal(this->axis) += value;
   }
 
   /* ---------------------------------------------------------------------- */
   inline void Increment::operator()(__attribute__((unused)) UInt node,
                                     Vector<bool> & flags, Vector<Real> & primal,
                                     __attribute__((unused))
                                     const Vector<Real> & coord) const {
     DIRICHLET_SANITY_CHECK;
     flags.set(true);
     primal += value;
   }
 
   /* ------------------------------------------------------------------------ */
   /* Neumann */
   /* ------------------------------------------------------------------------ */
-  // inline void FreeBoundary::
-  // operator()(__attribute__((unused)) const IntegrationPoint & quad_point,
-  //            Vector<Real> & dual,
-  //            __attribute__((unused)) const Vector<Real> & coord,
-  //            __attribute__((unused)) const Vector<Real> & normals) const {
-  //   for (UInt i(0); i < dual.size(); ++i) {
-  //     dual(i) = 0.0;
-  //   }
-  // }
+  inline void FreeBoundary::
+  operator()(__attribute__((unused)) const IntegrationPoint & quad_point,
+             Vector<Real> & dual,
+             __attribute__((unused)) const Vector<Real> & coord,
+             __attribute__((unused)) const Vector<Real> & normals) const {
+    for (UInt i(0); i < dual.size(); ++i) {
+      dual(i) = 0.0;
+    }
+  }
 
   /* ---------------------------------------------------------------------- */
   inline void FromHigherDim::operator()(__attribute__((unused))
                                         const IntegrationPoint & quad_point,
                                         Vector<Real> & dual,
                                         __attribute__((unused))
                                         const Vector<Real> & coord,
                                         const Vector<Real> & normals) const {
     dual.mul<false>(this->bc_data, normals);
   }
 
   /* ---------------------------------------------------------------------- */
   inline void FromSameDim::
   operator()(__attribute__((unused)) const IntegrationPoint & quad_point,
              Vector<Real> & dual,
              __attribute__((unused)) const Vector<Real> & coord,
              __attribute__((unused)) const Vector<Real> & normals) const {
     dual = this->bc_data;
   }
 } // namespace BC
 } // namespace akantu
 
 #endif /* __AKANTU_BOUNDARY_CONDITION_FUNCTOR_INLINE_IMPL_CC__ */
diff --git a/src/model/solid_mechanics/solid_mechanics_model.hh b/src/model/solid_mechanics/solid_mechanics_model.hh
index 0edab13c7..8d53cb1b7 100644
--- a/src/model/solid_mechanics/solid_mechanics_model.hh
+++ b/src/model/solid_mechanics/solid_mechanics_model.hh
@@ -1,568 +1,566 @@
 /**
  * @file   solid_mechanics_model.hh
  *
  * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
  * @author Daniel Pino Muñoz <daniel.pinomunoz@epfl.ch>
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Tue Jul 27 2010
  * @date last modification: Wed Feb 21 2018
  *
  * @brief  Model of Solid Mechanics
  *
  * @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 "boundary_condition.hh"
 #include "data_accessor.hh"
 #include "fe_engine.hh"
 #include "model.hh"
 #include "non_local_manager.hh"
 #include "solid_mechanics_model_event_handler.hh"
 /* -------------------------------------------------------------------------- */
 
 #ifndef __AKANTU_SOLID_MECHANICS_MODEL_HH__
 #define __AKANTU_SOLID_MECHANICS_MODEL_HH__
 
 namespace akantu {
 class Material;
 class MaterialSelector;
 class DumperIOHelper;
 class NonLocalManager;
 template <ElementKind kind, class IntegrationOrderFunctor>
 class IntegratorGauss;
 template <ElementKind kind> class ShapeLagrange;
 } // namespace akantu
 
 /* -------------------------------------------------------------------------- */
 namespace akantu {
 
 /* -------------------------------------------------------------------------- */
 class SolidMechanicsModel
     : public Model,
       public DataAccessor<Element>,
       public DataAccessor<UInt>,
       public BoundaryCondition<SolidMechanicsModel>,
       public NonLocalManagerCallback,
       public EventHandlerManager<SolidMechanicsModelEventHandler> {
 
   /* ------------------------------------------------------------------------ */
   /* Constructors/Destructors                                                 */
   /* ------------------------------------------------------------------------ */
 public:
   class NewMaterialElementsEvent : public NewElementsEvent {
   public:
     AKANTU_GET_MACRO_NOT_CONST(MaterialList, material, Array<UInt> &);
     AKANTU_GET_MACRO(MaterialList, material, const Array<UInt> &);
 
   protected:
     Array<UInt> material;
   };
 
   using MyFEEngineType = FEEngineTemplate<IntegratorGauss, ShapeLagrange>;
 
 protected:
   using EventManager = EventHandlerManager<SolidMechanicsModelEventHandler>;
 
 public:
   SolidMechanicsModel(
       Mesh & mesh, UInt spatial_dimension = _all_dimensions,
       const ID & id = "solid_mechanics_model", const MemoryID & memory_id = 0,
       const ModelType model_type = ModelType::_solid_mechanics_model);
 
   ~SolidMechanicsModel() override;
 
   /* ------------------------------------------------------------------------ */
   /* Methods                                                                  */
   /* ------------------------------------------------------------------------ */
 protected:
   /// initialize completely the model
   void initFullImpl(
       const ModelOptions & options = SolidMechanicsModelOptions()) override;
 
   /// initialize all internal arrays for materials
   virtual void initMaterials();
 
   /// initialize the model
   void initModel() override;
 
   /// function to print the containt of the class
   void printself(std::ostream & stream, int indent = 0) const override;
 
   /// get some default values for derived classes
   std::tuple<ID, TimeStepSolverType>
   getDefaultSolverID(const AnalysisMethod & method) override;
 
   /* ------------------------------------------------------------------------ */
   /* Solver interface                                                         */
   /* ------------------------------------------------------------------------ */
 public:
   /// assembles the stiffness matrix,
   virtual void assembleStiffnessMatrix();
   /// assembles the internal forces in the array internal_forces
   virtual void assembleInternalForces();
 
 protected:
   /// callback for the solver, this adds f_{ext} - f_{int} to the residual
   void assembleResidual() override;
 
   /// callback for the solver, this adds f_{ext} or  f_{int} to the residual
   void assembleResidual(const ID & residual_part) override;
   bool canSplitResidual() override { return true; }
 
   /// get the type of matrix needed
   MatrixType getMatrixType(const ID & matrix_id) override;
 
   /// callback for the solver, this assembles different matrices
   void assembleMatrix(const ID & matrix_id) override;
 
   /// callback for the solver, this assembles the stiffness matrix
   void assembleLumpedMatrix(const ID & matrix_id) override;
 
   /// callback for the solver, this is called at beginning of solve
   void predictor() override;
   /// callback for the solver, this is called at end of solve
   void corrector() override;
 
   /// callback for the solver, this is called at beginning of solve
   void beforeSolveStep() override;
   /// callback for the solver, this is called at end of solve
   void afterSolveStep() override;
 
   /// Callback for the model to instantiate the matricees when needed
   void initSolver(TimeStepSolverType time_step_solver_type,
                   NonLinearSolverType non_linear_solver_type) override;
 
 protected:
   /* ------------------------------------------------------------------------ */
   TimeStepSolverType getDefaultSolverType() const override;
   /* ------------------------------------------------------------------------ */
   ModelSolverOptions
   getDefaultSolverOptions(const TimeStepSolverType & type) const override;
 
 public:
   bool isDefaultSolverExplicit() {
     return method == _explicit_lumped_mass ||
            method == _explicit_consistent_mass;
   }
 
 protected:
   /// update the current position vector
   void updateCurrentPosition();
 
   /* ------------------------------------------------------------------------ */
   /* Materials (solid_mechanics_model_material.cc)                            */
   /* ------------------------------------------------------------------------ */
 public:
   /// register an empty material of a given type
   Material & registerNewMaterial(const ID & mat_name, const ID & mat_type,
                                  const ID & opt_param);
 
   /// reassigns materials depending on the material selector
   virtual void reassignMaterial();
 
   /// apply a constant eigen_grad_u on all quadrature points of a given material
   virtual void applyEigenGradU(const Matrix<Real> & prescribed_eigen_grad_u,
                                const ID & material_name,
                                const GhostType ghost_type = _not_ghost);
 
 protected:
   /// register a material in the dynamic database
   Material & registerNewMaterial(const ParserSection & mat_section);
 
   /// read the material files to instantiate all the materials
   void instantiateMaterials();
 
   /// set the element_id_by_material and add the elements to the good materials
   virtual void
   assignMaterialToElements(const ElementTypeMapArray<UInt> * filter = nullptr);
 
   /* ------------------------------------------------------------------------ */
   /* Mass (solid_mechanics_model_mass.cc)                                     */
   /* ------------------------------------------------------------------------ */
 public:
   /// assemble the lumped mass matrix
   void assembleMassLumped();
 
   /// assemble the mass matrix for consistent mass resolutions
   void assembleMass();
 
 protected:
   /// assemble the lumped mass matrix for local and ghost elements
   void assembleMassLumped(GhostType ghost_type);
 
   /// assemble the mass matrix for either _ghost or _not_ghost elements
   void assembleMass(GhostType ghost_type);
 
   /// fill a vector of rho
   void computeRho(Array<Real> & rho, ElementType type, GhostType ghost_type);
 
   /// compute the kinetic energy
   Real getKineticEnergy();
   Real getKineticEnergy(const ElementType & type, UInt index);
 
   /// compute the external work (for impose displacement, the velocity should be
   /// given too)
   Real getExternalWork();
 
   /* ------------------------------------------------------------------------ */
   /* NonLocalManager inherited members                                        */
   /* ------------------------------------------------------------------------ */
 protected:
   void initializeNonLocal() override;
 
   void updateDataForNonLocalCriterion(ElementTypeMapReal & criterion) override;
 
   void computeNonLocalStresses(const GhostType & ghost_type) override;
 
   void
   insertIntegrationPointsInNeighborhoods(const GhostType & ghost_type) override;
 
   /// update the values of the non local internal
   void updateLocalInternal(ElementTypeMapReal & internal_flat,
                            const GhostType & ghost_type,
                            const ElementKind & kind) override;
 
   /// copy the results of the averaging in the materials
   void updateNonLocalInternal(ElementTypeMapReal & internal_flat,
                               const GhostType & ghost_type,
                               const ElementKind & kind) override;
 
   /* ------------------------------------------------------------------------ */
   /* Data Accessor inherited members                                          */
   /* ------------------------------------------------------------------------ */
 public:
   UInt getNbData(const Array<Element> & elements,
                  const SynchronizationTag & tag) const override;
 
   void packData(CommunicationBuffer & buffer, const Array<Element> & elements,
                 const SynchronizationTag & tag) const override;
 
   void unpackData(CommunicationBuffer & buffer, const Array<Element> & elements,
                   const SynchronizationTag & tag) override;
 
   UInt getNbData(const Array<UInt> & dofs,
                  const SynchronizationTag & tag) const override;
 
   void packData(CommunicationBuffer & buffer, const Array<UInt> & dofs,
                 const SynchronizationTag & tag) const override;
 
   void unpackData(CommunicationBuffer & buffer, const Array<UInt> & dofs,
                   const SynchronizationTag & tag) override;
 
 protected:
   void
   splitElementByMaterial(const Array<Element> & elements,
                          std::vector<Array<Element>> & elements_per_mat) const;
 
   template <typename Operation>
   void splitByMaterial(const Array<Element> & elements, Operation && op) const;
 
   /* ------------------------------------------------------------------------ */
   /* Mesh Event Handler inherited members                                     */
   /* ------------------------------------------------------------------------ */
 protected:
   void onNodesAdded(const Array<UInt> & nodes_list,
                     const NewNodesEvent & event) override;
   void onNodesRemoved(const Array<UInt> & element_list,
                       const Array<UInt> & new_numbering,
                       const RemovedNodesEvent & event) override;
   void onElementsAdded(const Array<Element> & nodes_list,
                        const NewElementsEvent & event) override;
   void onElementsRemoved(const Array<Element> & element_list,
                          const ElementTypeMapArray<UInt> & new_numbering,
                          const RemovedElementsEvent & event) override;
   void onElementsChanged(const Array<Element> &, const Array<Element> &,
                          const ElementTypeMapArray<UInt> &,
                          const ChangedElementsEvent &) override{};
 
   /* ------------------------------------------------------------------------ */
   /* Dumpable interface (kept for convenience) and dumper relative functions  */
   /* ------------------------------------------------------------------------ */
 public:
   virtual void onDump();
 
   //! decide wether a field is a material internal or not
   bool isInternal(const std::string & field_name,
                   const ElementKind & element_kind);
 #ifndef SWIG
   //! give the amount of data per element
   virtual ElementTypeMap<UInt>
   getInternalDataPerElem(const std::string & field_name,
                          const ElementKind & kind);
 
   //! flatten a given material internal field
   ElementTypeMapArray<Real> &
   flattenInternal(const std::string & field_name, const ElementKind & kind,
                   const GhostType ghost_type = _not_ghost);
   //! flatten all the registered material internals
   void flattenAllRegisteredInternals(const ElementKind & kind);
 #endif
 
   dumper::Field * createNodalFieldReal(const std::string & field_name,
                                        const std::string & group_name,
                                        bool padding_flag) override;
 
   dumper::Field * createNodalFieldBool(const std::string & field_name,
                                        const std::string & group_name,
                                        bool padding_flag) override;
 
   dumper::Field * createElementalField(const std::string & field_name,
                                        const std::string & group_name,
                                        bool padding_flag,
                                        const UInt & spatial_dimension,
                                        const ElementKind & kind) override;
 
   virtual void dump(const std::string & dumper_name);
 
   virtual void dump(const std::string & dumper_name, UInt step);
 
   virtual void dump(const std::string & dumper_name, Real time, UInt step);
 
   void dump() override;
 
   virtual void dump(UInt step);
 
   virtual void dump(Real time, UInt step);
 
   /* ------------------------------------------------------------------------ */
   /* Accessors                                                                */
   /* ------------------------------------------------------------------------ */
 public:
   /// return the dimension of the system space
   AKANTU_GET_MACRO(SpatialDimension, Model::spatial_dimension, UInt);
 
   /// set the value of the time step
   void setTimeStep(Real time_step, const ID & solver_id = "") override;
 
   /// get the value of the conversion from forces/ mass to acceleration
   AKANTU_GET_MACRO(F_M2A, f_m2a, Real);
 
   /// set the value of the conversion from forces/ mass to acceleration
   AKANTU_SET_MACRO(F_M2A, f_m2a, Real);
 
   /// get the SolidMechanicsModel::displacement vector
   AKANTU_GET_MACRO(Displacement, *displacement, Array<Real> &);
 
   /// get the SolidMechanicsModel::previous_displacement vector
   AKANTU_GET_MACRO(PreviousDisplacement, *previous_displacement, Array<Real> &);
 
   /// get the SolidMechanicsModel::current_position vector \warn only consistent
   /// after a call to SolidMechanicsModel::updateCurrentPosition
   const Array<Real> & getCurrentPosition();
 
   /// get  the SolidMechanicsModel::increment  vector \warn  only  consistent if
   /// SolidMechanicsModel::setIncrementFlagOn has been called before
   AKANTU_GET_MACRO(Increment, *displacement_increment, Array<Real> &);
 
   /// get the lumped SolidMechanicsModel::mass vector
   AKANTU_GET_MACRO(Mass, *mass, Array<Real> &);
 
   /// get the SolidMechanicsModel::velocity vector
   AKANTU_GET_MACRO(Velocity, *velocity, Array<Real> &);
 
   /// get    the    SolidMechanicsModel::acceleration    vector,   updated    by
   /// SolidMechanicsModel::updateAcceleration
   AKANTU_GET_MACRO(Acceleration, *acceleration, Array<Real> &);
 
   /// get the SolidMechanicsModel::external_force vector (external forces)
   AKANTU_GET_MACRO(ExternalForce, *external_force, Array<Real> &);
 
   /// get the SolidMechanicsModel::force vector (external forces)
   Array<Real> & getForce() {
     AKANTU_DEBUG_WARNING("getForce was maintained for backward compatibility, "
                          "use getExternalForce instead");
     return *external_force;
   }
 
   /// get the SolidMechanicsModel::internal_force vector (internal forces)
   AKANTU_GET_MACRO(InternalForce, *internal_force, Array<Real> &);
 
   /// get the SolidMechanicsModel::blocked_dofs vector
   AKANTU_GET_MACRO(BlockedDOFs, *blocked_dofs, Array<bool> &);
 
   /// get the value of the SolidMechanicsModel::increment_flag
   AKANTU_GET_MACRO(IncrementFlag, increment_flag, bool);
 
   /// get a particular material (by material index)
   inline Material & getMaterial(UInt mat_index);
 
   /// get a particular material (by material index)
   inline const Material & getMaterial(UInt mat_index) const;
 
   /// get a particular material (by material name)
   inline Material & getMaterial(const std::string & name);
 
   /// get a particular material (by material name)
   inline const Material & getMaterial(const std::string & name) const;
 
   /// get a particular material id from is name
   inline UInt getMaterialIndex(const std::string & name) const;
 
   /// give the number of materials
   inline UInt getNbMaterials() const { return materials.size(); }
 
   /// give the material internal index from its id
   Int getInternalIndexFromID(const ID & id) const;
 
   /// compute the stable time step
   Real getStableTimeStep();
 
   /// get the energies
   Real getEnergy(const std::string & energy_id);
 
   /// compute the energy for energy
   Real getEnergy(const std::string & energy_id, const ElementType & type,
                  UInt index);
 
   AKANTU_GET_MACRO(MaterialByElement, material_index,
                    const ElementTypeMapArray<UInt> &);
   AKANTU_GET_MACRO(MaterialLocalNumbering, material_local_numbering,
                    const ElementTypeMapArray<UInt> &);
 
   /// vectors containing local material element index for each global element
   /// index
   AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(MaterialByElement, material_index,
                                          UInt);
   AKANTU_GET_MACRO_BY_ELEMENT_TYPE(MaterialByElement, material_index, UInt);
   AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(MaterialLocalNumbering,
                                          material_local_numbering, UInt);
   AKANTU_GET_MACRO_BY_ELEMENT_TYPE(MaterialLocalNumbering,
                                    material_local_numbering, UInt);
 
   AKANTU_GET_MACRO_NOT_CONST(MaterialSelector, *material_selector,
                              MaterialSelector &);
   AKANTU_SET_MACRO(MaterialSelector, material_selector,
                    std::shared_ptr<MaterialSelector>);
 
   /// Access the non_local_manager interface
   AKANTU_GET_MACRO(NonLocalManager, *non_local_manager, NonLocalManager &);
 
   /// get the FEEngine object to integrate or interpolate on the boundary
   FEEngine & getFEEngineBoundary(const ID & name = "") override;
 
 protected:
   friend class Material;
 
 protected:
   /// compute the stable time step
   Real getStableTimeStep(const GhostType & ghost_type);
 
   /* ------------------------------------------------------------------------ */
   /* Class Members                                                            */
   /* ------------------------------------------------------------------------ */
 protected:
   /// conversion coefficient form force/mass to acceleration
   Real f_m2a;
 
   /// displacements array
   Array<Real> * displacement;
   UInt displacement_release{0};
 
   /// displacements array at the previous time step (used in finite deformation)
   Array<Real> * previous_displacement{nullptr};
 
   /// increment of displacement
   Array<Real> * displacement_increment{nullptr};
 
   /// lumped mass array
   Array<Real> * mass{nullptr};
 
   /// Check if materials need to recompute the mass array
   bool need_to_reassemble_lumped_mass{true};
   /// Check if materials need to recompute the mass matrix
   bool need_to_reassemble_mass{true};
 
   /// velocities array
   Array<Real> * velocity{nullptr};
 
   /// accelerations array
   Array<Real> * acceleration{nullptr};
 
   /// accelerations array
   // Array<Real> * increment_acceleration;
 
   /// external forces array
   Array<Real> * external_force{nullptr};
 
   /// internal forces array
   Array<Real> * internal_force{nullptr};
 
   /// array specifing if a degree of freedom is blocked or not
   Array<bool> * blocked_dofs{nullptr};
 
   /// array of current position used during update residual
   Array<Real> * current_position{nullptr};
   UInt current_position_release{0};
 
   /// Arrays containing the material index for each element
   ElementTypeMapArray<UInt> material_index;
 
   /// Arrays containing the position in the element filter of the material
   /// (material's local numbering)
   ElementTypeMapArray<UInt> material_local_numbering;
 
   /// list of used materials
   std::vector<std::unique_ptr<Material>> materials;
 
   /// mapping between material name and material internal id
   std::map<std::string, UInt> materials_names_to_id;
 
   /// class defining of to choose a material
   std::shared_ptr<MaterialSelector> material_selector;
 
   /// flag defining if the increment must be computed or not
   bool increment_flag;
 
   /// tells if the material are instantiated
   bool are_materials_instantiated;
 
   using flatten_internal_map = std::map<std::pair<std::string, ElementKind>,
                                         ElementTypeMapArray<Real> *>;
 
   /// map a registered internals to be flattened for dump purposes
   flatten_internal_map registered_internals;
 
   /// non local manager
   std::unique_ptr<NonLocalManager> non_local_manager;
 };
 
 /* -------------------------------------------------------------------------- */
 namespace BC {
-  namespace Neumann {
-    using FromStress = FromHigherDim;
-    using FromTraction = FromSameDim;
-  } // namespace Neumann
+  using FromStress = FromHigherDim;
+  using FromTraction = FromSameDim;
 } // namespace BC
 
 } // namespace akantu
 
 /* -------------------------------------------------------------------------- */
 /* inline functions                                                           */
 /* -------------------------------------------------------------------------- */
 #include "material.hh"
 #include "parser.hh"
 
 #include "solid_mechanics_model_inline_impl.cc"
 #include "solid_mechanics_model_tmpl.hh"
 /* -------------------------------------------------------------------------- */
 
 #endif /* __AKANTU_SOLID_MECHANICS_MODEL_HH__ */
diff --git a/test/test_model/test_solid_mechanics_model/test_energies/test_solid_mechanics_model_work_dynamics.cc b/test/test_model/test_solid_mechanics_model/test_energies/test_solid_mechanics_model_work_dynamics.cc
index 13e33c967..8376713b6 100644
--- a/test/test_model/test_solid_mechanics_model/test_energies/test_solid_mechanics_model_work_dynamics.cc
+++ b/test/test_model/test_solid_mechanics_model/test_energies/test_solid_mechanics_model_work_dynamics.cc
@@ -1,157 +1,158 @@
 /**
  * @file   test_solid_mechanics_model_work_dynamics.cc
  *
  * @author Tobias Brink <tobias.brink@epfl.ch>
  *
  * @date creation: Fri Dec 15 2017
  * @date last modification: Fri Jan 26 2018
  *
  * @brief  test work in dynamic simulations
  *
  * @section LICENSE
  *
  * Copyright (©) 2016-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/>.
  *
  * @section description
  *
  * Assuming that the kinetic energy and the potential energy of the
  * linear elastic material are bug free, the work in a dynamic
  * simulation must equal the change in internal energy (first law of
  * thermodynamics). Work in dynamics is an infinitesimal work Fds,
  * thus we need to integrate it and compare at the end. In this test,
  * we use one Dirichlet boundary condition (with u = 0.0, 0.01, and
  * -0.01) and one Neumann boundary condition for F on the opposite
  * side. Then we do a few steps to get reference energies for work and
  * internal energy. After more steps, the change in both work and
  * internal energy must be equal.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 #include "../test_solid_mechanics_model_fixture.hh"
 /* -------------------------------------------------------------------------- */
 
 using namespace akantu;
 
 namespace {
 
 template <typename type_>
 class TestSMMFixtureWorkDynamic : public TestSMMFixture<type_> {
 public:
   void SetUp() override {
-    this->mesh_file = "../../patch_tests/data/bar" + aka::to_string(this->type) + ".msh";
+    this->mesh_file =
+        "../../patch_tests/data/bar" + aka::to_string(this->type) + ".msh";
     TestSMMFixture<type_>::SetUp();
 
     getStaticParser().parse("test_solid_mechanics_model_"
                             "work_material.dat");
 
     /// model initialization
     this->model->initFull();
 
     /// Create a node group for Neumann BCs.
     auto & apply_force_grp = this->mesh->createNodeGroup("apply_force");
     auto & fixed_grp = this->mesh->createNodeGroup("fixed");
 
     const auto & pos = this->mesh->getNodes();
     const auto & lower = this->mesh->getLowerBounds();
     const auto & upper = this->mesh->getUpperBounds();
 
     UInt i = 0;
     for (auto && posv : make_view(pos, this->spatial_dimension)) {
       if (posv(_x) > upper(_x) - 1e-6) {
         apply_force_grp.add(i);
       } else if (posv(_x) < lower(_x) + 1e-6) {
         fixed_grp.add(i);
       }
       ++i;
     }
 
     this->mesh->createElementGroupFromNodeGroup("el_apply_force", "apply_force",
                                                 this->spatial_dimension - 1);
     this->mesh->createElementGroupFromNodeGroup("el_fixed", "fixed",
                                                 this->spatial_dimension - 1);
 
     Vector<Real> surface_traction(this->spatial_dimension);
     surface_traction(_x) = 0.5;
     if (this->spatial_dimension == 1) {
       // TODO: this is a hack to work
       //      around non-implemented
       //      BC::Neumann::FromTraction for 1D
       auto & force = this->model->getForce();
 
       for (auto && pair : zip(make_view(pos, this->spatial_dimension),
                               make_view(force, this->spatial_dimension))) {
         auto & posv = std::get<0>(pair);
         auto & forcev = std::get<1>(pair);
         if (posv(_x) > upper(_x) - 1e-6) {
           forcev(_x) = surface_traction(_x);
         }
       }
     } else {
-      this->model->applyBC(BC::Neumann::FromTraction(surface_traction),
+      this->model->applyBC(BC::FromTraction(surface_traction),
                            "el_apply_force");
     }
 
     /// set up timestep
     auto time_step = this->model->getStableTimeStep() * 0.1;
     this->model->setTimeStep(time_step);
   }
 };
 
 TYPED_TEST_CASE(TestSMMFixtureWorkDynamic, gtest_element_types);
 
 /* TODO: this is currently disabled for terrible results and performance
 TYPED_TEST(TestSMMFixtureBar, WorkImplicit) {
   test_body(*(this->model), *(this->mesh), _implicit_dynamic, 500);
 }
 */
 // model.assembleMassLumped();
 TYPED_TEST(TestSMMFixtureWorkDynamic, WorkExplicit) {
   /// Do the sim
   std::vector<Real> displacements{0.00, 0.01, -0.01};
 
   for (auto && u : displacements) {
-    this->model->applyBC(BC::Dirichlet::FixedValue(u, _x), "el_fixed");
+    this->model->applyBC(BC::FixedValue(u, _x), "el_fixed");
 
     // First, "equilibrate" a bit to get a reference state of total
     // energy and work. This is needed when we have a Dirichlet with
     // finite displacement on one side.
     for (UInt i = 0; i < 25; ++i) {
       this->model->solveStep();
     }
     // Again, work reported by Akantu is infinitesimal (dW) and we
     // need to integrate a while to get a decent value.
     double Etot0 =
         this->model->getEnergy("potential") + this->model->getEnergy("kinetic");
     double W = 0.0;
     for (UInt i = 0; i < 200; ++i) {
       /// Solve.
       this->model->solveStep();
 
       const auto dW = this->model->getEnergy("external work");
 
       W += dW;
     }
     // Finally check.
     const auto Epot = this->model->getEnergy("potential");
     const auto Ekin = this->model->getEnergy("kinetic");
 
     EXPECT_NEAR(W, Ekin + Epot - Etot0, 5e-2);
     // Sadly not very exact for such a coarse mesh.
   }
 }
-}
+} // namespace
diff --git a/test/test_model/test_solid_mechanics_model/test_energies/test_solid_mechanics_model_work_quasistatic.cc b/test/test_model/test_solid_mechanics_model/test_energies/test_solid_mechanics_model_work_quasistatic.cc
index e4e48cadc..c45346f01 100644
--- a/test/test_model/test_solid_mechanics_model/test_energies/test_solid_mechanics_model_work_quasistatic.cc
+++ b/test/test_model/test_solid_mechanics_model/test_energies/test_solid_mechanics_model_work_quasistatic.cc
@@ -1,147 +1,147 @@
 /**
  * @file   test_solid_mechanics_model_work_quasistatic.cc
  *
  * @author Tobias Brink <tobias.brink@epfl.ch>
  *
  * @date creation: Wed Nov 29 2017
  * @date last modification: Fri Jan 26 2018
  *
  * @brief  test work in quasistatic
  *
  * @section LICENSE
  *
  * Copyright (©) 2016-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/>.
  *
  * @section description
  *
  * Assuming that the potential energy of a linear elastic material
  * works correctly, the work in a static simulation must equal the
  * potential energy of the material. Since the work in static is an
  * infinitesimal work Fds, we need to integrate by increasing F from 0
  * to F_final in steps. This test uses one Dirichlet boundary
  * condition (with u = 0.0, 0.1, and -0.1) and one Neumann boundary
  * condition for F on the opposite side. The final work must be the
  * same for all u.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 #include "../test_solid_mechanics_model_fixture.hh"
 #include "mesh_utils.hh"
 /* -------------------------------------------------------------------------- */
 
 using namespace akantu;
 
 namespace {
 
 TYPED_TEST(TestSMMFixture, WorkQuasistatic) {
   const auto spatial_dimension = this->spatial_dimension;
 
   getStaticParser().parse("test_solid_mechanics_model_"
                           "work_material.dat");
 
   /// model initialization
   MeshUtils::buildFacets(*this->mesh);
   this->model->initFull(_analysis_method = _static);
 
   /// Create a node group for Neumann BCs.
   auto & apply_force_grp = this->mesh->createNodeGroup("apply_force");
   auto & fixed_grp = this->mesh->createNodeGroup("fixed");
 
   const auto & pos = this->mesh->getNodes();
   auto & flags = this->model->getBlockedDOFs();
   auto & lower = this->mesh->getLowerBounds();
   auto & upper = this->mesh->getUpperBounds();
 
   UInt i = 0;
   for (auto && data : zip(make_view(pos, spatial_dimension),
                           make_view(flags, spatial_dimension))) {
     const auto & posv = std::get<0>(data);
     auto & flag = std::get<1>(data);
     if (posv(_x) > upper(_x) - 1e-6) {
       apply_force_grp.add(i);
     } else if (posv(_x) < lower(_x) + 1e-6) {
       fixed_grp.add(i);
 
       if ((spatial_dimension > 1) and (posv(_y) < lower(_y) + 1e-6)) {
         flag(_y) = true;
         if ((spatial_dimension > 2) and (posv(_z) < lower(_z) + 1e-6)) {
           flag(_z) = true;
         }
       }
     }
     ++i;
   }
 
   this->mesh->createElementGroupFromNodeGroup("el_apply_force", "apply_force",
                                               spatial_dimension - 1);
   this->mesh->createElementGroupFromNodeGroup("el_fixed", "fixed",
                                               spatial_dimension - 1);
 
   std::vector<Real> displacements{0.0, 0.1, -0.1};
   for (auto && u : displacements) {
-    this->model->applyBC(BC::Dirichlet::FixedValue(u, _x), "el_fixed");
+    this->model->applyBC(BC::FixedValue(u, _x), "el_fixed");
 
     Vector<Real> surface_traction(spatial_dimension);
     Real work = 0.0;
     Real Epot;
     static const UInt N = 100;
     for (UInt i = 0; i <= N; ++i) {
       this->model->getForce().clear(); // reset external forces to zero
 
       surface_traction(_x) = (1.0 * i) / N;
 
       if (spatial_dimension == 1) {
         // \TODO: this is a hack to work
         //      around non-implemented
         //      BC::Neumann::FromTraction for 1D
         auto & force = this->model->getForce();
         for (auto && pair : zip(make_view(pos, spatial_dimension),
                                 make_view(force, spatial_dimension))) {
           auto & posv = std::get<0>(pair);
           auto & forcev = std::get<1>(pair);
           if (posv(_x) > upper(_x) - 1e-6) {
             forcev(_x) = surface_traction(_x);
           }
         }
       } else {
-        this->model->applyBC(BC::Neumann::FromTraction(surface_traction),
+        this->model->applyBC(BC::FromTraction(surface_traction),
                              "el_apply_force");
       }
 
       /// Solve.
       this->model->solveStep();
 
       Epot = this->model->getEnergy("potential");
       // In static, this is infinitesimal work!
       auto Fds = this->model->getEnergy("external work");
 
       work += Fds; // integrate
 
       /// Check that no work was done for zero force.
       if (i == 0) {
         EXPECT_NEAR(work, 0.0, 1e-12);
       }
     }
 
     // Due to the finite integration steps, we make a rather large error
     // in our work integration, thus the allowed delta is 1e-2.
     EXPECT_NEAR(work, Epot, 1e-2);
   }
 }
 
 } // namespace
diff --git a/test/test_model/test_solid_mechanics_model/test_materials/test_local_material.cc b/test/test_model/test_solid_mechanics_model/test_materials/test_local_material.cc
index 5f085aa7a..d31133f04 100644
--- a/test/test_model/test_solid_mechanics_model/test_materials/test_local_material.cc
+++ b/test/test_model/test_solid_mechanics_model/test_materials/test_local_material.cc
@@ -1,119 +1,119 @@
 /**
  * @file   test_local_material.cc
  *
  * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
  * @author Marion Estelle Chambart <marion.chambart@epfl.ch>
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  * @author Clement Roux <clement.roux@epfl.ch>
  *
  * @date creation: Wed Aug 04 2010
  * @date last modification: Thu Dec 14 2017
  *
  * @brief  test of the class SolidMechanicsModel with custom local damage on a
  * notched plate
  *
  * @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 <iostream>
 /* -------------------------------------------------------------------------- */
 #include "local_material_damage.hh"
 #include "solid_mechanics_model.hh"
 /* -------------------------------------------------------------------------- */
 
 using namespace akantu;
 
 int main(int argc, char * argv[]) {
   akantu::initialize("material.dat", argc, argv);
   UInt max_steps = 1100;
 
   const UInt spatial_dimension = 2;
   Mesh mesh(spatial_dimension);
   mesh.read("mesh_section_gap.msh");
 
   /// model initialization
   MaterialFactory::getInstance().registerAllocator(
       "local_damage", [](UInt, const ID &, SolidMechanicsModel & model,
                          const ID & id) -> std::unique_ptr<Material> {
         return std::make_unique<LocalMaterialDamage>(model, id);
       });
 
   SolidMechanicsModel model(mesh);
   model.initFull();
 
   std::cout << model.getMaterial(0) << std::endl;
 
   model.addDumpField("damage");
 
   model.addDumpField("strain");
   model.addDumpField("stress");
   model.addDumpFieldVector("displacement");
   model.addDumpFieldVector("external_force");
   model.addDumpFieldVector("internal_force");
   model.dump();
 
   Real time_step = model.getStableTimeStep();
   model.setTimeStep(time_step / 2.5);
 
   /// Dirichlet boundary conditions
-  model.applyBC(BC::Dirichlet::FixedValue(0.0, _x), "Fixed");
+  model.applyBC(BC::FixedValue(0.0, _x), "Fixed");
   // model.applyBC(BC::Dirichlet::FixedValue(0.0, _y), "Fixed");
   Matrix<Real> stress(2, 2);
   stress.eye(5e7);
-  model.applyBC(BC::Neumann::FromHigherDim(stress), "Traction");
+  model.applyBC(BC::FromHigherDim(stress), "Traction");
 
   for (UInt s = 0; s < max_steps; ++s)
     model.solveStep();
 
   model.dump();
   const auto & lower_bounds = mesh.getLowerBounds();
   const auto & upper_bounds = mesh.getUpperBounds();
   Real L = upper_bounds(_x) - lower_bounds(_x);
   Real H = upper_bounds(_y) - lower_bounds(_y);
 
   const auto & filter = model.getMaterial("concrete").getElementFilter();
 
   Vector<Real> barycenter(spatial_dimension);
 
   for (auto & type : filter.elementTypes(spatial_dimension)) {
     UInt nb_elem = mesh.getNbElement(type);
     const UInt nb_gp = model.getFEEngine().getNbIntegrationPoints(type);
     const auto & material_damage_array =
         model.getMaterial(0).getArray<Real>("damage", type);
     UInt cpt = 0;
     for (auto nel : arange(nb_elem)) {
       mesh.getBarycenter({type, nel, _not_ghost}, barycenter);
       if ((std::abs(barycenter(_x) - (L / 2) + 0.025) < 0.025) &&
           (std::abs(barycenter(_y) - (H / 2) + 0.045) < 0.045)) {
         if (material_damage_array(cpt, 0) < 0.9) {
           std::terminate();
         } else {
           std::cout << "element " << nel << " is correctly broken" << std::endl;
         }
       }
 
       cpt += nb_gp;
     }
   }
 
   akantu::finalize();
   return 0;
 }
diff --git a/test/test_model/test_solid_mechanics_model/test_materials/test_material_elasto_plastic_linear_isotropic_hardening/test_material_elasto_plastic_linear_isotropic_hardening.cc b/test/test_model/test_solid_mechanics_model/test_materials/test_material_elasto_plastic_linear_isotropic_hardening/test_material_elasto_plastic_linear_isotropic_hardening.cc
index 90fc154a9..4b710dff2 100644
--- a/test/test_model/test_solid_mechanics_model/test_materials/test_material_elasto_plastic_linear_isotropic_hardening/test_material_elasto_plastic_linear_isotropic_hardening.cc
+++ b/test/test_model/test_solid_mechanics_model/test_materials/test_material_elasto_plastic_linear_isotropic_hardening/test_material_elasto_plastic_linear_isotropic_hardening.cc
@@ -1,91 +1,91 @@
 /**
  * @file   test_material_elasto_plastic_linear_isotropic_hardening.cc
  *
  * @author Jaehyun Cho <jaehyun.cho@epfl.ch>
  * @author Lucas Frerot <lucas.frerot@epfl.ch>
  *
  * @date creation: Sun Oct 19 2014
  * @date last modification: Mon Sep 11 2017
  *
  * @brief  test for material type elasto plastic linear isotropic hardening
  * using tension-compression test
  *
  * @section LICENSE
  *
  * Copyright (©) 2015-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 "non_linear_solver.hh"
 #include "solid_mechanics_model.hh"
 /* -------------------------------------------------------------------------- */
 #include <iostream>
 /* -------------------------------------------------------------------------- */
 
 using namespace akantu;
 
 /* -------------------------------------------------------------------------- */
 
 int main(int argc, char * argv[]) {
   initialize("test_material_elasto_plastic_linear_isotropic_hardening.dat",
              argc, argv);
 
   const UInt spatial_dimension = 2;
   const Real u_increment = 0.1;
   const UInt steps = 20;
 
   Mesh mesh(spatial_dimension);
   mesh.read("test_material_elasto_plastic_linear_isotropic_hardening.msh");
 
   SolidMechanicsModel model(mesh);
   model.initFull(_analysis_method = _static);
 
   auto & solver = model.getNonLinearSolver("static");
   solver.set("max_iterations", 300);
   solver.set("threshold", 1e-5);
 
-  model.applyBC(BC::Dirichlet::FixedValue(0.0, _x), "left");
-  model.applyBC(BC::Dirichlet::FixedValue(0.0, _y), "bottom");
+  model.applyBC(BC::FixedValue(0.0, _x), "left");
+  model.applyBC(BC::FixedValue(0.0, _y), "bottom");
 
   std::cout.precision(4);
   for (UInt i = 0; i < steps; ++i) {
 
-    model.applyBC(BC::Dirichlet::FixedValue(i * u_increment, _x), "right");
+    model.applyBC(BC::FixedValue(i * u_increment, _x), "right");
 
     try {
       model.solveStep();
     } catch (debug::NLSNotConvergedException & e) {
       std::cout << e.niter << " " << e.error << std::endl;
       throw;
     }
     Real strainxx = i * u_increment / 10.;
 
     const Array<UInt> & edge_nodes = mesh.getElementGroup("right").getNodes();
     Array<Real> & residual = model.getInternalForce();
     Real reaction = 0;
 
     for (UInt n = 0; n < edge_nodes.size(); n++) {
       reaction -= residual(edge_nodes(n), 0);
     }
 
     std::cout << strainxx << "," << reaction << std::endl;
   }
 
   finalize();
   return 0;
 }
diff --git a/test/test_model/test_solid_mechanics_model/test_materials/test_multi_material_elastic.cc b/test/test_model/test_solid_mechanics_model/test_materials/test_multi_material_elastic.cc
index 554c706c1..cf0c0b98a 100644
--- a/test/test_model/test_solid_mechanics_model/test_materials/test_multi_material_elastic.cc
+++ b/test/test_model/test_solid_mechanics_model/test_materials/test_multi_material_elastic.cc
@@ -1,127 +1,127 @@
 /**
  * @file   test_multi_material_elastic.cc
  *
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Fri Mar 03 2017
  * @date last modification: Tue Feb 20 2018
  *
  * @brief  Test with 2 elastic materials
  *
  * @section LICENSE
  *
  * Copyright (©) 2016-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 "non_linear_solver.hh"
 #include <solid_mechanics_model.hh>
 
 using namespace akantu;
 
 int main(int argc, char * argv[]) {
   initialize("test_multi_material_elastic.dat", argc, argv);
 
   UInt spatial_dimension = 2;
   Mesh mesh(spatial_dimension);
 
   mesh.read("test_multi_material_elastic.msh");
 
   SolidMechanicsModel model(mesh);
 
   auto && mat_sel = std::make_shared<MeshDataMaterialSelector<std::string>>(
       "physical_names", model);
   model.setMaterialSelector(mat_sel);
 
   model.initFull(_analysis_method = _static);
 
-  model.applyBC(BC::Dirichlet::FlagOnly(_y), "ground");
-  model.applyBC(BC::Dirichlet::FlagOnly(_x), "corner");
+  model.applyBC(BC::FlagOnly(_y), "ground");
+  model.applyBC(BC::FlagOnly(_x), "corner");
 
   Vector<Real> trac(spatial_dimension, 0.);
   trac(_y) = 1.;
-  model.applyBC(BC::Neumann::FromTraction(trac), "air");
+  model.applyBC(BC::FromTraction(trac), "air");
 
   model.addDumpField("external_force");
   model.addDumpField("internal_force");
   model.addDumpField("blocked_dofs");
   model.addDumpField("displacement");
   model.addDumpField("stress");
   model.addDumpField("grad_u");
 
   // model.dump();
   auto & solver = model.getNonLinearSolver("static");
   solver.set("max_iterations", 1);
   solver.set("threshold", 1e-8);
   solver.set("convergence_type", _scc_residual);
 
   model.solveStep();
   // model.dump();
 
   std::map<std::string, Matrix<Real>> ref_strain;
   ref_strain["strong"] = Matrix<Real>(spatial_dimension, spatial_dimension, 0.);
   ref_strain["strong"](_y, _y) = .5;
 
   ref_strain["weak"] = Matrix<Real>(spatial_dimension, spatial_dimension, 0.);
   ref_strain["weak"](_y, _y) = 1;
 
   Matrix<Real> ref_stress(spatial_dimension, spatial_dimension, 0.);
   ref_stress(_y, _y) = 1.;
 
   std::vector<std::string> mats = {"strong", "weak"};
 
   typedef Array<Real>::const_matrix_iterator mat_it;
   auto check = [](mat_it it, mat_it end, const Matrix<Real> & ref) -> bool {
     for (; it != end; ++it) {
       Real dist = (*it - ref).norm<L_2>();
 
       // std::cout << *it << " " << dist << " " << (dist < 1e-10 ? "OK" : "Not
       // OK") << std::endl;
 
       if (dist > 1e-10)
         return false;
     }
 
     return true;
   };
 
   auto tit = mesh.firstType(spatial_dimension);
   auto tend = mesh.lastType(spatial_dimension);
   for (; tit != tend; ++tit) {
     auto & type = *tit;
 
     for (auto mat_id : mats) {
       auto & stress = model.getMaterial(mat_id).getStress(type);
       auto & grad_u = model.getMaterial(mat_id).getGradU(type);
 
       auto sit = stress.begin(spatial_dimension, spatial_dimension);
       auto send = stress.end(spatial_dimension, spatial_dimension);
 
       auto git = grad_u.begin(spatial_dimension, spatial_dimension);
       auto gend = grad_u.end(spatial_dimension, spatial_dimension);
 
       if (!check(sit, send, ref_stress))
         AKANTU_ERROR("The stresses are not correct");
       if (!check(git, gend, ref_strain[mat_id]))
         AKANTU_ERROR("The grad_u are not correct");
     }
   }
 
   finalize();
 
   return 0;
 }
diff --git a/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_dynamics.cc b/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_dynamics.cc
index 5e9c2fb79..6f1a6a360 100644
--- a/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_dynamics.cc
+++ b/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_dynamics.cc
@@ -1,313 +1,313 @@
 /**
  * @file   test_solid_mechanics_model_dynamics.cc
  *
  * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
  *
  * @date creation: Wed Nov 29 2017
  * @date last modification: Wed Feb 21 2018
  *
  * @brief  test of the class SolidMechanicsModel on the 3d cube
  *
  * @section LICENSE
  *
  * Copyright (©) 2016-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 "boundary_condition_functor.hh"
 #include "test_solid_mechanics_model_fixture.hh"
 /* -------------------------------------------------------------------------- */
 
 using namespace akantu;
 
 namespace {
 
 const Real A = 1e-1;
 const Real E = 1.;
 const Real poisson = 3. / 10;
 const Real lambda = E * poisson / (1 + poisson) / (1 - 2 * poisson);
 const Real mu = E / 2 / (1. + poisson);
 const Real rho = 1.;
 const Real cp = std::sqrt((lambda + 2 * mu) / rho);
 const Real cs = std::sqrt(mu / rho);
 const Real c = std::sqrt(E / rho);
 
 const Vector<Real> k = {.5, 0., 0.};
 const Vector<Real> psi1 = {0., 0., 1.};
 const Vector<Real> psi2 = {0., 1., 0.};
 const Real knorm = k.norm();
 
 /* -------------------------------------------------------------------------- */
 template <UInt dim> struct Verification {};
 
 /* -------------------------------------------------------------------------- */
 template <> struct Verification<1> {
   void displacement(Vector<Real> & disp, const Vector<Real> & coord,
                     Real current_time) {
     const auto & x = coord(_x);
     const Real omega = c * k[0];
     disp(_x) = A * cos(k[0] * x - omega * current_time);
   }
 
   void velocity(Vector<Real> & vel, const Vector<Real> & coord,
                 Real current_time) {
     const auto & x = coord(_x);
     const Real omega = c * k[0];
     vel(_x) = omega * A * sin(k[0] * x - omega * current_time);
   }
 };
 
 /* -------------------------------------------------------------------------- */
 template <> struct Verification<2> {
   void displacement(Vector<Real> & disp, const Vector<Real> & X,
                     Real current_time) {
     Vector<Real> kshear = {k[1], k[0]};
     Vector<Real> kpush = {k[0], k[1]};
 
     const Real omega_p = knorm * cp;
     const Real omega_s = knorm * cs;
 
     Real phase_p = X.dot(kpush) - omega_p * current_time;
     Real phase_s = X.dot(kpush) - omega_s * current_time;
 
     disp = A * (kpush * cos(phase_p) + kshear * cos(phase_s));
   }
 
   void velocity(Vector<Real> & vel, const Vector<Real> & X, Real current_time) {
     Vector<Real> kshear = {k[1], k[0]};
     Vector<Real> kpush = {k[0], k[1]};
 
     const Real omega_p = knorm * cp;
     const Real omega_s = knorm * cs;
 
     Real phase_p = X.dot(kpush) - omega_p * current_time;
     Real phase_s = X.dot(kpush) - omega_s * current_time;
 
     vel =
         A * (kpush * omega_p * cos(phase_p) + kshear * omega_s * cos(phase_s));
   }
 };
 
 /* -------------------------------------------------------------------------- */
 template <> struct Verification<3> {
   void displacement(Vector<Real> & disp, const Vector<Real> & coord,
                     Real current_time) {
     const auto & X = coord;
     Vector<Real> kpush = k;
     Vector<Real> kshear1(3);
     Vector<Real> kshear2(3);
     kshear1.crossProduct(k, psi1);
     kshear2.crossProduct(k, psi2);
 
     const Real omega_p = knorm * cp;
     const Real omega_s = knorm * cs;
 
     Real phase_p = X.dot(kpush) - omega_p * current_time;
     Real phase_s = X.dot(kpush) - omega_s * current_time;
 
     disp = A * (kpush * cos(phase_p) + kshear1 * cos(phase_s) +
                 kshear2 * cos(phase_s));
   }
 
   void velocity(Vector<Real> & vel, const Vector<Real> & coord,
                 Real current_time) {
     const auto & X = coord;
     Vector<Real> kpush = k;
     Vector<Real> kshear1(3);
     Vector<Real> kshear2(3);
     kshear1.crossProduct(k, psi1);
     kshear2.crossProduct(k, psi2);
 
     const Real omega_p = knorm * cp;
     const Real omega_s = knorm * cs;
 
     Real phase_p = X.dot(kpush) - omega_p * current_time;
     Real phase_s = X.dot(kpush) - omega_s * current_time;
 
     vel =
         A * (kpush * omega_p * cos(phase_p) + kshear1 * omega_s * cos(phase_s) +
              kshear2 * omega_s * cos(phase_s));
   }
 };
 
 /* -------------------------------------------------------------------------- */
 template <ElementType _type>
-class SolutionFunctor : public BC::Dirichlet::DirichletFunctor {
+class SolutionFunctor : public BC::DirichletFunctor {
 public:
   SolutionFunctor(Real current_time, SolidMechanicsModel & model)
       : current_time(current_time), model(model) {}
 
 public:
   static constexpr UInt dim = ElementClass<_type>::getSpatialDimension();
 
   inline void operator()(UInt node, Vector<bool> & flags, Vector<Real> & primal,
                          const Vector<Real> & coord) const {
     flags(0) = true;
     auto & vel = model.getVelocity();
     auto it = vel.begin(model.getSpatialDimension());
     Vector<Real> v = it[node];
 
     Verification<dim> verif;
     verif.displacement(primal, coord, current_time);
     verif.velocity(v, coord, current_time);
   }
 
 private:
   Real current_time;
   SolidMechanicsModel & model;
 };
 
 /* -------------------------------------------------------------------------- */
 // This fixture uses somewhat finer meshes representing bars.
 template <typename type_, typename analysis_method_>
 class TestSMMFixtureBar : public TestSMMFixture<type_> {
   using parent = TestSMMFixture<type_>;
 
 public:
   void SetUp() override {
     this->mesh_file =
         "../patch_tests/data/bar" + aka::to_string(this->type) + ".msh";
     parent::SetUp();
 
     auto analysis_method = analysis_method_::value;
     this->initModel("test_solid_mechanics_model_"
                     "dynamics_material.dat", analysis_method);
 
     const auto & position = this->mesh->getNodes();
     auto & displacement = this->model->getDisplacement();
     auto & velocity = this->model->getVelocity();
 
     constexpr auto dim = parent::spatial_dimension;
 
     Verification<dim> verif;
 
     for (auto && tuple :
          zip(make_view(position, dim), make_view(displacement, dim),
              make_view(velocity, dim))) {
       verif.displacement(std::get<1>(tuple), std::get<0>(tuple), 0.);
       verif.velocity(std::get<2>(tuple), std::get<0>(tuple), 0.);
     }
 
     if (this->dump_paraview)
       this->model->dump();
 
     /// boundary conditions
     this->model->applyBC(SolutionFunctor<parent::type>(0., *this->model), "BC");
 
     wave_velocity = 1.; // sqrt(E/rho) = sqrt(1/1) = 1
     simulation_time = 5 / wave_velocity;
     time_step = this->model->getTimeStep();
 
     max_steps = simulation_time / time_step; // 100
   }
 
   void solveStep() {
     constexpr auto dim = parent::spatial_dimension;
     Real current_time = 0.;
     const auto & position = this->mesh->getNodes();
     const auto & displacement = this->model->getDisplacement();
     UInt nb_nodes = this->mesh->getNbNodes();
     UInt nb_global_nodes = this->mesh->getNbGlobalNodes();
 
     max_error = -1.;
 
     Array<Real> displacement_solution(nb_nodes, dim);
     Verification<dim> verif;
 
     auto ndump = 50;
     auto dump_freq = max_steps / ndump;
 
     for (UInt s = 0; s < this->max_steps;
          ++s, current_time += this->time_step) {
       if (s % dump_freq == 0 && this->dump_paraview)
         this->model->dump();
 
       /// boundary conditions
       this->model->applyBC(
           SolutionFunctor<parent::type>(current_time, *this->model), "BC");
 
       // compute the disp solution
       for (auto && tuple : zip(make_view(position, dim),
                                make_view(displacement_solution, dim))) {
         verif.displacement(std::get<1>(tuple), std::get<0>(tuple),
                            current_time);
       }
 
       // compute the error solution
       Real disp_error = 0.;
       auto n = 0;
       for (auto && tuple : zip(make_view(displacement, dim),
                                make_view(displacement_solution, dim))) {
         if (this->mesh->isLocalOrMasterNode(n)) {
           auto diff = std::get<1>(tuple) - std::get<0>(tuple);
           disp_error += diff.dot(diff);
         }
         ++n;
       }
 
       this->mesh->getCommunicator().allReduce(disp_error,
                                               SynchronizerOperation::_sum);
 
       disp_error = sqrt(disp_error) / nb_global_nodes;
       max_error = std::max(disp_error, max_error);
 
       this->model->solveStep();
     }
   }
 
 protected:
   Real time_step;
   Real wave_velocity;
   Real simulation_time;
   UInt max_steps;
   Real max_error{-1};
 };
 
 template <AnalysisMethod t>
 using analysis_method_t = std::integral_constant<AnalysisMethod, t>;
 
 using TestTypes = gtest_list_t<TestElementTypes>;
 
 template <typename type_>
 using TestSMMFixtureBarExplicit =
     TestSMMFixtureBar<type_, analysis_method_t<_explicit_lumped_mass>>;
 
 TYPED_TEST_CASE(TestSMMFixtureBarExplicit, TestTypes);
 
 /* -------------------------------------------------------------------------- */
 TYPED_TEST(TestSMMFixtureBarExplicit, Dynamics) {
   this->solveStep();
   EXPECT_NEAR(this->max_error, 0., 2e-3);
   // std::cout << "max error: " << max_error << std::endl;
 }
 
 
 /* -------------------------------------------------------------------------- */
 template <typename type_>
 using TestSMMFixtureBarImplicit =
     TestSMMFixtureBar<type_, analysis_method_t<_implicit_dynamic>>;
 
 TYPED_TEST_CASE(TestSMMFixtureBarImplicit, TestTypes);
 
 TYPED_TEST(TestSMMFixtureBarImplicit, Dynamics) {
   this->solveStep();
   EXPECT_NEAR(this->max_error, 0., 2e-3);
 }
 
 
 }
diff --git a/test/test_python_interface/pybind11/CMakeLists.txt b/test/test_python_interface/pybind11/CMakeLists.txt
index 334aef6ce..aa11d71e0 100644
--- a/test/test_python_interface/pybind11/CMakeLists.txt
+++ b/test/test_python_interface/pybind11/CMakeLists.txt
@@ -1,10 +1,9 @@
 pybind11_add_module(py11_akantu_test_common test_common.cc)
 target_link_libraries(py11_akantu_test_common PRIVATE akantu)
 set_property(TARGET py11_akantu_test_common PROPERTY DEBUG_POSTFIX "")
 
 register_test(test_pybind11_common
   SCRIPT test_common.py
-  DEPENDS py11_akantu
   PACKAGE pybind11 python_interface
   PYTHON
   )
diff --git a/test/test_python_interface/pybind11/test_common.cc b/test/test_python_interface/pybind11/test_common.cc
index 13c20dd1e..3bf8976c0 100644
--- a/test/test_python_interface/pybind11/test_common.cc
+++ b/test/test_python_interface/pybind11/test_common.cc
@@ -1,13 +1,55 @@
-#include <pybind11/pybind11.h>
 #include <aka_array.hh>
+#include <pybind11/pybind11.h>
+
+#include "../../python/pybind11/py_aka_array.cc"
+#include "../../python/pybind11/py_aka_boundary_conditions.cc"
+#include "../../python/pybind11/py_aka_common.cc"
+#include "../../python/pybind11/py_aka_solid_mechanics_model.cc"
 
 namespace py = pybind11;
-using namespace akantu;
+namespace _aka = akantu;
+
+std::map<long, std::shared_ptr<_aka::Array<_aka::Real>>> arrays;
 
 PYBIND11_MODULE(py11_akantu_test_common, mod) {
   mod.doc() = "Akantu Test function for common ";
 
-  Array<Real> data{1000, 3};
+  register_enums(mod);
+  register_boundary_conditions(mod);
+  register_solid_mechanics_models(mod);
+
+  mod.def("createData",
+          [&](_aka::UInt size, _aka::UInt nb_components) {
+            auto ptr =
+                std::make_shared<_aka::Array<_aka::Real>>(size, nb_components);
+            ptr->clear();
+            long addr = (long)ptr->storage();
+            py::print("initial pointer: " + std::to_string(addr));
+            arrays[addr] = ptr;
+            return std::tuple<long, _aka::Array<_aka::Real> &>(addr, *ptr);
+          },
+          py::return_value_policy::reference);
+  mod.def("getData",
+          [&](long addr) -> _aka::Array<_aka::Real> & {
+            auto & array = *arrays[addr];
+            py::print("gotten pointer: " +
+                      std::to_string((long)array.storage()));
+            return array;
+          },
+          py::return_value_policy::reference);
+
+  mod.def("getCopyData",
+          [&](long addr) -> _aka::Array<_aka::Real> {
+            auto & array = *arrays[addr];
+            py::print("gotten pointer: " +
+                      std::to_string((long)array.storage()));
+            return array;
+          },
+          py::return_value_policy::copy);
 
-  mod.def("getData", [&data]() -> Array<Real>& { return data; });
+  mod.def("getRawPointer", [](_aka::Array<_aka::Real> & _data) {
+    py::print("received proxy: " + std::to_string((long)&_data));
+    py::print("raw pointer: " + std::to_string((long)_data.storage()));
+    return (long)_data.storage();
+  });
 } // Module akantu_test_common
diff --git a/test/test_python_interface/pybind11/test_common.py b/test/test_python_interface/pybind11/test_common.py
index 3f253b395..f8daf06f1 100644
--- a/test/test_python_interface/pybind11/test_common.py
+++ b/test/test_python_interface/pybind11/test_common.py
@@ -1,8 +1,34 @@
 # import pybind11
-import py11_akantu as pyaka
+import numpy as np
 import py11_akantu_test_common as pyaka_test
 
 
-def test_array():
-    array = pyaka_test.getData()
+def test_array_size():
+    ptr, array = pyaka_test.createData(1000, 3)
     assert array.shape == (1000, 3)
+
+
+def test_array_nocopy():
+    ptr, array = pyaka_test.createData(1000, 3)
+    through_python = pyaka_test.getRawPointer(array)
+    assert(ptr == through_python)
+
+
+def test_modify_values():
+    ptr, array = pyaka_test.createData(3, 3)
+    array[0, :] = (1., 2., 3.)
+    array2 = pyaka_test.getData(ptr)
+    assert(np.linalg.norm(array-array2) < 1e-15)
+
+    for i in [1, 2, 3]:
+        ptr, array = pyaka_test.createData(10000, i)
+        array[:, :] = np.random.random((10000, i))
+        array2 = pyaka_test.getData(ptr)
+        assert(np.linalg.norm(array-array2) < 1e-15)
+
+
+def test_array_copy():
+    ptr, array = pyaka_test.createData(1000, 3)
+    array2 = pyaka_test.getCopyData(ptr)
+    ptr2 = pyaka_test.getRawPointer(array2)
+    assert(ptr != ptr2)