diff --git a/src/common/aka_element_classes_info.hh.in b/src/common/aka_element_classes_info.hh.in
index 2b34cddb3..f3c201ba6 100644
--- a/src/common/aka_element_classes_info.hh.in
+++ b/src/common/aka_element_classes_info.hh.in
@@ -1,209 +1,141 @@
 /**
  * @file   aka_element_classes_info.hh.in
  *
  * @author Aurelia Isabel Cuba Ramos <aurelia.cubaramos@epfl.ch>
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Sun Jul 19 2015
  * @date last modification: Tue Feb 20 2018
  *
  * @brief  Declaration of the enums for the element classes
  *
  *
  * 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 "aka_safe_enum.hh"
 /* -------------------------------------------------------------------------- */
 #include <boost/preprocessor.hpp>
 /* -------------------------------------------------------------------------- */
 
 #ifndef AKANTU_AKA_ELEMENT_CLASSES_INFO_HH_
 #define AKANTU_AKA_ELEMENT_CLASSES_INFO_HH_
 
 namespace akantu {
 
 /* -------------------------------------------------------------------------- */
 /* Element Types                                                              */
 /* -------------------------------------------------------------------------- */
 
 /// @enum ElementType type of elements
 enum ElementType {
   _not_defined,
   @AKANTU_ELEMENT_TYPES_ENUM@
   _max_element_type
 };
 
 @AKANTU_ELEMENT_TYPES_BOOST_SEQ@
 
 @AKANTU_ALL_ELEMENT_BOOST_SEQ@
 
 /* -------------------------------------------------------------------------- */
 /* Element Kinds                                                              */
 /* -------------------------------------------------------------------------- */
 @AKANTU_ELEMENT_KINDS_BOOST_SEQ@
 
 @AKANTU_ELEMENT_KIND_BOOST_SEQ@
 
 enum ElementKind {
   BOOST_PP_SEQ_ENUM(AKANTU_ELEMENT_KIND),
   _ek_not_defined
 };
 
 
 /* -------------------------------------------------------------------------- */
 struct ElementKind_def {
   using type = ElementKind;
   static const type _begin_ = BOOST_PP_SEQ_HEAD(AKANTU_ELEMENT_KIND);
   static const type _end_   = _ek_not_defined;
 };
 
 using element_kind_t = safe_enum<ElementKind_def> ;
 
 /* -------------------------------------------------------------------------- */
 /// @enum GeometricalType type of element potentially contained in a Mesh
 enum GeometricalType {
   @AKANTU_GEOMETRICAL_TYPES_ENUM@
   _gt_not_defined
 };
 
 /* -------------------------------------------------------------------------- */
 /* Interpolation Types                                                        */
 /* -------------------------------------------------------------------------- */
 @AKANTU_INTERPOLATION_TYPES_BOOST_SEQ@
 
 /// @enum InterpolationType type of elements
 enum InterpolationType {
   BOOST_PP_SEQ_ENUM(AKANTU_INTERPOLATION_TYPES),
   _itp_not_defined
 };
 
 /* -------------------------------------------------------------------------- */
 /* Some sub types less probable to change                                     */
 /* -------------------------------------------------------------------------- */
 /// @enum GeometricalShapeType types of shapes to define the contains
 /// function in the element classes
 enum GeometricalShapeType {
   @AKANTU_GEOMETRICAL_SHAPES_ENUM@
   _gst_not_defined
 };
 
 /* -------------------------------------------------------------------------- */
 /// @enum GaussIntegrationType classes of types using common
 /// description of the gauss point position and weights
 enum GaussIntegrationType {
   @AKANTU_GAUSS_INTEGRATION_TYPES_ENUM@
   _git_not_defined
 };
 
 /* -------------------------------------------------------------------------- */
 /// @enum InterpolationKind the family of interpolation types
 enum InterpolationKind {
   @AKANTU_INTERPOLATION_KIND_ENUM@
   _itk_not_defined
 };
 
-/* -------------------------------------------------------------------------- */
-// BOOST PART: TOUCH ONLY IF YOU KNOW WHAT YOU ARE DOING
-#define AKANTU_BOOST_CASE_MACRO(r, macro, _type)                               \
-  case _type: {                                                                \
-    macro(_type);                                                              \
-    break;                                                                     \
-  }
-
-#define AKANTU_BOOST_LIST_SWITCH(macro1, list1, var)                           \
-  do {                                                                         \
-    switch (var) {                                                             \
-      BOOST_PP_SEQ_FOR_EACH(AKANTU_BOOST_CASE_MACRO, macro1, list1)            \
-    default: {                                                                 \
-      AKANTU_ERROR("Type ("                                                    \
-                   << var /* NOLINT */ << ") not handled by this function");   \
-    }                                                                          \
-    }                                                                          \
-  } while (0)
-
-#define AKANTU_BOOST_LIST_SWITCH_NO_DEFAULT(macro1, list1, var)                \
-  do {                                                                         \
-    switch (var) {                                                             \
-      BOOST_PP_SEQ_FOR_EACH(AKANTU_BOOST_CASE_MACRO, macro1, list1)            \
-    case _not_defined: /* FALLTHRU */                                          \
-    case _max_element_type:                                                    \
-      break;                                                                   \
-    }                                                                          \
-  } while (0)
-
-#define AKANTU_BOOST_ELEMENT_SWITCH(macro1, list1)                             \
-  AKANTU_BOOST_LIST_SWITCH(macro1, list1, type)
-
-#define AKANTU_BOOST_ELEMENT_SWITCH_NO_DEFAULT(macro1, list1)                  \
-  AKANTU_BOOST_LIST_SWITCH_NO_DEFAULT(macro1, list1, type)
-
-#define AKANTU_BOOST_ALL_ELEMENT_SWITCH(macro)                                 \
-  AKANTU_BOOST_ELEMENT_SWITCH(macro, AKANTU_ALL_ELEMENT_TYPE)
-
-#define AKANTU_BOOST_ALL_ELEMENT_SWITCH_NO_DEFAULT(macro)                      \
-  AKANTU_BOOST_ELEMENT_SWITCH_NO_DEFAULT(macro, AKANTU_ALL_ELEMENT_TYPE)
-
-#define AKANTU_BOOST_LIST_MACRO(r, macro, type) macro(type)
-
-#define AKANTU_BOOST_APPLY_ON_LIST(macro, list)                                \
-  BOOST_PP_SEQ_FOR_EACH(AKANTU_BOOST_LIST_MACRO, macro, list)
-
-#define AKANTU_BOOST_ALL_ELEMENT_LIST(macro)                                   \
-  AKANTU_BOOST_APPLY_ON_LIST(macro, AKANTU_ALL_ELEMENT_TYPE)
-
-#define AKANTU_GET_ELEMENT_LIST(kind) AKANTU##kind##_ELEMENT_TYPE
-
-#define AKANTU_BOOST_KIND_ELEMENT_SWITCH(macro, kind)                          \
-  AKANTU_BOOST_ELEMENT_SWITCH(macro, AKANTU_GET_ELEMENT_LIST(kind))
-
-// BOOST_PP_SEQ_TO_LIST does not exists in Boost < 1.49
-#define AKANTU_GENERATE_KIND_LIST(seq)                                         \
-  BOOST_PP_TUPLE_TO_LIST(BOOST_PP_SEQ_SIZE(seq), BOOST_PP_SEQ_TO_TUPLE(seq))
-
-#define AKANTU_ELEMENT_KIND_BOOST_LIST                                         \
-  AKANTU_GENERATE_KIND_LIST(AKANTU_ELEMENT_KIND)
-
-#define AKANTU_BOOST_ALL_KIND_LIST(macro, list)                                \
-  BOOST_PP_LIST_FOR_EACH(AKANTU_BOOST_LIST_MACRO, macro, list)
-
-#define AKANTU_BOOST_ALL_KIND(macro)                                           \
-  AKANTU_BOOST_ALL_KIND_LIST(macro, AKANTU_ELEMENT_KIND_BOOST_LIST)
-
-#define AKANTU_BOOST_ALL_KIND_SWITCH(macro)                                    \
-  AKANTU_BOOST_LIST_SWITCH(macro, AKANTU_ELEMENT_KIND, kind)
 
 @AKANTU_ELEMENT_KINDS_BOOST_MACROS@
 
 // /// define kept for compatibility reasons (they are most probably not needed
 // /// anymore) \todo check if they can be removed
 // #define AKANTU_REGULAR_ELEMENT_TYPE	AKANTU_ek_regular_ELEMENT_TYPE
 // #define AKANTU_COHESIVE_ELEMENT_TYPE	AKANTU_ek_cohesive_ELEMENT_TYPE
 // #define AKANTU_STRUCTURAL_ELEMENT_TYPE  AKANTU_ek_structural_ELEMENT_TYPE
 // #define AKANTU_IGFEM_ELEMENT_TYPE       AKANTU_ek_igfem_ELEMENT_TYPE
 
+
 /* -------------------------------------------------------------------------- */
 /* Lists of interests for FEEngineTemplate functions                          */
 /* -------------------------------------------------------------------------- */
 @AKANTU_FE_ENGINE_LISTS@
 
 } // akantu
 
 #endif /* AKANTU_AKA_ELEMENT_CLASSES_INFO_HH_ */
 
 #include "aka_element_classes_info_inline_impl.hh"
diff --git a/src/common/aka_element_classes_info_inline_impl.hh b/src/common/aka_element_classes_info_inline_impl.hh
index e5f889756..1407320ff 100644
--- a/src/common/aka_element_classes_info_inline_impl.hh
+++ b/src/common/aka_element_classes_info_inline_impl.hh
@@ -1,55 +1,226 @@
 /**
  * @file   aka_element_classes_info_inline_impl.hh
  *
  * @author Aurelia Isabel Cuba Ramos <aurelia.cubaramos@epfl.ch>
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Thu Jun 18 2015
  * @date last modification: Tue Sep 29 2020
  *
  * @brief  Implementation of the streaming fonction for the element classes
  * enums
  *
  *
  * @section LICENSE
  *
  * Copyright (©) 2015-2021 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_config.hh"
+#include "aka_tuple_tools.hh"
+/* -------------------------------------------------------------------------- */
+#include <tuple>
 #include <unordered_map>
 /* -------------------------------------------------------------------------- */
 
 #ifndef AKANTU_AKA_ELEMENT_CLASSES_INFO_INLINE_IMPL_HH_
 #define AKANTU_AKA_ELEMENT_CLASSES_INFO_INLINE_IMPL_HH_
 
 namespace akantu {
 
+/* -------------------------------------------------------------------------- */
+// BOOST PART: TOUCH ONLY IF YOU KNOW WHAT YOU ARE DOING
+#define AKANTU_BOOST_CASE_MACRO(r, macro, _type)                               \
+  case _type: {                                                                \
+    macro(_type);                                                              \
+    break;                                                                     \
+  }
+
+#define AKANTU_BOOST_LIST_SWITCH(macro1, list1, var)                           \
+  do {                                                                         \
+    switch (var) {                                                             \
+      BOOST_PP_SEQ_FOR_EACH(AKANTU_BOOST_CASE_MACRO, macro1, list1)            \
+    default: {                                                                 \
+      AKANTU_ERROR("Type (" << var << ") not handled by this function"); \
+    }                                                                          \
+    }                                                                          \
+  } while (0)
+
+#define AKANTU_BOOST_LIST_SWITCH_NO_DEFAULT(macro1, list1, var)                \
+  do {                                                                         \
+    switch (var) {                                                             \
+      BOOST_PP_SEQ_FOR_EACH(AKANTU_BOOST_CASE_MACRO, macro1, list1)            \
+    case _not_defined:                                                         \
+      break;                                                                   \
+    case _max_element_type:                                                    \
+      break;                                                                   \
+    }                                                                          \
+  } while (0)
+
+#define AKANTU_BOOST_LIST_SWITCH_CONSTEXPR(macro1, list1, var)                 \
+  do {                                                                         \
+    switch (var) {                                                             \
+      BOOST_PP_SEQ_FOR_EACH(AKANTU_BOOST_CASE_MACRO, macro1, list1)            \
+    default: {                                                                 \
+        macro1(_not_defined);                                                  \
+    }                                                                          \
+    }                                                                          \
+  } while (0)
+
+#define AKANTU_BOOST_ELEMENT_SWITCH(macro1, list1)                             \
+  AKANTU_BOOST_LIST_SWITCH(macro1, list1, type)
+
+#define AKANTU_BOOST_ELEMENT_SWITCH_NO_DEFAULT(macro1, list1)                  \
+  AKANTU_BOOST_LIST_SWITCH_NO_DEFAULT(macro1, list1, type)
+
+#define AKANTU_BOOST_ELEMENT_SWITCH_CONSTEXPR(macro1, list1)                   \
+  AKANTU_BOOST_LIST_SWITCH_CONSTEXPR(macro1, list1, type)
+
+#define AKANTU_BOOST_ALL_ELEMENT_SWITCH(macro)                                 \
+  AKANTU_BOOST_ELEMENT_SWITCH(macro, AKANTU_ALL_ELEMENT_TYPE)
+
+#define AKANTU_BOOST_ALL_ELEMENT_SWITCH_NO_DEFAULT(macro)                      \
+  AKANTU_BOOST_ELEMENT_SWITCH_NO_DEFAULT(macro, AKANTU_ALL_ELEMENT_TYPE)
+
+#define AKANTU_BOOST_ALL_ELEMENT_SWITCH_CONSTEXPR(macro)                       \
+  AKANTU_BOOST_ELEMENT_SWITCH_CONSTEXPR(macro, AKANTU_ALL_ELEMENT_TYPE)
+
+#define AKANTU_BOOST_LIST_MACRO(r, macro, type) macro(type)
+
+#define AKANTU_BOOST_APPLY_ON_LIST(macro, list)                                \
+  BOOST_PP_SEQ_FOR_EACH(AKANTU_BOOST_LIST_MACRO, macro, list)
+
+#define AKANTU_BOOST_ALL_ELEMENT_LIST(macro)                                   \
+  AKANTU_BOOST_APPLY_ON_LIST(macro, AKANTU_ALL_ELEMENT_TYPE)
+
+#define AKANTU_GET_ELEMENT_LIST(kind) AKANTU##kind##_ELEMENT_TYPE
+
+#define AKANTU_BOOST_KIND_ELEMENT_SWITCH(macro, kind)                          \
+  AKANTU_BOOST_ELEMENT_SWITCH(macro, AKANTU_GET_ELEMENT_LIST(kind))
+
+// BOOST_PP_SEQ_TO_LIST does not exists in Boost < 1.49
+#define AKANTU_GENERATE_KIND_LIST(seq)                                         \
+  BOOST_PP_TUPLE_TO_LIST(BOOST_PP_SEQ_SIZE(seq), BOOST_PP_SEQ_TO_TUPLE(seq))
+
+#define AKANTU_ELEMENT_KIND_BOOST_LIST                                         \
+  AKANTU_GENERATE_KIND_LIST(AKANTU_ELEMENT_KIND)
+
+#define AKANTU_BOOST_ALL_KIND_LIST(macro, list)                                \
+  BOOST_PP_LIST_FOR_EACH(AKANTU_BOOST_LIST_MACRO, macro, list)
+
+#define AKANTU_BOOST_ALL_KIND(macro)                                           \
+  AKANTU_BOOST_ALL_KIND_LIST(macro, AKANTU_ELEMENT_KIND_BOOST_LIST)
+
+#define AKANTU_BOOST_ALL_KIND_SWITCH(macro)                                    \
+  AKANTU_BOOST_LIST_SWITCH(macro, AKANTU_ELEMENT_KIND, kind)
+
+/* -------------------------------------------------------------------------- */
+/* -------------------------------------------------------------------------- */
+
 AKANTU_ENUM_OUTPUT_STREAM(
     ElementType, AKANTU_ALL_ELEMENT_TYPE(_not_defined)(_max_element_type))
 AKANTU_ENUM_INPUT_STREAM(ElementType, AKANTU_ALL_ELEMENT_TYPE)
 
 AKANTU_ENUM_OUTPUT_STREAM(InterpolationType, AKANTU_INTERPOLATION_TYPES)
 AKANTU_ENUM_INPUT_STREAM(InterpolationType, AKANTU_INTERPOLATION_TYPES)
 
 AKANTU_ENUM_OUTPUT_STREAM(ElementKind, AKANTU_ELEMENT_KIND)
 AKANTU_ENUM_INPUT_STREAM(ElementKind, AKANTU_ELEMENT_KIND)
 
+template <::akantu::ElementType t>
+using element_type_t = std::integral_constant<::akantu::ElementType, t>;
+
+#define OP_CAT(s, data, elem) BOOST_PP_CAT(_element_type, elem)
+
+// creating a type instead of a using helps to debug
+#define AKANTU_DECLARE_ELEMENT_TYPE_STRUCT(r, data, elem)                      \
+  struct BOOST_PP_CAT(_element_type, elem)                                     \
+      : public element_type_t<::akantu::elem> {};
+
+BOOST_PP_SEQ_FOR_EACH(AKANTU_DECLARE_ELEMENT_TYPE_STRUCT, _,
+                      AKANTU_ALL_ELEMENT_TYPE)
+
+#undef AKANTU_DECLARE_ELEMENT_TYPE_STRUCT
+
+template <ElementKind kind> struct ElementTypes;
+
+template <> struct ElementTypes<_ek_regular> {
+  using type = std::tuple<BOOST_PP_SEQ_ENUM(
+      BOOST_PP_SEQ_TRANSFORM(OP_CAT, _, AKANTU_ek_regular_ELEMENT_TYPE))>;
+};
+
+#if defined(AKANTU_COHESIVE_ELEMENT)
+template <> struct ElementTypes<_ek_cohesive> {
+  using type = std::tuple<BOOST_PP_SEQ_ENUM(
+      BOOST_PP_SEQ_TRANSFORM(OP_CAT, _, AKANTU_ek_cohesive_ELEMENT_TYPE))>;
+};
+#endif
+
+#if defined(AKANTU_STRUCTURAL_MECHANICS)
+template <> struct ElementTypes<_ek_structural> {
+  using type = std::tuple<BOOST_PP_SEQ_ENUM(
+      BOOST_PP_SEQ_TRANSFORM(OP_CAT, _, AKANTU_ek_structural_ELEMENT_TYPE))>;
+};
+#endif
+
+#undef OP_CAT
+
+template <ElementKind kind>
+using ElementTypes_t = typename ElementTypes<kind>::type;
+
+#define OP_CAT(s, data, elem) ElementTypes_t<elem>
+using AllElementTypes = tuple::cat_t<BOOST_PP_SEQ_ENUM(
+      BOOST_PP_SEQ_TRANSFORM(OP_CAT, _, AKANTU_ELEMENT_KIND))>;
+#undef OP_CAT
+    
+namespace details {
+  template <std::size_t S> struct visit_tuple_impl {
+    template <class Function, class DynamicType, class Tuple>
+    static constexpr decltype(auto) visit(const Tuple &, Function && function,
+                                          const DynamicType & type) {
+      using integral_type = std::tuple_element_t<S - 1, Tuple>;
+      if (integral_type::value == type) {
+        return function(integral_type{});
+      } else {
+        return visit_tuple_impl<S - 1>::visit(
+            Tuple{}, std::forward<Function>(function), type);
+      }
+    }
+  };
+
+  template <> struct visit_tuple_impl<0> {
+    template <class Function, class DynamicType, class Tuple>
+    static constexpr decltype(auto) visit(const Tuple &, Function && function,
+                                          const DynamicType & /*type*/) {
+      return function(element_type_t<_not_defined>{});
+    }
+  };
+} // namespace details
+
+template <class Tuple, class Function, class DynamicType>
+constexpr decltype(auto) tuple_dispatch(Function && function,
+                                        const DynamicType & type) {
+  return details::visit_tuple_impl<std::tuple_size<Tuple>::value>::visit(
+      Tuple{}, std::forward<Function>(function), type);
+}
+
+
 } // namespace akantu
 
 #endif /* AKANTU_AKA_ELEMENT_CLASSES_INFO_INLINE_IMPL_HH_ */
diff --git a/src/fe_engine/fe_engine_template.hh b/src/fe_engine/fe_engine_template.hh
index 7eaa35dba..b6c8bac0a 100644
--- a/src/fe_engine/fe_engine_template.hh
+++ b/src/fe_engine/fe_engine_template.hh
@@ -1,432 +1,450 @@
 /**
  * @file   fe_engine_template.hh
  *
  * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
  * @author Sébastien Hartmann <sebastien.hartmann@epfl.ch>
  * @author Mohit Pundir <mohit.pundir@epfl.ch>
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Fri Jun 18 2010
  * @date last modification: Fri May 14 2021
  *
  * @brief  templated class that calls integration and shape objects
  *
  *
  * @section LICENSE
  *
  * Copyright (©) 2010-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  * Akantu is free software: you can redistribute it and/or modify it under the
  * terms of the GNU Lesser General Public License as published by the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  * 
  * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
  * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
  * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
  * details.
  * 
  * You should have received a copy of the GNU Lesser General Public License
  * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 #include "fe_engine.hh"
 //#include "integrator.hh"
 //#include "shape_functions.hh"
 /* -------------------------------------------------------------------------- */
 #include <type_traits>
 /* -------------------------------------------------------------------------- */
 
 /* -------------------------------------------------------------------------- */
 #ifndef AKANTU_FE_ENGINE_TEMPLATE_HH_
 #define AKANTU_FE_ENGINE_TEMPLATE_HH_
 
 namespace akantu {
 class Integrator;
 class ShapeFunctions;
 } // namespace akantu
 
 namespace akantu {
 class DOFManager;
 namespace fe_engine {
   namespace details {
     template <ElementKind> struct AssembleLumpedTemplateHelper;
     template <ElementKind> struct AssembleFieldMatrixHelper;
   } // namespace details
 } // namespace fe_engine
 
 template <ElementKind, typename> struct AssembleFieldMatrixStructHelper;
 
 struct DefaultIntegrationOrderFunctor {
   template <ElementType type> static inline constexpr int getOrder() {
     return ElementClassProperty<type>::polynomial_degree;
   }
 };
 
 /* -------------------------------------------------------------------------- */
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind = _ek_regular,
           class IntegrationOrderFunctor = DefaultIntegrationOrderFunctor>
 class FEEngineTemplate : public FEEngine {
   /* ------------------------------------------------------------------------ */
   /* Constructors/Destructors                                                 */
   /* ------------------------------------------------------------------------ */
 public:
   using Integ = I<kind, IntegrationOrderFunctor>;
   using Shape = S<kind>;
 
   FEEngineTemplate(Mesh & mesh, Int spatial_dimension = _all_dimensions,
                    const ID & id = "fem");
 
   ~FEEngineTemplate() override;
 
   /* ------------------------------------------------------------------------ */
   /* Methods                                                                  */
   /* ------------------------------------------------------------------------ */
 public:
   /// pre-compute all the shape functions, their derivatives and the jacobians
   void initShapeFunctions(GhostType ghost_type = _not_ghost) override;
   void initShapeFunctions(const Array<Real> & nodes,
                           GhostType ghost_type = _not_ghost);
 
   /* ------------------------------------------------------------------------ */
   /* Integration method bridges                                               */
   /* ------------------------------------------------------------------------ */
   /// integrate f for all elements of type "type"
   void
   integrate(const Array<Real> & f, Array<Real> & intf, Int nb_degree_of_freedom,
             const ElementType & type, const GhostType & ghost_type = _not_ghost,
             const Array<Idx> & filter_elements = empty_filter) const override;
 
   /// integrate a scalar value on all elements of type "type"
   Real
   integrate(const Array<Real> & f, const ElementType & type,
             const GhostType & ghost_type = _not_ghost,
             const Array<Idx> & filter_elements = empty_filter) const override;
 
   /// integrate one element scalar value on all elements of type "type"
   Real integrate(const Ref<const VectorXr> & f, const ElementType & type,
                  Int index,
                  const GhostType & ghost_type = _not_ghost) const override;
 
   /// integrate partially around an integration point (@f$ intf_q = f_q * J_q *
   /// w_q @f$)
   void integrateOnIntegrationPoints(
       const Array<Real> & f, Array<Real> & intf, Int nb_degree_of_freedom,
       const ElementType & type, const GhostType & ghost_type = _not_ghost,
       const Array<Idx> & filter_elements = empty_filter) const override;
 
   /// interpolate on a phyiscal point inside an element
   void interpolate(const Ref<const VectorXr> & real_coords,
                    const Ref<const MatrixXr> & nodal_values,
                    Ref<VectorXr> interpolated,
                    const Element & element) const override;
 
   /// get the number of integration points
   Int getNbIntegrationPoints(
       const ElementType & type,
       const GhostType & ghost_type = _not_ghost) const override;
 
   /// get shapes precomputed
   const Array<Real> & getShapes(const ElementType & type,
                                 const GhostType & ghost_type = _not_ghost,
                                 Int id = 0) const override;
 
   /// get the derivatives of shapes
   const Array<Real> &
   getShapesDerivatives(const ElementType & type,
                        const GhostType & ghost_type = _not_ghost,
                        Int id = 0) const override;
 
   /// get integration points
   const inline Matrix<Real> &
   getIntegrationPoints(ElementType type,
                        GhostType ghost_type = _not_ghost) const override;
 
   /* ------------------------------------------------------------------------ */
   /* Shape method bridges                                                     */
   /* ------------------------------------------------------------------------ */
 
   /// compute the gradient of a nodal field on the integration points
   void gradientOnIntegrationPoints(
       const Array<Real> & u, Array<Real> & nablauq,
       const Int nb_degree_of_freedom, const ElementType & type,
       const GhostType & ghost_type = _not_ghost,
       const Array<Idx> & filter_elements = empty_filter) const override;
 
   /// interpolate a nodal field on the integration points
   void interpolateOnIntegrationPoints(
       const Array<Real> & u, Array<Real> & uq, Int nb_degree_of_freedom,
       const ElementType & type, const GhostType & ghost_type = _not_ghost,
       const Array<Idx> & filter_elements = empty_filter) const override;
 
   /// interpolate a nodal field on the integration points given a
   /// by_element_type
   void interpolateOnIntegrationPoints(
       const Array<Real> & u, ElementTypeMapArray<Real> & uq,
       const ElementTypeMapArray<Idx> * filter_elements =
           nullptr) const override;
 
   /// pre multiplies a tensor by the shapes derivaties
   void
   computeBtD(const Array<Real> & Ds, Array<Real> & BtDs,
              const ElementType & type, const GhostType & ghost_type,
              const Array<Idx> & filter_elements = empty_filter) const override;
 
   /// left and right  multiplies a tensor by the shapes derivaties
   void
   computeBtDB(const Array<Real> & Ds, Array<Real> & BtDBs, Int order_d,
               const ElementType & type, const GhostType & ghost_type,
               const Array<Idx> & filter_elements = empty_filter) const override;
 
   /// left multiples a vector by the shape functions
   void computeNtb(const Array<Real> & bs, Array<Real> & Ntbs,
                   const ElementType & type, const GhostType & ghost_type,
                   const Array<Idx> & filter_elements) const override;
 
   /// compute the position of integration points given by an element_type_map
   /// from nodes position
   inline void computeIntegrationPointsCoordinates(
       ElementTypeMapArray<Real> & quadrature_points_coordinates,
       const ElementTypeMapArray<Idx> * filter_elements =
           nullptr) const override;
 
   /// compute the position of integration points from nodes position
   inline void computeIntegrationPointsCoordinates(
       Array<Real> & quadrature_points_coordinates, const ElementType & type,
       const GhostType & ghost_type = _not_ghost,
       const Array<Idx> & filter_elements = empty_filter) const override;
 
   /// interpolate field at given position (interpolation_points) from given
   /// values of this field at integration points (field)
   inline void interpolateElementalFieldFromIntegrationPoints(
       const ElementTypeMapArray<Real> & field,
       const ElementTypeMapArray<Real> & interpolation_points_coordinates,
       ElementTypeMapArray<Real> & result, const GhostType ghost_type,
       const ElementTypeMapArray<Idx> * element_filter) const override;
 
   /// Interpolate field at given position from given values of this field at
   /// integration points (field)
   /// using matrices precomputed with
   /// initElementalFieldInterplationFromIntegrationPoints
   inline void interpolateElementalFieldFromIntegrationPoints(
       const ElementTypeMapArray<Real> & field,
       const ElementTypeMapArray<Real> &
           interpolation_points_coordinates_matrices,
       const ElementTypeMapArray<Real> & quad_points_coordinates_inv_matrices,
       ElementTypeMapArray<Real> & result, const GhostType ghost_type,
       const ElementTypeMapArray<Idx> * element_filter) const override;
 
   /// Build pre-computed matrices for interpolation of field form integration
   /// points at other given positions (interpolation_points)
   inline void initElementalFieldInterpolationFromIntegrationPoints(
       const ElementTypeMapArray<Real> & interpolation_points_coordinates,
       ElementTypeMapArray<Real> & interpolation_points_coordinates_matrices,
       ElementTypeMapArray<Real> & quad_points_coordinates_inv_matrices,
       const ElementTypeMapArray<Idx> * element_filter = nullptr) const override;
 
   /// find natural coords from real coords provided an element
   void inverseMap(const Vector<Real> & real_coords, Int element,
                   const ElementType & type, Vector<Real> & natural_coords,
                   const GhostType & ghost_type = _not_ghost) const;
 
   /// return true if the coordinates provided are inside the element, false
   /// otherwise
   inline bool contains(const Vector<Real> & real_coords, Int element,
                        const ElementType & type,
                        const GhostType & ghost_type = _not_ghost) const;
 
   /// compute the shape on a provided point
   inline void
   computeShapes(const Ref<const VectorXr> & real_coords, Int element,
                 const ElementType & type, Ref<VectorXr> shapes,
                 const GhostType & ghost_type = _not_ghost) const override;
 
   /// compute the shape derivatives on a provided point
   inline void computeShapeDerivatives(
       const Ref<const VectorXr> & real__coords, Int element,
       const ElementType & type, Ref<MatrixXr> shape_derivatives,
       const GhostType & ghost_type = _not_ghost) const override;
 
   /* ------------------------------------------------------------------------ */
   /* Other methods                                                            */
   /* ------------------------------------------------------------------------ */
   /// pre-compute normals on integration points
   void
   computeNormalsOnIntegrationPoints(GhostType ghost_type = _not_ghost) override;
   void
   computeNormalsOnIntegrationPoints(const Array<Real> & field,
                                     GhostType ghost_type = _not_ghost) override;
   void computeNormalsOnIntegrationPoints(
-      const Array<Real> & field, Array<Real> & normal, ElementType type,
-      GhostType ghost_type = _not_ghost) const override;
-  template <ElementType type>
+      const GhostType & ghost_type = _not_ghost) override;
+  void computeNormalsOnIntegrationPoints(
+      const Array<Real> & field,
+      const GhostType & ghost_type = _not_ghost) override;
+  void computeNormalsOnIntegrationPoints(
+      const Array<Real> & field, Array<Real> & normal, const ElementType & type,
+      const GhostType & ghost_type = _not_ghost) const override;
+
+  template <
+      ElementType type, ElementKind kind_ = kind,
+      std::enable_if_t<kind_ != _ek_regular> * = nullptr>
+  void computeNormalsOnIntegrationPoints(const Array<Real> & /*field*/,
+                                         Array<Real> & /*normal*/,
+                                         const GhostType & /*ghost_type*/) const {
+    AKANTU_TO_IMPLEMENT();
+  }
+
+
+  template <
+      ElementType type, ElementKind kind_ = kind,
+      std::enable_if_t<kind_ == _ek_regular and type != _point_1> * = nullptr>
   void computeNormalsOnIntegrationPoints(const Array<Real> & field,
                                          Array<Real> & normal,
-                                         GhostType ghost_type) const;
+                                         const GhostType & ghost_type) const;
 
 private:
   // To avoid a weird full specialization of a method in a non specalized class
   void computeNormalsOnIntegrationPointsPoint1(const Array<Real> & /*unused*/,
                                                Array<Real> & normal,
                                                GhostType ghost_type) const;
 
 public:
   /// function to print the contain of the class
   void printself(std::ostream & stream, int indent = 0) const override;
 
   void assembleFieldLumped(
       const std::function<void(Matrix<Real> &, const Element &)> & field_funct,
       const ID & matrix_id, const ID & dof_id, DOFManager & dof_manager,
       ElementType type, GhostType ghost_type) const override;
 
   /// assemble a field as a matrix (ex. rho to mass matrix)
   void assembleFieldMatrix(
       const std::function<void(Matrix<Real> &, const Element &)> & field_funct,
       const ID & matrix_id, const ID & dof_id, DOFManager & dof_manager,
       ElementType type, GhostType ghost_type) const override;
 
   /// assemble a field as a lumped matrix (ex. rho in lumped mass)
   // template <class Functor>
   // void assembleFieldLumped(const Functor & field_funct, const ID & matrix_id,
   //                          const ID & dof_id, DOFManager & dof_manager,
   //                          ElementType type,
   //                          GhostType ghost_type) const;
 
   // /// assemble a field as a matrix (ex. rho to mass matrix)
   // template <class Functor>
   // void assembleFieldMatrix(const Functor & field_funct, const ID & matrix_id,
   //                          const ID & dof_id, DOFManager & dof_manager,
   //                          ElementType type,
   //                          GhostType ghost_type) const;
 
   // #ifdef AKANTU_STRUCTURAL_MECHANICS
   //   /// assemble a field as a matrix (ex. rho to mass matrix)
   //   void assembleFieldMatrix(const Array<Real> & field_1,
   //                            UInt nb_degree_of_freedom, SparseMatrix & M,
   //                            Array<Real> * n,
   //                            ElementTypeMapArray<Real> & rotation_mat,
   //                            ElementType type,
   //                            GhostType ghost_type = _not_ghost)
   //                            const;
 
   //   /// compute shapes function in a matrix for structural elements
   //   void
   //   computeShapesMatrix(ElementType type, UInt nb_degree_of_freedom,
   //                       UInt nb_nodes_per_element, Array<Real> * n, UInt id,
   //                       UInt degree_to_interpolate, UInt degree_interpolated,
   //                       const bool sign,
   //                       GhostType ghost_type = _not_ghost) const
   //                       override;
   // #endif
 
 private:
   friend struct fe_engine::details::AssembleLumpedTemplateHelper<kind>;
   friend struct fe_engine::details::AssembleFieldMatrixHelper<kind>;
   friend struct AssembleFieldMatrixStructHelper<kind, void>;
 
   /// templated function to compute the scaling to assemble a lumped matrix
   template <ElementType type>
   void assembleFieldLumped(
       const std::function<void(Matrix<Real> &, const Element &)> & field_funct,
       const ID & matrix_id, const ID & dof_id, DOFManager & dof_manager,
       GhostType ghost_type) const;
 
   /// @f$ \tilde{M}_{i} = \sum_j M_{ij} = \sum_j \int \rho \varphi_i \varphi_j
   /// dV = \int \rho \varphi_i dV @f$
   template <ElementType type>
   void assembleLumpedRowSum(const Array<Real> & field, const ID & matrix_id,
                             const ID & dof_id, DOFManager & dof_manager,
                             GhostType ghost_type) const;
 
   /// @f$ \tilde{M}_{i} = c * M_{ii} = \int_{V_e} \rho dV @f$
   template <ElementType type>
   void assembleLumpedDiagonalScaling(const Array<Real> & field,
                                      const ID & matrix_id, const ID & dof_id,
                                      DOFManager & dof_manager,
                                      GhostType ghost_type) const;
 
   /// assemble a field as a matrix (ex. rho to mass matrix)
   template <ElementType type>
   void assembleFieldMatrix(
       const std::function<void(Matrix<Real> &, const Element &)> & field_funct,
       const ID & matrix_id, const ID & dof_id, DOFManager & dof_manager,
       GhostType ghost_type) const;
 
 #ifdef AKANTU_STRUCTURAL_MECHANICS
 
   /// assemble a field as a matrix for structural elements (ex. rho to mass
   /// matrix)
   template <ElementType type>
   void assembleFieldMatrix(const Array<Real> & field_1,
                            Int nb_degree_of_freedom, SparseMatrix & M,
                            Array<Real> * n,
                            ElementTypeMapArray<Real> & rotation_mat,
                            __attribute__((unused)) GhostType ghost_type) const;
 
 #endif
 
   /* ------------------------------------------------------------------------ */
   /* Mesh Event Handler interface                                             */
   /* ------------------------------------------------------------------------ */
 public:
   void onElementsAdded(const Array<Element> & /*new_elements*/,
                        const NewElementsEvent & /*unused*/) override;
   void onElementsRemoved(const Array<Element> & /*unused*/,
                          const ElementTypeMapArray<UInt> & /*unused*/,
                          const RemovedElementsEvent & /*unused*/) override;
   void onElementsChanged(const Array<Element> & /*unused*/,
                          const Array<Element> & /*unused*/,
                          const ElementTypeMapArray<UInt> & /*unused*/,
                          const ChangedElementsEvent & /*unused*/) override;
 
   /* ------------------------------------------------------------------------ */
   /* Accessors                                                                */
   /* ------------------------------------------------------------------------ */
 public:
   /// get the shape class (probably useless: see getShapeFunction)
   const ShapeFunctions & getShapeFunctionsInterface() const override {
     return shape_functions;
   };
   /// get the shape class
   const Shape & getShapeFunctions() const { return shape_functions; };
 
   /// get the integrator class (probably useless: see getIntegrator)
   const Integrator & getIntegratorInterface() const override {
     return integrator;
   };
   /// get the integrator class
   const Integ & getIntegrator() const { return integrator; };
 
   /* ------------------------------------------------------------------------ */
   /* Class Members                                                            */
   /* ------------------------------------------------------------------------ */
 private:
   Integ integrator;
   Shape shape_functions;
 };
 
 } // namespace akantu
 
 /* -------------------------------------------------------------------------- */
 /* inline functions                                                           */
 /* -------------------------------------------------------------------------- */
 #include "fe_engine_template_tmpl.hh"
 #include "fe_engine_template_tmpl_field.hh"
 /* -------------------------------------------------------------------------- */
 /* Shape Linked specialization                                                */
 /* -------------------------------------------------------------------------- */
 #if defined(AKANTU_STRUCTURAL_MECHANICS)
 #include "fe_engine_template_tmpl_struct.hh"
 #endif
 /* -------------------------------------------------------------------------- */
 /* Shape IGFEM specialization                                                 */
 /* -------------------------------------------------------------------------- */
 #if defined(AKANTU_IGFEM)
 #include "fe_engine_template_tmpl_igfem.hh"
 #endif
 
 #endif /* AKANTU_FE_ENGINE_TEMPLATE_HH_ */
diff --git a/src/fe_engine/fe_engine_template_tmpl.hh b/src/fe_engine/fe_engine_template_tmpl.hh
index fa2216887..b44b6522a 100644
--- a/src/fe_engine/fe_engine_template_tmpl.hh
+++ b/src/fe_engine/fe_engine_template_tmpl.hh
@@ -1,1446 +1,1134 @@
 /**
  * @file   fe_engine_template_tmpl.hh
  *
  * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
  * @author Dana Christen <dana.christen@gmail.com>
  * @author Mauro Corrado <mauro.corrado@epfl.ch>
  * @author Aurelia Isabel Cuba Ramos <aurelia.cubaramos@epfl.ch>
  * @author Sébastien Hartmann <sebastien.hartmann@epfl.ch>
  * @author Mohit Pundir <mohit.pundir@epfl.ch>
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  * @author Marco Vocialta <marco.vocialta@epfl.ch>
  *
  * @date creation: Tue Feb 15 2011
  * @date last modification: Fri May 14 2021
  *
  * @brief  Template implementation of FEEngineTemplate
  *
  *
  * @section LICENSE
  *
  * Copyright (©) 2010-2021 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 "dof_manager.hh"
 #include "fe_engine_template.hh"
 /* -------------------------------------------------------------------------- */
 
+#define DISPATCH_HELPER(Func, ...)                                             \
+  auto && call = [&](auto && enum_type) {                                      \
+    return this->Func<ElementType(enum_type)>(__VA_ARGS__);                    \
+  };                                                                           \
+  tuple_dispatch<ElementTypes_t<kind>>(call, type);
+
 namespace akantu {
 
 /* -------------------------------------------------------------------------- */
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::FEEngineTemplate(
     Mesh & mesh, Int spatial_dimension, const ID & id)
     : FEEngine(mesh, spatial_dimension, id),
       integrator(mesh, spatial_dimension, id),
       shape_functions(mesh, spatial_dimension, id) {}
 
 /* -------------------------------------------------------------------------- */
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::~FEEngineTemplate() =
     default;
 
 /* -------------------------------------------------------------------------- */
-/**
- * Helper class to be able to write a partial specialization on the element kind
- */
-namespace fe_engine {
-  namespace details {
-    template <ElementKind kind> struct GradientOnIntegrationPointsHelper {
-      template <class S>
-      static void call(const S &, Mesh &, const Array<Real> &, Array<Real> &,
-                       Int, const ElementType &, const GhostType &,
-                       const Array<Idx> &) {
-        AKANTU_TO_IMPLEMENT();
-      }
-    };
-
-#define COMPUTE_GRADIENT(type)                                                 \
-  if (element_dimension == ElementClass<type>::getSpatialDimension())          \
-    shape_functions.template gradientOnIntegrationPoints<type>(                \
-        u, nablauq, nb_degree_of_freedom, ghost_type, filter_elements);
-
-#define AKANTU_SPECIALIZE_GRADIENT_ON_INTEGRATION_POINTS_HELPER(kind)          \
-  template <> struct GradientOnIntegrationPointsHelper<kind> {                 \
-    template <class S>                                                         \
-    static void call(const S & shape_functions, Mesh & mesh,                   \
-                     const Array<Real> & u, Array<Real> & nablauq,             \
-                     Int nb_degree_of_freedom, const ElementType & type,       \
-                     const GhostType & ghost_type,                             \
-                     const Array<Idx> & filter_elements) {                     \
-      auto element_dimension = mesh.getSpatialDimension(type);                 \
-      AKANTU_BOOST_KIND_ELEMENT_SWITCH(COMPUTE_GRADIENT, kind);                \
-    }                                                                          \
-  };
-
-    AKANTU_BOOST_ALL_KIND_LIST(
-        AKANTU_SPECIALIZE_GRADIENT_ON_INTEGRATION_POINTS_HELPER,
-        AKANTU_FE_ENGINE_LIST_GRADIENT_ON_INTEGRATION_POINTS)
-
-#undef AKANTU_SPECIALIZE_GRADIENT_ON_INTEGRATION_POINTS_HELPER
-#undef COMPUTE_GRADIENT
-  } // namespace details
-} // namespace fe_engine
-
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::
     gradientOnIntegrationPoints(const Array<Real> & u, Array<Real> & nablauq,
                                 Int nb_degree_of_freedom,
                                 const ElementType & type,
                                 const GhostType & ghost_type,
                                 const Array<Int> & filter_elements) const {
   AKANTU_DEBUG_IN();
 
   auto nb_element = mesh.getNbElement(type, ghost_type);
   if (filter_elements != empty_filter) {
     nb_element = filter_elements.size();
   }
   auto nb_points =
       shape_functions.getIntegrationPoints(type, ghost_type).cols();
 
 #ifndef AKANTU_NDEBUG
 
   auto element_dimension = mesh.getSpatialDimension(type);
 
   AKANTU_DEBUG_ASSERT(u.size() == mesh.getNbNodes(),
                       "The vector u(" << u.getID()
                                       << ") has not the good size.");
   AKANTU_DEBUG_ASSERT(u.getNbComponent() == nb_degree_of_freedom,
                       "The vector u("
                           << u.getID()
                           << ") has not the good number of component.");
 
   AKANTU_DEBUG_ASSERT(
       nablauq.getNbComponent() == nb_degree_of_freedom * element_dimension,
       "The vector nablauq(" << nablauq.getID()
                             << ") has not the good number of component.");
 
 // AKANTU_DEBUG_ASSERT(nablauq.size() == nb_element * nb_points,
 //                  "The vector nablauq(" << nablauq.getID()
 //                  << ") has not the good size.");
 #endif
 
   nablauq.resize(nb_element * nb_points);
 
-  fe_engine::details::GradientOnIntegrationPointsHelper<kind>::call(
-      shape_functions, mesh, u, nablauq, nb_degree_of_freedom, type, ghost_type,
-      filter_elements);
+  auto call = [&](auto && integral_type) {
+    constexpr ElementType type = integral_type;
+    if (element_dimension == ElementClass<type>::getSpatialDimension())
+      shape_functions.template gradientOnIntegrationPoints<type>(
+          u, nablauq, nb_degree_of_freedom, ghost_type, filter_elements);
+  };
+
+  tuple_dispatch<ElementTypes_t<kind>>(call, type);
+  // fe_engine::details::GradientOnIntegrationPointsHelper<kind>::call(
+  //     shape_functions, mesh, u, nablauq, nb_degree_of_freedom, type,
+  //     ghost_type, filter_elements);
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::initShapeFunctions(
     GhostType ghost_type) {
   initShapeFunctions(mesh.getNodes(), ghost_type);
 }
 
 /* -------------------------------------------------------------------------- */
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::initShapeFunctions(
     const Array<Real> & nodes, GhostType ghost_type) {
   AKANTU_DEBUG_IN();
 
   for (auto & type : mesh.elementTypes(element_dimension, ghost_type, kind)) {
     integrator.initIntegrator(nodes, type, ghost_type);
     const auto & control_points = getIntegrationPoints(type, ghost_type);
     shape_functions.initShapeFunctions(nodes, control_points, type, ghost_type);
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
-/**
- * Helper class to be able to write a partial specialization on the element kind
- */
-namespace fe_engine {
-  namespace details {
-    template <ElementKind kind> struct IntegrateHelper {};
-
-#define INTEGRATE(type)                                                        \
-  integrator.template integrate<type>(f, intf, nb_degree_of_freedom,           \
-                                      ghost_type, filter_elements);
-
-#define AKANTU_SPECIALIZE_INTEGRATE_HELPER(kind)                               \
-  template <> struct IntegrateHelper<kind> {                                   \
-    template <class I>                                                         \
-    static void call(const I & integrator, const Array<Real> & f,              \
-                     Array<Real> & intf, Int nb_degree_of_freedom,             \
-                     const ElementType & type, const GhostType & ghost_type,   \
-                     const Array<Idx> & filter_elements) {                     \
-      AKANTU_BOOST_KIND_ELEMENT_SWITCH(INTEGRATE, kind);                       \
-    }                                                                          \
-  };
-
-    AKANTU_BOOST_ALL_KIND(AKANTU_SPECIALIZE_INTEGRATE_HELPER)
-
-#undef AKANTU_SPECIALIZE_INTEGRATE_HELPER
-#undef INTEGRATE
-  } // namespace details
-} // namespace fe_engine
-
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::integrate(
     const Array<Real> & f, Array<Real> & intf, Int nb_degree_of_freedom,
     const ElementType & type, const GhostType & ghost_type,
     const Array<Idx> & filter_elements) const {
 
   auto nb_element = mesh.getNbElement(type, ghost_type);
   if (filter_elements != empty_filter) {
     nb_element = filter_elements.size();
   }
 #ifndef AKANTU_NDEBUG
 
   auto nb_quadrature_points = getNbIntegrationPoints(type);
 
   AKANTU_DEBUG_ASSERT(f.size() == nb_element * nb_quadrature_points,
                       "The vector f(" << f.getID() << " size " << f.size()
                                       << ") has not the good size ("
                                       << nb_element << ").");
   AKANTU_DEBUG_ASSERT(f.getNbComponent() == nb_degree_of_freedom,
                       "The vector f("
                           << f.getID()
                           << ") has not the good number of component.");
   AKANTU_DEBUG_ASSERT(intf.getNbComponent() == nb_degree_of_freedom,
                       "The vector intf("
                           << intf.getID()
                           << ") has not the good number of component.");
 #endif
 
   intf.resize(nb_element);
 
-  fe_engine::details::IntegrateHelper<kind>::call(integrator, f, intf,
-                                                  nb_degree_of_freedom, type,
-                                                  ghost_type, filter_elements);
-}
-
-/* -------------------------------------------------------------------------- */
-/**
- * Helper class to be able to write a partial specialization on the element kind
- */
-namespace fe_engine {
-  namespace details {
-    template <ElementKind kind> struct IntegrateScalarHelper {};
-
-#define INTEGRATE(type)                                                        \
-  integral =                                                                   \
-      integrator.template integrate<type>(f, ghost_type, filter_elements);
-
-#define AKANTU_SPECIALIZE_INTEGRATE_SCALAR_HELPER(kind)                        \
-  template <> struct IntegrateScalarHelper<kind> {                             \
-    template <class I>                                                         \
-    static Real call(const I & integrator, const Array<Real> & f,              \
-                     const ElementType & type, const GhostType & ghost_type,   \
-                     const Array<Int> & filter_elements) {                     \
-      Real integral = 0.;                                                      \
-      AKANTU_BOOST_KIND_ELEMENT_SWITCH(INTEGRATE, kind);                       \
-      return integral;                                                         \
-    }                                                                          \
+  auto call = [&](auto && integral_type) {
+    integrator.template integrate<ElementType(integral_type)>(
+        f, intf, nb_degree_of_freedom, ghost_type, filter_elements);
   };
 
-    AKANTU_BOOST_ALL_KIND(AKANTU_SPECIALIZE_INTEGRATE_SCALAR_HELPER)
-
-#undef AKANTU_SPECIALIZE_INTEGRATE_SCALAR_HELPER
-#undef INTEGRATE
-  } // namespace details
-} // namespace fe_engine
+  tuple_dispatch<ElementTypes_t<kind>>(call, type);
+}
 
+/* -------------------------------------------------------------------------- */
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 Real FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::integrate(
     const Array<Real> & f, const ElementType & type,
     const GhostType & ghost_type, const Array<Int> & filter_elements) const {
   AKANTU_DEBUG_IN();
 
 #ifndef AKANTU_NDEBUG
-  //   std::stringstream sstr; sstr << ghost_type;
-  //   AKANTU_DEBUG_ASSERT(sstr.str() == nablauq.getTag(),
-  //                  "The vector " << nablauq.getID() << " is not taged " <<
-  //                  ghost_type);
   auto nb_element = mesh.getNbElement(type, ghost_type);
   if (filter_elements != empty_filter) {
     nb_element = filter_elements.size();
   }
 
   auto nb_quadrature_points = getNbIntegrationPoints(type, ghost_type);
 
   AKANTU_DEBUG_ASSERT(
       f.size() == nb_element * nb_quadrature_points,
       "The vector f(" << f.getID() << ") has not the good size. (" << f.size()
                       << "!=" << nb_quadrature_points * nb_element << ")");
   AKANTU_DEBUG_ASSERT(f.getNbComponent() == 1,
                       "The vector f("
                           << f.getID()
                           << ") has not the good number of component.");
 #endif
 
-  Real integral = fe_engine::details::IntegrateScalarHelper<kind>::call(
-      integrator, f, type, ghost_type, filter_elements);
+  auto call = [&](auto && integral_type) {
+    return integrator.template integrate<ElementType(integral_type)>(
+        f, ghost_type, filter_elements);
+  };
+
+  auto integral = tuple_dispatch<ElementTypes_t<kind>>(call, type);
+
   AKANTU_DEBUG_OUT();
   return integral;
 }
 
 /* -------------------------------------------------------------------------- */
-/**
- * Helper class to be able to write a partial specialization on the element kind
- */
-namespace fe_engine {
-  namespace details {
-    template <ElementKind kind> struct IntegrateScalarOnOneElementHelper {};
-
-#define INTEGRATE(type)                                                        \
-  res = integrator.template integrate<type>(f, index, ghost_type);
-
-#define AKANTU_SPECIALIZE_INTEGRATE_SCALAR_ON_ONE_ELEMENT_HELPER(kind)         \
-  template <> struct IntegrateScalarOnOneElementHelper<kind> {                 \
-    template <class I>                                                         \
-    static Real call(const I & integrator, const Ref<const VectorXr> & f,      \
-                     const ElementType & type, Int index,                      \
-                     const GhostType & ghost_type) {                           \
-      Real res = 0.;                                                           \
-      AKANTU_BOOST_KIND_ELEMENT_SWITCH(INTEGRATE, kind);                       \
-      return res;                                                              \
-    }                                                                          \
-  };
-
-    AKANTU_BOOST_ALL_KIND(
-        AKANTU_SPECIALIZE_INTEGRATE_SCALAR_ON_ONE_ELEMENT_HELPER)
-
-#undef AKANTU_SPECIALIZE_INTEGRATE_SCALAR_ON_ONE_ELEMENT_HELPER
-#undef INTEGRATE
-  } // namespace details
-} // namespace fe_engine
-
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 Real FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::integrate(
     const Ref<const VectorXr> & f, const ElementType & type, Int index,
     const GhostType & ghost_type) const {
 
-  Real res = fe_engine::details::IntegrateScalarOnOneElementHelper<kind>::call(
-      integrator, f, type, index, ghost_type);
-  return res;
+  auto call = [&](auto && integral_type) {
+    return integrator.template integrate<ElementType(type)>(f, index,
+                                                            ghost_type);
+  };
+  return tuple_dispatch<ElementTypes_t<kind>>(call, type);
 }
 
 /* -------------------------------------------------------------------------- */
-/**
- * Helper class to be able to write a partial specialization on the element kind
- */
-namespace fe_engine {
-  namespace details {
-    template <ElementKind kind> struct IntegrateOnIntegrationPointsHelper {};
-
-#define INTEGRATE(type)                                                        \
-  integrator.template integrateOnIntegrationPoints<type>(                      \
-      f, intf, nb_degree_of_freedom, ghost_type, filter_elements);
-
-#define AKANTU_SPECIALIZE_INTEGRATE_ON_INTEGRATION_POINTS_HELPER(kind)         \
-  template <> struct IntegrateOnIntegrationPointsHelper<kind> {                \
-    template <class I>                                                         \
-    static void call(const I & integrator, const Array<Real> & f,              \
-                     Array<Real> & intf, Int nb_degree_of_freedom,             \
-                     const ElementType & type, const GhostType & ghost_type,   \
-                     const Array<Int> & filter_elements) {                     \
-      AKANTU_BOOST_KIND_ELEMENT_SWITCH(INTEGRATE, kind);                       \
-    }                                                                          \
-  };
-
-    AKANTU_BOOST_ALL_KIND(
-        AKANTU_SPECIALIZE_INTEGRATE_ON_INTEGRATION_POINTS_HELPER)
-
-#undef AKANTU_SPECIALIZE_INTEGRATE_ON_INTEGRATION_POINTS_HELPER
-#undef INTEGRATE
-  } // namespace details
-} // namespace fe_engine
-
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::
     integrateOnIntegrationPoints(const Array<Real> & f, Array<Real> & intf,
                                  Int nb_degree_of_freedom,
                                  const ElementType & type,
                                  const GhostType & ghost_type,
                                  const Array<Int> & filter_elements) const {
 
   auto nb_element = mesh.getNbElement(type, ghost_type);
   if (filter_elements != empty_filter)
     nb_element = filter_elements.size();
   auto nb_quadrature_points = getNbIntegrationPoints(type);
 #ifndef AKANTU_NDEBUG
   //   std::stringstream sstr; sstr << ghost_type;
   //   AKANTU_DEBUG_ASSERT(sstr.str() == nablauq.getTag(),
   //                  "The vector " << nablauq.getID() << " is not taged "
   //                  << ghost_type);
 
   AKANTU_DEBUG_ASSERT(f.size() == nb_element * nb_quadrature_points,
                       "The vector f(" << f.getID() << " size " << f.size()
                                       << ") has not the good size ("
                                       << nb_element << ").");
   AKANTU_DEBUG_ASSERT(f.getNbComponent() == nb_degree_of_freedom,
                       "The vector f("
                           << f.getID()
                           << ") has not the good number of component.");
   AKANTU_DEBUG_ASSERT(intf.getNbComponent() == nb_degree_of_freedom,
                       "The vector intf("
                           << intf.getID()
                           << ") has not the good number of component.");
 #endif
 
   intf.resize(nb_element * nb_quadrature_points);
-  fe_engine::details::IntegrateOnIntegrationPointsHelper<kind>::call(
-      integrator, f, intf, nb_degree_of_freedom, type, ghost_type,
-      filter_elements);
-}
-
-/* -------------------------------------------------------------------------- */
-/**
- * Helper class to be able to write a partial specialization on the element kind
- */
-namespace fe_engine {
-  namespace details {
-    template <ElementKind kind> struct InterpolateOnIntegrationPointsHelper {
-      template <class S>
-      static void call(const S &, const Array<Real> &, Array<Real> &, Int,
-                       const ElementType &, const GhostType &,
-                       const Array<Idx> &) {
-        AKANTU_TO_IMPLEMENT();
-      }
-    };
-
-#define INTERPOLATE(type)                                                      \
-  shape_functions.template interpolateOnIntegrationPoints<type>(               \
-      u, uq, nb_degree_of_freedom, ghost_type, filter_elements);
 
-#define AKANTU_SPECIALIZE_INTERPOLATE_ON_INTEGRATION_POINTS_HELPER(kind)       \
-  template <> struct InterpolateOnIntegrationPointsHelper<kind> {              \
-    template <class S>                                                         \
-    static void call(const S & shape_functions, const Array<Real> & u,         \
-                     Array<Real> & uq, Int nb_degree_of_freedom,               \
-                     const ElementType & type, const GhostType & ghost_type,   \
-                     const Array<Int> & filter_elements) {                     \
-      AKANTU_BOOST_KIND_ELEMENT_SWITCH(INTERPOLATE, kind);                     \
-    }                                                                          \
+  auto call = [&](auto && integral_type) {
+    integrator
+        .template integrateOnIntegrationPoints<ElementType(integral_type)>(
+            f, intf, nb_degree_of_freedom, ghost_type, filter_elements);
   };
-    AKANTU_BOOST_ALL_KIND_LIST(
-        AKANTU_SPECIALIZE_INTERPOLATE_ON_INTEGRATION_POINTS_HELPER,
-        AKANTU_FE_ENGINE_LIST_INTERPOLATE_ON_INTEGRATION_POINTS)
+  tuple_dispatch<ElementTypes_t<kind>>(call, type);
+}
 
-#undef AKANTU_SPECIALIZE_INTERPOLATE_ON_INTEGRATION_POINTS_HELPER
-#undef INTERPOLATE
-  } // namespace details
-} // namespace fe_engine
+#define DISPATCH_SHAPE_FUNCTIONS_HELPER(Func, ...)                             \
+  auto && call = [&](auto && enum_type) {                                      \
+    return shape_functions.template Func<ElementType(enum_type)>(__VA_ARGS__); \
+  };                                                                           \
+  tuple_dispatch<ElementTypes_t<kind>>(call, type);
 
+/* -------------------------------------------------------------------------- */
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::
     interpolateOnIntegrationPoints(const Array<Real> & u, Array<Real> & uq,
                                    Int nb_degree_of_freedom,
                                    const ElementType & type,
                                    const GhostType & ghost_type,
                                    const Array<Int> & filter_elements) const {
   AKANTU_DEBUG_IN();
 
   auto nb_points =
       shape_functions.getIntegrationPoints(type, ghost_type).cols();
   auto nb_element = mesh.getNbElement(type, ghost_type);
   if (filter_elements != empty_filter) {
     nb_element = filter_elements.size();
   }
 #ifndef AKANTU_NDEBUG
 
   AKANTU_DEBUG_ASSERT(u.size() == mesh.getNbNodes(),
                       "The vector u(" << u.getID()
                                       << ") has not the good size.");
   AKANTU_DEBUG_ASSERT(u.getNbComponent() == nb_degree_of_freedom,
                       "The vector u("
                           << u.getID()
                           << ") has not the good number of component.");
   AKANTU_DEBUG_ASSERT(uq.getNbComponent() == nb_degree_of_freedom,
                       "The vector uq("
                           << uq.getID()
                           << ") has not the good number of component.");
 #endif
 
   uq.resize(nb_element * nb_points);
 
-  fe_engine::details::InterpolateOnIntegrationPointsHelper<kind>::call(
-      shape_functions, u, uq, nb_degree_of_freedom, type, ghost_type,
-      filter_elements);
+  DISPATCH_SHAPE_FUNCTIONS_HELPER(interpolateOnIntegrationPoints, u, uq,
+                                  nb_degree_of_freedom, ghost_type,
+                                  filter_elements);
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::
     interpolateOnIntegrationPoints(
         const Array<Real> & u, ElementTypeMapArray<Real> & uq,
         const ElementTypeMapArray<Idx> * filter_elements) const {
   AKANTU_DEBUG_IN();
 
   const Array<Idx> * filter = nullptr;
 
   for (auto ghost_type : ghost_types) {
     for (auto & type : uq.elementTypes(_all_dimensions, ghost_type, kind)) {
       auto nb_quad_per_element = getNbIntegrationPoints(type, ghost_type);
 
       Int nb_element = 0;
 
       if (filter_elements != nullptr) {
         filter = &((*filter_elements)(type, ghost_type));
         nb_element = filter->size();
       } else {
         filter = &empty_filter;
         nb_element = mesh.getNbElement(type, ghost_type);
       }
 
       auto nb_tot_quad = nb_quad_per_element * nb_element;
 
       Array<Real> & quad = uq(type, ghost_type);
       quad.resize(nb_tot_quad);
 
       interpolateOnIntegrationPoints(u, quad, quad.getNbComponent(), type,
                                      ghost_type, *filter);
     }
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
-namespace fe_engine {
-  namespace details {
-    template <ElementKind kind> struct ComputeBtDHelper {};
-
-#define COMPUTE_BTD(type)                                                      \
-  shape_functions.template computeBtD<type>(Ds, BtDs, ghost_type,              \
-                                            filter_elements);
-
-#define AKANTU_SPECIALIZE_COMPUTE_BtD_HELPER(kind)                             \
-  template <> struct ComputeBtDHelper<kind> {                                  \
-    template <class S>                                                         \
-    static void call(const S & shape_functions, const Array<Real> & Ds,        \
-                     Array<Real> & BtDs, const ElementType & type,             \
-                     const GhostType & ghost_type,                             \
-                     const Array<Idx> & filter_elements) {                     \
-      AKANTU_BOOST_KIND_ELEMENT_SWITCH(COMPUTE_BTD, kind);                     \
-    }                                                                          \
-  };
-
-    AKANTU_BOOST_ALL_KIND(AKANTU_SPECIALIZE_COMPUTE_BtD_HELPER)
-
-#undef AKANTU_SPECIALIZE_COMPUTE_BtD_HELPER
-#undef COMPUTE_BTD
-  } // namespace details
-} // namespace fe_engine
-
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 inline void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::computeBtD(
     const Array<Real> & Ds, Array<Real> & BtDs, const ElementType & type,
     const GhostType & ghost_type, const Array<Idx> & filter_elements) const {
-  fe_engine::details::ComputeBtDHelper<kind>::call(
-      shape_functions, Ds, BtDs, type, ghost_type, filter_elements);
+
+  DISPATCH_SHAPE_FUNCTIONS_HELPER(computeBtD, Ds, BtDs, ghost_type,
+                                  filter_elements);
 }
 
 /* -------------------------------------------------------------------------- */
-namespace fe_engine {
-  namespace details {
-    template <ElementKind kind> struct ComputeBtDBHelper {};
-
-#define COMPUTE_BTDB(type)                                                     \
-  shape_functions.template computeBtDB<type>(Ds, BtDBs, order_d, ghost_type,   \
-                                             filter_elements);
-
-#define AKANTU_SPECIALIZE_COMPUTE_BtDB_HELPER(kind)                            \
-  template <> struct ComputeBtDBHelper<kind> {                                 \
-    template <class S>                                                         \
-    static void call(const S & shape_functions, const Array<Real> & Ds,        \
-                     Array<Real> & BtDBs, Int order_d,                         \
-                     const ElementType & type, const GhostType & ghost_type,   \
-                     const Array<Idx> & filter_elements) {                     \
-      AKANTU_BOOST_KIND_ELEMENT_SWITCH(COMPUTE_BTDB, kind);                    \
-    }                                                                          \
-  };
-
-    AKANTU_BOOST_ALL_KIND(AKANTU_SPECIALIZE_COMPUTE_BtDB_HELPER)
-
-#undef AKANTU_SPECIALIZE_COMPUTE_BtDB_HELPER
-#undef COMPUTE_BTDB
-  } // namespace details
-} // namespace fe_engine
-
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 inline void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::computeBtDB(
     const Array<Real> & Ds, Array<Real> & BtDBs, Int order_d,
     const ElementType & type, const GhostType & ghost_type,
     const Array<Idx> & filter_elements) const {
-  fe_engine::details::ComputeBtDBHelper<kind>::call(
-      shape_functions, Ds, BtDBs, order_d, type, ghost_type, filter_elements);
+  DISPATCH_SHAPE_FUNCTIONS_HELPER(computeBtDB, Ds, BtDBs, order_d, ghost_type,
+                                  filter_elements);
 }
 
 /* -------------------------------------------------------------------------- */
-namespace fe_engine {
-  namespace details {
-    template <ElementKind kind> struct ComputeNtbNHelper {};
-
-#define COMPUTE_NtbN(type)                                                     \
-  shape_functions.template computeNtbN<type>(bs, NtbNs, ghost_type,            \
-                                             filter_elements);
-
-#define AKANTU_SPECIALIZE_COMPUTE_NtbN_HELPER(kind)                            \
-  template <> struct ComputeNtbNHelper<kind> {                                 \
-    template <class S>                                                         \
-    static void call(const S & shape_functions, const Array<Real> & bs,        \
-                     Array<Real> & NtbNs, ElementType type,                    \
-                     GhostType ghost_type,                                     \
-                     const Array<UInt> & filter_elements) {                    \
-      AKANTU_BOOST_KIND_ELEMENT_SWITCH(COMPUTE_NtbN, kind);                    \
-    }                                                                          \
-  };
-
-    AKANTU_BOOST_ALL_KIND(AKANTU_SPECIALIZE_COMPUTE_NtbN_HELPER)
-
-#undef AKANTU_SPECIALIZE_COMPUTE_NtbN_HELPER
-#undef COMPUTE_NtbN
-  } // namespace details
-} // namespace fe_engine
-
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 inline void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::computeNtbN(
     const Array<Real> & bs, Array<Real> & NtbNs, ElementType type,
-    GhostType ghost_type, const Array<UInt> & filter_elements) const {
-  fe_engine::details::ComputeNtbNHelper<kind>::call(
-      shape_functions, bs, NtbNs, type, ghost_type, filter_elements);
+    GhostType ghost_type, const Array<Idx> & filter_elements) const {
+
+  DISPATCH_SHAPE_FUNCTIONS_HELPER(computeNtbN, bs, NtbNs, type, ghost_type,
+                                  filter_elements);
 }
 
 /* -------------------------------------------------------------------------- */
-namespace fe_engine {
-  namespace details {
-    template <ElementKind kind> struct ComputeNtbHelper {};
-
-#define COMPUTE_Ntb(type)                                                      \
-  shape_functions.template computeNtb<type>(bs, Ntbs, ghost_type,              \
-                                            filter_elements);
-
-#define AKANTU_SPECIALIZE_COMPUTE_Ntb_HELPER(kind)                             \
-  template <> struct ComputeNtbHelper<kind> {                                  \
-    template <class S>                                                         \
-    static void call(const S & shape_functions, const Array<Real> & bs,        \
-                     Array<Real> & Ntbs, const ElementType & type,             \
-                     const GhostType & ghost_type,                             \
-                     const Array<Idx> & filter_elements) {                     \
-      AKANTU_BOOST_KIND_ELEMENT_SWITCH(COMPUTE_Ntb, kind);                     \
-    }                                                                          \
-  };
-
-    AKANTU_BOOST_ALL_KIND(AKANTU_SPECIALIZE_COMPUTE_Ntb_HELPER)
-
-#undef AKANTU_SPECIALIZE_COMPUTE_Ntb_HELPER
-#undef COMPUTE_Ntb
-  } // namespace details
-} // namespace fe_engine
-
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 inline void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::computeNtb(
     const Array<Real> & bs, Array<Real> & Ntbs, const ElementType & type,
     const GhostType & ghost_type, const Array<Idx> & filter_elements) const {
-  fe_engine::details::ComputeNtbHelper<kind>::call(
-      shape_functions, bs, Ntbs, type, ghost_type, filter_elements);
+
+  DISPATCH_SHAPE_FUNCTIONS_HELPER(computeNtb, bs, Ntbs, ghost_type,
+                                  filter_elements);
 }
 
 /* -------------------------------------------------------------------------- */
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 inline void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::
     computeIntegrationPointsCoordinates(
         ElementTypeMapArray<Real> & quadrature_points_coordinates,
         const ElementTypeMapArray<Idx> * filter_elements) const {
 
   const Array<Real> & nodes_coordinates = mesh.getNodes();
 
   interpolateOnIntegrationPoints(
       nodes_coordinates, quadrature_points_coordinates, filter_elements);
 }
 
 /* -------------------------------------------------------------------------- */
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 inline void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::
     computeIntegrationPointsCoordinates(
         Array<Real> & quadrature_points_coordinates, const ElementType & type,
         const GhostType & ghost_type,
         const Array<Idx> & filter_elements) const {
 
   const Array<Real> & nodes_coordinates = mesh.getNodes();
 
   auto spatial_dimension = mesh.getSpatialDimension();
 
   interpolateOnIntegrationPoints(
       nodes_coordinates, quadrature_points_coordinates, spatial_dimension, type,
       ghost_type, filter_elements);
 }
 
 /* -------------------------------------------------------------------------- */
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 inline void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::
     initElementalFieldInterpolationFromIntegrationPoints(
         const ElementTypeMapArray<Real> & interpolation_points_coordinates,
         ElementTypeMapArray<Real> & interpolation_points_coordinates_matrices,
         ElementTypeMapArray<Real> & quad_points_coordinates_inv_matrices,
         const ElementTypeMapArray<Idx> * element_filter) const {
 
   AKANTU_DEBUG_IN();
 
   auto spatial_dimension = this->mesh.getSpatialDimension();
 
   ElementTypeMapArray<Real> quadrature_points_coordinates(
       "quadrature_points_coordinates_for_interpolation", getID());
 
   quadrature_points_coordinates.initialize(*this,
                                            _nb_component = spatial_dimension);
 
   computeIntegrationPointsCoordinates(quadrature_points_coordinates,
                                       element_filter);
   shape_functions.initElementalFieldInterpolationFromIntegrationPoints(
       interpolation_points_coordinates,
       interpolation_points_coordinates_matrices,
       quad_points_coordinates_inv_matrices, quadrature_points_coordinates,
       element_filter);
 }
 
 /* -------------------------------------------------------------------------- */
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 inline void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::
     interpolateElementalFieldFromIntegrationPoints(
         const ElementTypeMapArray<Real> & field,
         const ElementTypeMapArray<Real> & interpolation_points_coordinates,
         ElementTypeMapArray<Real> & result, const GhostType ghost_type,
         const ElementTypeMapArray<Idx> * element_filter) const {
 
   ElementTypeMapArray<Real> interpolation_points_coordinates_matrices(
       "interpolation_points_coordinates_matrices", id);
   ElementTypeMapArray<Real> quad_points_coordinates_inv_matrices(
       "quad_points_coordinates_inv_matrices", id);
 
   initElementalFieldInterpolationFromIntegrationPoints(
       interpolation_points_coordinates,
       interpolation_points_coordinates_matrices,
       quad_points_coordinates_inv_matrices, element_filter);
 
   interpolateElementalFieldFromIntegrationPoints(
       field, interpolation_points_coordinates_matrices,
       quad_points_coordinates_inv_matrices, result, ghost_type, element_filter);
 }
 
 /* -------------------------------------------------------------------------- */
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 inline void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::
     interpolateElementalFieldFromIntegrationPoints(
         const ElementTypeMapArray<Real> & field,
         const ElementTypeMapArray<Real> &
             interpolation_points_coordinates_matrices,
         const ElementTypeMapArray<Real> & quad_points_coordinates_inv_matrices,
         ElementTypeMapArray<Real> & result, const GhostType ghost_type,
         const ElementTypeMapArray<Idx> * element_filter) const {
 
   shape_functions.interpolateElementalFieldFromIntegrationPoints(
       field, interpolation_points_coordinates_matrices,
       quad_points_coordinates_inv_matrices, result, ghost_type, element_filter);
 }
 
 /* -------------------------------------------------------------------------- */
 /**
  * Helper class to be able to write a partial specialization on the element kind
  */
 namespace fe_engine {
   namespace details {
     template <ElementKind kind> struct InterpolateHelper {
       template <class S>
       static void call(const S &, const Ref<const VectorXr> &, Int,
                        const Ref<const MatrixXr> &, Ref<VectorXr>,
                        const ElementType &, const GhostType &) {
         AKANTU_TO_IMPLEMENT();
       }
     };
 
 #define INTERPOLATE(type)                                                      \
   shape_functions.template interpolate<type>(                                  \
       real_coords, element, nodal_values, interpolated, ghost_type);
 
 #define AKANTU_SPECIALIZE_INTERPOLATE_HELPER(kind)                             \
   template <> struct InterpolateHelper<kind> {                                 \
     template <class S>                                                         \
     static void call(const S & shape_functions,                                \
                      const Ref<const VectorXr> & real_coords, Int element,     \
                      const Ref<const MatrixXr> & nodal_values,                 \
                      Ref<VectorXr> interpolated, const ElementType & type,     \
                      const GhostType & ghost_type) {                           \
       AKANTU_BOOST_KIND_ELEMENT_SWITCH(INTERPOLATE, kind);                     \
     }                                                                          \
   };
 
     AKANTU_BOOST_ALL_KIND_LIST(AKANTU_SPECIALIZE_INTERPOLATE_HELPER,
                                AKANTU_FE_ENGINE_LIST_INTERPOLATE)
 
 #undef AKANTU_SPECIALIZE_INTERPOLATE_HELPER
 #undef INTERPOLATE
   } // namespace details
 } // namespace fe_engine
 
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 inline void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::interpolate(
     const Ref<const VectorXr> & real_coords,
     const Ref<const MatrixXr> & nodal_values, Ref<VectorXr> interpolated,
     const Element & element) const {
 
   AKANTU_DEBUG_IN();
 
-  fe_engine::details::InterpolateHelper<kind>::call(
-      shape_functions, real_coords, element.element, nodal_values, interpolated,
-      element.type, element.ghost_type);
+  auto type = element.type;
+  DISPATCH_SHAPE_FUNCTIONS_HELPER(interpolate, real_coords, element,
+                                  nodal_values, interpolated,
+                                  element.ghost_type);
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::
     computeNormalsOnIntegrationPoints(GhostType ghost_type) {
   AKANTU_DEBUG_IN();
 
   computeNormalsOnIntegrationPoints(mesh.getNodes(), ghost_type);
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::
     computeNormalsOnIntegrationPoints(const Array<Real> & field,
                                       GhostType ghost_type) {
   AKANTU_DEBUG_IN();
 
   //  Real * coord = mesh.getNodes().data();
   auto spatial_dimension = mesh.getSpatialDimension();
 
   // allocate the normal arrays
   normals_on_integration_points.initialize(
       *this, _nb_component = spatial_dimension,
       _spatial_dimension = element_dimension, _ghost_type = ghost_type,
       _element_kind = kind);
 
   // loop over the type to build the normals
   for (auto & type : mesh.elementTypes(element_dimension, ghost_type, kind)) {
     auto & normals_on_quad = normals_on_integration_points(type, ghost_type);
     computeNormalsOnIntegrationPoints(field, normals_on_quad, type, ghost_type);
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
-/**
- * Helper class to be able to write a partial specialization on the element kind
- */
-namespace fe_engine {
-  namespace details {
-    template <ElementKind kind> struct ComputeNormalsOnIntegrationPoints {
-      template <template <ElementKind, class> class I,
-                template <ElementKind> class S, ElementKind k, class IOF>
-      static void call(const FEEngineTemplate<I, S, k, IOF> & /*unused*/,
-                       const Array<Real> & /*unused*/, Array<Real> & /*unused*/,
-                       ElementType /*unused*/, GhostType /*unused*/) {
-        AKANTU_TO_IMPLEMENT();
-      }
-    };
+template <template <ElementKind, class> class I, template <ElementKind> class S,
+          ElementKind kind, class IntegrationOrderFunctor>
+template <ElementType type, ElementKind kind_,
+          std::enable_if_t<kind_ == _ek_regular and type != _point_1> *>
+void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::
+    computeNormalsOnIntegrationPoints(const Array<Real> & field,
+                                      Array<Real> & normal,
+                                      const GhostType & ghost_type) const {
+  AKANTU_DEBUG_IN();
 
-#define COMPUTE_NORMALS_ON_INTEGRATION_POINTS(type)                            \
-  fem.template computeNormalsOnIntegrationPoints<type>(field, normal,          \
-                                                       ghost_type);
-
-#define AKANTU_SPECIALIZE_COMPUTE_NORMALS_ON_INTEGRATION_POINTS(kind)          \
-  template <> struct ComputeNormalsOnIntegrationPoints<kind> {                 \
-    template <template <ElementKind, class> class I,                           \
-              template <ElementKind> class S, ElementKind k, class IOF>        \
-    static void call(const FEEngineTemplate<I, S, k, IOF> & fem,               \
-                     const Array<Real> & field, Array<Real> & normal,          \
-                     ElementType type, GhostType ghost_type) {                 \
-      AKANTU_BOOST_KIND_ELEMENT_SWITCH(COMPUTE_NORMALS_ON_INTEGRATION_POINTS,  \
-                                       kind);                                  \
-    }                                                                          \
-  };
+  auto spatial_dimension = mesh.getSpatialDimension();
+  constexpr auto nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
+  auto nb_points = getNbIntegrationPoints(type, ghost_type);
 
-    AKANTU_BOOST_ALL_KIND_LIST(
-        AKANTU_SPECIALIZE_COMPUTE_NORMALS_ON_INTEGRATION_POINTS,
-        AKANTU_FE_ENGINE_LIST_COMPUTE_NORMALS_ON_INTEGRATION_POINTS)
+  auto nb_element = mesh.getConnectivity(type, ghost_type).size();
+  normal.resize(nb_element * nb_points);
 
-#undef AKANTU_SPECIALIZE_COMPUTE_NORMALS_ON_INTEGRATION_POINTS
-#undef COMPUTE_NORMALS_ON_INTEGRATION_POINTS
-  } // namespace details
-} // namespace fe_engine
+  Array<Real> f_el(0, spatial_dimension * nb_nodes_per_element);
+  FEEngine::extractNodalToElementField(mesh, field, f_el, type, ghost_type);
 
+  const auto & quads =
+      integrator.template getIntegrationPoints<type>(ghost_type);
+
+  for (auto && data : zip(make_view(normal, spatial_dimension, nb_points),
+                          make_view<Eigen::Dynamic, nb_nodes_per_element>(
+                              f_el, spatial_dimension, nb_nodes_per_element))) {
+    ElementClass<type>::computeNormalsOnNaturalCoordinates(
+        quads, std::get<1>(data), std::get<0>(data));
+  }
+
+  AKANTU_DEBUG_OUT();
+}
+
+/* -------------------------------------------------------------------------- */
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
-void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::
-    computeNormalsOnIntegrationPoints(const Array<Real> & field,
-                                      Array<Real> & normal, ElementType type,
-                                      GhostType ghost_type) const {
-  fe_engine::details::ComputeNormalsOnIntegrationPoints<kind>::call(
-      *this, field, normal, type, ghost_type);
+template <ElementType type, ElementKind kind_,
+          std::enable_if_t<kind_ == _ek_regular and type == _point_1> *>
+inline void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::
+    computeNormalsOnIntegrationPoints(const Array<Real> &, Array<Real> & normal,
+                                      const GhostType & ghost_type) const {
+  AKANTU_DEBUG_IN();
+
+  AKANTU_DEBUG_ASSERT(mesh.getSpatialDimension() == 1,
+                      "Mesh dimension must be 1 to compute normals on points!");
+
+  auto spatial_dimension = mesh.getSpatialDimension();
+  // Int nb_nodes_per_element  = Mesh::getNbNodesPerElement(type);
+  auto nb_points = getNbIntegrationPoints(type, ghost_type);
+  const auto & connectivity = mesh.getConnectivity(type, ghost_type);
+  auto nb_element = connectivity.size();
+
+  normal.resize(nb_element * nb_points);
+  auto normals_on_quad =
+      make_view(normal, spatial_dimension, nb_points).begin();
+  const auto & segments = mesh.getElementToSubelement(type, ghost_type);
+  const auto & coords = mesh.getNodes();
+
+  const Mesh * mesh_segment;
+  if (mesh.isMeshFacets())
+    mesh_segment = &(mesh.getMeshParent());
+  else
+    mesh_segment = &mesh;
+
+  for (Idx elem = 0; elem < nb_element; ++elem) {
+    auto nb_segment = segments(elem).size();
+    AKANTU_DEBUG_ASSERT(
+        nb_segment > 0,
+        "Impossible to compute a normal on a point connected to 0 segments");
+
+    Real normal_value = 1;
+    if (nb_segment == 1) {
+      auto point = connectivity(elem);
+      const auto segment = segments(elem)[0];
+      const auto & segment_connectivity =
+          mesh_segment->getConnectivity(segment.type, segment.ghost_type);
+      Vector<Idx> segment_points = segment_connectivity.begin(
+          Mesh::getNbNodesPerElement(segment.type))[segment.element];
+      Real difference;
+      if (segment_points(0) == point) {
+        difference = coords(elem) - coords(segment_points(1));
+      } else {
+        difference = coords(elem) - coords(segment_points(0));
+      }
+
+      normal_value = difference / std::abs(difference);
+    }
+
+    for (Idx n(0); n < nb_points; ++n) {
+      (*normals_on_quad)(0, n) = normal_value;
+    }
+    ++normals_on_quad;
+  }
+
+  AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
-template <ElementType type>
 void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::
     computeNormalsOnIntegrationPoints(const Array<Real> & field,
                                       Array<Real> & normal,
-                                      GhostType ghost_type) const {
-  AKANTU_DEBUG_IN();
+                                      const ElementType & type,
+                                      const GhostType & ghost_type) const {
 
-  static_if(aka::bool_constant<type == _point_1>{})
-      .then([&](auto && field) {
-        computeNormalsOnIntegrationPointsPoint1(field, normal, ghost_type);
-      })
-      .else_([&](auto && field) {
-        auto spatial_dimension = mesh.getSpatialDimension();
-        auto nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
-        auto nb_points = getNbIntegrationPoints(type, ghost_type);
-
-        auto nb_element = mesh.getConnectivity(type, ghost_type).size();
-        normal.resize(nb_element * nb_points);
-
-        Array<Real> f_el(0, spatial_dimension * nb_nodes_per_element);
-        FEEngine::extractNodalToElementField(mesh, field, f_el, type,
-                                             ghost_type);
-
-        const auto & quads =
-            integrator.template getIntegrationPoints<type>(ghost_type);
-
-        for (auto && data :
-             zip(make_view(normal, spatial_dimension, nb_points),
-                 make_view(f_el, spatial_dimension, nb_nodes_per_element))) {
-          ElementClass<type>::computeNormalsOnNaturalCoordinates(
-              quads, std::get<1>(data), std::get<0>(data));
-        }
-      })(field);
-
-  AKANTU_DEBUG_OUT();
+  DISPATCH_HELPER(computeNormalsOnIntegrationPoints, field, normal, ghost_type);
 }
 
 /* -------------------------------------------------------------------------- */
 /**
  * Helper class to be able to write a partial specialization on the element kind
  */
 template <ElementKind kind> struct InverseMapHelper {
   template <class S>
   static void
   call(const S & /*shape_functions*/, const Vector<Real> & /*real_coords*/,
        Int /*element*/, const ElementType & /*type*/,
        Vector<Real> & /*natural_coords*/, const GhostType & /*ghost_type*/) {
     AKANTU_TO_IMPLEMENT();
   }
 };
 
 #define INVERSE_MAP(type)                                                      \
   shape_functions.template inverseMap<type>(real_coords, element,              \
                                             natural_coords, ghost_type);
 
 #define AKANTU_SPECIALIZE_INVERSE_MAP_HELPER(kind)                             \
   template <> struct InverseMapHelper<kind> {                                  \
     template <class S>                                                         \
     static void call(const S & shape_functions,                                \
                      const Vector<Real> & real_coords, Int element,            \
                      const ElementType & type, Vector<Real> & natural_coords,  \
                      const GhostType & ghost_type) {                           \
       AKANTU_BOOST_KIND_ELEMENT_SWITCH(INVERSE_MAP, kind);                     \
     }                                                                          \
   };
 
 AKANTU_BOOST_ALL_KIND_LIST(AKANTU_SPECIALIZE_INVERSE_MAP_HELPER,
                            AKANTU_FE_ENGINE_LIST_INVERSE_MAP)
 
 #undef AKANTU_SPECIALIZE_INVERSE_MAP_HELPER
 #undef INVERSE_MAP
 
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 inline void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::inverseMap(
     const Vector<Real> & real_coords, Int element, const ElementType & type,
     Vector<Real> & natural_coords, const GhostType & ghost_type) const {
 
   AKANTU_DEBUG_IN();
 
   InverseMapHelper<kind>::call(shape_functions, real_coords, element, type,
                                natural_coords, ghost_type);
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 /**
  * Helper class to be able to write a partial specialization on the element kind
  */
 namespace fe_engine {
   namespace details {
     template <ElementKind kind> struct ContainsHelper {
       template <class S>
       static void call(const S &, const Vector<Real> &, Int,
                        const ElementType &, const GhostType &) {
         AKANTU_TO_IMPLEMENT();
       }
     };
 
 #define CONTAINS(type)                                                         \
   contain = shape_functions.template contains<type>(real_coords, element,      \
                                                     ghost_type);
 
 #define AKANTU_SPECIALIZE_CONTAINS_HELPER(kind)                                \
   template <> struct ContainsHelper<kind> {                                    \
     template <template <ElementKind> class S, ElementKind k>                   \
     static bool call(const S<k> & shape_functions,                             \
                      const Vector<Real> & real_coords, Int element,            \
                      const ElementType & type, const GhostType & ghost_type) { \
       bool contain = false;                                                    \
       AKANTU_BOOST_KIND_ELEMENT_SWITCH(CONTAINS, kind);                        \
       return contain;                                                          \
     }                                                                          \
   };
 
     AKANTU_BOOST_ALL_KIND_LIST(AKANTU_SPECIALIZE_CONTAINS_HELPER,
                                AKANTU_FE_ENGINE_LIST_CONTAINS)
 
 #undef AKANTU_SPECIALIZE_CONTAINS_HELPER
 #undef CONTAINS
   } // namespace details
 } // namespace fe_engine
 
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 inline bool FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::contains(
     const Vector<Real> & real_coords, Int element, const ElementType & type,
     const GhostType & ghost_type) const {
   return fe_engine::details::ContainsHelper<kind>::call(
       shape_functions, real_coords, element, type, ghost_type);
 }
 
 /* -------------------------------------------------------------------------- */
 /**
  * Helper class to be able to write a partial specialization on the element kind
  */
 namespace fe_engine {
   namespace details {
     template <ElementKind kind> struct ComputeShapesHelper {
       template <class S>
       static void call(const S &, const Ref<const VectorXr> &, Int,
                        const ElementType, Ref<VectorXr>, const GhostType &) {
         AKANTU_TO_IMPLEMENT();
       }
     };
 
 #define COMPUTE_SHAPES(type)                                                   \
   shape_functions.template computeShapes<type>(real_coords, element, shapes,   \
                                                ghost_type);
 
 #define AKANTU_SPECIALIZE_COMPUTE_SHAPES_HELPER(kind)                          \
   template <> struct ComputeShapesHelper<kind> {                               \
     template <class S>                                                         \
     static void call(const S & shape_functions,                                \
                      const Ref<const VectorXr> & real_coords, Int element,     \
                      const ElementType type, Ref<VectorXr> & shapes,           \
                      const GhostType & ghost_type) {                           \
       AKANTU_BOOST_KIND_ELEMENT_SWITCH(COMPUTE_SHAPES, kind);                  \
     }                                                                          \
   };
 
     AKANTU_BOOST_ALL_KIND_LIST(AKANTU_SPECIALIZE_COMPUTE_SHAPES_HELPER,
                                AKANTU_FE_ENGINE_LIST_COMPUTE_SHAPES)
 
 #undef AKANTU_SPECIALIZE_COMPUTE_SHAPES_HELPER
 #undef COMPUTE_SHAPES
   } // namespace details
 } // namespace fe_engine
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 inline void
 FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::computeShapes(
     const Ref<const VectorXr> & real_coords, Int element,
     const ElementType & type, Ref<VectorXr> shapes,
     const GhostType & ghost_type) const {
 
   AKANTU_DEBUG_IN();
 
   fe_engine::details::ComputeShapesHelper<kind>::call(
       shape_functions, real_coords, element, type, shapes, ghost_type);
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 /**
  * Helper class to be able to write a partial specialization on the element kind
  */
 namespace fe_engine {
   namespace details {
     template <ElementKind kind> struct ComputeShapeDerivativesHelper {
       template <class S>
       static void call(const S &, const Ref<const VectorXr> &, Int,
                        const ElementType, Ref<MatrixXr>, const GhostType &) {
         AKANTU_TO_IMPLEMENT();
       }
     };
 
 #define COMPUTE_SHAPE_DERIVATIVES(type)                                        \
   VectorProxy<const Real> coords_mat(real_coords.data(),                       \
                                      shape_derivatives.rows(), 1);             \
   Tensor3Proxy<Real> shapesd_tensor(shape_derivatives.data(),                  \
                                     shape_derivatives.rows(),                  \
                                     shape_derivatives.cols(), 1);              \
   shape_functions.template computeShapeDerivatives<type>(                      \
       coords_mat, element, shapesd_tensor, ghost_type);
 
 #define AKANTU_SPECIALIZE_COMPUTE_SHAPE_DERIVATIVES_HELPER(kind)               \
   template <> struct ComputeShapeDerivativesHelper<kind> {                     \
     template <class S>                                                         \
     static void call(const S & shape_functions,                                \
                      const Ref<const VectorXr> & real_coords, Int element,     \
                      const ElementType type, Ref<MatrixXr> shape_derivatives,  \
                      const GhostType & ghost_type) {                           \
       AKANTU_BOOST_KIND_ELEMENT_SWITCH(COMPUTE_SHAPE_DERIVATIVES, kind);       \
     }                                                                          \
   };
 
     AKANTU_BOOST_ALL_KIND_LIST(
         AKANTU_SPECIALIZE_COMPUTE_SHAPE_DERIVATIVES_HELPER,
         AKANTU_FE_ENGINE_LIST_COMPUTE_SHAPES_DERIVATIVES)
 
 #undef AKANTU_SPECIALIZE_COMPUTE_SHAPE_DERIVATIVES_HELPER
 #undef COMPUTE_SHAPE_DERIVATIVES
   } // namespace details
 } // namespace fe_engine
 
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 inline void
 FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::computeShapeDerivatives(
     const Ref<const VectorXr> & real_coords, Int element,
     const ElementType & type, Ref<MatrixXr> shape_derivatives,
     const GhostType & ghost_type) const {
   AKANTU_DEBUG_IN();
 
   fe_engine::details::ComputeShapeDerivativesHelper<kind>::call(
       shape_functions, real_coords, element, type, shape_derivatives,
       ghost_type);
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 /**
  * Helper class to be able to write a partial specialization on the element kind
  */
 namespace fe_engine {
   namespace details {
     template <ElementKind kind> struct GetNbIntegrationPointsHelper {};
 
 #define GET_NB_INTEGRATION_POINTS(type)                                        \
   nb_quad_points = integrator.template getNbIntegrationPoints<type>(ghost_type);
 
 #define AKANTU_SPECIALIZE_GET_NB_INTEGRATION_POINTS_HELPER(kind)               \
   template <> struct GetNbIntegrationPointsHelper<kind> {                      \
     template <template <ElementKind, class> class I, ElementKind k, class IOF> \
     static Int call(const I<k, IOF> & integrator, const ElementType type,      \
                     const GhostType & ghost_type) {                            \
       Int nb_quad_points = 0;                                                  \
       AKANTU_BOOST_KIND_ELEMENT_SWITCH(GET_NB_INTEGRATION_POINTS, kind);       \
       return nb_quad_points;                                                   \
     }                                                                          \
   };
 
     AKANTU_BOOST_ALL_KIND(AKANTU_SPECIALIZE_GET_NB_INTEGRATION_POINTS_HELPER)
 
 #undef AKANTU_SPECIALIZE_GET_NB_INTEGRATION_POINTS_HELPER
 #undef GET_NB_INTEGRATION
   } // namespace details
 } // namespace fe_engine
 
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 inline Int
 FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::getNbIntegrationPoints(
     ElementType type, GhostType ghost_type) const {
   return fe_engine::details::GetNbIntegrationPointsHelper<kind>::call(
       integrator, type, ghost_type);
 }
 
 /* -------------------------------------------------------------------------- */
 /**
  * Helper class to be able to write a partial specialization on the element kind
  */
 namespace fe_engine {
   namespace details {
     template <ElementKind kind> struct GetShapesHelper {};
 
 #define GET_SHAPES(type) ret = &(shape_functions.getShapes(type, ghost_type));
 
 #define AKANTU_SPECIALIZE_GET_SHAPES_HELPER(kind)                              \
   template <> struct GetShapesHelper<kind> {                                   \
     template <class S>                                                         \
     static const Array<Real> & call(const S & shape_functions,                 \
                                     const ElementType type,                    \
                                     GhostType ghost_type) {                    \
       const Array<Real> * ret = NULL;                                          \
       AKANTU_BOOST_KIND_ELEMENT_SWITCH(GET_SHAPES, kind);                      \
       return *ret;                                                             \
     }                                                                          \
   };
 
     AKANTU_BOOST_ALL_KIND(AKANTU_SPECIALIZE_GET_SHAPES_HELPER)
 
 #undef AKANTU_SPECIALIZE_GET_SHAPES_HELPER
 #undef GET_SHAPES
   } // namespace details
 } // namespace fe_engine
 
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 inline const Array<Real> &
 FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::getShapes(
     const ElementType & type, const GhostType & ghost_type,
     __attribute__((unused)) Int id) const {
   return fe_engine::details::GetShapesHelper<kind>::call(shape_functions, type,
                                                          ghost_type);
 }
 
 /* -------------------------------------------------------------------------- */
 /**
  * Helper class to be able to write a partial specialization on the element kind
  */
 namespace fe_engine {
   namespace details {
     template <ElementKind kind> struct GetShapesDerivativesHelper {
       template <template <ElementKind> class S, ElementKind k>
       static const Array<Real> & call(const S<k> &, const ElementType &,
                                       const GhostType &, Int) {
         AKANTU_TO_IMPLEMENT();
       }
     };
 
 #define GET_SHAPES_DERIVATIVES(type)                                           \
   ret = &(shape_functions.getShapesDerivatives(type, ghost_type));
 
 #define AKANTU_SPECIALIZE_GET_SHAPES_DERIVATIVES_HELPER(kind)                  \
   template <> struct GetShapesDerivativesHelper<kind> {                        \
     template <template <ElementKind> class S, ElementKind k>                   \
     static const Array<Real> &                                                 \
     call(const S<k> & shape_functions, const ElementType type,                 \
          const GhostType & ghost_type, __attribute__((unused)) Int id) {       \
       const Array<Real> * ret = NULL;                                          \
       AKANTU_BOOST_KIND_ELEMENT_SWITCH(GET_SHAPES_DERIVATIVES, kind);          \
       return *ret;                                                             \
     }                                                                          \
   };
 
     AKANTU_BOOST_ALL_KIND_LIST(AKANTU_SPECIALIZE_GET_SHAPES_DERIVATIVES_HELPER,
                                AKANTU_FE_ENGINE_LIST_GET_SHAPES_DERIVATIVES)
 
 #undef AKANTU_SPECIALIZE_GET_SHAPE_DERIVATIVES_HELPER
 #undef GET_SHAPES_DERIVATIVES
   } // namespace details
 } // namespace fe_engine
 
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 inline const Array<Real> &
 FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::getShapesDerivatives(
     const ElementType & type, const GhostType & ghost_type,
     __attribute__((unused)) Int id) const {
   return fe_engine::details::GetShapesDerivativesHelper<kind>::call(
       shape_functions, type, ghost_type, id);
 }
 
 /* -------------------------------------------------------------------------- */
 /**
  * Helper class to be able to write a partial specialization on the element kind
  */
 namespace fe_engine {
   namespace details {
     template <ElementKind kind> struct GetIntegrationPointsHelper {};
 
 #define GET_INTEGRATION_POINTS(type)                                           \
   ret = &(integrator.template getIntegrationPoints<type>(ghost_type));
 
 #define AKANTU_SPECIALIZE_GET_INTEGRATION_POINTS_HELPER(kind)                  \
   template <> struct GetIntegrationPointsHelper<kind> {                        \
     template <template <ElementKind, class> class I, ElementKind k, class IOF> \
     static const Matrix<Real> & call(const I<k, IOF> & integrator,             \
                                      const ElementType type,                   \
                                      GhostType ghost_type) {                   \
       const Matrix<Real> * ret = NULL;                                         \
       AKANTU_BOOST_KIND_ELEMENT_SWITCH(GET_INTEGRATION_POINTS, kind);          \
       return *ret;                                                             \
     }                                                                          \
   };
 
     AKANTU_BOOST_ALL_KIND(AKANTU_SPECIALIZE_GET_INTEGRATION_POINTS_HELPER)
 
 #undef AKANTU_SPECIALIZE_GET_INTEGRATION_POINTS_HELPER
 #undef GET_INTEGRATION_POINTS
   } // namespace details
 } // namespace fe_engine
 
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 inline const Matrix<Real> &
 FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::getIntegrationPoints(
     ElementType type, GhostType ghost_type) const {
   return fe_engine::details::GetIntegrationPointsHelper<kind>::call(
       integrator, type, ghost_type);
 }
 
 /* -------------------------------------------------------------------------- */
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::printself(
     std::ostream & stream, int indent) const {
   std::string space(indent, AKANTU_INDENT);
 
   stream << space << "FEEngineTemplate [" << std::endl;
   stream << space << " + parent [" << std::endl;
   FEEngine::printself(stream, indent + 3);
   stream << space << "   ]" << std::endl;
   stream << space << " + shape functions [" << std::endl;
   shape_functions.printself(stream, indent + 3);
   stream << space << "   ]" << std::endl;
   stream << space << " + integrator [" << std::endl;
   integrator.printself(stream, indent + 3);
   stream << space << "   ]" << std::endl;
   stream << space << "]" << std::endl;
 }
 
 /* -------------------------------------------------------------------------- */
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::onElementsAdded(
     const Array<Element> & new_elements, const NewElementsEvent & /*unused*/) {
   integrator.onElementsAdded(new_elements);
   shape_functions.onElementsAdded(new_elements);
 }
 
 /* -------------------------------------------------------------------------- */
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::onElementsRemoved(
     const Array<Element> & /*unused*/,
     const ElementTypeMapArray<UInt> & /*unused*/,
     const RemovedElementsEvent & /*unused*/) {}
 
 /* -------------------------------------------------------------------------- */
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::onElementsChanged(
     const Array<Element> & /*unused*/, const Array<Element> & /*unused*/,
     const ElementTypeMapArray<UInt> & /*unused*/,
     const ChangedElementsEvent & /*unused*/) {}
 
-/* -------------------------------------------------------------------------- */
-template <template <ElementKind, class> class I, template <ElementKind> class S,
-          ElementKind kind, class IntegrationOrderFunctor>
-inline void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::
-    computeNormalsOnIntegrationPointsPoint1(const Array<Real> & /*unused*/,
-                                            Array<Real> & normal,
-                                            GhostType ghost_type) const {
-  AKANTU_DEBUG_IN();
-
-  AKANTU_DEBUG_ASSERT(mesh.getSpatialDimension() == 1,
-                      "Mesh dimension must be 1 to compute normals on points!");
-  const auto type = _point_1;
-  auto spatial_dimension = mesh.getSpatialDimension();
-  // Int nb_nodes_per_element  = Mesh::getNbNodesPerElement(type);
-  auto nb_points = getNbIntegrationPoints(type, ghost_type);
-  const auto & connectivity = mesh.getConnectivity(type, ghost_type);
-  auto nb_element = connectivity.size();
-
-  normal.resize(nb_element * nb_points);
-  auto normals_on_quad =
-      make_view(normal, spatial_dimension, nb_points).begin();
-  const auto & segments = mesh.getElementToSubelement(type, ghost_type);
-  const auto & coords = mesh.getNodes();
-
-  const Mesh * mesh_segment;
-  if (mesh.isMeshFacets()) {
-    mesh_segment = &(mesh.getMeshParent());
-  } else {
-    mesh_segment = &mesh;
-  }
-
-  for (Idx elem = 0; elem < nb_element; ++elem) {
-    auto nb_segment = segments(elem).size();
-    AKANTU_DEBUG_ASSERT(
-        nb_segment > 0,
-        "Impossible to compute a normal on a point connected to 0 segments");
-
-    Real normal_value = 1;
-    if (nb_segment == 1) {
-      auto point = connectivity(elem);
-      const auto segment = segments(elem)[0];
-      const auto & segment_connectivity =
-          mesh_segment->getConnectivity(segment.type, segment.ghost_type);
-      Vector<Idx> segment_points = segment_connectivity.begin(
-          Mesh::getNbNodesPerElement(segment.type))[segment.element];
-      Real difference;
-      if (segment_points(0) == point) {
-        difference = coords(elem) - coords(segment_points(1));
-      } else {
-        difference = coords(elem) - coords(segment_points(0));
-      }
-
-      normal_value = difference / std::abs(difference);
-    }
-
-    for (Idx n(0); n < nb_points; ++n) {
-      (*normals_on_quad)(0, n) = normal_value;
-    }
-    ++normals_on_quad;
-  }
-
-  AKANTU_DEBUG_OUT();
-}
-
 } // namespace akantu
diff --git a/test/test_gtest_utils.hh b/test/test_gtest_utils.hh
index 0f5b6269d..2e016690e 100644
--- a/test/test_gtest_utils.hh
+++ b/test/test_gtest_utils.hh
@@ -1,260 +1,230 @@
 /**
  * @file   test_gtest_utils.hh
  *
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Tue Nov 14 2017
  * @date last modification:  Fri Jan 10 2020
  *
  * @brief  Utils to help write tests
  *
  *
  * @section LICENSE
  *
  * Copyright (©) 2016-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  * Akantu is free software: you can redistribute it and/or modify it under the
  * terms of the GNU Lesser General Public License as published by the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  * 
  * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
  * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
  * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
  * details.
  * 
  * You should have received a copy of the GNU Lesser General Public License
  * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 #include "aka_common.hh"
 #include "aka_iterators.hh"
 /* -------------------------------------------------------------------------- */
 #include <boost/preprocessor.hpp>
 #include <gtest/gtest.h>
 /* -------------------------------------------------------------------------- */
 
 #ifndef AKANTU_TEST_GTEST_UTILS_HH_
 #define AKANTU_TEST_GTEST_UTILS_HH_
 
 #if !defined(TYPED_TEST_SUITE)
 #define TYPED_TEST_SUITE(...) TYPED_TEST_CASE(__VA_ARGS__)
 #endif
 
 #if !defined(TYPED_TEST_SUITE_P)
 #define TYPED_TEST_SUITE_P(...) TYPED_TEST_CASE_P(__VA_ARGS__)
 #endif
 
 #if !defined(REGISTER_TYPED_TEST_SUITE_P)
 #define REGISTER_TYPED_TEST_SUITE_P(...) REGISTER_TYPED_TEST_CASE_P(__VA_ARGS__)
 #endif
 
 #if !defined(INSTANTIATE_TYPED_TEST_SUITE_P)
 #define INSTANTIATE_TYPED_TEST_SUITE_P(...)                                    \
   INSTANTIATE_TYPED_TEST_CASE_P(__VA_ARGS__)
 #endif
 
 namespace {
 
-/* -------------------------------------------------------------------------- */
-template <::akantu::ElementType t>
-using element_type_t = std::integral_constant<::akantu::ElementType, t>;
-
 /* -------------------------------------------------------------------------- */
 template <typename... T> struct gtest_list {};
 
 template <typename... Ts> struct gtest_list<std::tuple<Ts...>> {
   using type = ::testing::Types<Ts...>;
 };
 
 template <typename... T> using gtest_list_t = typename gtest_list<T...>::type;
 
 /* -------------------------------------------------------------------------- */
 //template <typename... T> struct tuple_concat {};
 
-template <typename... Ts>
-struct tuple_concat {
-  using type = decltype(std::tuple_cat(std::declval<Ts>()...));
-};
-
-template <typename... T>
-using tuple_concat_t = typename tuple_concat<T...>::type;
-
-/* -------------------------------------------------------------------------- */
-template <template <typename> class Pred, typename... Ts>
-struct tuple_filter {};
-
-template <template <typename> class Pred, typename T>
-struct tuple_filter<Pred, std::tuple<T>> {
-  using type = std::conditional_t<Pred<T>::value, std::tuple<T>, std::tuple<>>;
-};
-
-template <template <typename> class Pred, typename T, typename... Ts>
-struct tuple_filter<Pred, std::tuple<T, Ts...>> {
-  using type =
-      tuple_concat_t<typename tuple_filter<Pred, std::tuple<T>>::type,
-                     typename tuple_filter<Pred, std::tuple<Ts...>>::type>;
-};
-
-template <template <typename> class Pred, typename... Ts>
-using tuple_filter_t = typename tuple_filter<Pred, Ts...>::type;
 
 /* -------------------------------------------------------------------------- */
 template <size_t N, typename... Ts> struct tuple_split {};
 
 template <size_t N, typename T, typename... Ts>
 struct tuple_split<N, std::tuple<T, Ts...>> {
 protected:
   using split = tuple_split<N - 1, std::tuple<Ts...>>;
 
 public:
-  using type = tuple_concat_t<std::tuple<T>, typename split::type>;
+  using type = tuple::cat_t<std::tuple<T>, typename split::type>;
   using type_tail = typename split::type_tail;
 };
 
 template <typename T, typename... Ts>
 struct tuple_split<1, std::tuple<T, Ts...>> {
   using type = std::tuple<T>;
   using type_tail = std::tuple<Ts...>;
 };
 
 template <size_t N, typename... T>
 using tuple_split_t = typename tuple_split<N, T...>::type;
 
 template <size_t N, typename... T>
 using tuple_split_tail_t = typename tuple_split<N, T...>::type_tail;
 
 /* -------------------------------------------------------------------------- */
 template <typename... T> struct cross_product {};
 
 template <typename... T2s>
 struct cross_product<std::tuple<>, std::tuple<T2s...>> {
   using type = std::tuple<>;
 };
 
 template <typename T1, typename... T1s, typename... T2s>
 struct cross_product<std::tuple<T1, T1s...>, std::tuple<T2s...>> {
-  using type = tuple_concat_t<
+  using type = tuple::cat_t<
       std::tuple<std::tuple<T1, T2s>...>,
       typename cross_product<std::tuple<T1s...>, std::tuple<T2s...>>::type>;
 };
 
 template <typename... T>
 using cross_product_t = typename cross_product<T...>::type;
 /* -------------------------------------------------------------------------- */
 
 } // namespace
 
 #define OP_CAT(s, data, elem) BOOST_PP_CAT(_element_type, elem)
 
 // creating a type instead of a using helps to debug
 #define AKANTU_DECLARE_ELEMENT_TYPE_STRUCT(r, data, elem)                      \
   struct BOOST_PP_CAT(_element_type, elem)                                     \
       : public element_type_t<::akantu::elem> {};
 
 BOOST_PP_SEQ_FOR_EACH(AKANTU_DECLARE_ELEMENT_TYPE_STRUCT, _,
                       AKANTU_ALL_ELEMENT_TYPE)
 
 #undef AKANTU_DECLARE_ELEMENT_TYPE_STRUCT
 
 using TestElementTypesAll = std::tuple<BOOST_PP_SEQ_ENUM(
     BOOST_PP_SEQ_TRANSFORM(OP_CAT, _, AKANTU_ek_regular_ELEMENT_TYPE))>;
 
 #if defined(AKANTU_COHESIVE_ELEMENT)
 using TestCohesiveElementTypes = std::tuple<BOOST_PP_SEQ_ENUM(
     BOOST_PP_SEQ_TRANSFORM(OP_CAT, _, AKANTU_ek_cohesive_ELEMENT_TYPE))>;
 #endif
 
 #if defined(AKANTU_STRUCTURAL_MECHANICS)
 using TestElementTypesStructural = std::tuple<BOOST_PP_SEQ_ENUM(
     BOOST_PP_SEQ_TRANSFORM(OP_CAT, _, AKANTU_ek_structural_ELEMENT_TYPE))>;
 #endif
 
 using TestAllDimensions = std::tuple<std::integral_constant<unsigned int, 1>,
                                      std::integral_constant<unsigned int, 2>,
                                      std::integral_constant<unsigned int, 3>>;
 
 template <typename T, ::akantu::ElementType type>
 using is_element = aka::bool_constant<T::value == type>;
 
 template <typename T>
 using not_is_point_1 = aka::negation<is_element<T, ::akantu::_point_1>>;
 
-using TestElementTypes = tuple_filter_t<not_is_point_1, TestElementTypesAll>;
+using TestElementTypes = tuple::filter_t<not_is_point_1, TestElementTypesAll>;
 
 #if defined(AKANTU_STRUCTURAL_MECHANICS)
 using StructuralTestElementTypes =
-    tuple_filter_t<not_is_point_1, TestElementTypesStructural>;
+    tuple::filter_t<not_is_point_1, TestElementTypesStructural>;
 #endif
 
 /* -------------------------------------------------------------------------- */
 /* -------------------------------------------------------------------------- */
 template <size_t degree> class Polynomial {
 public:
   Polynomial() = default;
 
   Polynomial(std::initializer_list<double> && init) {
     for (auto && pair : akantu::zip(init, constants))
       std::get<1>(pair) = std::get<0>(pair);
   }
 
   double operator()(double x) {
     double res = 0.;
     for (auto && vals : akantu::enumerate(constants)) {
       double a;
       int k;
       std::tie(k, a) = vals;
       res += a * std::pow(x, k);
     }
     return res;
   }
 
   Polynomial extract(size_t pdegree) {
     Polynomial<degree> extract(*this);
     for (size_t d = pdegree + 1; d < degree + 1; ++d)
       extract.constants[d] = 0;
     return extract;
   }
 
   auto integral() {
     Polynomial<degree + 1> integral_;
     integral_.set(0, 0.);
     ;
     for (size_t d = 0; d < degree + 1; ++d) {
       integral_.set(1 + d, get(d) / double(d + 1));
     }
     return integral_;
   }
 
   auto integrate(double a, double b) {
     auto primitive = integral();
     return (primitive(b) - primitive(a));
   }
 
   double get(int i) const { return constants[i]; }
 
   void set(int i, double a) { constants[i] = a; }
 
 protected:
   std::array<double, degree + 1> constants;
 };
 
 template <size_t degree>
 std::ostream & operator<<(std::ostream & stream, const Polynomial<degree> & p) {
   for (size_t d = 0; d < degree + 1; ++d) {
     if (d != 0)
       stream << " + ";
 
     stream << p.get(degree - d);
     if (d != degree)
       stream << "x ^ " << degree - d;
   }
   return stream;
 }
 
 /* -------------------------------------------------------------------------- */
 
 #endif /* AKANTU_TEST_GTEST_UTILS_HH_ */
diff --git a/third-party/akantu_iterators/include/aka_tuple_tools.hh b/third-party/akantu_iterators/include/aka_tuple_tools.hh
index 38eac8e3c..b2932f71a 100644
--- a/third-party/akantu_iterators/include/aka_tuple_tools.hh
+++ b/third-party/akantu_iterators/include/aka_tuple_tools.hh
@@ -1,407 +1,431 @@
 /**
  * @file   aka_tuple_tools.hh
  *
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Fri Aug 11 2017
  * @date last modification: Mon Jan 29 2018
  *
  * @brief  iterator interfaces
  *
  *
  * Copyright 2019 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_compatibilty_with_cpp_standard.hh"
 #include "aka_str_hash.hh"
 /* -------------------------------------------------------------------------- */
 #include <tuple>
 /* -------------------------------------------------------------------------- */
 
 #ifndef AKANTU_AKA_TUPLE_TOOLS_HH
 #define AKANTU_AKA_TUPLE_TOOLS_HH
 
 #ifndef AKANTU_ITERATORS_NAMESPACE
 #define AKANTU_ITERATORS_NAMESPACE akantu
 #endif
 
 namespace AKANTU_ITERATORS_NAMESPACE {
 
 namespace tuple {
   /* ---------------------------------------------------------------------- */
   template <typename tag, typename type> struct named_tag {
     using _tag = tag;   ///< key
     using _type = type; ///< value type
 
     template <
         typename T,
         std::enable_if_t<not std::is_same<named_tag, T>::value> * = nullptr>
     explicit named_tag(T && value) // NOLINT
         : _value(std::forward<T>(value)) {}
 
     type _value;
   };
 
   namespace details {
 /* ---------------------------------------------------------------------- */
 #if (defined(__GNUC__) || defined(__GNUG__))
 #pragma GCC diagnostic push
 #pragma GCC diagnostic ignored "-Weffc++"
 #endif
     template <typename tag> struct named_tag_proxy {
       using _tag = tag;
 
       template <typename T> decltype(auto) operator=(T && value) {
         return named_tag<_tag, T>{std::forward<T>(value)};
       }
     };
 #if (defined(__GNUC__) || defined(__GNUG__))
 #pragma GCC diagnostic pop
 #endif
   } // namespace details
 
   /* ---------------------------------------------------------------------- */
   template <typename T> struct is_named_tag : public std::false_type {};
   template <typename tag>
   struct is_named_tag<details::named_tag_proxy<tag>> : public std::true_type {};
   template <typename tag, typename type>
   struct is_named_tag<named_tag<tag, type>> : public std::true_type {};
   /* ---------------------------------------------------------------------- */
 
   template <class... Params>
   struct named_tuple : public std::tuple<typename Params::_type...> {
     using Names_t = std::tuple<typename Params::_tag...>;
     using parent = std::tuple<typename Params::_type...>;
 
     named_tuple(Params &&... params)
         : parent(std::forward<typename Params::_type>(params._value)...) {}
 
     named_tuple(typename Params::_type &&... args)
         : parent(std::forward<typename Params::_type>(args)...) {}
 
   private:
     template <typename tag, std::size_t Idx,
               std::enable_if_t<Idx == sizeof...(Params)> * = nullptr>
     static constexpr std::size_t get_element_index() noexcept {
       return -1;
     }
 
     template <typename tag, std::size_t Idx,
               std::enable_if_t<(Idx < sizeof...(Params))> * = nullptr>
     static constexpr std::size_t get_element_index() noexcept {
       using _tag = std::tuple_element_t<Idx, Names_t>;
       return (std::is_same<_tag, tag>::value)
                  ? Idx
                  : get_element_index<tag, Idx + 1>();
     }
 
   public:
     template <typename NT,
               std::enable_if_t<is_named_tag<NT>::value> * = nullptr>
     constexpr decltype(auto) get(NT && /*unused*/) noexcept {
       const auto index = get_element_index<typename NT::_tag, 0>();
       static_assert((index != -1), "wrong named_tag");
       return (std::get<index>(*this));
     }
 
     template <typename NT,
               std::enable_if_t<is_named_tag<NT>::value> * = nullptr>
     constexpr decltype(auto) get(NT && /*unused*/) const noexcept {
       const auto index = get_element_index<typename NT::_tag, 0>();
       static_assert((index != -1), "wrong named_tag");
       return std::get<index>(*this);
     }
   };
 
   /* ---------------------------------------------------------------------- */
   template <typename T> struct is_named_tuple : public std::false_type {};
   template <typename... Params>
   struct is_named_tuple<named_tuple<Params...>> : public std::true_type {};
   /* ---------------------------------------------------------------------- */
 
   template <typename... Params>
   constexpr decltype(auto) make_named_tuple(Params &&... params) noexcept {
     return named_tuple<Params...>(std::forward<Params>(params)...);
   }
 
   template <typename tag> constexpr decltype(auto) make_named_tag() noexcept {
     return details::named_tag_proxy<tag>{};
   }
 
   template <size_t HashCode> constexpr decltype(auto) get() {
     return make_named_tag<std::integral_constant<size_t, HashCode>>();
   }
 
   template <size_t HashCode, class Tuple>
   constexpr decltype(auto) get(Tuple && tuple) {
     return tuple.get(get<HashCode>());
   }
 
   template <typename Param, typename Tuple,
             std::enable_if_t<is_named_tag<Param>::value> * = nullptr>
   constexpr decltype(auto) get(Tuple && tuple) noexcept {
     return tuple.template get<typename Param::hash>();
   }
 
 #if !defined(__INTEL_COMPILER)
-#  if defined(__clang__)
+#if defined(__clang__)
 // clang warnings here
-#    pragma clang diagnostic push
-#    pragma clang diagnostic ignored "-Wgnu-string-literal-operator-template"
-#  elif (defined(__GNUC__) || defined(__GNUG__))
+#pragma clang diagnostic push
+#pragma clang diagnostic ignored "-Wgnu-string-literal-operator-template"
+#elif (defined(__GNUC__) || defined(__GNUG__))
 // gcc warnings here
-#    pragma GCC diagnostic push
-#    pragma GCC diagnostic ignored "-Wpedantic"
-#  endif
+#pragma GCC diagnostic push
+#pragma GCC diagnostic ignored "-Wpedantic"
+#endif
   /// this is a GNU exstension
   /// http://www.open-std.org/jtc1/sc22/wg21/docs/papers/2013/n3599.html
   template <class CharT, CharT... chars>
   constexpr decltype(auto) operator"" _n() {
     return make_named_tag<std::integral_constant<
         std::size_t, string_literal<CharT, chars...>::hash>>();
   }
-#  if defined(__clang__)
-#    pragma clang diagnostic pop
-#  elif (defined(__GNUC__) || defined(__GNUG__))
-#    pragma GCC diagnostic pop
-#  endif
+#if defined(__clang__)
+#pragma clang diagnostic pop
+#elif (defined(__GNUC__) || defined(__GNUG__))
+#pragma GCC diagnostic pop
+#endif
 #endif
   /* ------------------------------------------------------------------------ */
   namespace details {
     template <std::size_t N> struct Foreach {
       template <class Tuple>
 
       static constexpr inline auto not_equal(Tuple && a, Tuple && b) -> bool {
         if (std::get<N - 1>(std::forward<Tuple>(a)) ==
             std::get<N - 1>(std::forward<Tuple>(b))) {
           return false;
         }
         return Foreach<N - 1>::not_equal(std::forward<Tuple>(a),
                                          std::forward<Tuple>(b));
       }
 
       template <class Tuple, class V>
       static constexpr inline auto find(Tuple && tuple, V && value)
           -> std::size_t {
         constexpr auto size = std::tuple_size<std::decay_t<Tuple>>::value;
         if (std::get<size - N>(std::forward<Tuple>(tuple)) == value) {
           return size - N;
         }
 
         return Foreach<N - 1>::find(std::forward<Tuple>(tuple),
                                     std::forward<V>(value));
       }
     };
 
     /* ---------------------------------------------------------------------- */
     template <> struct Foreach<0> {
       template <class Tuple>
       static constexpr inline auto not_equal(Tuple && a, Tuple && b) -> bool {
         return std::get<0>(std::forward<Tuple>(a)) !=
                std::get<0>(std::forward<Tuple>(b));
       }
 
       template <class Tuple, class V>
       static constexpr inline auto find(Tuple && tuple, V && value)
           -> std::size_t {
         constexpr auto size = std::tuple_size<std::decay_t<Tuple>>::value;
         if (std::get<size - 1>(std::forward<Tuple>(tuple)) == value) {
           return size - 1;
         }
 
         return size;
       }
     };
 
     template <typename... Ts>
     decltype(auto) make_tuple_no_decay(Ts &&... args) {
       return std::tuple<Ts...>(std::forward<Ts>(args)...);
     }
 
     template <typename... Names, typename... Ts>
-    constexpr decltype(auto)
-    make_named_tuple_no_decay(std::tuple<Names...> /*unused*/, Ts &&... args) {
+    decltype(auto) make_named_tuple_no_decay(std::tuple<Names...> /*unused*/,
+                                             Ts &&... args) {
       return named_tuple<named_tag<Names, Ts>...>(std::forward<Ts>(args)...);
     }
 
     template <class F, class Tuple, std::size_t... Is>
     constexpr void foreach_impl(F && func, Tuple && tuple,
                                 std::index_sequence<Is...> && /*unused*/) {
       (void)std::initializer_list<int>{
           (std::forward<F>(func)(std::get<Is>(std::forward<Tuple>(tuple))),
            0)...};
     }
 
     template <class F, class Tuple, std::size_t... Is>
     constexpr decltype(auto)
     transform_impl(F && func, Tuple && tuple,
                    std::index_sequence<Is...> && /*unused*/) {
       return make_tuple_no_decay(
           std::forward<F>(func)(std::get<Is>(std::forward<Tuple>(tuple)))...);
     }
 
     template <class F, class Tuple, std::size_t... Is>
     constexpr decltype(auto)
     transform_impl(F && func, Tuple && tuple1, Tuple && tuple2,
                    std::index_sequence<Is...> && /*unused*/) {
       return make_tuple_no_decay(
           std::forward<F>(func)(std::get<Is>(std::forward<Tuple>(tuple1)),
                                 std::get<Is>(std::forward<Tuple>(tuple2)))...);
     }
 
     template <class F, class Tuple, std::size_t... Is>
     constexpr decltype(auto)
     transform_named_impl(F && func, Tuple && tuple,
                          std::index_sequence<Is...> && /*unused*/) {
       return make_named_tuple_no_decay(
           typename std::decay_t<Tuple>::Names_t{},
           std::forward<F>(func)(std::get<Is>(std::forward<Tuple>(tuple)))...);
     }
 
   } // namespace details
 
   /* ------------------------------------------------------------------------
    */
   template <class Tuple,
             std::enable_if_t<not is_named_tuple<std::decay_t<Tuple>>::value> * =
                 nullptr>
   constexpr auto are_not_equal(Tuple && a, Tuple && b) -> bool {
     return details::Foreach<std::tuple_size<std::decay_t<Tuple>>::value>::
         not_equal(std::forward<Tuple>(a), std::forward<Tuple>(b));
   }
 
   template <
       class Tuple,
       std::enable_if_t<is_named_tuple<std::decay_t<Tuple>>::value> * = nullptr>
   constexpr auto are_not_equal(Tuple && a, Tuple && b) -> bool {
     return details::Foreach<
         std::tuple_size<typename std::decay_t<Tuple>::parent>::value>::
         not_equal(std::forward<Tuple>(a), std::forward<Tuple>(b));
   }
 
   template <class Tuple, class V>
   constexpr decltype(auto) find(Tuple && tuple, V && value) {
     return details::Foreach<std::tuple_size<std::decay_t<Tuple>>::value>::find(
         std::forward<Tuple>(tuple), std::forward<V>(value));
   }
 
   template <class F, class Tuple,
             std::enable_if_t<not is_named_tuple<std::decay_t<Tuple>>::value> * =
                 nullptr>
   constexpr void foreach (F && func, Tuple && tuple) {
     return details::foreach_impl(
         std::forward<F>(func), std::forward<Tuple>(tuple),
         std::make_index_sequence<
             std::tuple_size<std::decay_t<Tuple>>::value>{});
   }
 
   template <
       class F, class Tuple,
       std::enable_if_t<is_named_tuple<std::decay_t<Tuple>>::value> * = nullptr>
   constexpr void foreach (F && func, Tuple && tuple) {
     return details::foreach_impl(
         std::forward<F>(func), std::forward<Tuple>(tuple),
         std::make_index_sequence<
             std::tuple_size<typename std::decay_t<Tuple>::parent>::value>{});
   }
 
   template <class F, class Tuple,
             std::enable_if_t<not is_named_tuple<std::decay_t<Tuple>>::value> * =
                 nullptr>
   constexpr decltype(auto) transform(F && func, Tuple && tuple) {
     return details::transform_impl(
         std::forward<F>(func), std::forward<Tuple>(tuple),
         std::make_index_sequence<
             std::tuple_size<std::decay_t<Tuple>>::value>{});
   }
 
   template <class F, class Tuple,
             std::enable_if_t<not is_named_tuple<std::decay_t<Tuple>>::value> * =
                 nullptr>
   constexpr decltype(auto) transform(F && func, Tuple && tuple1,
                                      Tuple && tuple2) {
     return details::transform_impl(
         std::forward<F>(func), std::forward<Tuple>(tuple1),
         std::forward<Tuple>(tuple2),
         std::make_index_sequence<
             std::tuple_size<std::decay_t<Tuple>>::value>{});
   }
 
   template <
       class F, class Tuple,
       std::enable_if_t<is_named_tuple<std::decay_t<Tuple>>::value> * = nullptr>
   constexpr decltype(auto) transform(F && func, Tuple && tuple) {
     return details::transform_named_impl(
         std::forward<F>(func), std::forward<Tuple>(tuple),
         std::make_index_sequence<
             std::tuple_size<typename std::decay_t<Tuple>::parent>::value>{});
   }
 
   namespace details {
     template <class Tuple, std::size_t... Is>
     constexpr decltype(auto) flatten(Tuple && tuples,
                                      std::index_sequence<Is...> /*unused*/) {
       return std::tuple_cat(std::get<Is>(tuples)...);
     }
   } // namespace details
 
   template <class Tuple> constexpr decltype(auto) flatten(Tuple && tuples) {
     return details::flatten(std::forward<Tuple>(tuples),
                             std::make_index_sequence<
                                 std::tuple_size<std::decay_t<Tuple>>::value>());
   }
 
   namespace details {
     template <size_t n, class Tuple>
     decltype(auto) dynamic_get_impl(size_t i, Tuple && tuple) {
       constexpr auto size = std::tuple_size<std::decay_t<Tuple>>::value;
       if (i == n) {
         return std::get<n>(tuple);
       } else if (n == size - 1) {
         throw std::out_of_range("Tuple element out of range.");
       } else {
         return dynamic_get_impl<(n < size - 1 ? n + 1 : 0)>(i, tuple);
       }
     }
   } // namespace details
 
+  template <template <typename> class Pred, typename... Ts>
+  using filter_t = typename filter<Pred, Ts...>::type;
   template <class Tuple>
   constexpr decltype(auto) dynamic_get(std::size_t i, Tuple && tuple) {
     return details::dynamic_get_impl<0>(i, std::forward<Tuple>(tuple));
   }
+
+  template <typename... Ts> struct cat {
+    using type = decltype(std::tuple_cat(std::declval<Ts>()...));
+  };
+
+  template <typename... T> using cat_t = typename cat<T...>::type;
+
+  /* --------------------------------------------------------------------------
+   */
+  template <template <typename> class Pred, typename... Ts> struct filter {};
+
+  template <template <typename> class Pred, typename T>
+  struct filter<Pred, std::tuple<T>> {
+    using type =
+        std::conditional_t<Pred<T>::value, std::tuple<T>, std::tuple<>>;
+  };
+
+  template <template <typename> class Pred, typename T, typename... Ts>
+  struct filter<Pred, std::tuple<T, Ts...>> {
+    using type = cat_t<typename filter<Pred, std::tuple<T>>::type,
+                       typename filter<Pred, std::tuple<Ts...>>::type>;
+  };
 } // namespace tuple
 
 } // namespace AKANTU_ITERATORS_NAMESPACE
 
 /* -------------------------------------------------------------------------- */
 #include <iterator>
 /* -------------------------------------------------------------------------- */
 
 namespace std {
 template <typename tag, typename type>
 struct iterator_traits<
     ::AKANTU_ITERATORS_NAMESPACE::tuple::named_tag<tag, type>> {
   using iterator_category = typename type::iterator_category;
   using value_type = typename type::value_type;
   using difference_type = typename type::difference_type;
   using pointer = typename type::pointer;
   using reference = typename type::reference;
 };
 } // namespace std
 #endif /* AKANTU_AKA_TUPLE_TOOLS_HH */