diff --git a/CMakeLists.txt b/CMakeLists.txt
index 508eed3b1..ce25d7b56 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,183 +1,187 @@
 #===============================================================================
 # @file   CMakeLists.txt
 #
 # @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
 # @author Nicolas Richart <nicolas.richart@epfl.ch>
 #
 # @date creation: Mon Jun 14 2010
 # @date last modification: Fri Jan 22 2016
 #
 # @brief  main configuration file
 #
 # @section LICENSE
 #
 # Copyright (©)  2010-2012, 2014,  2015 EPFL  (Ecole Polytechnique  Fédérale de
 # Lausanne)  Laboratory (LSMS  -  Laboratoire de  Simulation  en Mécanique  des
 # Solides)
 #
 # Akantu is free  software: you can redistribute it and/or  modify it under the
 # terms  of the  GNU Lesser  General Public  License as  published by  the Free
 # Software Foundation, either version 3 of the License, or (at your option) any
 # later version.
 #
 # Akantu is  distributed in the  hope that it  will be useful, but  WITHOUT ANY
 # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
 # A  PARTICULAR PURPOSE. See  the GNU  Lesser General  Public License  for more
 # details.
 #
 # You should  have received  a copy  of the GNU  Lesser General  Public License
 # along with Akantu. If not, see <http://www.gnu.org/licenses/>.
 #
 # @section DESCRIPTION
 #-------------------------------------------------------------------------------
 #                      _               _
 #                     | |             | |
 #                 __ _| | ____ _ _ __ | |_ _   _
 #                / _` | |/ / _` | '_ \| __| | | |
 #               | (_| |   < (_| | | | | |_| |_| |
 #                \__,_|_|\_\__,_|_| |_|\__|\__,_|
 #
 #===============================================================================
 
 #===============================================================================
 # CMake Project
 #===============================================================================
 cmake_minimum_required(VERSION 3.1.3)
 
 # add this options before PROJECT keyword
 set(CMAKE_DISABLE_SOURCE_CHANGES ON)
 set(CMAKE_DISABLE_IN_SOURCE_BUILD ON)
 
 project(Akantu)
 enable_language(CXX)
 
 #===============================================================================
 # Misc. config for cmake
 #===============================================================================
 set(AKANTU_CMAKE_DIR "${PROJECT_SOURCE_DIR}/cmake")
 set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${PROJECT_SOURCE_DIR}/cmake")
 set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${PROJECT_SOURCE_DIR}/cmake/Modules")
 set(BUILD_SHARED_LIBS ON CACHE BOOL "Build shared libraries.")
 mark_as_advanced(BUILD_SHARED_LIBS)
 
 if(NOT AKANTU_TARGETS_EXPORT)
   set(AKANTU_TARGETS_EXPORT AkantuLibraryDepends)
 endif()
 
 include(CMakeVersionGenerator)
 include(CMakePackagesSystem)
 include(CMakeFlagsHandling)
 
 include(AkantuPackagesSystem)
 include(AkantuMacros)
 include(AkantuCleaning)
 
 #cmake_activate_debug_message()
 
 #===============================================================================
 # Version Number
 #===============================================================================
 # AKANTU version number.  An even minor number corresponds to releases.
 set(AKANTU_MAJOR_VERSION 3)
 set(AKANTU_MINOR_VERSION 0)
 set(AKANTU_PATCH_VERSION 0)
 define_project_version()
 
 #===============================================================================
 # Options
 #===============================================================================
 # Debug
 set(CMAKE_CXX_FLAGS_RELEASE "-O3 -DNDEBUG -DAKANTU_NDEBUG"
   CACHE STRING "Flags used by the compiler during release builds" FORCE)
 
 #add_flags(cxx "-Wall -Wextra -pedantic -Werror")
 if (CMAKE_CXX_COMPILER_ID STREQUAL "GNU" OR CMAKE_CXX_COMPILER_ID STREQUAL "Clang")
   add_flags(cxx "-Wall -Wextra -pedantic") # -Weffc++
 else()
   add_flags(cxx "-Wall")
 endif()
 
 option(AKANTU_EXAMPLES "Activate examples" OFF)
 option(AKANTU_TESTS "Activate tests" OFF)
 
 include(AkantuExtraCompilationProfiles)
 #===============================================================================
 # Dependencies
 #===============================================================================
 declare_akantu_types()
 
 package_list_packages(${PROJECT_SOURCE_DIR}/packages
   EXTRA_PACKAGES_FOLDER ${PROJECT_SOURCE_DIR}/extra_packages
   NO_AUTO_COMPILE_FLAGS)
 
 ## meta option \todo better way to do it when multiple package give enable the
 ## same feature
 if(AKANTU_SCOTCH)
   set(AKANTU_PARTITIONER ON)
 else()
   set(AKANTU_PARTITIONER OFF)
 endif()
 
 if(AKANTU_MUMPS)
   set(AKANTU_SOLVER ON)
 else()
   set(AKANTU_SOLVER OFF)
 endif()
 
 #===============================================================================
 # Akantu library
 #===============================================================================
 add_subdirectory(src)
 
 #===============================================================================
 # Documentation
 #===============================================================================
 if(AKANTU_DOCUMENTATION_DOXYGEN OR AKANTU_DOCUMENTATION_MANUAL)
   add_subdirectory(doc)
 else()
   set(AKANTU_DOC_EXCLUDE_FILES "${PROJECT_SOURCE_DIR}/doc/manual" CACHE INTERNAL "")
 endif()
 
 #===============================================================================
 # Examples and tests
 #===============================================================================
 include(AkantuTestsMacros)
 include(AkantuExampleMacros)
 
 if(AKANTU_TESTS)
   option(AKANTU_BUILD_ALL_TESTS "Build all tests" ON)
   find_package(GMSH REQUIRED)
+  package_is_activated(pybind11 _pybind11_act)
+  if(_pybind11_act)
+    find_package(pybind11 CONFIG REQUIRED QUIET) # to get the pybind11_add_module macro
+  endif()
 endif()
 
 if(AKANTU_EXAMPLES)
   find_package(GMSH REQUIRED)
   add_subdirectory(examples)
 endif()
 
 # tests
 add_test_tree(test)
 
 #===============================================================================
 # Python interface
 #===============================================================================
 package_is_activated(python_interface _python_act)
 if(_python_act)
   if(IS_ABSOLUTE "${CMAKE_INSTALL_PREFIX}")
     set(AKANTU_PYTHON_INSTALL_PREFIX ${CMAKE_INSTALL_PREFIX})
   else()
     set(AKANTU_PYTHON_INSTALL_PREFIX "${PROJECT_BINARY_DIR}/${CMAKE_INSTALL_PREFIX}")
   endif()
   add_subdirectory(python)
 endif()
 
 #===============================================================================
 # Install and Packaging
 #===============================================================================
 include(AkantuInstall)
 
 option(AKANTU_DISABLE_CPACK
   "This option commands the generation of extra info for the \"make package\" target" ON)
 mark_as_advanced(AKANTU_DISABLE_CPACK)
 if(NOT AKANTU_DISABLE_CPACK)
   include(AkantuCPack)
 endif()
diff --git a/src/common/aka_array.hh b/src/common/aka_array.hh
index c580f85f1..7694df0b2 100644
--- a/src/common/aka_array.hh
+++ b/src/common/aka_array.hh
@@ -1,371 +1,371 @@
 /**
  * @file   aka_array.hh
  *
  * @author Till Junge <till.junge@epfl.ch>
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Fri Jun 18 2010
  * @date last modification: Fri Jan 22 2016
  *
  * @brief  Array container for Akantu
  * This container differs from the std::vector from the fact it as 2 dimensions
  * a main dimension and the size stored per entries
  *
  * @section LICENSE
  *
  * Copyright (©)  2010-2012, 2014,  2015 EPFL  (Ecole Polytechnique  Fédérale de
  * Lausanne)  Laboratory (LSMS  -  Laboratoire de  Simulation  en Mécanique  des
  * Solides)
  *
  * Akantu is free  software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as  published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Akantu is  distributed in the  hope that it  will be useful, but  WITHOUT ANY
  * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
  * A  PARTICULAR PURPOSE. See  the GNU  Lesser General  Public License  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 
 #ifndef __AKANTU_VECTOR_HH__
 #define __AKANTU_VECTOR_HH__
 
 /* -------------------------------------------------------------------------- */
 #include "aka_common.hh"
 /* -------------------------------------------------------------------------- */
 #include <typeinfo>
 #include <vector>
 /* -------------------------------------------------------------------------- */
 
 namespace akantu {
 
 /// class that afford to store vectors in static memory
 class ArrayBase {
   /* ------------------------------------------------------------------------ */
   /* Constructors/Destructors                                                 */
   /* ------------------------------------------------------------------------ */
 public:
   explicit ArrayBase(ID id = "");
 
   virtual ~ArrayBase();
 
   /* ------------------------------------------------------------------------ */
   /* Methods                                                                  */
   /* ------------------------------------------------------------------------ */
 public:
   /// get the amount of space allocated in bytes
   inline UInt getMemorySize() const;
 
   /// set the size to zero without freeing the allocated space
   inline void empty();
 
   /// function to print the containt of the class
   virtual void printself(std::ostream & stream, int indent = 0) const;
 
   /* ------------------------------------------------------------------------ */
   /* Accessors */
   /* ------------------------------------------------------------------------ */
 public:
   /// Get the real size allocated in memory
   AKANTU_GET_MACRO(AllocatedSize, allocated_size, UInt);
   /// Get the Size of the Array
   UInt getSize() const __attribute__((deprecated)) { return size_; }
   UInt size() const { return size_; }
   /// Get the number of components
   AKANTU_GET_MACRO(NbComponent, nb_component, UInt);
   /// Get the name of th array
   AKANTU_GET_MACRO(ID, id, const ID &);
   /// Set the name of th array
   AKANTU_SET_MACRO(ID, id, const ID &);
 
   /* ------------------------------------------------------------------------ */
   /* Class Members                                                            */
   /* ------------------------------------------------------------------------ */
 protected:
   /// id of the vector
   ID id;
 
   /// the size allocated
   UInt allocated_size{0};
 
   /// the size used
   UInt size_{0};
 
   /// number of components
   UInt nb_component{1};
 
   /// size of the stored type
   UInt size_of_type{0};
 };
 
 /* -------------------------------------------------------------------------- */
 namespace {
   template <std::size_t dim, typename T> struct IteratorHelper {};
 
   template <typename T> struct IteratorHelper<0, T> { using type = T; };
   template <typename T> struct IteratorHelper<1, T> { using type = Vector<T>; };
   template <typename T> struct IteratorHelper<2, T> { using type = Matrix<T>; };
   template <typename T> struct IteratorHelper<3, T> {
     using type = Tensor3<T>;
   };
 
   template <std::size_t dim, typename T>
   using IteratorHelper_t = typename IteratorHelper<dim, T>::type;
 } // namespace
 
 /* -------------------------------------------------------------------------- */
 template <typename T, bool is_scal> class Array : public ArrayBase {
   /* ------------------------------------------------------------------------ */
   /* Constructors/Destructors                                                 */
   /* ------------------------------------------------------------------------ */
 public:
   using value_type = T;
   using reference = value_type &;
   using pointer_type = value_type *;
   using const_reference = const value_type &;
 
   /// Allocation of a new vector
   explicit inline Array(UInt size = 0, UInt nb_component = 1,
                         const ID & id = "");
 
   /// Allocation of a new vector with a default value
   Array(UInt size, UInt nb_component, const value_type def_values[],
         const ID & id = "");
 
   /// Allocation of a new vector with a default value
   Array(UInt size, UInt nb_component, const_reference value,
         const ID & id = "");
 
   /// Copy constructor (deep copy if deep=true)
   Array(const Array<value_type, is_scal> & vect, bool deep = true,
         const ID & id = "");
 
 #ifndef SWIG
   /// Copy constructor (deep copy)
   explicit Array(const std::vector<value_type> & vect);
 #endif
 
   inline ~Array() override;
 
   Array & operator=(const Array & a) {
     /// this is to let STL allocate and copy arrays in the case of
     /// std::vector::resize
-    AKANTU_DEBUG_ASSERT(this->size == 0, "Cannot copy akantu::Array");
+    AKANTU_DEBUG_ASSERT(this->size_ == 0, "Cannot copy akantu::Array");
     return const_cast<Array &>(a);
   }
 
 #ifndef SWIG
   /* ------------------------------------------------------------------------ */
   /* Iterator                                                                 */
   /* ------------------------------------------------------------------------ */
   /// \todo protected: does not compile with intel  check why
 public:
   template <class R, class it, class IR = R, bool is_tensor_ = is_tensor<R>{}>
   class iterator_internal;
 
 public:
   /* ------------------------------------------------------------------------ */
 
   /* ------------------------------------------------------------------------ */
   template <typename R = T> class const_iterator;
   template <typename R = T> class iterator;
 
   /* ------------------------------------------------------------------------ */
 
   /// iterator for Array of nb_component = 1
   using scalar_iterator = iterator<T>;
   /// const_iterator for Array of nb_component = 1
   using const_scalar_iterator = const_iterator<T>;
 
   /// iterator returning Vectors of size n  on entries of Array with
   /// nb_component = n
   using vector_iterator = iterator<Vector<T>>;
   /// const_iterator returning Vectors of n size on entries of Array with
   /// nb_component = n
   using const_vector_iterator = const_iterator<Vector<T>>;
 
   /// iterator returning Matrices of size (m, n) on entries of Array with
   /// nb_component = m*n
   using matrix_iterator = iterator<Matrix<T>>;
   /// const iterator returning Matrices of size (m, n) on entries of Array with
   /// nb_component = m*n
   using const_matrix_iterator = const_iterator<Matrix<T>>;
 
   /// iterator returning Tensor3 of size (m, n, k) on entries of Array with
   /// nb_component = m*n*k
   using tensor3_iterator = iterator<Tensor3<T>>;
   /// const iterator returning Tensor3 of size (m, n, k) on entries of Array
   /// with nb_component = m*n*k
   using const_tensor3_iterator = const_iterator<Tensor3<T>>;
 
   /* ------------------------------------------------------------------------ */
   template <typename... Ns> inline decltype(auto) begin(Ns&&... n);
   template <typename... Ns> inline decltype(auto) end(Ns&&... n);
 
   template <typename... Ns> inline decltype(auto) begin(Ns&&... n) const;
   template <typename... Ns> inline decltype(auto) end(Ns&&... n) const;
 
   template <typename... Ns> inline decltype(auto) begin_reinterpret(Ns&&... n);
   template <typename... Ns> inline decltype(auto) end_reinterpret(Ns&&... n);
 
   template <typename... Ns>
   inline decltype(auto) begin_reinterpret(Ns&&... n) const;
   template <typename... Ns>
   inline decltype(auto) end_reinterpret(Ns&&... n) const;
 #endif // SWIG
 
   /* ------------------------------------------------------------------------ */
   /* Methods                                                                  */
   /* ------------------------------------------------------------------------ */
 public:
   /// append a tuple of size nb_component containing value
   inline void push_back(const_reference value);
   /// append a vector
   // inline void push_back(const value_type new_elem[]);
 
 #ifndef SWIG
   /// append a Vector or a Matrix
   template <template <typename> class C,
             typename = std::enable_if_t<is_tensor<C<T>>{}>>
   inline void push_back(const C<T> & new_elem);
   /// append the value of the iterator
   template <typename Ret> inline void push_back(const iterator<Ret> & it);
 
   /// erase the value at position i
   inline void erase(UInt i);
   /// ask Nico, clarify
   template <typename R> inline iterator<R> erase(const iterator<R> & it);
 #endif
 
   /// changes the allocated size but not the size
   virtual void reserve(UInt size);
 
   /// change the size of the Array
   virtual void resize(UInt size);
 
   /// change the size of the Array and initialize the values
   virtual void resize(UInt size, const T & val);
 
   /// change the number of components by interlacing data
   /// @param multiplicator number of interlaced components add
   /// @param block_size blocks of data in the array
   /// Examaple for block_size = 2, multiplicator = 2
   /// array = oo oo oo -> new array = oo nn nn oo nn nn oo nn nn
   void extendComponentsInterlaced(UInt multiplicator, UInt stride);
 
   /// search elem in the vector, return  the position of the first occurrence or
   /// -1 if not found
   UInt find(const_reference elem) const;
 
   /// @see Array::find(const_reference elem) const
   UInt find(T elem[]) const;
 
 #ifndef SWIG
   /// @see Array::find(const_reference elem) const
   template <template <typename> class C,
             typename = std::enable_if_t<is_tensor<C<T>>{}>>
   inline UInt find(const C<T> & elem);
 #endif
 
   /// set all entries of the array to 0
   inline void clear() { std::fill_n(values, size_ * nb_component, T()); }
 
   /// set all entries of the array to the value t
   /// @param t value to fill the array with
   inline void set(T t) { std::fill_n(values, size_ * nb_component, t); }
 
 #ifndef SWIG
   /// set all tuples of the array to a given vector or matrix
   /// @param vm Matrix or Vector to fill the array with
   template <template <typename> class C,
             typename = std::enable_if_t<is_tensor<C<T>>{}>>
   inline void set(const C<T> & vm);
 #endif
 
   /// Append the content of the other array to the current one
   void append(const Array<T> & other);
 
   /// copy another Array in the current Array, the no_sanity_check allows you to
   /// force the copy in cases where you know what you do with two non matching
   /// Arrays in terms of n
   void copy(const Array<T, is_scal> & other, bool no_sanity_check = false);
 
   /// give the address of the memory allocated for this vector
   T * storage() const { return values; };
 
   /// function to print the containt of the class
   void printself(std::ostream & stream, int indent = 0) const override;
 
 protected:
   /// perform the allocation for the constructors
   void allocate(UInt size, UInt nb_component = 1);
 
   /// resize initializing with uninitialized_fill if fill is set
   void resizeUnitialized(UInt new_size, bool fill, const T & val = T());
 
   /* ------------------------------------------------------------------------ */
   /* Operators                                                                */
   /* ------------------------------------------------------------------------ */
 public:
   /// substraction entry-wise
   Array<T, is_scal> & operator-=(const Array<T, is_scal> & other);
   /// addition entry-wise
   Array<T, is_scal> & operator+=(const Array<T, is_scal> & other);
   /// multiply evry entry by alpha
   Array<T, is_scal> & operator*=(const T & alpha);
 
   /// check if the array are identical entry-wise
   bool operator==(const Array<T, is_scal> & other) const;
   /// @see Array::operator==(const Array<T, is_scal> & other) const
   bool operator!=(const Array<T, is_scal> & other) const;
 
   /// return a reference to the j-th entry of the i-th tuple
   inline reference operator()(UInt i, UInt j = 0);
   /// return a const reference to the j-th entry of the i-th tuple
   inline const_reference operator()(UInt i, UInt j = 0) const;
 
   /// return a reference to the ith component of the 1D array
   inline reference operator[](UInt i);
   /// return a const reference to the ith component of the 1D array
   inline const_reference operator[](UInt i) const;
 
   /* ------------------------------------------------------------------------ */
   /* Class Members                                                            */
   /* ------------------------------------------------------------------------ */
 protected:
   /// array of values
   T * values; // /!\ very dangerous
 };
 
 /* -------------------------------------------------------------------------- */
 /* Inline Functions Array<T, is_scal>                                         */
 /* -------------------------------------------------------------------------- */
 template <typename T, bool is_scal>
 inline std::ostream & operator<<(std::ostream & stream,
                                  const Array<T, is_scal> & _this) {
   _this.printself(stream);
   return stream;
 }
 
 /* -------------------------------------------------------------------------- */
 /* Inline Functions ArrayBase                                                 */
 /* -------------------------------------------------------------------------- */
 inline std::ostream & operator<<(std::ostream & stream,
                                  const ArrayBase & _this) {
   _this.printself(stream);
   return stream;
 }
 
 } // namespace akantu
 
 #include "aka_array_tmpl.hh"
 #include "aka_types.hh"
 
 #endif /* __AKANTU_VECTOR_HH__ */
diff --git a/src/fe_engine/element_classes/element_class_pentahedron_15_inline_impl.cc b/src/fe_engine/element_classes/element_class_pentahedron_15_inline_impl.cc
index b69d2d6b1..55b8b9d92 100644
--- a/src/fe_engine/element_classes/element_class_pentahedron_15_inline_impl.cc
+++ b/src/fe_engine/element_classes/element_class_pentahedron_15_inline_impl.cc
@@ -1,170 +1,173 @@
 /**
  * @file   element_class_pentahedron_15_inline_impl.cc
  *
  * @author Mauro Corrado <mauro.corrado@epfl.ch>
  * @author Sacha Laffely <sacha.laffely@epfl.ch>
  * @author Damien Scantamburlo <damien.scantamburlo@epfl.ch>
  *
  * @date creation: Tue Mar 31 2015
  * @date last modification: Wed Sep 16 2015
  *
- * @brief  Specialization of the element_class class for the type _pentahedron_15
+ * @brief  Specialization of the element_class class for the type
+_pentahedron_15
  *
  * @section LICENSE
  *
  * Copyright (©) 2015 EPFL (Ecole Polytechnique Fédérale de Lausanne) Laboratory
  * (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  * Akantu is free  software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as  published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Akantu is  distributed in the  hope that it  will be useful, but  WITHOUT ANY
  * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
  * A  PARTICULAR PURPOSE. See  the GNU  Lesser General  Public License  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
  *
  * @section DESCRIPTION
  *
  * @verbatim
              z
              ^
              |
              |
              |  1
              | /|\
              |/ | \
              10 7  6
             /   |   \
            /    |    \
            4    2--8--0
           | \  /      /
           |  \11     /
           13  12    9---------->y
           | /  \   /
           |/    \ /
           5--14--3
          /
         /
        /
       v
      x
        x   y    z
 * N0  -1   1    0
 * N1  -1   0    1
 * N2  -1   0    0
 * N3   1   1    0
 * N4   1   0    1
 * N5   1   0    0
 * N6  -1   0.5  0.5
 * N7  -1   0    0.5
 * N8  -1   0.5  0
 * N9   0   1    0
 * N10  0   0    1
 * N11  0   0    0
 * N12  1   0.5  0.5
 * N13  1   0    0.5
 * N14  1   0.5  0
  */
 
 /* -------------------------------------------------------------------------- */
-AKANTU_DEFINE_ELEMENT_CLASS_PROPERTY(_pentahedron_15,
-                                     _gt_pentahedron_15,
-                                     _itp_lagrange_pentahedron_15,
-                                     _ek_regular,
-                                     3,
-                                     _git_pentahedron,
-                                     2);
+AKANTU_DEFINE_ELEMENT_CLASS_PROPERTY(_pentahedron_15, _gt_pentahedron_15,
+                                     _itp_lagrange_pentahedron_15, _ek_regular,
+                                     3, _git_pentahedron, 2);
 
 /* -------------------------------------------------------------------------- */
 template <>
 template <class vector_type>
-inline void
-InterpolationElement<_itp_lagrange_pentahedron_15>::computeShapes(const vector_type & c,
-                                                                 vector_type & N) {
+inline void InterpolationElement<_itp_lagrange_pentahedron_15>::computeShapes(
+    const vector_type & c, vector_type & N) {
+  auto & x = c(0);
+  auto & y = c(1);
+  auto & z = c(2);
+
   // Shape Functions, Natural coordinates
-  N( 0) = 0.5 * c(1) * (1 - c(0)) * (2 * c(1) - 2 - c(0));
-  N( 1) = 0.5 * c(2) * (1 - c(0)) * (2 * c(2) - 2 - c(0));
-  N( 2) = 0.5 * (c(0) - 1) * (1 - c(1) - c(2)) * (c(0) + 2 * c(1) + 2 * c(2));
-  N( 3) = 0.5 * c(1) * (1 + c(0)) * (2 * c(1) - 2 + c(0));
-  N( 4) = 0.5 * c(2) * (1 + c(0)) * (2 * c(2) - 2 + c(0));
-  N( 5) = 0.5 * (-c(0) - 1) * (1 - c(1) - c(2)) * (-c(0) + 2 * c(1) + 2 * c(2));
-  N( 6) = 2.0 * c(1) * c(2) * (1 - c(0));
-  N( 7) = 2.0 * c(2) * (1 - c(1) - c(2)) * (1 - c(0));
-  N( 8) = 2.0 * c(1) * (1 - c(0)) * (1 - c(1) - c(2));
-  N( 9) = c(1) * (1 - c(0) * c(0));
-  N(10) = c(2) * (1 - c(0) * c(0));
-  N(11) = (1 - c(1) - c(2)) * (1 - c(0) * c(0));
-  N(12) = 2.0 * c(1) * c(2) * (1 + c(0));
-  N(13) = 2.0 * c(2) * (1 - c(1) - c(2)) * (1 + c(0));
-  N(14) = 2.0 * c(1) * (1 - c(1) - c(2)) * (1 + c(0));
+  N(0) = 0.5 * y * (1 - x) * (2 * y - 2 - x);
+  N(1) = 0.5 * z * (1 - x) * (2 * z - 2 - x);
+  N(2) = 0.5 * (x - 1) * (1 - y - z) * (x + 2 * y + 2 * z);
+  N(3) = 0.5 * y * (1 + x) * (2 * y - 2 + x);
+  N(4) = 0.5 * z * (1 + x) * (2 * z - 2 + x);
+  N(5) = 0.5 * (-x - 1) * (1 - y - z) * (-x + 2 * y + 2 * z);
+  N(6) = 2.0 * y * z * (1 - x);
+  N(7) = 2.0 * z * (1 - y - z) * (1 - x);
+  N(8) = 2.0 * y * (1 - x) * (1 - y - z);
+  N(9) = y * (1 - x * x);
+  N(10) = z * (1 - x * x);
+  N(11) = (1 - y - z) * (1 - x * x);
+  N(12) = 2.0 * y * z * (1 + x);
+  N(13) = 2.0 * z * (1 - y - z) * (1 + x);
+  N(14) = 2.0 * y * (1 - y - z) * (1 + x);
 }
 
 /* -------------------------------------------------------------------------- */
 template <>
 template <class vector_type, class matrix_type>
-inline void
-InterpolationElement<_itp_lagrange_pentahedron_15>::computeDNDS(const vector_type & c,
-                                                              matrix_type & dnds) {
-  //ddx
-  dnds(0, 0) =  0.5 * c(1) * (2 * c(0) - 2 * c(1) + 1);
-  dnds(0, 1) =  0.5 * c(2) * (2 * c(0) - 2 * c(2) + 1);
-  dnds(0, 2) = -0.5 * (2 * c(0) + 2 * c(1) + 2 * c(2) - 1) * (c(1) + c(2) -1);
-  dnds(0, 3) =  0.5 * c(1) * (2 * c(0) + 2 * c(1) - 1);
-  dnds(0, 4) =  0.5 * c(2) * (2 * c(0) + 2 * c(2) - 1);
-  dnds(0, 5) = -0.5 * (c(1) + c(2) - 1) * (2 * c(0) - 2 * c(1) - 2 * c(2) + 1);
-  dnds(0, 6) = -2.0 * c(1) * c(2);
-  dnds(0, 7) =  2.0 * c(2) * (c(1) + c(2) - 1);
-  dnds(0, 8) =  2.0 * c(1) * (c(1) + c(2) - 1);
-  dnds(0, 9) = -2.0 * c(0) * c(1);
-  dnds(0,10) = -2.0 * c(0) * c(2);
-  dnds(0,11) =  2.0 * c(0) * (c(1) + c(2) - 1);
-  dnds(0,12) =  2.0 * c(1) * c(2);
-  dnds(0,13) = -2.0 * c(2) * (c(1) + c(2) - 1);
-  dnds(0,14) = -2.0 * c(1) * (c(1) + c(2) - 1);
+inline void InterpolationElement<_itp_lagrange_pentahedron_15>::computeDNDS(
+    const vector_type & c, matrix_type & dnds) {
+  auto & x = c(0);
+  auto & y = c(1);
+  auto & z = c(2);
+
+  // ddx
+  dnds(0, 0) = 0.5 * y * (2 * x - 2 * y + 1);
+  dnds(0, 1) = 0.5 * z * (2 * x - 2 * z + 1);
+  dnds(0, 2) = -0.5 * (2 * x + 2 * y + 2 * z - 1) * (y + z - 1);
+  dnds(0, 3) = 0.5 * y * (2 * x + 2 * y - 1);
+  dnds(0, 4) = 0.5 * z * (2 * x + 2 * z - 1);
+  dnds(0, 5) = -0.5 * (y + z - 1) * (2 * x - 2 * y - 2 * z + 1);
+  dnds(0, 6) = -2.0 * y * z;
+  dnds(0, 7) = 2.0 * z * (y + z - 1);
+  dnds(0, 8) = 2.0 * y * (y + z - 1);
+  dnds(0, 9) = -2.0 * x * y;
+  dnds(0, 10) = -2.0 * x * z;
+  dnds(0, 11) = 2.0 * x * (y + z - 1);
+  dnds(0, 12) = 2.0 * y * z;
+  dnds(0, 13) = -2.0 * z * (y + z - 1);
+  dnds(0, 14) = -2.0 * y * (y + z - 1);
 
-  //ddy
-  dnds(1, 0) = -0.5 * (c(0) - 1) * (4 * c(1) - c(0) - 2);
-  dnds(1, 1) =  0.0;
-  dnds(1, 2) = -0.5 * (c(0) - 1) * (4 * c(1) + c(0) + 2 * (2 * c(2) - 1));
-  dnds(1, 3) =  0.5 * (c(0) + 1) * (4 * c(1) + c(0) - 2);
-  dnds(1, 4) =  0.0;
-  dnds(1, 5) =  0.5 * (c(0) + 1) * (4 * c(1) - c(0) + 2 * (2 * c(2) - 1));
-  dnds(1, 6) = -2.0 * (c(0) - 1) * c(2);
-  dnds(1, 7) =  2.0 * c(2) * (c(0) - 1);
-  dnds(1, 8) =  2.0 * (2 * c(1) + c(2) - 1) * (c(0) - 1);
-  dnds(1, 9) = -(c(0) * c(0) - 1);
-  dnds(1,10) =  0.0;
-  dnds(1,11) =  (c(0) * c(0) - 1);
-  dnds(1,12) =  2.0 * c(2) * (c(0) + 1);
-  dnds(1,13) = -2.0 * c(2) * (c(0) + 1);
-  dnds(1,14) = -2.0 * (2 * c(1) + c(2) - 1) * (c(0) + 1);
+  // ddy
+  dnds(1, 0) = -0.5 * (x - 1) * (4 * y - x - 2);
+  dnds(1, 1) = 0.0;
+  dnds(1, 2) = -0.5 * (x - 1) * (4 * y + x + 2 * (2 * z - 1));
+  dnds(1, 3) = 0.5 * (x + 1) * (4 * y + x - 2);
+  dnds(1, 4) = 0.0;
+  dnds(1, 5) = 0.5 * (x + 1) * (4 * y - x + 2 * (2 * z - 1));
+  dnds(1, 6) = -2.0 * (x - 1) * z;
+  dnds(1, 7) = 2.0 * z * (x - 1);
+  dnds(1, 8) = 2.0 * (2 * y + z - 1) * (x - 1);
+  dnds(1, 9) = -(x * x - 1);
+  dnds(1, 10) = 0.0;
+  dnds(1, 11) = (x * x - 1);
+  dnds(1, 12) = 2.0 * z * (x + 1);
+  dnds(1, 13) = -2.0 * z * (x + 1);
+  dnds(1, 14) = -2.0 * (2 * y + z - 1) * (x + 1);
 
-  //ddz
-  dnds(2, 0) =  0.0;
-  dnds(2, 1) = -0.5 * (c(0) - 1) * (4 * c(2) - c(0) - 2);
-  dnds(2, 2) = -0.5 * (c(0) - 1) * (4 * c(2) + c(0) + 2 * (2 * c(1) - 1));
-  dnds(2, 3) =  0.0;
-  dnds(2, 4) =  0.5 * (c(0) + 1) * (4 * c(2) + c(0) - 2);
-  dnds(2, 5) =  0.5 * (c(0) + 1) * (4 * c(2) - c(0) + 2 * (2 * c(1) - 1));
-  dnds(2, 6) = -2.0 * (c(0) - 1) * c(1);
-  dnds(2, 7) =  2.0 * (c(0) - 1) * (2 * c(2) + c(1) - 1);
-  dnds(2, 8) =  2.0 * c(1) * (c(0) - 1);
-  dnds(2, 9) =  0.0;
-  dnds(2,10) = -(c(0) * c(0) - 1);
-  dnds(2,11) =  (c(0) * c(0) - 1);
-  dnds(2,12) =  2.0 * (c(0) + 1) * c(1);
-  dnds(2,13) = -2.0 * (c(0) + 1) * (2 * c(2) + c(1) - 1);
-  dnds(2,14) = -2.0 * (c(0) + 1) * c(1);
+  // ddz
+  dnds(2, 0) = 0.0;
+  dnds(2, 1) = -0.5 * (x - 1) * (4 * z - x - 2);
+  dnds(2, 2) = -0.5 * (x - 1) * (4 * z + x + 2 * (2 * y - 1));
+  dnds(2, 3) = 0.0;
+  dnds(2, 4) = 0.5 * (x + 1) * (4 * z + x - 2);
+  dnds(2, 5) = 0.5 * (x + 1) * (4 * z - x + 2 * (2 * y - 1));
+  dnds(2, 6) = -2.0 * (x - 1) * y;
+  dnds(2, 7) = 2.0 * (x - 1) * (2 * z + y - 1);
+  dnds(2, 8) = 2.0 * y * (x - 1);
+  dnds(2, 9) = 0.0;
+  dnds(2, 10) = -(x * x - 1);
+  dnds(2, 11) = (x * x - 1);
+  dnds(2, 12) = 2.0 * (x + 1) * y;
+  dnds(2, 13) = -2.0 * (x + 1) * (2 * z + y - 1);
+  dnds(2, 14) = -2.0 * (x + 1) * y;
 }
 
 /* -------------------------------------------------------------------------- */
-template<>
-inline Real
-GeometricalElement<_gt_pentahedron_15>::getInradius(const Matrix<Real> & coord) {
+template <>
+inline Real GeometricalElement<_gt_pentahedron_15>::getInradius(
+    const Matrix<Real> & coord) {
   return GeometricalElement<_gt_pentahedron_6>::getInradius(coord);
 }
diff --git a/test/test_fe_engine/CMakeLists.txt b/test/test_fe_engine/CMakeLists.txt
index ff7cbf343..90475bf34 100644
--- a/test/test_fe_engine/CMakeLists.txt
+++ b/test/test_fe_engine/CMakeLists.txt
@@ -1,119 +1,121 @@
 #===============================================================================
 # @file   CMakeLists.txt
 #
 # @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
 #
 # @date creation: Fri Sep 03 2010
 # @date last modification: Mon Dec 07 2015
 #
 # @brief  configuration for FEM tests
 #
 # @section LICENSE
 #
 # Copyright (©)  2010-2012, 2014,  2015 EPFL  (Ecole Polytechnique  Fédérale de
 # Lausanne)  Laboratory (LSMS  -  Laboratoire de  Simulation  en Mécanique  des
 # Solides)
 #
 # Akantu is free  software: you can redistribute it and/or  modify it under the
 # terms  of the  GNU Lesser  General Public  License as  published by  the Free
 # Software Foundation, either version 3 of the License, or (at your option) any
 # later version.
 #
 # Akantu is  distributed in the  hope that it  will be useful, but  WITHOUT ANY
 # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
 # A  PARTICULAR PURPOSE. See  the GNU  Lesser General  Public License  for more
 # details.
 #
 # You should  have received  a copy  of the GNU  Lesser General  Public License
 # along with Akantu. If not, see <http://www.gnu.org/licenses/>.
 #
 # @section DESCRIPTION
 #
 #===============================================================================
 
 #===============================================================================
 function(register_fem_test operation type)
   set(_target test_${operation}${type})
 
   register_test(${_target}
     SOURCES test_${operation}.cc
     FILES_TO_COPY ${type}.msh
     COMPILE_OPTIONS TYPE=${type}
     PACKAGE core
     )
 endfunction()
 
 #===============================================================================
 package_get_element_types(core _types)
 
 foreach(_type ${_types})
   if(EXISTS ${CMAKE_CURRENT_SOURCE_DIR}/${_type}.msh)
     list(APPEND _meshes ${_type}.msh)
     #register_fem_test(fe_engine_precomputation ${_type})
   else()
     if(NOT ${_type} STREQUAL _point_1)
       message("The mesh ${_type}.msh is missing, the fe_engine test cannot be activated without it")
     endif()
   endif()
 endforeach()
 
 
 register_test(test_fe_engine_precomupation_bernoulli_beam_2 PACKAGE core SOURCES test_fe_engine_precomputation_bernoulli_2.cc)
 #add_mesh(test_fem_circle_1_mesh circle.geo 2 1 OUTPUT circle1.msh)
 #add_mesh(test_fem_circle_2_mesh circle.geo 2 2 OUTPUT circle2.msh)
 
 # Tests for class MeshData
 macro(register_typed_test test_name type value1 value2)
   set(target test_${test_name}_${type})
   register_test(${target}
     SOURCES test_${test_name}.cc
     COMPILE_OPTIONS "TYPE=${type};VALUE1=${value1};VALUE2=${value2}"
     PACKAGE core
     )
 endmacro()
 
 register_typed_test(mesh_data string \"5\" \"10\")
 register_typed_test(mesh_data UInt 5 10)
 
 add_mesh(test_boundary_msh cube.geo 3 1)
 add_mesh(test_boundary_msh_physical_names cube_physical_names.geo 3 1)
 
 register_test(test_mesh_boundary
   SOURCES test_mesh_boundary.cc
   DEPENDS test_boundary_msh test_boundary_msh_physical_names
   PACKAGE core)
 
 register_test(test_facet_element_mapping
   SOURCES test_facet_element_mapping.cc
   DEPENDS test_boundary_msh_physical_names
   PACKAGE core)
 
+pybind11_add_module(aka_test MODULE pybind11_akantu.cc)
+
 register_gtest_sources(
   SOURCES test_fe_engine_precomputation.cc
   PACKAGE core pybind11
   DEPENDS aka_test
   )
 
 
 register_gtest_sources(
   SOURCES test_fe_engine_gauss_integration.cc
   PACKAGE core
   )
 
 register_gtest_sources(
   SOURCES test_gradient.cc
   PACKAGE core
   )
 
 register_gtest_sources(
   SOURCES test_integrate.cc
   PACKAGE core
   )
 
 register_gtest_sources(
   SOURCES test_inverse_map.cc
   PACKAGE core
   )
 
 register_gtest_test(test_fe_engine
   FILES_TO_COPY ${_meshes})
diff --git a/test/test_fe_engine/py_engine/py_engine.py b/test/test_fe_engine/py_engine/py_engine.py
index 5ecb388df..943fd2353 100644
--- a/test/test_fe_engine/py_engine/py_engine.py
+++ b/test/test_fe_engine/py_engine/py_engine.py
@@ -1,188 +1,344 @@
 #!/usr/bin/env python3
 
 __all__ = ['Shapes']
 
 import numpy as np
 import numpy.polynomial.polynomial as poly
 import aka_test
 
 class Shapes(object):
     NATURAL_COORDS = {
-        (1, 'quadrangle'): np.array([[-1.], [1.], [0.]]),
-        (2, 'quadrangle'): np.array([[-1., -1.], [ 1., -1.], [ 1.,  1.], [-1.,  1.],
-                                     [ 0., -1.], [ 1.,  0.], [ 0.,  1.], [-1.,  0.]]),
-        (3, 'quadrangle'): np.array([[-1., -1., -1.], [ 1., -1., -1.], [ 1.,  1., -1.], [-1.,  1., -1.],
-                                     [-1., -1.,  1.], [ 1., -1.,  1.], [ 1.,  1.,  1.], [-1.,  1.,  1.],
-                                     [ 0., -1., -1.], [ 1.,  0., -1.], [ 0.,  1., -1.], [-1.,  0., -1.],
-                                     [-1., -1.,  0.], [ 1., -1.,  0.], [ 1.,  1.,  0.], [-1.,  1.,  0.],
-                                     [ 0., -1.,  1.], [ 1.,  0.,  1.], [ 0.,  1.,  1.], [-1.,  0.,  1.]]),
-        (2, 'triangle'):   np.array([[0., 0.], [1., 0.], [0., 1,], [.5, 0.], [.5, .5], [0., .5]]),
-        (3, 'triangle'):   np.array([[0., 0., 0.], [1., 0., 0.], [0., 1., 0.], [0., 0., 1.],
-                                     [.5, 0., 0.], [.5, .5, 0.], [0., .5, 0.],
-                                     [0., 0., .5], [.5, 0., .5], [0., .5, .5]]),
+        (1, 'quadrangle'):  np.array([[-1.], [1.], [0.]]),
+        (2, 'quadrangle'):  np.array([[-1., -1.], [ 1., -1.], [ 1.,  1.], [-1.,  1.],
+                                      [ 0., -1.], [ 1.,  0.], [ 0.,  1.], [-1.,  0.]]),
+        (3, 'quadrangle'):  np.array([[-1., -1., -1.], [ 1., -1., -1.], [ 1.,  1., -1.], [-1.,  1., -1.],
+                                      [-1., -1.,  1.], [ 1., -1.,  1.], [ 1.,  1.,  1.], [-1.,  1.,  1.],
+                                      [ 0., -1., -1.], [ 1.,  0., -1.], [ 0.,  1., -1.], [-1.,  0., -1.],
+                                      [-1., -1.,  0.], [ 1., -1.,  0.], [ 1.,  1.,  0.], [-1.,  1.,  0.],
+                                      [ 0., -1.,  1.], [ 1.,  0.,  1.], [ 0.,  1.,  1.], [-1.,  0.,  1.]]),
+        (2, 'triangle'):    np.array([[0., 0.], [1., 0.], [0., 1,], [.5, 0.], [.5, .5], [0., .5]]),
+        (3, 'triangle'):    np.array([[0., 0., 0.], [1., 0., 0.], [0., 1., 0.], [0., 0., 1.],
+                                      [.5, 0., 0.], [.5, .5, 0.], [0., .5, 0.],
+                                      [0., 0., .5], [.5, 0., .5], [0., .5, .5]]),
+        (3, 'pentahedron'): np.array([[-1.,  1.,  0.], [-1.,  0.,  1.], [-1.,  0.,  0.],
+                                      [ 1.,  1.,  0.], [ 1.,  0.,  1.], [ 1.,  0.,  0.],
+                                      [-1.,  .5,  .5], [-1.,  0.,  .5], [-1.,  .5,  0.],
+                                      [ 0.,  1.,  0.], [ 0.,  0.,  1.], [ 0.,  0.,  0.],
+                                      [ 1.,  .5,  .5], [ 1.,  0.,  .5], [ 1.,  .5,  0.],
+                                      [ 0.,  .5,  .5], [ 0.,  0.,  .5], [ 0.,  .5,  0.]]),
     }
 
     QUADRATURE_W = {
         (1, 'quadrangle', 1): np.array([2.]),
         (1, 'quadrangle', 2): np.array([1., 1.]),
         (2, 'triangle', 1): np.array([1./2.]),
         (2, 'triangle', 2): np.array([1., 1., 1.])/6.,
-        (3, 'triangle', 1): np.array([]),
-        (3, 'triangle', 2): np.array([]),
+        (3, 'triangle', 1): np.array([1./6.]),
+        (3, 'triangle', 2): np.array([1., 1., 1., 1.])/24.,
         (2, 'quadrangle', 1): np.array([1., 1., 1., 1.]),
-        (2, 'quadrangle', 2): np.array([]),
+        (2, 'quadrangle', 2): np.array([1., 1., 1., 1.]),
         (3, 'quadrangle', 1): np.array([1., 1., 1., 1.,
                                         1., 1., 1., 1.]),
-        (3, 'quadrangle', 2): np.array([]),
-
+        (3, 'quadrangle', 2): np.array([1., 1., 1., 1.,
+                                        1., 1., 1., 1.]),
+        (3, 'pentahedron', 1): np.array([1., 1., 1., 1., 1., 1.])/6.,
+        (3, 'pentahedron', 2): np.array([1., 1., 1., 1., 1., 1.])/6.,
     }
+
+    _tet_a = (5. - np.sqrt(5.))/20.
+    _tet_b = (5. + 3.*np.sqrt(5.))/20.
     QUADRATURE_G = {
         (1, 'quadrangle', 1): np.array([[0.]]),
         (1, 'quadrangle', 2): np.array([[-1.], [1.]])/np.sqrt(3.),
         (2, 'triangle', 1): np.array([[1., 1.]])/3.,
         (2, 'triangle', 2): np.array([[1./6., 1./6.], [2./3, 1./6], [1./6., 2./3.]]),
-        (2, 'quadrangle', 1): np.array([[-1., -1.], [-1.,  1.],
-                                        [ 1.,  1.], [ 1., -1.]])/np.sqrt(3.),
-        (2, 'quadrangle', 2): np.array([]),
-        (3, 'quadrangle', 1): np.array([[-1., -1., -1.], [-1.,  1., -1.],
-                                        [ 1.,  1., -1.], [ 1., -1., -1.],
-                                        [-1., -1.,  1.], [-1.,  1.,  1.],
-                                        [ 1.,  1.,  1.], [ 1., -1.,  1.]])/np.sqrt(3.),
-        (3, 'quadrangle', 2): np.array([]),
+        (3, 'triangle', 1): np.array([[1., 1., 1.]])/4.,
+        (3, 'triangle', 2): np.array([[_tet_a, _tet_a, _tet_a],
+                                      [_tet_b, _tet_a, _tet_a],
+                                      [_tet_a, _tet_b, _tet_a],
+                                      [_tet_a, _tet_a, _tet_b]]),
+        (2, 'quadrangle', 1): np.array([[-1., -1.], [ 1., -1.],
+                                        [-1.,  1.], [ 1.,  1.]])/np.sqrt(3.),
+        (2, 'quadrangle', 2): np.array([[-1., -1.], [ 1., -1.],
+                                        [-1.,  1.], [ 1.,  1.]])/np.sqrt(3.),
+        (3, 'quadrangle', 1): np.array([[-1., -1., -1.],
+                                        [ 1., -1., -1.],
+                                        [-1.,  1., -1.],
+                                        [ 1.,  1., -1.],
+                                        [-1., -1.,  1.],
+                                        [ 1., -1.,  1.],
+                                        [-1.,  1.,  1.],
+                                        [ 1.,  1.,  1.]])/np.sqrt(3.),
+        (3, 'quadrangle', 2): np.array([[-1., -1., -1.],
+                                        [ 1., -1., -1.],
+                                        [-1.,  1., -1.],
+                                        [ 1.,  1., -1.],
+                                        [-1., -1.,  1.],
+                                        [ 1., -1.,  1.],
+                                        [-1.,  1.,  1.],
+                                        [ 1.,  1.,  1.]])/np.sqrt(3.),
+ #       (3, 'pentahedron', 1): np.array([[-1./np.sqrt(3.), 0.5, 0.5],
+ #                                        [-1./np.sqrt(3.), 0. , 0.5],
+ #                                        [-1./np.sqrt(3.), 0.5, 0. ],
+ #                                        [ 1./np.sqrt(3.), 0.5, 0.5],
+ #                                        [ 1./np.sqrt(3.), 0. , 0.5],
+ #                                        [ 1./np.sqrt(3.), 0.5 ,0. ]]),
+        (3, 'pentahedron', 1): np.array([[-1./np.sqrt(3.), 1./6., 1./6.],
+                                         [-1./np.sqrt(3.), 2./3., 1./6.],
+                                         [-1./np.sqrt(3.), 1./6., 2./3.],
+                                         [ 1./np.sqrt(3.), 1./6., 1./6.],
+                                         [ 1./np.sqrt(3.), 2./3., 1./6.],
+                                         [ 1./np.sqrt(3.), 1./6., 2./3.]]),
+        (3, 'pentahedron', 2): np.array([[-1./np.sqrt(3.), 1./6., 1./6.],
+                                         [-1./np.sqrt(3.), 2./3., 1./6.],
+                                         [-1./np.sqrt(3.), 1./6., 2./3.],
+                                         [ 1./np.sqrt(3.), 1./6., 1./6.],
+                                         [ 1./np.sqrt(3.), 2./3., 1./6.],
+                                         [ 1./np.sqrt(3.), 1./6., 2./3.]]),
     }
 
     ELEMENT_TYPES = {
         '_segment_2':      ('quadrangle', 1, 'lagrange', 1, 2),
         '_segment_3':      ('quadrangle', 2, 'lagrange', 1, 3),
         '_triangle_3':     ('triangle', 1, 'lagrange', 2, 3),
         '_triangle_6':     ('triangle', 2, 'lagrange', 2, 6),
         '_quadrangle_4':   ('quadrangle', 1, 'serendip', 2, 4),
         '_quadrangle_8':   ('quadrangle', 2, 'serendip', 2, 8),
         '_tetrahedron_4':  ('triangle', 1, 'lagrange', 3, 4),
         '_tetrahedron_10': ('triangle', 2, 'lagrange', 3, 10),
         '_pentahedron_6':  ('pentahedron', 1, 'lagrange', 3, 6),
         '_pentahedron_15': ('pentahedron', 2, 'lagrange', 3, 15),
         '_hexahedron_8':   ('quadrangle', 1, 'serendip', 3, 8),
         '_hexahedron_20':  ('quadrangle', 2, 'serendip', 3, 20),
     }
 
     MONOMES = {(1, 'quadrangle'): np.array([[0], [1], [2], [3], [4], [5]]),
                (2, 'triangle'): np.array([[0, 0],                     #        1
                                           [1, 0], [0, 1],             #     x     y
                                           [2, 0], [1, 1], [0, 2]]),   # x^2   x.y   y^2
                (2, 'quadrangle'): np.array([[0, 0],
                                             [1, 0], [1, 1], [0, 1],
                                             [2, 0], [2, 1], [1, 2], [0, 2]]),
                (3, 'triangle'): np.array([[0, 0, 0],
                                           [1, 0, 0], [0, 1, 0], [0, 0, 1],
                                           [2, 0, 0], [1, 1, 0], [0, 2, 0], [0, 1, 1], [0, 0, 2],
                                           [1, 0, 1]]),
                (3, 'quadrangle'): np.array([[0, 0, 0],
                                             [1, 0, 0], [0, 1, 0], [0, 0, 1],
-                                            [1, 1, 0], [0, 1, 1], [1, 0, 1], [1, 1, 1],
+                                            [1, 1, 0], [1, 0, 1],
+                                            [0, 1, 1], [1, 1, 1],
                                             [2, 0, 0], [0, 2, 0], [0, 0, 2],
                                             [2, 1, 0], [2, 0, 1], [2, 1, 1],
                                             [1, 2, 0], [0, 2, 1], [1, 2, 1],
                                             [1, 0, 2], [0, 1, 2], [1, 1, 2]]),
     }
 
+    SHAPES = {
+        (3, 'pentahedron', 1): np.array([
+            [[[ 0.,  0.], [ 1.,  0.]], [[ 0.,  0.], [-1.,  0.]]],
+            [[[ 0.,  1.], [ 0.,  0.]], [[ 0., -1.], [ 0.,  0.]]],
+            [[[ 1., -1.], [-1.,  0.]], [[-1.,  1.], [ 1.,  0.]]],
+            [[[ 0.,  0.], [ 1.,  0.]], [[ 0.,  0.], [ 1.,  0.]]],
+            [[[ 0.,  1.], [ 0.,  0.]], [[ 0.,  1.], [ 0.,  0.]]],
+            [[[ 1., -1.], [-1.,  0.]], [[ 1., -1.], [-1.,  0.]]]
+        ])/2.,
+        (3, 'pentahedron', 2): np.array([
+            # 0
+            [[[ 0. ,  0. ,  0. ], [-1. ,  0. ,  0. ], [ 1. ,  0. ,  0. ]],
+             [[ 0. ,  0. ,  0. ], [ 0.5,  0. ,  0. ], [-1. ,  0. ,  0. ]],
+             [[ 0. ,  0. ,  0. ], [ 0.5,  0. ,  0. ], [ 0. ,  0. ,  0. ]]],
+
+            # 1
+            [[[ 0. , -1. ,  1. ], [ 0. ,  0. ,  0. ], [ 0. ,  0. ,  0. ]],
+             [[ 0. ,  0.5, -1. ], [ 0. ,  0. ,  0. ], [ 0. ,  0. ,  0. ]],
+             [[ 0. ,  0.5,  0. ], [ 0. ,  0. ,  0. ], [ 0. ,  0. ,  0. ]]],
+
+            # 2
+            [[[ 0. , -1. ,  1. ], [-1. ,  2. ,  0. ], [ 1. ,  0. ,  0. ]],
+             [[-0.5,  1.5, -1. ], [ 1.5, -2. ,  0. ], [-1. ,  0. ,  0. ]],
+             [[ 0.5, -0.5,  0. ], [-0.5,  0. ,  0. ], [ 0. ,  0. ,  0. ]]],
+
+            # 3
+            [[[ 0. ,  0. ,  0. ], [-1. ,  0. ,  0. ], [ 1. ,  0. ,  0. ]],
+             [[ 0. ,  0. ,  0. ], [-0.5,  0. ,  0. ], [ 1. ,  0. ,  0. ]],
+             [[ 0. ,  0. ,  0. ], [ 0.5,  0. ,  0. ], [ 0. ,  0. ,  0. ]]],
+
+            # 4
+            [[[ 0. , -1. ,  1. ], [ 0. ,  0. ,  0. ], [ 0. ,  0. ,  0. ]],
+             [[ 0. , -0.5,  1. ], [ 0. ,  0. ,  0. ], [ 0. ,  0. ,  0. ]],
+             [[ 0. ,  0.5,  0. ], [ 0. ,  0. ,  0. ], [ 0. ,  0. ,  0. ]]],
+
+            # 5
+            [[[ 0. , -1. ,  1. ], [-1. ,  2. ,  0. ], [ 1. ,  0. ,  0. ]],
+             [[ 0.5, -1.5,  1. ], [-1.5,  2. ,  0. ], [ 1. ,  0. ,  0. ]],
+             [[ 0.5, -0.5,  0. ], [-0.5,  0. ,  0. ], [ 0. ,  0. ,  0. ]]],
+
+            # 6
+            [[[ 0. ,  0. ,  0. ], [ 0. ,  2. ,  0. ], [ 0. ,  0. ,  0. ]],
+             [[ 0. ,  0. ,  0. ], [ 0. , -2. ,  0. ], [ 0. ,  0. ,  0. ]],
+             [[ 0. ,  0. ,  0. ], [ 0. ,  0. ,  0. ], [ 0. ,  0. ,  0. ]]],
+
+            # 7
+            [[[ 0. ,  2. , -2. ], [ 0. , -2. ,  0. ], [ 0. ,  0. ,  0. ]],
+             [[ 0. , -2. ,  2. ], [ 0. ,  2. ,  0. ], [ 0. ,  0. ,  0. ]],
+             [[ 0. ,  0. ,  0. ], [ 0. ,  0. ,  0. ], [ 0. ,  0. ,  0. ]]],
+
+            # 8
+            [[[ 0. ,  0. ,  0. ], [ 2. , -2. ,  0. ], [-2. ,  0. ,  0. ]],
+             [[ 0. ,  0. ,  0. ], [-2. ,  2. ,  0. ], [ 2. ,  0. ,  0. ]],
+             [[ 0. ,  0. ,  0. ], [ 0. ,  0. ,  0. ], [ 0. ,  0. ,  0. ]]],
+
+            # 9
+            [[[ 0. ,  0. ,  0. ], [ 1. ,  0. ,  0. ], [ 0. ,  0. ,  0. ]],
+             [[ 0. ,  0. ,  0. ], [ 0. ,  0. ,  0. ], [ 0. ,  0. ,  0. ]],
+             [[ 0. ,  0. ,  0. ], [-1. ,  0. ,  0. ], [ 0. ,  0. ,  0. ]]],
+
+            # 10
+            [[[ 0. ,  1. ,  0. ], [ 0. ,  0. ,  0. ], [ 0. ,  0. ,  0. ]],
+             [[ 0. ,  0. ,  0. ], [ 0. ,  0. ,  0. ], [ 0. ,  0. ,  0. ]],
+             [[ 0. , -1. ,  0. ], [ 0. ,  0. ,  0. ], [ 0. ,  0. ,  0. ]]],
+
+            # 11
+            [[[ 1. , -1. ,  0. ], [-1. ,  0. ,  0. ], [ 0. ,  0. ,  0. ]],
+             [[ 0. ,  0. ,  0. ], [ 0. ,  0. ,  0. ], [ 0. ,  0. ,  0. ]],
+             [[-1. ,  1. ,  0. ], [ 1. ,  0. ,  0. ], [ 0. ,  0. ,  0. ]]],
+
+            # 12
+            [[[ 0. ,  0. ,  0. ], [ 0. ,  2. ,  0. ], [ 0. ,  0. ,  0. ]],
+             [[ 0. ,  0. ,  0. ], [ 0. ,  2. ,  0. ], [ 0. ,  0. ,  0. ]],
+             [[ 0. ,  0. ,  0. ], [ 0. ,  0. ,  0. ], [ 0. ,  0. ,  0. ]]],
+
+            # 13
+            [[[ 0. ,  2. , -2. ], [ 0. , -2. ,  0. ], [ 0. ,  0. ,  0. ]],
+             [[ 0. ,  2. , -2. ], [ 0. , -2. ,  0. ], [ 0. ,  0. ,  0. ]],
+             [[ 0. ,  0. ,  0. ], [ 0. ,  0. ,  0. ], [ 0. ,  0. ,  0. ]]],
+
+            # 14
+            [[[ 0. ,  0. ,  0. ], [ 2. , -2. ,  0. ], [-2. ,  0. ,  0. ]],
+             [[ 0. ,  0. ,  0. ], [ 2. , -2. ,  0. ], [-2. ,  0. ,  0. ]],
+             [[ 0. ,  0. ,  0. ], [ 0. ,  0. ,  0. ], [ 0. ,  0. ,  0. ]]],
+        ])}
+
     def __init__(self,  element):
         self._shape, self._order, self._inter_poly, self._dim, self._nnodes = self.ELEMENT_TYPES[element]
         self._ksi = self.NATURAL_COORDS[(self._dim, self._shape)][:self._nnodes]
         self._g = self.QUADRATURE_G[(self._dim, self._shape, self._order)]
         self._w = self.QUADRATURE_W[(self._dim, self._shape, self._order)]
+
+    def polyval(self, x, p):
+        if 1 == self._dim:
+            return poly.polyval(x[0], p)
+        if 2 == self._dim:
+            return poly.polyval2d(x[0], x[1], p)
+        if 3 == self._dim:
+            return poly.polyval3d(x[0], x[1], x[2], p)
+
+    def shape_from_monomes(self):
         momo = self.MONOMES[(self._dim, self._shape)][:self._nnodes]
 
         _shape = list(momo[0])
 
         for s in range(len(_shape)):
             _shape[s] = max(momo[:,s])+1
         self._poly_shape = tuple(_shape)
 
         self._monome = []
         for m in momo:
             p = np.zeros(self._poly_shape)
             p[tuple(m)] = 1
             self._monome.append(p)
 
-    def polyval(self, x, p):
-        if 1 == self._dim:
-            return poly.polyval(x[0], p)
-        if 2 == self._dim:
-            return poly.polyval2d(x[0], x[1], p)
-        if 3 == self._dim:
-            return poly.polyval3d(x[0], x[1], x[2], p)
-
-    def precompute(self, **kwargs):
-        X = np.array(kwargs["X"],  copy=False)
-        nb_element = X.shape[0]
-        X = X.reshape(nb_element, self._nnodes, self._dim)
-
         # evaluate polynomial constant for shapes
         _x = self._ksi
         _xe = np.zeros((self._nnodes, self._nnodes))
         for n in range(self._nnodes):
             _xe[:,n] = [self.polyval(_x[n], m) for m in self._monome]
 
         _a = np.linalg.inv(_xe)
         _n = np.zeros((self._nnodes,) + self._monome[0].shape)
         # set shapes polynomials
         for n in range(self._nnodes):
             for m in range(len(self._monome)):
                 _n[n] += _a[n, m] * self._monome[m]
 
+        return _n
+
+    def compute_shapes(self):
+        if (self._dim, self._shape) in self.MONOMES:
+            return self.shape_from_monomes()
+        else:
+            _n = self.SHAPES[(self._dim, self._shape, self._order)]
+            self._poly_shape = _n[0].shape
+            return _n
+
+    def precompute(self, **kwargs):
+        X = np.array(kwargs["X"],  copy=False)
+        nb_element = X.shape[0]
+        X = X.reshape(nb_element, self._nnodes, self._dim)
+
+        _x = self._ksi
+        _n = self.compute_shapes()
+
         # sanity check on shapes
         for n in range(self._nnodes):
             for m in range(self._nnodes):
                 v = self.polyval(_x[n], _n[m])
                 ve = 1. if n == m else 0.
                 test = np.isclose(v, ve)
                 if not test:
-                    raise("Most probably an error in the shapes evaluation")
+                    raise Exception("Most probably an error in the shapes evaluation")
 
         # compute shapes derivatives
         _b = np.zeros((self._dim, self._nnodes,) + self._poly_shape)
         for d in range(self._dim):
             for n in range(self._nnodes):
                 _der = poly.polyder(_n[n], axis=d)
                 _mshape = np.array(self._poly_shape)
                 _mshape[d] = _mshape[d] - _der.shape[d]
                 _mshape = tuple(_mshape)
                 _comp = np.zeros(_mshape)
 
                 if 1 == self._dim:
                     _bt = np.hstack((_der, _comp))
                 else:
                     if 0 == d:
                         _bt = np.vstack((_der, _comp))
                     if 1 == d:
                         _bt = np.hstack((_der, _comp))
                     if 2 == d:
                         _bt = np.dstack((_der, _comp))
 
                 _b[d, n] = _bt
 
         _nb_quads = len(self._g)
         _nq = np.zeros((_nb_quads, self._nnodes))
         _bq = np.zeros((_nb_quads, self._dim, self._nnodes))
 
         # evaluate shapes and shapes derivatives on gauss points
         for q in range(_nb_quads):
             _g = self._g[q]
             for n in range(self._nnodes):
                 _nq[q, n] = self.polyval(_g, _n[n])
-
                 for d in range(self._dim):
                     _bq[q, d, n] = self.polyval(_g, _b[d, n])
 
         _j = np.array(kwargs['j'], copy=False).reshape((nb_element, _nb_quads))
-        _B = np.array(kwargs['B'], copy=False).reshape((nb_element, _nb_quads, self._dim, self._nnodes))
+        _B = np.array(kwargs['B'], copy=False).reshape((nb_element, _nb_quads, self._nnodes, self._dim))
         _N = np.array(kwargs['N'], copy=False).reshape((nb_element, _nb_quads, self._nnodes))
+        _Q = np.array(kwargs['Q'], copy=False)
+        if np.linalg.norm(_Q - self._g.T) > 1e-15:
+            raise Exception('Not using the same quadrature points')
 
         for e in range(nb_element):
             for q in range(_nb_quads):
                 _J = np.matmul(_bq[q], X[e])
 
+                if(np.linalg.norm(_N[e, q] - _nq[q]) > 1e-10):
+                    print(f"{e},{q}")
+                    print(_N[e, q])
+                    print(_nq[q])
                 _N[e, q] = _nq[q]
-                _B[e, q] = np.matmul(np.linalg.inv(_J), _bq[q])
+                _tmp = np.matmul(np.linalg.inv(_J), _bq[q])
+                _B[e, q] = _tmp.T
+
                 _j[e, q] = np.linalg.det(_J) * self._w[q]
diff --git a/test/test_fe_engine/pybind11_akantu.cc b/test/test_fe_engine/pybind11_akantu.cc
index 1c86b55bf..594469b10 100644
--- a/test/test_fe_engine/pybind11_akantu.cc
+++ b/test/test_fe_engine/pybind11_akantu.cc
@@ -1,44 +1,72 @@
 /* -------------------------------------------------------------------------- */
 #include "pybind11_akantu.hh"
 /* -------------------------------------------------------------------------- */
 #include <pybind11/numpy.h>
 #include <pybind11/pybind11.h>
 /* -------------------------------------------------------------------------- */
 
 namespace py = pybind11;
 
 namespace akantu {
 
 PYBIND11_MODULE(aka_test, m) {
   m.doc() = "Module for the tests of akantu";
 
   py::class_<ArrayProxy<Real>>(m, "ArrayProxy", py::buffer_protocol())
       .def_buffer([](ArrayProxy<Real> & a) {
         return py::buffer_info(
             a.storage(),                           /* Pointer to buffer */
             sizeof(Real),                          /* Size of one scalar */
             py::format_descriptor<Real>::format(), /* Python struct-style
                                                        format descriptor */
             2,                                     /* Number of dimensions */
             {a.size(), a.getNbComponent()},        /* Buffer dimensions */
             {sizeof(Real) *
                  a.getNbComponent(), /* Strides (in bytes) for each index */
              sizeof(Real)});
       })
       .def(py::init([](py::array & n) {
         /* Request a buffer descriptor from Python */
         py::buffer_info info = n.request();
 
         /* Some sanity checks ... */
         if (info.format != py::format_descriptor<Real>::format())
           throw std::runtime_error(
               "Incompatible format: expected a double array!");
 
         if (info.ndim != 2)
           throw std::runtime_error("Incompatible buffer dimension!");
 
         return std::make_unique<ArrayProxy<Real>>(static_cast<Real *>(info.ptr),
                                                   info.shape[0], info.shape[1]);
       }));
-}
+
+  py::class_<MatrixProxy<Real>>(m, "Matrix", py::buffer_protocol())
+      .def_buffer([](MatrixProxy<Real> & a) {
+        return py::buffer_info(
+            a.storage(),                           /* Pointer to buffer */
+            sizeof(Real),                          /* Size of one scalar */
+            py::format_descriptor<Real>::format(), /* Python struct-style
+                                                       format descriptor */
+            2,                                     /* Number of dimensions */
+            {a.size(0), a.size(1)},        /* Buffer dimensions */
+            {sizeof(Real), a.size(0)*sizeof(Real)}  /* Strides (in bytes) for each index */
+          );
+      })
+      .def(py::init([](py::array & n) {
+        /* Request a buffer descriptor from Python */
+        py::buffer_info info = n.request();
+
+        /* Some sanity checks ... */
+        if (info.format != py::format_descriptor<Real>::format())
+          throw std::runtime_error(
+              "Incompatible format: expected a double array!");
+
+        if (info.ndim != 2)
+          throw std::runtime_error("Incompatible buffer dimension!");
+
+        return std::make_unique<MatrixProxy<Real>>(static_cast<Real *>(info.ptr),
+                                                   info.shape[0], info.shape[1]);
+      }));
+} // Module aka test
 } // namespace akantu
diff --git a/test/test_fe_engine/pybind11_akantu.hh b/test/test_fe_engine/pybind11_akantu.hh
index 080dc21f4..98bacea8a 100644
--- a/test/test_fe_engine/pybind11_akantu.hh
+++ b/test/test_fe_engine/pybind11_akantu.hh
@@ -1,104 +1,111 @@
 #include "aka_array.hh"
+#include "aka_types.hh"
 /* -------------------------------------------------------------------------- */
 #include <pybind11/cast.h>
 #include <pybind11/numpy.h>
 #include <pybind11/pybind11.h>
 /* -------------------------------------------------------------------------- */
 
 #ifndef __AKANTU_PYBIND11_AKANTU_HH__
 #define __AKANTU_PYBIND11_AKANTU_HH__
 
 namespace akantu {
 
 template <typename T> class ArrayProxy : public Array<T> {
 public:
   ArrayProxy(T * wrapped_memory, UInt size = 0, UInt nb_component = 1,
              const ID & id = "")
       : Array<T>(0, nb_component, id) {
     this->values = wrapped_memory;
     this->size_ = size;
   }
 
-  ArrayProxy(const ArrayProxy<T> & array) : Array<T>(0, array.nb_component, array.id){
+  ArrayProxy(const ArrayProxy<T> & array)
+      : Array<T>(0, array.nb_component, array.id) {
     this->values = array.values;
     this->size_ = array.size_;
   }
 
   ArrayProxy(ArrayProxy<T> && array) {
     this->nb_component = std::move(array.nb_component);
     this->values = std::move(array.values);
     this->size_ = std::move(array.size_);
     this->id = std::move(array.id);
   }
 
   ArrayProxy(Array<T> & array)
       : Array<T>(0, array.getNbComponent(), array.getID()) {
     this->values = array.storage();
     this->size_ = array.size();
   }
 
   ~ArrayProxy() {
     this->values = nullptr;
     this->size_ = 0;
   }
 
   void setNbComponent(UInt nb_component) {
     UInt new_size = this->size_ / nb_component;
     AKANTU_DEBUG_ASSERT(
         nb_component * new_size == this->nb_component * this->size_,
         nb_component
             << " is not valid as a new number of component for this array");
 
     this->nb_component = nb_component;
     this->size_ = new_size;
   }
 
   void resize(UInt new_size) {
     AKANTU_DEBUG_ASSERT(this->size_ == new_size,
                         "cannot resize a temporary vector");
   }
 };
 
 template <typename T> decltype(auto) make_proxy(Array<T> & array) {
   return ArrayProxy<T>(array);
 }
 
+template <typename T>
+decltype(auto) make_proxy(const Matrix<T> & array) {
+  return MatrixProxy<T>(array);
+}
+
 } // namespace akantu
 
 // namespace pybind11 {
 // namespace detail {
 //   template <> struct type_caster<akantu::ArrayProxy<akantu::Real>> {
 //   public:
 //     PYBIND11_TYPE_CASTER(akantu::ArrayProxy<akantu::Real>, _("ArrayProxy"));
 //     bool load(handle, bool) {
 //       // /* Extract PyObject from handle */
 //       // PyObject *source = src.ptr();
 //       // /* Try converting into a Python integer value */
 //       // PyObject *tmp = PyNumber_Long(source);
 //       // if (!tmp)
 //       //   return false;
 //       // /* Now try to convert into a C++ int */
 //       // value.long_value = PyLong_AsLong(tmp);
 //       // Py_DECREF(tmp);
 //       // /* Ensure return code was OK (to avoid out-of-range errors etc) */
 //       // return !(value.long_value == -1 && !PyErr_Occurred());
 //       return false;
 //     }
 
 //     static handle cast(akantu::ArrayProxy<akantu::Real> & src,
 //                        return_value_policy /* policy */, handle parent) {
 //       constexpr ssize_t elem_size = sizeof(akantu::Real);
 //       ssize_t nb_comp = src.getNbComponent();
 //       ssize_t size = src.size();
 //       std::cout << "<memory at " << reinterpret_cast<void *>(src.storage())
 //                 << ">\n";
 
 //       auto a = array({size, nb_comp}, {elem_size * nb_comp, elem_size},
 //                      src.storage(), handle());
 //       return a.release();
 //     }
 //   };
 // } // namespace detail
 // } //  namespace pybind11
 
 #endif /* __AKANTU_PYBIND11_AKANTU_HH__ */
diff --git a/test/test_fe_engine/test_fe_engine_fixture.hh b/test/test_fe_engine/test_fe_engine_fixture.hh
index 443897734..a637f8bbf 100644
--- a/test/test_fe_engine/test_fe_engine_fixture.hh
+++ b/test/test_fe_engine/test_fe_engine_fixture.hh
@@ -1,71 +1,69 @@
 /* -------------------------------------------------------------------------- */
 #include "element_class.hh"
 #include "fe_engine.hh"
 #include "integrator_gauss.hh"
 #include "shape_lagrange.hh"
 #include "test_gtest_utils.hh"
 /* -------------------------------------------------------------------------- */
 #include <gtest/gtest.h>
 /* -------------------------------------------------------------------------- */
 
 #ifndef __AKANTU_TEST_FE_ENGINE_FIXTURE_HH__
 #define __AKANTU_TEST_FE_ENGINE_FIXTURE_HH__
 
 using namespace akantu;
 
 template <typename type_> class TestFEMFixture : public ::testing::Test {
 public:
   static constexpr const ElementType type = type_::value;
   static constexpr const size_t dim = ElementClass<type>::getSpatialDimension();
   using FEM = FEEngineTemplate<IntegratorGauss, ShapeLagrange>;
 
   void SetUp() override {
     const auto dim = this->dim;
     const auto type = this->type;
     mesh = std::make_unique<Mesh>(dim);
 
     std::stringstream meshfilename;
     meshfilename << type << ".msh";
     mesh->read(meshfilename.str());
 
     lower = mesh->getLowerBounds();
     upper = mesh->getUpperBounds();
 
     nb_element = this->mesh->getNbElement(type);
 
     fem = std::make_unique<FEM>(*mesh, dim, "my_fem");
     fem->initShapeFunctions();
 
     nb_quadrature_points_total = fem->getNbIntegrationPoints(type) * nb_element;
 
-    std::stringstream sstr;
-    sstr << type;
-    SCOPED_TRACE(sstr.str().c_str());
+    SCOPED_TRACE(aka::to_string(type));
   }
 
   void TearDown() override {
     fem.reset(nullptr);
     mesh.reset(nullptr);
   }
 
 protected:
   std::unique_ptr<FEM> fem;
   std::unique_ptr<Mesh> mesh;
   UInt nb_element;
   UInt nb_quadrature_points_total;
   Vector<Real> lower;
   Vector<Real> upper;
 };
 
 template <typename type_>
 constexpr const ElementType TestFEMFixture<type_>::type;
 
 template <typename type_>
 constexpr const size_t TestFEMFixture<type_>::dim;
 
 
 using types = gtest_list_t<TestElementTypes>;
 
 TYPED_TEST_CASE(TestFEMFixture, types);
 
 #endif /* __AKANTU_TEST_FE_ENGINE_FIXTURE_HH__ */
diff --git a/test/test_fe_engine/test_fe_engine_precomputation.cc b/test/test_fe_engine/test_fe_engine_precomputation.cc
index c9efe534b..c66054cc6 100644
--- a/test/test_fe_engine/test_fe_engine_precomputation.cc
+++ b/test/test_fe_engine/test_fe_engine_precomputation.cc
@@ -1,113 +1,112 @@
 /**
  * @file   test_fe_engine_precomputation.cc
  *
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Mon Jun 14 2010
  * @date last modification: Mon Jul 13 2015
  *
  * @brief  test of the fem class
  *
  * @section LICENSE
  *
  * Copyright (©)  2010-2012, 2014,  2015 EPFL  (Ecole Polytechnique  Fédérale de
  * Lausanne)  Laboratory (LSMS  -  Laboratoire de  Simulation  en Mécanique  des
  * Solides)
  *
  * Akantu is free  software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as  published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Akantu is  distributed in the  hope that it  will be useful, but  WITHOUT ANY
  * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
  * A  PARTICULAR PURPOSE. See  the GNU  Lesser General  Public License  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 #include "pybind11_akantu.hh"
 #include "test_fe_engine_fixture.hh"
 /* -------------------------------------------------------------------------- */
 #include <pybind11/embed.h>
 #include <pybind11/numpy.h>
 /* -------------------------------------------------------------------------- */
 using namespace akantu;
 
 namespace py = pybind11;
 using namespace py::literals;
 
 template <typename type_>
 class TestFEMPyFixture : public TestFEMFixture<type_> {
   using parent = TestFEMFixture<type_>;
 
 public:
   void SetUp() override {
     parent::SetUp();
 
     const auto & connectivities = this->mesh->getConnectivity(this->type);
     const auto & nodes = this->mesh->getNodes().begin(this->dim);
     coordinates = std::make_unique<Array<Real>>(
         connectivities.size(), connectivities.getNbComponent() * this->dim);
 
     for (auto && tuple :
          zip(make_view(connectivities, connectivities.getNbComponent()),
              make_view(*coordinates, this->dim,
                        connectivities.getNbComponent()))) {
       const auto & conn = std::get<0>(tuple);
       const auto & X = std::get<1>(tuple);
       for (auto s : arange(conn.size())) {
         Vector<Real>(X(s)) = Vector<Real>(nodes[conn(s)]);
       }
     }
   }
 
   void TearDown() override {
     parent::TearDown();
     coordinates.reset(nullptr);
   }
 
 protected:
   std::unique_ptr<Array<Real>> coordinates;
 };
 
 TYPED_TEST_CASE(TestFEMPyFixture, types);
 
-TYPED_TEST(TestFEMPyFixture, DISABLED_Precompute) {
+TYPED_TEST(TestFEMPyFixture, Precompute) {
+  SCOPED_TRACE(aka::to_string(this->type));
   const auto & N = this->fem->getShapeFunctions().getShapes(this->type);
   const auto & B =
       this->fem->getShapeFunctions().getShapesDerivatives(this->type);
   const auto & j = this->fem->getIntegrator().getJacobians(this->type);
 
-  Array<Real> ref_N(this->nb_quadrature_points_total, N.getNbComponent());
-  Array<Real> ref_B(this->nb_quadrature_points_total, B.getNbComponent());
+  // Array<Real> ref_N(this->nb_quadrature_points_total, N.getNbComponent());
+  // Array<Real> ref_B(this->nb_quadrature_points_total, B.getNbComponent());
   Array<Real> ref_j(this->nb_quadrature_points_total, j.getNbComponent());
-
+  auto ref_N(N);
+  auto ref_B(B);
   py::module py_engine = py::module::import("py_engine");
   auto py_shape = py_engine.attr("Shapes")(py::str(aka::to_string(this->type)));
-  auto kwargs = py::dict("N"_a = make_proxy(ref_N), "B"_a = make_proxy(ref_B),
-                         "j"_a = make_proxy(ref_j),
-                         "X"_a = make_proxy(*this->coordinates));
+  auto kwargs = py::dict(
+      "N"_a = make_proxy(ref_N), "B"_a = make_proxy(ref_B),
+      "j"_a = make_proxy(ref_j), "X"_a = make_proxy(*this->coordinates),
+      "Q"_a = make_proxy(
+          this->fem->getIntegrationPoints(this->type)));
 
   auto ret = py_shape.attr("precompute")(**kwargs);
-
-  auto check = [&](auto & ref_A, auto & A) {
-    for (auto && n :
-         zip(make_view(ref_A,
-                       this->nb_quadrature_points_total / this->nb_element,
-                       ref_A.getNbComponent()),
-             make_view(A, this->nb_quadrature_points_total / this->nb_element,
-                       A.getNbComponent()))) {
+  auto check = [&](auto & ref_A, auto & A, const auto & id) {
+    SCOPED_TRACE(aka::to_string(this->type) + " " + id);
+    for (auto && n : zip(make_view(ref_A, ref_A.getNbComponent()),
+                         make_view(A, A.getNbComponent()))) {
       auto diff = (std::get<0>(n) - std::get<1>(n)).template norm<L_inf>();
-      //std::cout << "||" << std::get<0>(n) << " - " << std::get<1>(n) << "|| = " << diff << std::endl;
-      EXPECT_NEAR(0., diff, 1e-14);
+      EXPECT_NEAR(0., diff, 1e-10);
     }
   };
-  check(ref_N, N);
-  check(ref_B, B);
-  check(ref_j, j);
+  check(ref_N, N, "N");
+  check(ref_B, B, "B");
+  check(ref_j, j, "j");
 }