diff --git a/include/TFEL/Material/FiniteStrainBehaviourTangentOperator.hxx b/include/TFEL/Material/FiniteStrainBehaviourTangentOperator.hxx
index fc7fafe65..017170025 100644
--- a/include/TFEL/Material/FiniteStrainBehaviourTangentOperator.hxx
+++ b/include/TFEL/Material/FiniteStrainBehaviourTangentOperator.hxx
@@ -17,9 +17,7 @@
 #include <vector>
 #include <string>
 #include <utility>
-
 #include "TFEL/Config/TFELConfig.hxx"
-
 #include "TFEL/Metaprogramming/Implements.hxx"
 #include "TFEL/TypeTraits/BaseType.hxx"
 #include "TFEL/Math/Forward/st2tost2.hxx"
@@ -28,7 +26,6 @@
 #include "TFEL/Math/T2toST2/T2toST2Concept.hxx"
 #include "TFEL/Math/Forward/t2tot2.hxx"
 #include "TFEL/Math/T2toT2/T2toT2Concept.hxx"
-
 #include "TFEL/Utilities/GenTypeBase.hxx"
 #include "TFEL/Material/FiniteStrainBehaviourTangentOperatorBase.hxx"
 
@@ -115,28 +112,22 @@ namespace tfel::material {
       GenType::operator=(src);
       return *this;
     }  // end of FiniteStrainBehaviourTangentOperator& operator=
-       /*!
-        * \param const T1&, the value affected to the GenType.
-        * \pre   T1 must be a type that the GenType can hold.
-        */
+    /*!
+     * \param const T1&, the value affected to the GenType.
+     * \pre   T1 must be a type that the GenType can hold.
+     */
     template <typename T1>
-    TFEL_INLINE
-        typename std::enable_if<tfel::meta::TLCountNbrOfT<T1, TOTypes>::value ==
-                                    1,
-                                void>::type
-        set_uninitialised() {
+    TFEL_INLINE void set_uninitialised() requires(
+        tfel::meta::TLCountNbrOfT<T1, TOTypes>::value == 1) {
       GenType::template set_uninitialised<T1>();
     }
     /*!
-     * assignement operator
+     * \brief assignement operator
      */
-    template <typename T>
-    std::enable_if_t<
-        tfel::math::implementsT2toST2Concept<T>() &&
-            tfel::math::getSpaceDimension<T>() == N &&
-            std::is_same_v<tfel::math::numeric_type<T>, StressType>,
-        FiniteStrainBehaviourTangentOperator&>
-    operator=(const T& e) {
+    template <tfel::math::T2toST2Concept T>
+    FiniteStrainBehaviourTangentOperator& operator=(const T& e) requires(
+        tfel::math::getSpaceDimension<T>() == N &&
+        std::is_same_v<tfel::math::numeric_type<T>, StressType>) {
       using namespace tfel::math;
       if (this->template is<t2tost2<N, StressType>*>()) {
         *(this->template get<t2tost2<N, StressType>*>()) = e;
@@ -152,11 +143,9 @@ namespace tfel::material {
      * \brief assignement operator
      */
     template <tfel::math::ST2toST2Concept T>
-    std::enable_if_t<
+    FiniteStrainBehaviourTangentOperator& operator=(const T& e) requires(
         tfel::math::getSpaceDimension<T>() == N &&
-            std::is_same_v<tfel::math::numeric_type<T>, StressType>,
-        FiniteStrainBehaviourTangentOperator&>
-    operator=(const T& e) {
+        std::is_same_v<tfel::math::numeric_type<T>, StressType>) {
       using namespace tfel::math;
       if (this->template is<st2tost2<N, StressType>*>()) {
         *(this->template get<st2tost2<N, StressType>*>()) = e;
diff --git a/include/TFEL/Math/ST2toST2/ConvertT2toST2ToST2toST2Expr.hxx b/include/TFEL/Math/ST2toST2/ConvertT2toST2ToST2toST2Expr.hxx
index 0a4110636..5399e44ae 100644
--- a/include/TFEL/Math/ST2toST2/ConvertT2toST2ToST2toST2Expr.hxx
+++ b/include/TFEL/Math/ST2toST2/ConvertT2toST2ToST2toST2Expr.hxx
@@ -26,10 +26,8 @@ namespace tfel::math {
   template <unsigned short N>
   struct ConvertT2toST2ToST2toST2Expr;
 
-  /*!
-   * Partial specialisation for 1D stensor
-   */
-  template <typename ST2toST2ResultType>
+  //! \brief partial specialisation in 1D
+  template <ST2toST2Concept ST2toST2ResultType>
   struct Expr<ST2toST2ResultType, ConvertT2toST2ToST2toST2Expr<1u>>
       : public ST2toST2ConceptBase<
             Expr<ST2toST2ResultType, ConvertT2toST2ToST2toST2Expr<1u>>>,
@@ -42,23 +40,22 @@ namespace tfel::math {
     /*!
      * \param[in] A : t2tost2 to be converted
      */
-    template <typename T2toST2Type>
-    Expr(const T2toST2Type& A) {
+    template <T2toST2Concept T2toST2Type>
+    TFEL_HOST_DEVICE constexpr Expr(const T2toST2Type& A) noexcept {
       using tfel::fsalgo::copy;
-      static_assert(implementsT2toST2Concept<T2toST2Type>());
       static_assert(getSpaceDimension<ST2toST2ResultType>() ==
                     getSpaceDimension<T2toST2Type>());
       static_assert(isAssignableTo<numeric_type<T2toST2Type>,
                                    numeric_type<ST2toST2ResultType>>());
       copy<9u>::exe(&A(0, 0), this->v);
     }  // end of Expr
-    /*!
-     * \brief access operator
-     * \param[in] i : line   index
-     * \param[in] j : column index
-     */
-    const value_type& operator()(const unsigned short i,
-                                 const unsigned short j) const {
+       /*!
+        * \brief access operator
+        * \param[in] i : line   index
+        * \param[in] j : column index
+        */
+    TFEL_HOST_DEVICE constexpr const value_type& operator()(
+        const unsigned short i, const unsigned short j) const noexcept {
       return this->v[i * 3 + j];
     }  // end of operator()
     /*!
@@ -66,16 +63,13 @@ namespace tfel::math {
      * In this case, the number of lines and columns
      * are deduced from the template parameter
      */
-    RunTimeProperties getRunTimeProperties() const {
+    TFEL_HOST_DEVICE constexpr auto getRunTimeProperties() const noexcept {
       return RunTimeProperties();
     }
-  };  // end of struct Expr<ST2toST2ResultType,ConvertT2toST2ToST2toST2Expr<1u>
-      // >
+  };  // end of Expr<ST2toST2ResultType,ConvertT2toST2ToST2toST2Expr<1u>>
 
-  /*!
-   * Partial specialisation for 2D stensor
-   */
-  template <typename ST2toST2ResultType>
+  //! \brief partial specialisation in 2D
+  template <ST2toST2Concept ST2toST2ResultType>
   struct Expr<ST2toST2ResultType, ConvertT2toST2ToST2toST2Expr<2u>>
       : public ST2toST2ConceptBase<
             Expr<ST2toST2ResultType, ConvertT2toST2ToST2toST2Expr<2u>>>,
@@ -89,9 +83,8 @@ namespace tfel::math {
     /*!
      * \param[in] A : t2tost2 to be converted
      */
-    template <typename T2toST2Type>
-    Expr(const T2toST2Type& A) {
-      static_assert(implementsT2toST2Concept<T2toST2Type>());
+    template <T2toST2Concept T2toST2Type>
+    TFEL_HOST_DEVICE constexpr Expr(const T2toST2Type& A) noexcept {
       static_assert(getSpaceDimension<ST2toST2ResultType>() ==
                     getSpaceDimension<T2toST2Type>());
       static_assert(isAssignableTo<numeric_type<T2toST2Type>,
@@ -120,8 +113,8 @@ namespace tfel::math {
      * \param[in] i : line   index
      * \param[in] j : column index
      */
-    const value_type& operator()(const unsigned short i,
-                                 const unsigned short j) const {
+    TFEL_HOST_DEVICE constexpr const value_type& operator()(
+        const unsigned short i, const unsigned short j) const noexcept {
       return this->v[i * 4 + j];
     }  // end of operator()
     /*!
@@ -129,16 +122,13 @@ namespace tfel::math {
      * In this case, the number of lines and columns
      * are deduced from the template parameter
      */
-    RunTimeProperties getRunTimeProperties() const {
+    TFEL_HOST_DEVICE constexpr auto getRunTimeProperties() const noexcept {
       return RunTimeProperties();
     }
-  };  // end of struct Expr<ST2toST2ResultType,ConvertT2toST2ToST2toST2Expr<2u>
-      // >
+  };  // end of Expr<ST2toST2ResultType,ConvertT2toST2ToST2toST2Expr<2u>>
 
-  /*!
-   * Partial specialisation for 3D stensor
-   */
-  template <typename ST2toST2ResultType>
+  //! \brief partial specialisation in 3D
+  template <ST2toST2Concept ST2toST2ResultType>
   struct Expr<ST2toST2ResultType, ConvertT2toST2ToST2toST2Expr<3u>>
       : public ST2toST2ConceptBase<
             Expr<ST2toST2ResultType, ConvertT2toST2ToST2toST2Expr<3u>>>,
@@ -151,9 +141,8 @@ namespace tfel::math {
     /*!
      * \param[in] A : t2tost2 to be converted
      */
-    template <typename T2toST2Type>
-    Expr(const T2toST2Type& A) {
-      static_assert(implementsT2toST2Concept<T2toST2Type>());
+    template <T2toST2Concept T2toST2Type>
+    TFEL_HOST_DEVICE constexpr Expr(const T2toST2Type& A) noexcept {
       static_assert(getSpaceDimension<ST2toST2ResultType>() ==
                     getSpaceDimension<T2toST2Type>());
       static_assert(isAssignableTo<numeric_type<T2toST2Type>,
@@ -202,8 +191,8 @@ namespace tfel::math {
      * \param[in] i : line   index
      * \param[in] j : column index
      */
-    const value_type& operator()(const unsigned short i,
-                                 const unsigned short j) const {
+    TFEL_HOST_DEVICE constexpr const value_type& operator()(
+        const unsigned short i, const unsigned short j) const noexcept {
       return this->v[i * 6 + j];
     }  // end of operator()
     /*!
@@ -211,11 +200,10 @@ namespace tfel::math {
      * In this case, the number of lines and columns
      * are deduced from the template parameter
      */
-    RunTimeProperties getRunTimeProperties() const {
+    TFEL_HOST_DEVICE constexpr auto getRunTimeProperties() const noexcept {
       return RunTimeProperties();
     }
-  };  // end of struct Expr<ST2toST2ResultType,ConvertT2toST2ToST2toST2Expr<3u>
-      // >
+  };  // end of Expr<ST2toST2ResultType,ConvertT2toST2ToST2toST2Expr<3u>>
 
 }  // end of namespace tfel::math
 
diff --git a/include/TFEL/Math/ST2toST2/ConvertToTangentModuli.hxx b/include/TFEL/Math/ST2toST2/ConvertToTangentModuli.hxx
index a9909838f..ca6e1706e 100644
--- a/include/TFEL/Math/ST2toST2/ConvertToTangentModuli.hxx
+++ b/include/TFEL/Math/ST2toST2/ConvertToTangentModuli.hxx
@@ -24,158 +24,129 @@
 
 namespace tfel::math {
 
-  template <typename T2toST2Type, TensorConcept TensorType>
-  typename std::enable_if<((getSpaceDimension<T2toST2Type>() == 1u) &&
-                           (getSpaceDimension<TensorType>() == 1u) &&
-                           implementsT2toST2Concept<T2toST2Type>()),
-                          st2tost2<1u,
-                                   result_type<numeric_type<T2toST2Type>,
-                                               numeric_type<TensorType>,
-                                               OpMult>>>::type
-  convertToTangentModuli(const T2toST2Type& K, const TensorType& F) {
-    st2tost2<1u, result_type<numeric_type<T2toST2Type>,
-                             numeric_type<TensorType>, OpMult>>
-        C;
-    C(0, 0) = K(0, 0) * F(0);
-    C(0, 1) = K(0, 1) * F(1);
-    C(0, 2) = K(0, 2) * F(2);
-    C(1, 0) = K(1, 0) * F(0);
-    C(1, 1) = K(1, 1) * F(1);
-    C(1, 2) = K(1, 2) * F(2);
-    C(2, 0) = K(2, 0) * F(0);
-    C(2, 1) = K(2, 1) * F(1);
-    C(2, 2) = K(2, 2) * F(2);
+  template <T2toST2Concept T2toST2Type, TensorConcept TensorType>
+  TFEL_HOST_DEVICE constexpr auto convertToTangentModuli(
+      const T2toST2Type& K, const TensorType& F) noexcept  //
+      requires(getSpaceDimension<T2toST2Type>() ==
+               getSpaceDimension<TensorType>()) {
+    constexpr auto N = getSpaceDimension<T2toST2Type>();
+    using value_type = result_type<numeric_type<T2toST2Type>,
+                                   numeric_type<TensorType>, OpMult>;
+    auto C = st2tost2<N, value_type>{};
+    if constexpr (N == 1) {
+      C(0, 0) = K(0, 0) * F(0);
+      C(0, 1) = K(0, 1) * F(1);
+      C(0, 2) = K(0, 2) * F(2);
+      C(1, 0) = K(1, 0) * F(0);
+      C(1, 1) = K(1, 1) * F(1);
+      C(1, 2) = K(1, 2) * F(2);
+      C(2, 0) = K(2, 0) * F(0);
+      C(2, 1) = K(2, 1) * F(1);
+      C(2, 2) = K(2, 2) * F(2);
+    } else if constexpr (N == 2) {
+      constexpr const auto icste2 = Cste<value_type>::isqrt2 / 2;
+      C(0, 0) = K(0, 0) * F(0) + K(0, 3) * F(3);
+      C(0, 3) =
+          (K(0, 0) * F(4) + K(0, 3) * F(1) + K(0, 4) * F(0) + K(0, 1) * F(3)) *
+          icste2;
+      C(0, 1) = K(0, 4) * F(4) + K(0, 1) * F(1);
+      C(0, 2) = K(0, 2) * F(2);
+      C(3, 0) = K(3, 0) * F(0) + K(3, 3) * F(3);
+      C(3, 3) =
+          (K(3, 0) * F(4) + K(3, 3) * F(1) + K(3, 4) * F(0) + K(3, 1) * F(3)) *
+          icste2;
+      C(3, 1) = K(3, 4) * F(4) + K(3, 1) * F(1);
+      C(3, 2) = K(3, 2) * F(2);
+      C(1, 0) = K(1, 0) * F(0) + K(1, 3) * F(3);
+      C(1, 3) =
+          (K(1, 0) * F(4) + K(1, 3) * F(1) + K(1, 4) * F(0) + K(1, 1) * F(3)) *
+          icste2;
+      C(1, 1) = K(1, 4) * F(4) + K(1, 1) * F(1);
+      C(1, 2) = K(1, 2) * F(2);
+      C(2, 0) = K(2, 0) * F(0) + K(2, 3) * F(3);
+      C(2, 3) =
+          (K(2, 0) * F(4) + K(2, 3) * F(1) + K(2, 4) * F(0) + K(2, 1) * F(3)) *
+          icste2;
+      C(2, 1) = K(2, 4) * F(4) + K(2, 1) * F(1);
+      C(2, 2) = K(2, 2) * F(2);
+    } else {
+      constexpr const auto icste2 = Cste<value_type>::isqrt2 / 2;
+      C(0, 0) = K(0, 0) * F(0) + K(0, 3) * F(3) + K(0, 5) * F(5);
+      C(0, 3) = (K(0, 0) * F(4) + K(0, 3) * F(1) + K(0, 5) * F(7) +
+                 K(0, 4) * F(0) + K(0, 1) * F(3) + K(0, 7) * F(5)) *
+                icste2;
+      C(0, 4) = (K(0, 0) * F(6) + K(0, 3) * F(8) + K(0, 5) * F(2) +
+                 K(0, 6) * F(0) + K(0, 8) * F(3) + K(0, 2) * F(5)) *
+                icste2;
+      C(0, 1) = K(0, 4) * F(4) + K(0, 1) * F(1) + K(0, 7) * F(7);
+      C(0, 5) = (K(0, 4) * F(6) + K(0, 1) * F(8) + K(0, 7) * F(2) +
+                 K(0, 6) * F(4) + K(0, 8) * F(1) + K(0, 2) * F(7)) *
+                icste2;
+      C(0, 2) = K(0, 6) * F(6) + K(0, 8) * F(8) + K(0, 2) * F(2);
+      C(3, 0) = K(3, 0) * F(0) + K(3, 3) * F(3) + K(3, 5) * F(5);
+      C(3, 3) = (K(3, 0) * F(4) + K(3, 3) * F(1) + K(3, 5) * F(7) +
+                 K(3, 4) * F(0) + K(3, 1) * F(3) + K(3, 7) * F(5)) *
+                icste2;
+      C(3, 4) = (K(3, 0) * F(6) + K(3, 3) * F(8) + K(3, 5) * F(2) +
+                 K(3, 6) * F(0) + K(3, 8) * F(3) + K(3, 2) * F(5)) *
+                icste2;
+      C(3, 1) = K(3, 4) * F(4) + K(3, 1) * F(1) + K(3, 7) * F(7);
+      C(3, 5) = (K(3, 4) * F(6) + K(3, 1) * F(8) + K(3, 7) * F(2) +
+                 K(3, 6) * F(4) + K(3, 8) * F(1) + K(3, 2) * F(7)) *
+                icste2;
+      C(3, 2) = K(3, 6) * F(6) + K(3, 8) * F(8) + K(3, 2) * F(2);
+      C(4, 0) = K(4, 0) * F(0) + K(4, 3) * F(3) + K(4, 5) * F(5);
+      C(4, 3) = (K(4, 0) * F(4) + K(4, 3) * F(1) + K(4, 5) * F(7) +
+                 K(4, 4) * F(0) + K(4, 1) * F(3) + K(4, 7) * F(5)) *
+                icste2;
+      C(4, 4) = (K(4, 0) * F(6) + K(4, 3) * F(8) + K(4, 5) * F(2) +
+                 K(4, 6) * F(0) + K(4, 8) * F(3) + K(4, 2) * F(5)) *
+                icste2;
+      C(4, 1) = K(4, 4) * F(4) + K(4, 1) * F(1) + K(4, 7) * F(7);
+      C(4, 5) = (K(4, 4) * F(6) + K(4, 1) * F(8) + K(4, 7) * F(2) +
+                 K(4, 6) * F(4) + K(4, 8) * F(1) + K(4, 2) * F(7)) *
+                icste2;
+      C(4, 2) = K(4, 6) * F(6) + K(4, 8) * F(8) + K(4, 2) * F(2);
+      C(1, 0) = K(1, 0) * F(0) + K(1, 3) * F(3) + K(1, 5) * F(5);
+      C(1, 3) = (K(1, 0) * F(4) + K(1, 3) * F(1) + K(1, 5) * F(7) +
+                 K(1, 4) * F(0) + K(1, 1) * F(3) + K(1, 7) * F(5)) *
+                icste2;
+      C(1, 4) = (K(1, 0) * F(6) + K(1, 3) * F(8) + K(1, 5) * F(2) +
+                 K(1, 6) * F(0) + K(1, 8) * F(3) + K(1, 2) * F(5)) *
+                icste2;
+      C(1, 1) = K(1, 4) * F(4) + K(1, 1) * F(1) + K(1, 7) * F(7);
+      C(1, 5) = (K(1, 4) * F(6) + K(1, 1) * F(8) + K(1, 7) * F(2) +
+                 K(1, 6) * F(4) + K(1, 8) * F(1) + K(1, 2) * F(7)) *
+                icste2;
+      C(1, 2) = K(1, 6) * F(6) + K(1, 8) * F(8) + K(1, 2) * F(2);
+      C(5, 0) = K(5, 0) * F(0) + K(5, 3) * F(3) + K(5, 5) * F(5);
+      C(5, 3) = (K(5, 0) * F(4) + K(5, 3) * F(1) + K(5, 5) * F(7) +
+                 K(5, 4) * F(0) + K(5, 1) * F(3) + K(5, 7) * F(5)) *
+                icste2;
+      C(5, 4) = (K(5, 0) * F(6) + K(5, 3) * F(8) + K(5, 5) * F(2) +
+                 K(5, 6) * F(0) + K(5, 8) * F(3) + K(5, 2) * F(5)) *
+                icste2;
+      C(5, 1) = K(5, 4) * F(4) + K(5, 1) * F(1) + K(5, 7) * F(7);
+      C(5, 5) = (K(5, 4) * F(6) + K(5, 1) * F(8) + K(5, 7) * F(2) +
+                 K(5, 6) * F(4) + K(5, 8) * F(1) + K(5, 2) * F(7)) *
+                icste2;
+      C(5, 2) = K(5, 6) * F(6) + K(5, 8) * F(8) + K(5, 2) * F(2);
+      C(2, 0) = K(2, 0) * F(0) + K(2, 3) * F(3) + K(2, 5) * F(5);
+      C(2, 3) = (K(2, 0) * F(4) + K(2, 3) * F(1) + K(2, 5) * F(7) +
+                 K(2, 4) * F(0) + K(2, 1) * F(3) + K(2, 7) * F(5)) *
+                icste2;
+      C(2, 4) = (K(2, 0) * F(6) + K(2, 3) * F(8) + K(2, 5) * F(2) +
+                 K(2, 6) * F(0) + K(2, 8) * F(3) + K(2, 2) * F(5)) *
+                icste2;
+      C(2, 1) = K(2, 4) * F(4) + K(2, 1) * F(1) + K(2, 7) * F(7);
+      C(2, 5) = (K(2, 4) * F(6) + K(2, 1) * F(8) + K(2, 7) * F(2) +
+                 K(2, 6) * F(4) + K(2, 8) * F(1) + K(2, 2) * F(7)) *
+                icste2;
+      C(2, 2) = K(2, 6) * F(6) + K(2, 8) * F(8) + K(2, 2) * F(2);
+    }
     return C;
   }  // end of convertToTangentModuli
 
-  template <typename T2toST2Type, TensorConcept TensorType>
-  typename std::enable_if<((getSpaceDimension<T2toST2Type>() == 2u) &&
-                           (getSpaceDimension<TensorType>() == 2u) &&
-                           implementsT2toST2Concept<T2toST2Type>()),
-                          st2tost2<2u,
-                                   result_type<numeric_type<T2toST2Type>,
-                                               numeric_type<TensorType>,
-                                               OpMult>>>::type
-  convertToTangentModuli(const T2toST2Type& K, const TensorType& F) {
-    using res = result_type<numeric_type<T2toST2Type>, numeric_type<TensorType>,
-                            OpMult>;
-    TFEL_HOST_DEVICE constexpr const auto icste2 = Cste<res>::isqrt2 / 2;
-    st2tost2<2u, res> C;
-    C(0, 0) = K(0, 0) * F(0) + K(0, 3) * F(3);
-    C(0, 3) =
-        (K(0, 0) * F(4) + K(0, 3) * F(1) + K(0, 4) * F(0) + K(0, 1) * F(3)) *
-        icste2;
-    C(0, 1) = K(0, 4) * F(4) + K(0, 1) * F(1);
-    C(0, 2) = K(0, 2) * F(2);
-    C(3, 0) = K(3, 0) * F(0) + K(3, 3) * F(3);
-    C(3, 3) =
-        (K(3, 0) * F(4) + K(3, 3) * F(1) + K(3, 4) * F(0) + K(3, 1) * F(3)) *
-        icste2;
-    C(3, 1) = K(3, 4) * F(4) + K(3, 1) * F(1);
-    C(3, 2) = K(3, 2) * F(2);
-    C(1, 0) = K(1, 0) * F(0) + K(1, 3) * F(3);
-    C(1, 3) =
-        (K(1, 0) * F(4) + K(1, 3) * F(1) + K(1, 4) * F(0) + K(1, 1) * F(3)) *
-        icste2;
-    C(1, 1) = K(1, 4) * F(4) + K(1, 1) * F(1);
-    C(1, 2) = K(1, 2) * F(2);
-    C(2, 0) = K(2, 0) * F(0) + K(2, 3) * F(3);
-    C(2, 3) =
-        (K(2, 0) * F(4) + K(2, 3) * F(1) + K(2, 4) * F(0) + K(2, 1) * F(3)) *
-        icste2;
-    C(2, 1) = K(2, 4) * F(4) + K(2, 1) * F(1);
-    C(2, 2) = K(2, 2) * F(2);
-    return C;
-  }  // end of convertToTangentModuli
-
-  template <typename T2toST2Type, TensorConcept TensorType>
-  typename std::enable_if<((getSpaceDimension<T2toST2Type>() == 3u) &&
-                           (getSpaceDimension<TensorType>() == 3u) &&
-                           implementsT2toST2Concept<T2toST2Type>()),
-                          st2tost2<3u,
-                                   result_type<numeric_type<T2toST2Type>,
-                                               numeric_type<TensorType>,
-                                               OpMult>>>::type
-  convertToTangentModuli(const T2toST2Type& K, const TensorType& F) {
-    using res = result_type<numeric_type<T2toST2Type>, numeric_type<TensorType>,
-                            OpMult>;
-    TFEL_HOST_DEVICE constexpr const auto icste2 = Cste<res>::isqrt2 / 2;
-    st2tost2<3u, res> C;
-    C(0, 0) = K(0, 0) * F(0) + K(0, 3) * F(3) + K(0, 5) * F(5);
-    C(0, 3) = (K(0, 0) * F(4) + K(0, 3) * F(1) + K(0, 5) * F(7) +
-               K(0, 4) * F(0) + K(0, 1) * F(3) + K(0, 7) * F(5)) *
-              icste2;
-    C(0, 4) = (K(0, 0) * F(6) + K(0, 3) * F(8) + K(0, 5) * F(2) +
-               K(0, 6) * F(0) + K(0, 8) * F(3) + K(0, 2) * F(5)) *
-              icste2;
-    C(0, 1) = K(0, 4) * F(4) + K(0, 1) * F(1) + K(0, 7) * F(7);
-    C(0, 5) = (K(0, 4) * F(6) + K(0, 1) * F(8) + K(0, 7) * F(2) +
-               K(0, 6) * F(4) + K(0, 8) * F(1) + K(0, 2) * F(7)) *
-              icste2;
-    C(0, 2) = K(0, 6) * F(6) + K(0, 8) * F(8) + K(0, 2) * F(2);
-    C(3, 0) = K(3, 0) * F(0) + K(3, 3) * F(3) + K(3, 5) * F(5);
-    C(3, 3) = (K(3, 0) * F(4) + K(3, 3) * F(1) + K(3, 5) * F(7) +
-               K(3, 4) * F(0) + K(3, 1) * F(3) + K(3, 7) * F(5)) *
-              icste2;
-    C(3, 4) = (K(3, 0) * F(6) + K(3, 3) * F(8) + K(3, 5) * F(2) +
-               K(3, 6) * F(0) + K(3, 8) * F(3) + K(3, 2) * F(5)) *
-              icste2;
-    C(3, 1) = K(3, 4) * F(4) + K(3, 1) * F(1) + K(3, 7) * F(7);
-    C(3, 5) = (K(3, 4) * F(6) + K(3, 1) * F(8) + K(3, 7) * F(2) +
-               K(3, 6) * F(4) + K(3, 8) * F(1) + K(3, 2) * F(7)) *
-              icste2;
-    C(3, 2) = K(3, 6) * F(6) + K(3, 8) * F(8) + K(3, 2) * F(2);
-    C(4, 0) = K(4, 0) * F(0) + K(4, 3) * F(3) + K(4, 5) * F(5);
-    C(4, 3) = (K(4, 0) * F(4) + K(4, 3) * F(1) + K(4, 5) * F(7) +
-               K(4, 4) * F(0) + K(4, 1) * F(3) + K(4, 7) * F(5)) *
-              icste2;
-    C(4, 4) = (K(4, 0) * F(6) + K(4, 3) * F(8) + K(4, 5) * F(2) +
-               K(4, 6) * F(0) + K(4, 8) * F(3) + K(4, 2) * F(5)) *
-              icste2;
-    C(4, 1) = K(4, 4) * F(4) + K(4, 1) * F(1) + K(4, 7) * F(7);
-    C(4, 5) = (K(4, 4) * F(6) + K(4, 1) * F(8) + K(4, 7) * F(2) +
-               K(4, 6) * F(4) + K(4, 8) * F(1) + K(4, 2) * F(7)) *
-              icste2;
-    C(4, 2) = K(4, 6) * F(6) + K(4, 8) * F(8) + K(4, 2) * F(2);
-    C(1, 0) = K(1, 0) * F(0) + K(1, 3) * F(3) + K(1, 5) * F(5);
-    C(1, 3) = (K(1, 0) * F(4) + K(1, 3) * F(1) + K(1, 5) * F(7) +
-               K(1, 4) * F(0) + K(1, 1) * F(3) + K(1, 7) * F(5)) *
-              icste2;
-    C(1, 4) = (K(1, 0) * F(6) + K(1, 3) * F(8) + K(1, 5) * F(2) +
-               K(1, 6) * F(0) + K(1, 8) * F(3) + K(1, 2) * F(5)) *
-              icste2;
-    C(1, 1) = K(1, 4) * F(4) + K(1, 1) * F(1) + K(1, 7) * F(7);
-    C(1, 5) = (K(1, 4) * F(6) + K(1, 1) * F(8) + K(1, 7) * F(2) +
-               K(1, 6) * F(4) + K(1, 8) * F(1) + K(1, 2) * F(7)) *
-              icste2;
-    C(1, 2) = K(1, 6) * F(6) + K(1, 8) * F(8) + K(1, 2) * F(2);
-    C(5, 0) = K(5, 0) * F(0) + K(5, 3) * F(3) + K(5, 5) * F(5);
-    C(5, 3) = (K(5, 0) * F(4) + K(5, 3) * F(1) + K(5, 5) * F(7) +
-               K(5, 4) * F(0) + K(5, 1) * F(3) + K(5, 7) * F(5)) *
-              icste2;
-    C(5, 4) = (K(5, 0) * F(6) + K(5, 3) * F(8) + K(5, 5) * F(2) +
-               K(5, 6) * F(0) + K(5, 8) * F(3) + K(5, 2) * F(5)) *
-              icste2;
-    C(5, 1) = K(5, 4) * F(4) + K(5, 1) * F(1) + K(5, 7) * F(7);
-    C(5, 5) = (K(5, 4) * F(6) + K(5, 1) * F(8) + K(5, 7) * F(2) +
-               K(5, 6) * F(4) + K(5, 8) * F(1) + K(5, 2) * F(7)) *
-              icste2;
-    C(5, 2) = K(5, 6) * F(6) + K(5, 8) * F(8) + K(5, 2) * F(2);
-    C(2, 0) = K(2, 0) * F(0) + K(2, 3) * F(3) + K(2, 5) * F(5);
-    C(2, 3) = (K(2, 0) * F(4) + K(2, 3) * F(1) + K(2, 5) * F(7) +
-               K(2, 4) * F(0) + K(2, 1) * F(3) + K(2, 7) * F(5)) *
-              icste2;
-    C(2, 4) = (K(2, 0) * F(6) + K(2, 3) * F(8) + K(2, 5) * F(2) +
-               K(2, 6) * F(0) + K(2, 8) * F(3) + K(2, 2) * F(5)) *
-              icste2;
-    C(2, 1) = K(2, 4) * F(4) + K(2, 1) * F(1) + K(2, 7) * F(7);
-    C(2, 5) = (K(2, 4) * F(6) + K(2, 1) * F(8) + K(2, 7) * F(2) +
-               K(2, 6) * F(4) + K(2, 8) * F(1) + K(2, 2) * F(7)) *
-              icste2;
-    C(2, 2) = K(2, 6) * F(6) + K(2, 8) * F(8) + K(2, 2) * F(2);
-    return C;
-  }
-
 }  // end of namespace tfel::math
 
 #endif /* LIB_TFEL_MATH_CONVERTTOTANGENTMODULI_HXX */
diff --git a/include/TFEL/Math/ST2toST2/ST2toST2Concept.hxx b/include/TFEL/Math/ST2toST2/ST2toST2Concept.hxx
index a563e96b1..b9fc26262 100644
--- a/include/TFEL/Math/ST2toST2/ST2toST2Concept.hxx
+++ b/include/TFEL/Math/ST2toST2/ST2toST2Concept.hxx
@@ -54,7 +54,10 @@ namespace tfel::math {
       (requires(const T t, const unsigned short i, const unsigned short j) {
         t(i, j);
       });
-  //! \brief partial specialisation for symmetric tensors
+  /*!
+   * \brief partial specialisation for fourth order tensor transforming a
+   * symmetric tensor in another symmetric tensor
+   */
   template <typename T>
   struct ConceptRebind<ST2toST2Tag, T> {
     //! \brief a simple alias
diff --git a/include/TFEL/Math/ST2toST2/st2tost2.ixx b/include/TFEL/Math/ST2toST2/st2tost2.ixx
index 1bd809d32..bfdb26884 100644
--- a/include/TFEL/Math/ST2toST2/st2tost2.ixx
+++ b/include/TFEL/Math/ST2toST2/st2tost2.ixx
@@ -53,26 +53,25 @@ namespace tfel::math {
 
   template <unsigned short N, typename T>
   template <StensorConcept StensorType>
-  std::enable_if_t<getSpaceDimension<StensorType>() == N &&
-                       isAssignableTo<numeric_type<StensorType>, T>(),
-                   tfel::math::st2tost2<N, T>>
-  st2tost2<N, T>::stpd(const StensorType& s) {
+  TFEL_HOST_DEVICE constexpr auto
+  st2tost2<N, T>::stpd(const StensorType& s) noexcept requires(
+      getSpaceDimension<StensorType>() == N &&
+      isAssignableTo<numeric_type<StensorType>, T>()) {
     return StensorSymmetricProductDerivative<N, T>::exe(s);
   }
 
   template <unsigned short N, typename T>
-  template <typename T2toST2Type>
-  std::enable_if_t<implementsT2toST2Concept<T2toST2Type>() &&
-                       getSpaceDimension<T2toST2Type>() == N &&
-                       isAssignableTo<numeric_type<T2toST2Type>, T>(),
-                   Expr<st2tost2<N, T>, ConvertT2toST2ToST2toST2Expr<N>>>
-  st2tost2<N, T>::convert(const T2toST2Type& src) {
+  template <T2toST2Concept T2toST2Type>
+  TFEL_HOST_DEVICE constexpr auto
+  st2tost2<N, T>::convert(const T2toST2Type& src) noexcept requires(
+      getSpaceDimension<T2toST2Type>() == N &&
+      isAssignableTo<numeric_type<T2toST2Type>, T>()) {
     return Expr<st2tost2<N, T>, ConvertT2toST2ToST2toST2Expr<N>>(src);
-  }
+  }  // end of convert
 
   template <unsigned short N, typename T>
-  st2tost2<N, base_type<T>> st2tost2<N, T>::fromRotationMatrix(
-      const rotation_matrix<T>& r) {
+  TFEL_HOST_DEVICE constexpr auto st2tost2<N, T>::fromRotationMatrix(
+      const rotation_matrix<T>& r) noexcept {
     return st2tost2_internals::BuildFromRotationMatrix<N, T>::exe(r);
   }  // end of st2tost2<N,T>::fromRotationMatrix
 
@@ -245,21 +244,21 @@ namespace tfel::math {
   }  // end of invert
 
   template <ST2toST2Concept ST2toST2Type>
-  TFEL_HOST_DEVICE constexpr st2tost2<getSpaceDimension<ST2toST2Type>(),
-                                      numeric_type<ST2toST2Type>>
-  change_basis(const ST2toST2Type& s,
-               const rotation_matrix<numeric_type<ST2toST2Type>>& r) noexcept {
+  TFEL_HOST_DEVICE constexpr auto change_basis(
+      const ST2toST2Type& s,
+      const rotation_matrix<numeric_type<ST2toST2Type>>& r) noexcept {
     constexpr auto N = getSpaceDimension<ST2toST2Type>();
+    using Result = st2tost2<N, numeric_type<ST2toST2Type>>;
     if constexpr (N == 1) {
-      return s;
+      return Result{s};
     } else {
       using st2tost2 =
           tfel::math::st2tost2<N, tfel::math::numeric_type<ST2toST2Type>>;
       const auto sr = st2tost2::fromRotationMatrix(r);
       const auto sir = st2tost2::fromRotationMatrix(transpose(r));
-      return sr * s * sir;
+      return Result{sr * s * sir};
     }
-  }  // end of change_basie
+  }  // end of change_basis
 
   template <ST2toST2Concept ST2toST2Type, TensorConcept TensorType>
   TFEL_HOST_DEVICE constexpr auto
diff --git a/include/TFEL/Math/ST2toT2/ST2toT2Concept.hxx b/include/TFEL/Math/ST2toT2/ST2toT2Concept.hxx
index 2133b121d..6ac2571c5 100644
--- a/include/TFEL/Math/ST2toT2/ST2toT2Concept.hxx
+++ b/include/TFEL/Math/ST2toT2/ST2toT2Concept.hxx
@@ -29,41 +29,52 @@ namespace tfel::math {
    * \brief Helper class to characterise st2tot2.
    */
   struct ST2toT2Tag {};
-
+  /*!
+   * \brief an helper class that simply exposes publically a member named
+   * ConceptTag as an alias to T2toT2Tag.
+   *
+   * The main reason for this alias is to properly implement the `ConceptRebind`
+   * metafunction.
+   */
   template <class T>
-  struct ST2toT2Concept {
-    typedef ST2toT2Tag ConceptTag;
-
-   protected:
-    ST2toT2Concept() = default;
-    ST2toT2Concept(ST2toT2Concept&&) = default;
-    ST2toT2Concept(const ST2toT2Concept&) = default;
-    ST2toT2Concept& operator=(const ST2toT2Concept&) = default;
-    ~ST2toT2Concept() = default;
+  struct ST2toT2ConceptBase {
+    using ConceptTag = ST2toT2Tag;
   };
-
+  /*!
+   * \brief definition of the ST2toT2Concept
+   *
+   * A class matching this concept must expose the `ST2toT2Tag` and have
+   * appropriate access operators.
+   */
+  template <typename T>
+  concept ST2toT2Concept =
+      (std::is_same_v<typename std::decay_t<T>::ConceptTag, ST2toT2Tag>)&&  //
+      (requires(const T t, const unsigned short i, const unsigned short j) {
+        t(i, j);
+      });
+  //! \brief partial specialisation for fourth order tensor
+  template <typename Type>
+  struct ConceptRebind<ST2toT2Tag, Type> {
+    using type = ST2toT2ConceptBase<Type>;
+  };
+  /*!
+   * \return the sum of the absolute values of all components of an
+   * linear application transforming a tensor in a tensor
+   * \param[in] t: fourth order tensor
+   */
+  constexpr auto abs(const ST2toT2Concept auto&) noexcept;
   /*!
    * \brief an helper function which returns if the given type implements the
    * `ST2toT2Concept`.
    * \tparam ST2toT2Type: type tested
+   * \note function given for backward compatibility with versions prior
+   * to 5.0
    */
   template <typename ST2toT2Type>
-  TFEL_HOST_DEVICE constexpr bool implementsST2toT2Concept() {
-    return tfel::meta::implements<ST2toT2Type, ST2toT2Concept>();
+  [[deprecated]] TFEL_HOST_DEVICE constexpr bool implementsST2toT2Concept() {
+    return ST2toT2Concept<ST2toT2Type>;
   }  // end of implementsST2toT2Concept
 
-  //! paratial specialisation for symmetric tensors
-  template <typename Type>
-  struct ConceptRebind<ST2toT2Tag, Type> {
-    using type = ST2toT2Concept<Type>;
-  };
-
-  template <typename ST2toT2Type>
-  typename std::enable_if<
-      implementsST2toT2Concept<ST2toT2Type>(),
-      typename tfel::typetraits::AbsType<numeric_type<ST2toT2Type>>::type>::type
-  abs(const ST2toT2Type&);
-
 }  // end of namespace tfel::math
 
 #include "TFEL/Math/ST2toT2/ST2toT2Concept.ixx"
diff --git a/include/TFEL/Math/ST2toT2/ST2toT2Concept.ixx b/include/TFEL/Math/ST2toT2/ST2toT2Concept.ixx
index 3ea84ecce..27d57413f 100644
--- a/include/TFEL/Math/ST2toT2/ST2toT2Concept.ixx
+++ b/include/TFEL/Math/ST2toT2/ST2toT2Concept.ixx
@@ -19,21 +19,18 @@
 
 namespace tfel::math {
 
-  template <typename ST2toT2Type>
-  typename std::enable_if<
-      implementsST2toT2Concept<ST2toT2Type>(),
-      typename tfel::typetraits::AbsType<numeric_type<ST2toT2Type>>::type>::type
-  abs(const ST2toT2Type& v) {
+  TFEL_HOST_DEVICE constexpr auto abs(const ST2toT2Concept auto& t) noexcept {
+    using ST2toT2Type = decltype(t);
     using NumType = numeric_type<ST2toT2Type>;
     using AbsNumType = typename tfel::typetraits::AbsType<NumType>::type;
     constexpr auto tsize =
         TensorDimeToSize<getSpaceDimension<ST2toT2Type>()>::value;
     constexpr auto ssize =
         StensorDimeToSize<getSpaceDimension<ST2toT2Type>()>::value;
-    AbsNumType a(0);
+    auto a = AbsNumType{};
     for (unsigned short i = 0; i < tsize; ++i) {
       for (unsigned short j = 0; j < ssize; ++j) {
-        a += abs(v(i, j));
+        a += abs(t(i, j));
       }
     }
     return a;
diff --git a/include/TFEL/Math/ST2toT2/ST2toT2ConceptIO.hxx b/include/TFEL/Math/ST2toT2/ST2toT2ConceptIO.hxx
index ea77f9b27..9022333d9 100644
--- a/include/TFEL/Math/ST2toT2/ST2toT2ConceptIO.hxx
+++ b/include/TFEL/Math/ST2toT2/ST2toT2ConceptIO.hxx
@@ -20,9 +20,9 @@
 namespace tfel::math {
 
   // Serialisation operator
-  template <typename ST2toT2Type>
-  std::enable_if_t<implementsST2toT2Concept<ST2toT2Type>(), std::ostream&>
-  operator<<(std::ostream& os, const ST2toT2Type& s) {
+  TFEL_HOST std::ostream& operator<<(std::ostream& os,
+                                     const ST2toT2Concept auto& t) {
+    using ST2toT2Type = decltype(t);
     constexpr auto N = getSpaceDimension<ST2toT2Type>();
     constexpr auto stensor_size = StensorDimeToSize<N>::value;
     constexpr auto tensor_size = TensorDimeToSize<N>::value;
@@ -34,7 +34,7 @@ namespace tfel::math {
         os << "[";
       }
       for (unsigned short j = 0; j < stensor_size;) {
-        os << s(i, j);
+        os << t(i, j);
         if (++j != stensor_size) {
           os << ",";
         }
diff --git a/include/TFEL/Math/ST2toT2/ST2toT2ConceptOperations.hxx b/include/TFEL/Math/ST2toT2/ST2toT2ConceptOperations.hxx
index 4d5032798..ade2d7175 100644
--- a/include/TFEL/Math/ST2toT2/ST2toT2ConceptOperations.hxx
+++ b/include/TFEL/Math/ST2toT2/ST2toT2ConceptOperations.hxx
@@ -177,15 +177,13 @@ namespace tfel::math {
                                       Expr<Result, UnaryOperation<A, OpNeg>>>;
   };
 
-  template <TensorConcept T1, typename T2>
-  TFEL_MATH_INLINE typename std::enable_if<
-      implementsST2toT2Concept<T2>() &&
-          !isInvalid<BinaryOperationResult<T1, T2, OpMult>>(),
-      BinaryOperationHandler<T1, T2, OpMult>>::type
-  operator|(const T1& a, const T2& b) {
+  template <TensorConcept T1, ST2toT2Concept T2>
+  TFEL_HOST_DEVICE constexpr auto
+  operator|(const T1& a, const T2& b) noexcept requires(
+      !isInvalid<BinaryOperationResult<T1, T2, OpMult>>()) {
     using Handle = BinaryOperationHandler<T1, T2, OpMult>;
     return Handle(a, b);
-  }
+  }  // end of operator|
 
 }  // end of namespace tfel::math
 
diff --git a/include/TFEL/Math/ST2toT2/ST2toT2ST2toST2ProductExpr.hxx b/include/TFEL/Math/ST2toT2/ST2toT2ST2toST2ProductExpr.hxx
index ce90d04cf..f888180bf 100644
--- a/include/TFEL/Math/ST2toT2/ST2toT2ST2toST2ProductExpr.hxx
+++ b/include/TFEL/Math/ST2toT2/ST2toT2ST2toST2ProductExpr.hxx
@@ -15,26 +15,25 @@
 #define LIB_TFEL_MATH_ST2TOT2ST2TOST2PRODUCTEXPR_HXX
 
 #include "TFEL/Config/TFELConfig.hxx"
-
 #include "TFEL/Math/General/EmptyRunTimeProperties.hxx"
 #include "TFEL/Math/ST2toT2/ST2toT2Concept.hxx"
 #include "TFEL/Math/ST2toST2/ST2toST2Concept.hxx"
 
 namespace tfel::math {
 
-  //! Empty structure used for partial specialisation of the
-  //! Expr class
+  /*!
+   * \brief empty structure used for partial specialisation of the
+   * Expr class
+   */
   template <unsigned short N>
   struct TFEL_VISIBILITY_LOCAL ST2toT2ST2toST2ProductExpr {
   };  // end of struct ST2toT2ST2toST2ProductExpr
 
-  /*!
-   * Partial specialisation
-   */
-  template <typename ST2toT2ResultType>
+  //! \brief partial specialisation in 1D
+  template <ST2toT2Concept ST2toT2ResultType>
   struct TFEL_VISIBILITY_LOCAL
       Expr<ST2toT2ResultType, ST2toT2ST2toST2ProductExpr<1u>>
-      : public ST2toT2Concept<
+      : public ST2toT2ConceptBase<
             Expr<ST2toT2ResultType, ST2toT2ST2toST2ProductExpr<1u>>>,
         public array_holder<
             TensorDimeToSize<getSpaceDimension<ST2toT2ResultType>()>::value *
@@ -50,10 +49,9 @@ namespace tfel::math {
      * \param[in] a : first term of the product
      * \param[in] b : second term of the product
      */
-    template <typename ST2toT2Type, ST2toST2Concept ST2toST2Type2>
+    template <ST2toT2Concept ST2toT2Type, ST2toST2Concept ST2toST2Type2>
     TFEL_HOST_DEVICE constexpr Expr(const ST2toT2Type& a,
                                     const ST2toST2Type2& b) noexcept {
-      static_assert(implementsST2toT2Concept<ST2toT2Type>());
       static_assert(getSpaceDimension<ST2toT2Type>() == 1u);
       static_assert(getSpaceDimension<ST2toST2Type2>() == 1u);
       this->v[0] = a(0, 0) * b(0, 0) + a(0, 1) * b(1, 0) + a(0, 2) * b(2, 0);
@@ -75,21 +73,17 @@ namespace tfel::math {
         const unsigned short i, const unsigned short j) const noexcept {
       return this->v[i * 3 + j];
     }  // end of operator()
-    /*!
-     * \return the runtime properties
-     */
+    //! \return the runtime properties
     TFEL_HOST_DEVICE constexpr auto getRunTimeProperties() const noexcept {
       return RunTimeProperties();
     }
   };  // end of struct Expr<ST2toT2ResultType,ST2toT2ST2toST2ProductExpr>
 
-  /*!
-   * Partial specialisation
-   */
-  template <typename ST2toT2ResultType>
+  //! \brief partial specialisation in 2D
+  template <ST2toT2Concept ST2toT2ResultType>
   struct TFEL_VISIBILITY_LOCAL
       Expr<ST2toT2ResultType, ST2toT2ST2toST2ProductExpr<2u>>
-      : public ST2toT2Concept<
+      : public ST2toT2ConceptBase<
             Expr<ST2toT2ResultType, ST2toT2ST2toST2ProductExpr<2u>>>,
         public array_holder<
             TensorDimeToSize<getSpaceDimension<ST2toT2ResultType>()>::value *
@@ -105,10 +99,9 @@ namespace tfel::math {
      * \param[in] a : first term of the product
      * \param[in] b : second term of the product
      */
-    template <typename ST2toT2Type, ST2toST2Concept ST2toST2Type2>
+    template <ST2toT2Concept ST2toT2Type, ST2toST2Concept ST2toST2Type2>
     TFEL_HOST_DEVICE constexpr Expr(const ST2toT2Type& a,
                                     const ST2toST2Type2& b) noexcept {
-      static_assert(implementsST2toT2Concept<ST2toT2Type>());
       static_assert(getSpaceDimension<ST2toT2Type>() == 2u);
       static_assert(getSpaceDimension<ST2toST2Type2>() == 2u);
       this->v[0] = a(0, 0) * b(0, 0) + a(0, 1) * b(1, 0) + a(0, 2) * b(2, 0) +
@@ -161,21 +154,17 @@ namespace tfel::math {
         const unsigned short i, const unsigned short j) const noexcept {
       return this->v[i * 4 + j];
     }  // end of operator()
-    /*!
-     * \return the runtime properties
-     */
+    //! \return the runtime properties
     TFEL_HOST_DEVICE constexpr auto getRunTimeProperties() const noexcept {
       return RunTimeProperties();
     }
   };  // end of struct Expr<ST2toT2ResultType,ST2toT2ST2toST2ProductExpr>
 
-  /*!
-   * Partial specialisation
-   */
-  template <typename ST2toT2ResultType>
+  //! \brief partial specialisation in 3D
+  template <ST2toT2Concept ST2toT2ResultType>
   struct TFEL_VISIBILITY_LOCAL
       Expr<ST2toT2ResultType, ST2toT2ST2toST2ProductExpr<3u>>
-      : public ST2toT2Concept<
+      : public ST2toT2ConceptBase<
             Expr<ST2toT2ResultType, ST2toT2ST2toST2ProductExpr<3u>>>,
         public array_holder<
             TensorDimeToSize<getSpaceDimension<ST2toT2ResultType>()>::value *
@@ -191,10 +180,9 @@ namespace tfel::math {
      * \param[in] a : first term of the product
      * \param[in] b : second term of the product
      */
-    template <typename ST2toT2Type, ST2toST2Concept ST2toST2Type2>
+    template <ST2toT2Concept ST2toT2Type, ST2toST2Concept ST2toST2Type2>
     TFEL_HOST_DEVICE constexpr Expr(const ST2toT2Type& a,
                                     const ST2toST2Type2& b) noexcept {
-      static_assert(implementsST2toT2Concept<ST2toT2Type>());
       static_assert(getSpaceDimension<ST2toT2Type>() == 3u);
       static_assert(getSpaceDimension<ST2toST2Type2>() == 3u);
       this->v[0] = a(0, 0) * b(0, 0) + a(0, 1) * b(1, 0) + a(0, 2) * b(2, 0) +
@@ -315,9 +303,7 @@ namespace tfel::math {
         const unsigned short i, const unsigned short j) const noexcept {
       return this->v[i * 6 + j];
     }  // end of operator()
-    /*!
-     * \return the runtime properties
-     */
+    //! \return the runtime properties
     TFEL_HOST_DEVICE constexpr auto getRunTimeProperties() const noexcept {
       return RunTimeProperties();
     }
diff --git a/include/TFEL/Math/ST2toT2/ST2toT2StensorProductExpr.hxx b/include/TFEL/Math/ST2toT2/ST2toT2StensorProductExpr.hxx
index 747a09d3f..e1395974d 100644
--- a/include/TFEL/Math/ST2toT2/ST2toT2StensorProductExpr.hxx
+++ b/include/TFEL/Math/ST2toT2/ST2toT2StensorProductExpr.hxx
@@ -28,9 +28,7 @@ namespace tfel::math {
   struct ST2toT2StensorProductExpr {
   };  // end of struct ST2toT2StensorProductExpr
 
-  /*!
-   * Partial specialisation
-   */
+  //! \brief partial specialisation in 1D
   template <TensorConcept TensorResultType>
   struct Expr<TensorResultType, ST2toT2StensorProductExpr<1u>>
       : public TensorConceptBase<
@@ -47,9 +45,9 @@ namespace tfel::math {
      * \param[in] a : first term of the product
      * \param[in] b : second term of the product
      */
-    template <typename ST2toT2Type, StensorConcept StensorType>
-    TFEL_MATH_INLINE Expr(const ST2toT2Type& a, const StensorType& b) {
-      static_assert(implementsST2toT2Concept<ST2toT2Type>());
+    template <ST2toT2Concept ST2toT2Type, StensorConcept StensorType>
+    TFEL_HOST_DEVICE constexpr Expr(const ST2toT2Type& a,
+                                    const StensorType& b) noexcept {
       static_assert(getSpaceDimension<ST2toT2Type>() == 1u);
       static_assert(getSpaceDimension<StensorType>() == 1u);
       this->v[0] = a(0, 0) * b(0) + a(0, 1) * b(1) + a(0, 2) * b(2);
@@ -60,21 +58,27 @@ namespace tfel::math {
      * \brief access operator
      * \param[in] i : index
      */
-    TFEL_MATH_INLINE const value_type& operator()(
-        const unsigned short i) const {
+    TFEL_HOST_DEVICE constexpr const value_type& operator[](
+        const unsigned short i) const noexcept {
+      return this->operator()(i);
+    }  // end of operator()
+    /*!
+     * \brief access operator
+     * \param[in] i : index
+     */
+    TFEL_HOST_DEVICE constexpr const value_type& operator()(
+        const unsigned short i) const noexcept {
       return this->v[i];
     }  // end of operator()
     /*!
      * \return the runtime properties
      */
-    TFEL_MATH_INLINE RunTimeProperties getRunTimeProperties() const {
+    TFEL_HOST_DEVICE constexpr auto getRunTimeProperties() const noexcept {
       return RunTimeProperties();
     }
   };  // end of struct Expr<TensorResultType,T2SST2toT2StensorProductExp>
 
-  /*!
-   * Partial specialisation
-   */
+  //! \brief partial specialisation in 2D
   template <TensorConcept TensorResultType>
   struct Expr<TensorResultType, ST2toT2StensorProductExpr<2u>>
       : public TensorConceptBase<
@@ -91,9 +95,9 @@ namespace tfel::math {
      * \param[in] a : first term of the product
      * \param[in] b : second term of the product
      */
-    template <typename ST2toT2Type, StensorConcept StensorType>
-    TFEL_MATH_INLINE Expr(const ST2toT2Type& a, const StensorType& b) {
-      static_assert(implementsST2toT2Concept<ST2toT2Type>());
+    template <ST2toT2Concept ST2toT2Type, StensorConcept StensorType>
+    TFEL_HOST_DEVICE constexpr Expr(const ST2toT2Type& a,
+                                    const StensorType& b) noexcept {
       static_assert(getSpaceDimension<ST2toT2Type>() == 2u);
       static_assert(getSpaceDimension<StensorType>() == 2u);
       this->v[0] =
@@ -111,21 +115,27 @@ namespace tfel::math {
      * \brief access operator
      * \param[in] i : index
      */
-    TFEL_MATH_INLINE const value_type& operator()(
+    TFEL_HOST_DEVICE constexpr const value_type& operator[](
+        const unsigned short i) const noexcept {
+      return this->operator()(i);
+    }  // end of operator()
+    /*!
+     * \brief access operator
+     * \param[in] i : index
+     */
+    TFEL_HOST_DEVICE constexpr const value_type& operator()(
         const unsigned short i) const {
       return this->v[i];
     }  // end of operator()
     /*!
      * \return the runtime properties
      */
-    TFEL_MATH_INLINE RunTimeProperties getRunTimeProperties() const {
+    TFEL_HOST_DEVICE constexpr auto getRunTimeProperties() const {
       return RunTimeProperties();
     }
   };  // end of struct Expr<TensorResultType,T2SST2toT2StensorProductExp>
 
-  /*!
-   * Partial specialisation
-   */
+  //! \brief partial specialisation in 3D
   template <TensorConcept TensorResultType>
   struct Expr<TensorResultType, ST2toT2StensorProductExpr<3u>>
       : public TensorConceptBase<
@@ -142,9 +152,9 @@ namespace tfel::math {
      * \param[in] a : first term of the product
      * \param[in] b : second term of the product
      */
-    template <typename ST2toT2Type, StensorConcept StensorType>
-    TFEL_MATH_INLINE Expr(const ST2toT2Type& a, const StensorType& b) {
-      static_assert(implementsST2toT2Concept<ST2toT2Type>());
+    template <ST2toT2Concept ST2toT2Type, StensorConcept StensorType>
+    TFEL_HOST_DEVICE constexpr Expr(const ST2toT2Type& a,
+                                    const StensorType& b) noexcept {
       static_assert(getSpaceDimension<ST2toT2Type>() == 3u);
       static_assert(getSpaceDimension<StensorType>() == 3u);
       this->v[0] = a(0, 0) * b(0) + a(0, 1) * b(1) + a(0, 2) * b(2) +
@@ -170,14 +180,22 @@ namespace tfel::math {
      * \brief access operator
      * \param[in] i : index
      */
-    TFEL_MATH_INLINE const value_type& operator()(
-        const unsigned short i) const {
+    TFEL_HOST_DEVICE constexpr const value_type& operator[](
+        const unsigned short i) const noexcept {
+      return this->operator()(i);
+    }  // end of operator()
+    /*!
+     * \brief access operator
+     * \param[in] i : index
+     */
+    TFEL_HOST_DEVICE constexpr const value_type& operator()(
+        const unsigned short i) const noexcept {
       return this->v[i];
     }  // end of operator()
     /*!
      * \return the runtime properties
      */
-    TFEL_MATH_INLINE RunTimeProperties getRunTimeProperties() const {
+    TFEL_HOST_DEVICE constexpr auto getRunTimeProperties() const noexcept {
       return RunTimeProperties();
     }
   };  // end of struct Expr<TensorResultType,T2SST2toT2StensorProductExp>
diff --git a/include/TFEL/Math/ST2toT2/ST2toT2T2toST2ProductExpr.hxx b/include/TFEL/Math/ST2toT2/ST2toT2T2toST2ProductExpr.hxx
index ef1c36b94..077cf83e5 100644
--- a/include/TFEL/Math/ST2toT2/ST2toT2T2toST2ProductExpr.hxx
+++ b/include/TFEL/Math/ST2toT2/ST2toT2T2toST2ProductExpr.hxx
@@ -28,9 +28,7 @@ namespace tfel::math {
   struct TFEL_VISIBILITY_LOCAL ST2toT2T2toST2ProductExpr {
   };  // end of struct ST2toT2T2toST2ProductExpr
 
-  /*!
-   * Partial specialisation
-   */
+  //! \brief partial specialisation in 1D
   template <typename T2toT2ResultType>
   struct TFEL_VISIBILITY_LOCAL
       Expr<T2toT2ResultType, ST2toT2T2toST2ProductExpr<1u>>
@@ -49,10 +47,9 @@ namespace tfel::math {
      * \param[in] a : first term of the product
      * \param[in] b : second term of the product
      */
-    template <typename ST2toT2Type, typename T2toST2Type2>
-    TFEL_MATH_INLINE Expr(const ST2toT2Type& a, const T2toST2Type2& b) {
-      static_assert(implementsST2toT2Concept<ST2toT2Type>());
-      static_assert(implementsT2toST2Concept<T2toST2Type2>());
+    template <typename ST2toT2Type, T2toST2Concept T2toST2Type2>
+    TFEL_HOST_DEVICE constexpr Expr(const ST2toT2Type& a,
+                                    const T2toST2Type2& b) noexcept {
       static_assert(getSpaceDimension<ST2toT2Type>() == 1u);
       static_assert(getSpaceDimension<T2toST2Type2>() == 1u);
       this->v[0] = a(0, 0) * b(0, 0) + a(0, 1) * b(1, 0) + a(0, 2) * b(2, 0);
@@ -70,21 +67,17 @@ namespace tfel::math {
      * \param[in] i : line index
      * \param[in] j : column index
      */
-    TFEL_MATH_INLINE const value_type& operator()(
-        const unsigned short i, const unsigned short j) const {
+    TFEL_HOST_DEVICE constexpr const value_type& operator()(
+        const unsigned short i, const unsigned short j) const noexcept {
       return this->v[i * 3 + j];
     }  // end of operator()
-    /*!
-     * \return the runtime properties
-     */
-    TFEL_MATH_INLINE RunTimeProperties getRunTimeProperties() const {
+    //! \return the runtime properties
+    TFEL_HOST_DEVICE constexpr auto getRunTimeProperties() const noexcept {
       return RunTimeProperties();
     }
-  };  // end of struct Expr<T2toT2ResultType,ST2toT2T2toST2ProductExpr>
+  };  // end of Expr<T2toT2ResultType,ST2toT2T2toST2ProductExpr>
 
-  /*!
-   * Partial specialisation
-   */
+  //! \brief partial specialisation in 2D
   template <typename T2toT2ResultType>
   struct TFEL_VISIBILITY_LOCAL
       Expr<T2toT2ResultType, ST2toT2T2toST2ProductExpr<2u>>
@@ -103,10 +96,9 @@ namespace tfel::math {
      * \param[in] a : first term of the product
      * \param[in] b : second term of the product
      */
-    template <typename ST2toT2Type, typename T2toST2Type2>
-    TFEL_MATH_INLINE Expr(const ST2toT2Type& a, const T2toST2Type2& b) {
-      static_assert(implementsST2toT2Concept<ST2toT2Type>());
-      static_assert(implementsT2toST2Concept<T2toST2Type2>());
+    template <typename ST2toT2Type, T2toST2Concept T2toST2Type2>
+    TFEL_HOST_DEVICE constexpr Expr(const ST2toT2Type& a,
+                                    const T2toST2Type2& b) noexcept {
       static_assert(getSpaceDimension<ST2toT2Type>() == 2u);
       static_assert(getSpaceDimension<T2toST2Type2>() == 2u);
       this->v[0] = a(0, 0) * b(0, 0) + a(0, 1) * b(1, 0) + a(0, 2) * b(2, 0) +
@@ -165,21 +157,17 @@ namespace tfel::math {
      * \param[in] i : line index
      * \param[in] j : column index
      */
-    TFEL_MATH_INLINE const value_type& operator()(
-        const unsigned short i, const unsigned short j) const {
+    TFEL_HOST_DEVICE constexpr const value_type& operator()(
+        const unsigned short i, const unsigned short j) const noexcept {
       return this->v[i * 5 + j];
     }  // end of operator()
-    /*!
-     * \return the runtime properties
-     */
-    TFEL_MATH_INLINE RunTimeProperties getRunTimeProperties() const {
+    //! \return the runtime properties
+    TFEL_HOST_DEVICE constexpr auto getRunTimeProperties() const noexcept {
       return RunTimeProperties();
     }
-  };  // end of struct Expr<T2toT2ResultType,ST2toT2T2toST2ProductExpr>
+  };  // end of Expr<T2toT2ResultType,ST2toT2T2toST2ProductExpr>
 
-  /*!
-   * Partial specialisation
-   */
+  //! \brief partial specialisation in 3D
   template <typename T2toT2ResultType>
   struct TFEL_VISIBILITY_LOCAL
       Expr<T2toT2ResultType, ST2toT2T2toST2ProductExpr<3u>>
@@ -198,10 +186,9 @@ namespace tfel::math {
      * \param[in] a : first term of the product
      * \param[in] b : second term of the product
      */
-    template <typename ST2toT2Type, typename T2toST2Type2>
-    TFEL_MATH_INLINE Expr(const ST2toT2Type& a, const T2toST2Type2& b) {
-      static_assert(implementsST2toT2Concept<ST2toT2Type>());
-      static_assert(implementsT2toST2Concept<T2toST2Type2>());
+    template <typename ST2toT2Type, T2toST2Concept T2toST2Type2>
+    TFEL_HOST_DEVICE constexpr Expr(const ST2toT2Type& a,
+                                    const T2toST2Type2& b) noexcept {
       static_assert(getSpaceDimension<ST2toT2Type>() == 3u);
       static_assert(getSpaceDimension<T2toST2Type2>() == 3u);
       this->v[0] = a(0, 0) * b(0, 0) + a(0, 1) * b(1, 0) + a(0, 2) * b(2, 0) +
@@ -372,17 +359,15 @@ namespace tfel::math {
      * \param[in] i : line index
      * \param[in] j : column index
      */
-    TFEL_MATH_INLINE const value_type& operator()(
-        const unsigned short i, const unsigned short j) const {
+    TFEL_HOST_DEVICE constexpr const value_type& operator()(
+        const unsigned short i, const unsigned short j) const noexcept {
       return this->v[i * 9 + j];
     }  // end of operator()
-    /*!
-     * \return the runtime properties
-     */
-    TFEL_MATH_INLINE RunTimeProperties getRunTimeProperties() const {
+    //! \return the runtime properties
+    TFEL_HOST_DEVICE constexpr auto getRunTimeProperties() const noexcept {
       return RunTimeProperties();
     }
-  };  // end of struct Expr<T2toT2ResultType,ST2toT2T2toST2ProductExpr>
+  };  // end of Expr<T2toT2ResultType,ST2toT2T2toST2ProductExpr>
 
 }  // end of namespace tfel::math
 
diff --git a/include/TFEL/Math/ST2toT2/StensorProductLeftDerivativeExpr.hxx b/include/TFEL/Math/ST2toT2/StensorProductLeftDerivativeExpr.hxx
index 24d0686a3..6aad7e04b 100644
--- a/include/TFEL/Math/ST2toT2/StensorProductLeftDerivativeExpr.hxx
+++ b/include/TFEL/Math/ST2toT2/StensorProductLeftDerivativeExpr.hxx
@@ -23,12 +23,10 @@ namespace tfel::math {
   struct StensorProductLeftDerivativeExpr {
   };  // end of struct StensorProductLeftDerivativeExpr
 
-  /*!
-   * Partial specialisation for 1D tensor
-   */
+  //! \brief partial specialisation in 2D
   template <typename ST2toT2ResultType>
   struct Expr<ST2toT2ResultType, StensorProductLeftDerivativeExpr<1u>>
-      : public ST2toT2Concept<
+      : public ST2toT2ConceptBase<
             Expr<ST2toT2ResultType, StensorProductLeftDerivativeExpr<1u>>>,
         public array_holder<9u, numeric_type<ST2toT2ResultType>> {
     static_assert(getSpaceDimension<ST2toT2ResultType>() == 1u);
@@ -96,15 +94,12 @@ namespace tfel::math {
     TFEL_HOST_DEVICE constexpr auto getRunTimeProperties() const noexcept {
       return RunTimeProperties();
     }
-  };  // end of struct
-      // Expr<ST2toT2ResultType,StensorProductLeftDerivativeExpr<1u> >
+  };  // end of Expr<ST2toT2ResultType,StensorProductLeftDerivativeExpr<1u>>
 
-  /*!
-   * Partial specialisation for 2D tensor
-   */
+  //! \brief partial specialisation in 2D
   template <typename ST2toT2ResultType>
   struct Expr<ST2toT2ResultType, StensorProductLeftDerivativeExpr<2u>>
-      : public ST2toT2Concept<
+      : public ST2toT2ConceptBase<
             Expr<ST2toT2ResultType, StensorProductLeftDerivativeExpr<2u>>>,
         public array_holder<20u, numeric_type<ST2toT2ResultType>> {
     static_assert(getSpaceDimension<ST2toT2ResultType>() == 2u);
@@ -203,15 +198,12 @@ namespace tfel::math {
     TFEL_HOST_DEVICE constexpr auto getRunTimeProperties() const noexcept {
       return RunTimeProperties();
     }
-  };  // end of struct
-      // Expr<ST2toT2ResultType,StensorProductLeftDerivativeExpr<2u> >
+  };  // end of Expr<ST2toT2ResultType,StensorProductLeftDerivativeExpr<2u>>
 
-  /*!
-   * Partial specialisation for 3D tensor
-   */
+  //! \brief partial specialisation in 3D
   template <typename ST2toT2ResultType>
   struct Expr<ST2toT2ResultType, StensorProductLeftDerivativeExpr<3u>>
-      : public ST2toT2Concept<
+      : public ST2toT2ConceptBase<
             Expr<ST2toT2ResultType, StensorProductLeftDerivativeExpr<3u>>>,
         public array_holder<54u, numeric_type<ST2toT2ResultType>> {
     static_assert(getSpaceDimension<ST2toT2ResultType>() == 3u);
@@ -431,8 +423,7 @@ namespace tfel::math {
     TFEL_HOST_DEVICE constexpr auto getRunTimeProperties() const noexcept {
       return RunTimeProperties();
     }
-  };  // end of struct
-      // Expr<ST2toT2ResultType,StensorProductLeftDerivativeExpr<2u> >
+  };  // end of Expr<ST2toT2ResultType,StensorProductLeftDerivativeExpr<3u>>
 
 }  // namespace tfel::math
 
diff --git a/include/TFEL/Math/ST2toT2/StensorProductRightDerivativeExpr.hxx b/include/TFEL/Math/ST2toT2/StensorProductRightDerivativeExpr.hxx
index b28847873..0fd914f3f 100644
--- a/include/TFEL/Math/ST2toT2/StensorProductRightDerivativeExpr.hxx
+++ b/include/TFEL/Math/ST2toT2/StensorProductRightDerivativeExpr.hxx
@@ -20,18 +20,17 @@
 namespace tfel::math {
 
   /*!
-   * Empty structure allowing partial specialisation
+   * \brief empty structure allowing partial specialisation
+   * of the Expr class
    */
   template <unsigned short N>
   struct StensorProductRightDerivativeExpr {
   };  // end of struct StensorProductRightDerivativeExpr
 
-  /*!
-   * Partial specialisation for 1D stensor
-   */
-  template <typename ST2toT2ResultType>
+  //! \brief partial specialisation in 1D
+  template <ST2toT2Concept ST2toT2ResultType>
   struct Expr<ST2toT2ResultType, StensorProductRightDerivativeExpr<1u>>
-      : public ST2toT2Concept<
+      : public ST2toT2ConceptBase<
             Expr<ST2toT2ResultType, StensorProductRightDerivativeExpr<1u>>>,
         public array_holder<9u, numeric_type<ST2toT2ResultType>> {
     static_assert(getSpaceDimension<ST2toT2ResultType>() == 1u);
@@ -99,15 +98,12 @@ namespace tfel::math {
     TFEL_HOST_DEVICE constexpr auto getRunTimeProperties() const noexcept {
       return RunTimeProperties();
     }
-  };  // end of struct
-      // Expr<ST2toT2ResultType,StensorProductRightDerivativeExpr<1u> >
+  };  // end of Expr<ST2toT2ResultType,StensorProductRightDerivativeExpr<1u>>
 
-  /*!
-   * Partial specialisation for 2D stensor
-   */
-  template <typename ST2toT2ResultType>
+  //! \brief partial specialisation in 3D
+  template <ST2toT2Concept ST2toT2ResultType>
   struct Expr<ST2toT2ResultType, StensorProductRightDerivativeExpr<2u>>
-      : public ST2toT2Concept<
+      : public ST2toT2ConceptBase<
             Expr<ST2toT2ResultType, StensorProductRightDerivativeExpr<2u>>>,
         public array_holder<20u, numeric_type<ST2toT2ResultType>> {
     static_assert(getSpaceDimension<ST2toT2ResultType>() == 2u);
@@ -206,15 +202,12 @@ namespace tfel::math {
     TFEL_HOST_DEVICE constexpr auto getRunTimeProperties() const noexcept {
       return RunTimeProperties();
     }
-  };  // end of struct
-      // Expr<ST2toT2ResultType,StensorProductRightDerivativeExpr<2u> >
+  };  // end of Expr<ST2toT2ResultType,StensorProductRightDerivativeExpr<2u>>
 
-  /*!
-   * Partial specialisation for 3D stensor
-   */
-  template <typename ST2toT2ResultType>
+  //! \brief partial specialisation in 3D
+  template <ST2toT2Concept ST2toT2ResultType>
   struct Expr<ST2toT2ResultType, StensorProductRightDerivativeExpr<3u>>
-      : public ST2toT2Concept<
+      : public ST2toT2ConceptBase<
             Expr<ST2toT2ResultType, StensorProductRightDerivativeExpr<3u>>>,
         public array_holder<54u, numeric_type<ST2toT2ResultType>> {
     static_assert(getSpaceDimension<ST2toT2ResultType>() == 3u);
@@ -435,8 +428,7 @@ namespace tfel::math {
     TFEL_HOST_DEVICE constexpr auto getRunTimeProperties() const noexcept {
       return RunTimeProperties();
     }
-  };  // end of struct
-      // Expr<ST2toT2ResultType,StensorProductRightDerivativeExpr<1u> >
+  };  // end of Expr<ST2toT2ResultType,StensorProductRightDerivativeExpr<3u>>
 
 }  // end of namespace tfel::math
 
diff --git a/include/TFEL/Math/ST2toT2/T2toST2ST2toT2ProductExpr.hxx b/include/TFEL/Math/ST2toT2/T2toST2ST2toT2ProductExpr.hxx
index b4a89626e..0297685f0 100644
--- a/include/TFEL/Math/ST2toT2/T2toST2ST2toT2ProductExpr.hxx
+++ b/include/TFEL/Math/ST2toT2/T2toST2ST2toT2ProductExpr.hxx
@@ -29,10 +29,8 @@ namespace tfel::math {
   struct TFEL_VISIBILITY_LOCAL T2toST2ST2toT2ProductExpr {
   };  // end of struct T2toST2ST2toT2ProductExpr
 
-  /*!
-   * Partial specialisation
-   */
-  template <typename ST2toST2ResultType>
+  //! \brief partial specialisation in 1D
+  template <ST2toST2Concept ST2toST2ResultType>
   struct TFEL_VISIBILITY_LOCAL
       Expr<ST2toST2ResultType, T2toST2ST2toT2ProductExpr<1u>>
       : public ST2toST2ConceptBase<
@@ -52,12 +50,9 @@ namespace tfel::math {
      * \param[in] a : first term of the product
      * \param[in] b : second term of the product
      */
-    template <typename T2toST2Type, typename ST2toT2Type2>
-    TFEL_MATH_INLINE Expr(const T2toST2Type& a, const ST2toT2Type2& b) {
-      //! a simple check
-      static_assert((implementsT2toST2Concept<T2toST2Type>()));
-      //! a simple check
-      static_assert((implementsST2toT2Concept<ST2toT2Type2>()));
+    template <T2toST2Concept T2toST2Type, ST2toT2Concept ST2toT2Type2>
+    TFEL_HOST_DEVICE constexpr Expr(const T2toST2Type& a,
+                                    const ST2toT2Type2& b) noexcept {
       //! a simple check
       static_assert((getSpaceDimension<T2toST2Type>() == 1u));
       //! a simple check
@@ -77,22 +72,18 @@ namespace tfel::math {
      * \param[in] i : line index
      * \param[in] j : column index
      */
-    TFEL_MATH_INLINE const value_type& operator()(
-        const unsigned short i, const unsigned short j) const {
+    TFEL_HOST_DEVICE constexpr const value_type& operator()(
+        const unsigned short i, const unsigned short j) const noexcept {
       return this->v[i * 3 + j];
     }  // end of operator()
-    /*!
-     * \return the runtime properties
-     */
-    TFEL_MATH_INLINE RunTimeProperties getRunTimeProperties() const {
+    //! \return the runtime properties
+    TFEL_HOST_DEVICE constexpr auto getRunTimeProperties() const noexcept {
       return RunTimeProperties();
     }
   };  // end of struct Expr<ST2toST2ResultType,T2toST2ST2toT2ProductExpr>
 
-  /*!
-   * Partial specialisation
-   */
-  template <typename ST2toST2ResultType>
+  //! \brief partial specialisation in 2D
+  template <ST2toST2Concept ST2toST2ResultType>
   struct TFEL_VISIBILITY_LOCAL
       Expr<ST2toST2ResultType, T2toST2ST2toT2ProductExpr<2u>>
       : public ST2toST2ConceptBase<
@@ -112,12 +103,9 @@ namespace tfel::math {
      * \param[in] a : first term of the product
      * \param[in] b : second term of the product
      */
-    template <typename T2toST2Type, typename ST2toT2Type2>
-    TFEL_MATH_INLINE Expr(const T2toST2Type& a, const ST2toT2Type2& b) {
-      //! a simple check
-      static_assert((implementsT2toST2Concept<T2toST2Type>()));
-      //! a simple check
-      static_assert((implementsST2toT2Concept<ST2toT2Type2>()));
+    template <T2toST2Concept T2toST2Type, ST2toT2Concept ST2toT2Type2>
+    TFEL_HOST_DEVICE constexpr Expr(const T2toST2Type& a,
+                                    const ST2toT2Type2& b) noexcept {
       //! a simple check
       static_assert((getSpaceDimension<T2toST2Type>() == 2u));
       //! a simple check
@@ -160,22 +148,18 @@ namespace tfel::math {
      * \param[in] i : line index
      * \param[in] j : column index
      */
-    TFEL_MATH_INLINE const value_type& operator()(
-        const unsigned short i, const unsigned short j) const {
+    TFEL_HOST_DEVICE constexpr const value_type& operator()(
+        const unsigned short i, const unsigned short j) const noexcept {
       return this->v[i * 4 + j];
     }  // end of operator()
-    /*!
-     * \return the runtime properties
-     */
-    TFEL_MATH_INLINE RunTimeProperties getRunTimeProperties() const {
+    //! \return the runtime properties
+    TFEL_HOST_DEVICE constexpr auto getRunTimeProperties() const noexcept {
       return RunTimeProperties();
     }
   };  // end of struct Expr<ST2toST2ResultType,T2toST2ST2toT2ProductExpr>
 
-  /*!
-   * Partial specialisation
-   */
-  template <typename ST2toST2ResultType>
+  //! \brief partial specialisation in 3D
+  template <ST2toST2Concept ST2toST2ResultType>
   struct TFEL_VISIBILITY_LOCAL
       Expr<ST2toST2ResultType, T2toST2ST2toT2ProductExpr<3u>>
       : public ST2toST2ConceptBase<
@@ -195,12 +179,9 @@ namespace tfel::math {
      * \param[in] a : first term of the product
      * \param[in] b : second term of the product
      */
-    template <typename T2toST2Type, typename ST2toT2Type2>
-    TFEL_MATH_INLINE Expr(const T2toST2Type& a, const ST2toT2Type2& b) {
-      //! a simple check
-      static_assert((implementsT2toST2Concept<T2toST2Type>()));
-      //! a simple check
-      static_assert((implementsST2toT2Concept<ST2toT2Type2>()));
+    template <T2toST2Concept T2toST2Type, ST2toT2Concept ST2toT2Type2>
+    TFEL_HOST_DEVICE constexpr Expr(const T2toST2Type& a,
+                                    const ST2toT2Type2& b) noexcept {
       //! a simple check
       static_assert((getSpaceDimension<T2toST2Type>() == 3u));
       //! a simple check
@@ -319,14 +300,12 @@ namespace tfel::math {
      * \param[in] i : line index
      * \param[in] j : column index
      */
-    TFEL_MATH_INLINE const value_type& operator()(
-        const unsigned short i, const unsigned short j) const {
+    TFEL_HOST_DEVICE constexpr const value_type& operator()(
+        const unsigned short i, const unsigned short j) const noexcept {
       return this->v[i * 6 + j];
     }  // end of operator()
-    /*!
-     * \return the runtime properties
-     */
-    TFEL_MATH_INLINE RunTimeProperties getRunTimeProperties() const {
+    //! \return the runtime properties
+    TFEL_HOST_DEVICE constexpr auto getRunTimeProperties() const noexcept {
       return RunTimeProperties();
     }
   };  // end of struct Expr<ST2toST2ResultType,T2toST2ST2toT2ProductExpr>
diff --git a/include/TFEL/Math/ST2toT2/T2toT2ST2toT2ProductExpr.hxx b/include/TFEL/Math/ST2toT2/T2toT2ST2toT2ProductExpr.hxx
index b2a506ac4..3b9b9885a 100644
--- a/include/TFEL/Math/ST2toT2/T2toT2ST2toT2ProductExpr.hxx
+++ b/include/TFEL/Math/ST2toT2/T2toT2ST2toT2ProductExpr.hxx
@@ -21,39 +21,38 @@
 
 namespace tfel::math {
 
-  //! Empty structure used for partial specialisation of the
-  //! Expr class
+  /*!
+   * \brief empty structure used for partial specialisation of the
+   * Expr class
+   */
   template <unsigned short N>
   struct TFEL_VISIBILITY_LOCAL T2toT2ST2toT2ProductExpr {
   };  // end of struct T2toT2ST2toT2ProductExpr
 
-  /*!
-   * Partial specialisation
-   */
-  template <typename ST2toT2ResultType>
+  //! \brief partial specialisation in 1D
+  template <ST2toT2Concept ST2toT2ResultType>
   struct TFEL_VISIBILITY_LOCAL
       Expr<ST2toT2ResultType, T2toT2ST2toT2ProductExpr<1u>>
-      : public ST2toT2Concept<
+      : public ST2toT2ConceptBase<
             Expr<ST2toT2ResultType, T2toT2ST2toT2ProductExpr<1u>>>,
         public array_holder<
             TensorDimeToSize<getSpaceDimension<ST2toT2ResultType>()>::value *
                 StensorDimeToSize<
                     getSpaceDimension<ST2toT2ResultType>()>::value,
             numeric_type<ST2toT2ResultType>> {
+    //! a simple check
     static_assert(getSpaceDimension<ST2toT2ResultType>() == 1u);
     //! a simple alias
     typedef EmptyRunTimeProperties RunTimeProperties;
     //! a simple alias
     typedef numeric_type<ST2toT2ResultType> value_type;
-    //! a simple check
     /*!
      * \param[in] a : first term of the product
      * \param[in] b : second term of the product
      */
-    template <T2toT2Concept T2toT2Type, typename ST2toT2Type2>
+    template <T2toT2Concept T2toT2Type, ST2toT2Concept ST2toT2Type2>
     TFEL_HOST_DEVICE constexpr Expr(const T2toT2Type& a,
                                     const ST2toT2Type2& b) noexcept {
-      static_assert(implementsST2toT2Concept<ST2toT2Type2>());
       static_assert(getSpaceDimension<T2toT2Type>() == 1u);
       static_assert(getSpaceDimension<ST2toT2Type2>() == 1u);
       this->v[0] = a(0, 0) * b(0, 0) + a(0, 1) * b(1, 0) + a(0, 2) * b(2, 0);
@@ -81,13 +80,11 @@ namespace tfel::math {
     }
   };  // end of struct Expr<ST2toT2ResultType,T2toT2ST2toT2ProductExpr>
 
-  /*!
-   * Partial specialisation
-   */
-  template <typename ST2toT2ResultType>
+  //! \brief partial specialisation in 2D
+  template <ST2toT2Concept ST2toT2ResultType>
   struct TFEL_VISIBILITY_LOCAL
       Expr<ST2toT2ResultType, T2toT2ST2toT2ProductExpr<2u>>
-      : public ST2toT2Concept<
+      : public ST2toT2ConceptBase<
             Expr<ST2toT2ResultType, T2toT2ST2toT2ProductExpr<2u>>>,
         public array_holder<
             TensorDimeToSize<getSpaceDimension<ST2toT2ResultType>()>::value *
@@ -103,10 +100,9 @@ namespace tfel::math {
      * \param[in] a : first term of the product
      * \param[in] b : second term of the product
      */
-    template <T2toT2Concept T2toT2Type, typename ST2toT2Type2>
+    template <T2toT2Concept T2toT2Type, ST2toT2Concept ST2toT2Type2>
     TFEL_HOST_DEVICE constexpr Expr(const T2toT2Type& a,
                                     const ST2toT2Type2& b) noexcept {
-      static_assert(implementsST2toT2Concept<ST2toT2Type2>());
       static_assert(getSpaceDimension<T2toT2Type>() == 2u);
       static_assert(getSpaceDimension<ST2toT2Type2>() == 2u);
       this->v[0] = a(0, 0) * b(0, 0) + a(0, 1) * b(1, 0) + a(0, 2) * b(2, 0) +
@@ -165,13 +161,11 @@ namespace tfel::math {
     }
   };  // end of struct Expr<ST2toT2ResultType,T2toT2ST2toT2ProductExpr>
 
-  /*!
-   * Partial specialisation
-   */
-  template <typename ST2toT2ResultType>
+  //! \brief partial specialisation in 3D
+  template <ST2toT2Concept ST2toT2ResultType>
   struct TFEL_VISIBILITY_LOCAL
       Expr<ST2toT2ResultType, T2toT2ST2toT2ProductExpr<3u>>
-      : public ST2toT2Concept<
+      : public ST2toT2ConceptBase<
             Expr<ST2toT2ResultType, T2toT2ST2toT2ProductExpr<3u>>>,
         public array_holder<
             TensorDimeToSize<getSpaceDimension<ST2toT2ResultType>()>::value *
@@ -188,10 +182,9 @@ namespace tfel::math {
      * \param[in] a : first term of the product
      * \param[in] b : second term of the product
      */
-    template <T2toT2Concept T2toT2Type, typename ST2toT2Type2>
+    template <T2toT2Concept T2toT2Type, ST2toT2Concept ST2toT2Type2>
     TFEL_HOST_DEVICE constexpr Expr(const T2toT2Type& a,
                                     const ST2toT2Type2& b) noexcept {
-      static_assert(implementsST2toT2Concept<ST2toT2Type2>());
       static_assert(getSpaceDimension<T2toT2Type>() == 3u);
       static_assert(getSpaceDimension<ST2toT2Type2>() == 3u);
       this->v[0] = a(0, 0) * b(0, 0) + a(0, 1) * b(1, 0) + a(0, 2) * b(2, 0) +
diff --git a/include/TFEL/Math/ST2toT2/TensorST2toT2ProductExpr.hxx b/include/TFEL/Math/ST2toT2/TensorST2toT2ProductExpr.hxx
index f902ec5e3..2135670ef 100644
--- a/include/TFEL/Math/ST2toT2/TensorST2toT2ProductExpr.hxx
+++ b/include/TFEL/Math/ST2toT2/TensorST2toT2ProductExpr.hxx
@@ -22,14 +22,14 @@
 
 namespace tfel::math {
 
-  //! Empty structure used for partial specialisation of the
-  //! Expr class
+  /*!
+   * \brief  empty structure used for partial specialisation of the
+   * Expr class
+   */
   template <unsigned short N>
   struct TensorST2toT2ProductExpr {};  // end of struct TensorST2toT2ProductExpr
 
-  /*!
-   * Partial specialisation
-   */
+  //! \brief partial specialisation in 1D
   template <StensorConcept StensorResultType>
   struct Expr<StensorResultType, TensorST2toT2ProductExpr<1u>>
       : public StensorConceptBase<
@@ -46,10 +46,9 @@ namespace tfel::math {
      * \param[in] a : first term of the product
      * \param[in] b : second term of the product
      */
-    template <TensorConcept TensorType, typename ST2toT2Type>
+    template <TensorConcept TensorType, ST2toT2Concept ST2toT2Type>
     TFEL_HOST_DEVICE constexpr Expr(const TensorType& a,
                                     const ST2toT2Type& b) noexcept {
-      static_assert(implementsST2toT2Concept<ST2toT2Type>());
       static_assert(getSpaceDimension<ST2toT2Type>() == 1u);
       static_assert(getSpaceDimension<TensorType>() == 1u);
       this->v[0] = b(0, 0) * a(0) + b(1, 0) * a(1) + b(2, 0) * a(2);
@@ -72,18 +71,13 @@ namespace tfel::math {
         const unsigned short i) const noexcept {
       return this->v[i];
     }  // end of operator()
-    /*!
-     * \return the runtime properties
-     */
-    TFEL_HOST_DEVICE constexpr RunTimeProperties getRunTimeProperties() const
-        noexcept {
+    //! \return the runtime properties
+    TFEL_HOST_DEVICE constexpr auto getRunTimeProperties() const noexcept {
       return RunTimeProperties();
     }
   };  // end of struct Expr<StensorResultType,T2STensorST2toT2ProductExp>
 
-  /*!
-   * Partial specialisation
-   */
+  //! \brief partial specialisation in 2D
   template <StensorConcept StensorResultType>
   struct Expr<StensorResultType, TensorST2toT2ProductExpr<2u>>
       : public StensorConceptBase<
@@ -100,10 +94,9 @@ namespace tfel::math {
      * \param[in] a : first term of the product
      * \param[in] b : second term of the product
      */
-    template <TensorConcept TensorType, typename ST2toT2Type>
+    template <TensorConcept TensorType, ST2toT2Concept ST2toT2Type>
     TFEL_HOST_DEVICE constexpr Expr(const TensorType& a,
                                     const ST2toT2Type& b) noexcept {
-      static_assert(implementsST2toT2Concept<ST2toT2Type>());
       static_assert(getSpaceDimension<ST2toT2Type>() == 2u);
       static_assert(getSpaceDimension<TensorType>() == 2u);
       this->v[0] = b(4, 0) * a(4) + b(3, 0) * a(3) + b(2, 0) * a(2) +
@@ -131,18 +124,13 @@ namespace tfel::math {
         const unsigned short i) const noexcept {
       return this->v[i];
     }  // end of operator()
-    /*!
-     * \return the runtime properties
-     */
-    TFEL_HOST_DEVICE constexpr RunTimeProperties getRunTimeProperties() const
-        noexcept {
+    //! \return the runtime properties
+    TFEL_HOST_DEVICE constexpr auto getRunTimeProperties() const noexcept {
       return RunTimeProperties();
     }
-  };  // end of struct Expr<StensorResultType,T2STensorST2toT2ProductExp>
+  };  // end of Expr<StensorResultType,T2STensorST2toT2ProductExp>
 
-  /*!
-   * Partial specialisation
-   */
+  //! \brief partial specialisation in 2D
   template <StensorConcept StensorResultType>
   struct Expr<StensorResultType, TensorST2toT2ProductExpr<3u>>
       : public StensorConceptBase<
@@ -159,10 +147,9 @@ namespace tfel::math {
      * \param[in] a : first term of the product
      * \param[in] b : second term of the product
      */
-    template <TensorConcept TensorType, typename ST2toT2Type>
+    template <TensorConcept TensorType, ST2toT2Concept ST2toT2Type>
     TFEL_HOST_DEVICE constexpr Expr(const TensorType& a,
                                     const ST2toT2Type& b) noexcept {
-      static_assert(implementsST2toT2Concept<ST2toT2Type>());
       static_assert(getSpaceDimension<ST2toT2Type>() == 3u);
       static_assert(getSpaceDimension<TensorType>() == 3u);
       this->v[0] = b(8, 0) * a(8) + b(7, 0) * a(7) + b(6, 0) * a(6) +
@@ -200,13 +187,11 @@ namespace tfel::math {
         const unsigned short i) const noexcept {
       return this->v[i];
     }  // end of operator()
-    /*!
-     * \return the runtime properties
-     */
+    //! \return the runtime properties
     TFEL_HOST_DEVICE constexpr auto getRunTimeProperties() const noexcept {
       return RunTimeProperties();
     }
-  };  // end of struct Expr<StensorResultType,T2STensorST2toT2ProductExp>
+  };  // end of Expr<StensorResultType,T2STensorST2toT2ProductExp>
 
 }  // end of namespace tfel::math
 
diff --git a/include/TFEL/Math/ST2toT2/st2tot2.ixx b/include/TFEL/Math/ST2toT2/st2tot2.ixx
index 245064d53..052723f11 100644
--- a/include/TFEL/Math/ST2toT2/st2tot2.ixx
+++ b/include/TFEL/Math/ST2toT2/st2tot2.ixx
@@ -18,7 +18,6 @@
 #include <iterator>
 #include <algorithm>
 #include <type_traits>
-
 #include "TFEL/FSAlgorithm/FSAlgorithm.hxx"
 #include "TFEL/Math/ST2toT2/StensorProductLeftDerivativeExpr.hxx"
 #include "TFEL/Math/ST2toT2/StensorProductRightDerivativeExpr.hxx"
@@ -27,52 +26,58 @@ namespace tfel::math {
 
   template <unsigned short N, typename T>
   template <StensorConcept StensorType>
-  std::enable_if_t<getSpaceDimension<StensorType>() == N &&
-                       isAssignableTo<numeric_type<StensorType>, T>(),
-                   Expr<st2tot2<N, T>, StensorProductLeftDerivativeExpr<N>>>
-  st2tot2<N, T>::tpld(const StensorType& b) {
+  TFEL_HOST_DEVICE constexpr auto
+  st2tot2<N, T>::tpld(const StensorType& b) noexcept requires(
+      getSpaceDimension<StensorType>() == N &&
+      isAssignableTo<numeric_type<StensorType>, T>()) {
     return Expr<st2tot2<N, T>, StensorProductLeftDerivativeExpr<N>>(b);
   }  // end of st2tot2<N,T>
 
   template <unsigned short N, typename ValueType>
   template <StensorConcept StensorType, ST2toST2Concept ST2toST2Type>
-  std::enable_if_t<
-      getSpaceDimension<StensorType>() == N &&
+  TFEL_HOST_DEVICE constexpr auto st2tot2<N, ValueType>::
+      tpld(const StensorType& b, const ST2toST2Type& C) noexcept requires(
+          getSpaceDimension<StensorType>() == N &&
           getSpaceDimension<ST2toST2Type>() == N &&
           isAssignableTo<BinaryOperationResult<numeric_type<StensorType>,
                                                numeric_type<ST2toST2Type>,
                                                OpMult>,
-                         ValueType>(),
-      Expr<st2tot2<N, ValueType>, StensorProductLeftDerivativeExpr<N>>>
-  st2tot2<N, ValueType>::tpld(const StensorType& b, const ST2toST2Type& C) {
+                         ValueType>()) {
     return Expr<st2tot2<N, ValueType>, StensorProductLeftDerivativeExpr<N>>(b,
                                                                             C);
   }
 
   template <unsigned short N, typename T>
   template <StensorConcept StensorType>
-  TFEL_HOST_DEVICE constexpr std::enable_if_t<
+  TFEL_HOST_DEVICE constexpr auto
+  st2tot2<N, T>::tprd(const StensorType& a) noexcept requires(
       getSpaceDimension<StensorType>() == N &&
-          isAssignableTo<numeric_type<StensorType>, T>(),
-      Expr<st2tot2<N, T>, StensorProductRightDerivativeExpr<N>>>
-  st2tot2<N, T>::tprd(const StensorType& a) noexcept {
+      isAssignableTo<numeric_type<StensorType>, T>()) {
     return Expr<st2tot2<N, T>, StensorProductRightDerivativeExpr<N>>(a);
   }
 
   template <unsigned short N, typename T>
   template <StensorConcept StensorType, ST2toST2Concept ST2toST2Type>
-  TFEL_HOST_DEVICE constexpr std::enable_if_t<
-      getSpaceDimension<StensorType>() == N &&
+  TFEL_HOST_DEVICE constexpr auto st2tot2<N, T>::
+      tprd(const StensorType& a, const ST2toST2Type& C) noexcept requires(
+          getSpaceDimension<StensorType>() == N &&
           getSpaceDimension<ST2toST2Type>() == N &&
           isAssignableTo<BinaryOperationResult<numeric_type<StensorType>,
                                                numeric_type<ST2toST2Type>,
                                                OpMult>,
-                         T>(),
-      Expr<st2tot2<N, T>, StensorProductRightDerivativeExpr<N>>>
-  st2tot2<N, T>::tprd(const StensorType& a, const ST2toST2Type& C) noexcept {
+                         T>()) {
     return Expr<st2tot2<N, T>, StensorProductRightDerivativeExpr<N>>(a, C);
   }
 
+  template <unsigned short N, typename T>
+  TFEL_HOST_DEVICE constexpr void st2tot2<N, T>::import(
+      const base_type<T>* const src) noexcept {
+    tfel::fsalgo::transform<
+        TensorDimeToSize<N>::value *
+        StensorDimeToSize<N>::value>::exe(src, this->begin(),
+                                          [](const auto& v) { return T(v); });
+  }  // end of import
+
   template <unsigned short N, typename T>
   TFEL_HOST_DEVICE constexpr void st2tot2<N, T>::copy(const auto p) noexcept {
     tfel::fsalgo::copy<TensorDimeToSize<N>::value *
diff --git a/include/TFEL/Math/Stensor/stensor.ixx b/include/TFEL/Math/Stensor/stensor.ixx
index 8d4a1020e..dda06605b 100644
--- a/include/TFEL/Math/Stensor/stensor.ixx
+++ b/include/TFEL/Math/Stensor/stensor.ixx
@@ -885,10 +885,9 @@ namespace tfel::math {
   }
 
   template <StensorConcept StensorType>
-  TFEL_HOST_DEVICE stensor<getSpaceDimension<StensorType>(),
-                           numeric_type<StensorType>>
-  change_basis(const StensorType& s,
-               const tmatrix<3u, 3u, base_type<numeric_type<StensorType>>>& r) {
+  TFEL_HOST_DEVICE constexpr auto change_basis(
+      const StensorType& s,
+      const tmatrix<3u, 3u, base_type<numeric_type<StensorType>>>& r) noexcept {
     stensor<getSpaceDimension<StensorType>(), numeric_type<StensorType>> ns{s};
     ns.changeBasis(r);
     return ns;
diff --git a/include/TFEL/Math/T2toST2/LeftCauchyGreenTensorDerivativeExpr.hxx b/include/TFEL/Math/T2toST2/LeftCauchyGreenTensorDerivativeExpr.hxx
index a4659a9de..565ddab85 100644
--- a/include/TFEL/Math/T2toST2/LeftCauchyGreenTensorDerivativeExpr.hxx
+++ b/include/TFEL/Math/T2toST2/LeftCauchyGreenTensorDerivativeExpr.hxx
@@ -34,7 +34,7 @@ namespace tfel::math {
    */
   template <typename T2toST2Type>
   struct Expr<T2toST2Type, LeftCauchyGreenTensorDerivativeExpr<1u>>
-      : public T2toST2Concept<
+      : public T2toST2ConceptBase<
             Expr<T2toST2Type, LeftCauchyGreenTensorDerivativeExpr<1u>>>,
         public array_holder<9u, numeric_type<T2toST2Type>> {
     static_assert(getSpaceDimension<T2toST2Type>() == 1u);
@@ -75,7 +75,7 @@ namespace tfel::math {
    */
   template <typename T2toST2Type>
   struct Expr<T2toST2Type, LeftCauchyGreenTensorDerivativeExpr<2u>>
-      : public T2toST2Concept<
+      : public T2toST2ConceptBase<
             Expr<T2toST2Type, LeftCauchyGreenTensorDerivativeExpr<2u>>>,
         public array_holder<20u, numeric_type<T2toST2Type>> {
     static_assert(getSpaceDimension<T2toST2Type>() == 2u);
@@ -128,7 +128,7 @@ namespace tfel::math {
    */
   template <typename T2toST2Type>
   struct Expr<T2toST2Type, LeftCauchyGreenTensorDerivativeExpr<3u>>
-      : public T2toST2Concept<
+      : public T2toST2ConceptBase<
             Expr<T2toST2Type, LeftCauchyGreenTensorDerivativeExpr<3u>>>,
         public array_holder<54u, numeric_type<T2toST2Type>> {
     static_assert(getSpaceDimension<T2toST2Type>() == 3u);
diff --git a/include/TFEL/Math/T2toST2/RightCauchyGreenTensorDerivativeExpr.hxx b/include/TFEL/Math/T2toST2/RightCauchyGreenTensorDerivativeExpr.hxx
index 6acd379cb..2fbccabe9 100644
--- a/include/TFEL/Math/T2toST2/RightCauchyGreenTensorDerivativeExpr.hxx
+++ b/include/TFEL/Math/T2toST2/RightCauchyGreenTensorDerivativeExpr.hxx
@@ -34,7 +34,7 @@ namespace tfel::math {
    */
   template <typename T2toST2Type>
   struct Expr<T2toST2Type, RightCauchyGreenTensorDerivativeExpr<1u>>
-      : public T2toST2Concept<
+      : public T2toST2ConceptBase<
             Expr<T2toST2Type, RightCauchyGreenTensorDerivativeExpr<1u>>>,
         public array_holder<9u, numeric_type<T2toST2Type>> {
     static_assert(getSpaceDimension<T2toST2Type>() == 1u);
@@ -81,7 +81,7 @@ namespace tfel::math {
    */
   template <typename T2toST2Type>
   struct Expr<T2toST2Type, RightCauchyGreenTensorDerivativeExpr<2u>>
-      : public T2toST2Concept<
+      : public T2toST2ConceptBase<
             Expr<T2toST2Type, RightCauchyGreenTensorDerivativeExpr<2u>>>,
         public array_holder<20u, numeric_type<T2toST2Type>> {
     static_assert(getSpaceDimension<T2toST2Type>() == 2u);
@@ -139,7 +139,7 @@ namespace tfel::math {
    */
   template <typename T2toST2Type>
   struct Expr<T2toST2Type, RightCauchyGreenTensorDerivativeExpr<3u>>
-      : public T2toST2Concept<
+      : public T2toST2ConceptBase<
             Expr<T2toST2Type, RightCauchyGreenTensorDerivativeExpr<3u>>>,
         public array_holder<54u, numeric_type<T2toST2Type>> {
     static_assert(getSpaceDimension<T2toST2Type>() == 3u);
diff --git a/include/TFEL/Math/T2toST2/ST2toST2T2toST2ProductExpr.hxx b/include/TFEL/Math/T2toST2/ST2toST2T2toST2ProductExpr.hxx
index e809b5fc0..d23e8a627 100644
--- a/include/TFEL/Math/T2toST2/ST2toST2T2toST2ProductExpr.hxx
+++ b/include/TFEL/Math/T2toST2/ST2toST2T2toST2ProductExpr.hxx
@@ -15,7 +15,6 @@
 #define LIB_TFEL_MATH_ST2TOST2T2TOST2PRODUCTEXPR_HXX
 
 #include "TFEL/Config/TFELConfig.hxx"
-
 #include "TFEL/Math/General/EmptyRunTimeProperties.hxx"
 #include "TFEL/Math/ST2toST2/ST2toST2Concept.hxx"
 #include "TFEL/Math/T2toST2/T2toST2Concept.hxx"
@@ -31,10 +30,10 @@ namespace tfel::math {
   };  // end of struct ST2toST2T2toST2ProductExpr
 
   //! \brief partial specialisation in 1D
-  template <typename T2toST2ResultType>
+  template <T2toST2Concept T2toST2ResultType>
   struct TFEL_VISIBILITY_LOCAL
       Expr<T2toST2ResultType, ST2toST2T2toST2ProductExpr<1u>>
-      : public T2toST2Concept<
+      : public T2toST2ConceptBase<
             Expr<T2toST2ResultType, ST2toST2T2toST2ProductExpr<1u>>>,
         public array_holder<
             StensorDimeToSize<getSpaceDimension<T2toST2ResultType>()>::value *
@@ -49,10 +48,9 @@ namespace tfel::math {
      * \param[in] a : first term of the product
      * \param[in] b : second term of the product
      */
-    template <ST2toST2Concept ST2toST2Type, typename T2toST2Type2>
+    template <ST2toST2Concept ST2toST2Type, T2toST2Concept T2toST2Type2>
     TFEL_HOST_DEVICE constexpr Expr(const ST2toST2Type& a,
                                     const T2toST2Type2& b) noexcept {
-      static_assert(implementsT2toST2Concept<T2toST2Type2>());
       static_assert(getSpaceDimension<ST2toST2Type>() == 1u);
       static_assert(getSpaceDimension<T2toST2Type2>() == 1u);
       this->v[0] = a(0, 0) * b(0, 0) + a(0, 1) * b(1, 0) + a(0, 2) * b(2, 0);
@@ -81,10 +79,10 @@ namespace tfel::math {
   };  // end of struct Expr<T2toST2ResultType,ST2toST2T2toST2ProductExpr>
 
   //! \brief partial specialisation in 2D
-  template <typename T2toST2ResultType>
+  template <T2toST2Concept T2toST2ResultType>
   struct TFEL_VISIBILITY_LOCAL
       Expr<T2toST2ResultType, ST2toST2T2toST2ProductExpr<2u>>
-      : public T2toST2Concept<
+      : public T2toST2ConceptBase<
             Expr<T2toST2ResultType, ST2toST2T2toST2ProductExpr<2u>>>,
         public array_holder<
             StensorDimeToSize<getSpaceDimension<T2toST2ResultType>()>::value *
@@ -99,10 +97,9 @@ namespace tfel::math {
      * \param[in] a : first term of the product
      * \param[in] b : second term of the product
      */
-    template <ST2toST2Concept ST2toST2Type, typename T2toST2Type2>
+    template <ST2toST2Concept ST2toST2Type, T2toST2Concept T2toST2Type2>
     TFEL_HOST_DEVICE constexpr Expr(const ST2toST2Type& a,
                                     const T2toST2Type2& b) noexcept {
-      static_assert(implementsT2toST2Concept<T2toST2Type2>());
       static_assert(getSpaceDimension<ST2toST2Type>() == 2u);
       static_assert(getSpaceDimension<T2toST2Type2>() == 2u);
       this->v[0] = a(0, 0) * b(0, 0) + a(0, 1) * b(1, 0) + a(0, 2) * b(2, 0) +
@@ -162,10 +159,10 @@ namespace tfel::math {
   };  // end of struct Expr<T2toST2ResultType,ST2toST2T2toST2ProductExpr>
 
   //! \brief partial specialisation in 3D
-  template <typename T2toST2ResultType>
+  template <T2toST2Concept T2toST2ResultType>
   struct TFEL_VISIBILITY_LOCAL
       Expr<T2toST2ResultType, ST2toST2T2toST2ProductExpr<3u>>
-      : public T2toST2Concept<
+      : public T2toST2ConceptBase<
             Expr<T2toST2ResultType, ST2toST2T2toST2ProductExpr<3u>>>,
         public array_holder<
             StensorDimeToSize<getSpaceDimension<T2toST2ResultType>()>::value *
@@ -180,10 +177,9 @@ namespace tfel::math {
      * \param[in] a : first term of the product
      * \param[in] b : second term of the product
      */
-    template <ST2toST2Concept ST2toST2Type, typename T2toST2Type2>
+    template <ST2toST2Concept ST2toST2Type, T2toST2Concept T2toST2Type2>
     TFEL_HOST_DEVICE constexpr Expr(const ST2toST2Type& a,
                                     const T2toST2Type2& b) noexcept {
-      static_assert(implementsT2toST2Concept<T2toST2Type2>());
       static_assert(getSpaceDimension<ST2toST2Type>() == 3u);
       static_assert(getSpaceDimension<T2toST2Type2>() == 3u);
       this->v[0] = a(0, 0) * b(0, 0) + a(0, 1) * b(1, 0) + a(0, 2) * b(2, 0) +
diff --git a/include/TFEL/Math/T2toST2/StensorT2toST2ProductExpr.hxx b/include/TFEL/Math/T2toST2/StensorT2toST2ProductExpr.hxx
index 866a1cdc6..d2cd93753 100644
--- a/include/TFEL/Math/T2toST2/StensorT2toST2ProductExpr.hxx
+++ b/include/TFEL/Math/T2toST2/StensorT2toST2ProductExpr.hxx
@@ -21,15 +21,15 @@
 
 namespace tfel::math {
 
-  //! Empty structure used for partial specialisation of the
-  //! Expr class
+  /*!
+   * \brief empty structure used for partial specialisation of the
+   * Expr class
+   */
   template <unsigned short N>
   struct StensorT2toST2ProductExpr {
   };  // end of struct StensorT2toST2ProductExpr
 
-  /*!
-   * Partial specialisation
-   */
+  //! \brief partial specialisation in 1D
   template <TensorConcept TensorResultType>
   struct Expr<TensorResultType, StensorT2toST2ProductExpr<1u>>
       : public TensorConceptBase<
@@ -46,9 +46,9 @@ namespace tfel::math {
      * \param[in] a : first term of the product
      * \param[in] b : second term of the product
      */
-    template <StensorConcept StensorType, typename T2toST2Type>
-    TFEL_MATH_INLINE Expr(const StensorType& a, const T2toST2Type& b) {
-      static_assert(implementsT2toST2Concept<T2toST2Type>());
+    template <StensorConcept StensorType, T2toST2Concept T2toST2Type>
+    TFEL_HOST_DEVICE constexpr Expr(const StensorType& a,
+                                    const T2toST2Type& b) noexcept {
       static_assert(getSpaceDimension<T2toST2Type>() == 1u);
       static_assert(getSpaceDimension<StensorType>() == 1u);
       this->v[0] = b(0, 0) * a(0) + b(1, 0) * a(1) + b(2, 0) * a(2);
@@ -59,21 +59,19 @@ namespace tfel::math {
      * \brief access operator
      * \param[in] i : index
      */
-    TFEL_MATH_INLINE const value_type& operator()(
-        const unsigned short i) const {
+    TFEL_HOST_DEVICE constexpr const value_type& operator()(
+        const unsigned short i) const noexcept {
       return this->v[i];
     }  // end of operator()
     /*!
      * \return the runtime properties
      */
-    TFEL_MATH_INLINE RunTimeProperties getRunTimeProperties() const {
+    TFEL_HOST_DEVICE constexpr auto getRunTimeProperties() const noexcept {
       return RunTimeProperties();
     }
   };  // end of struct Expr<TensorResultType,T2SStensorT2toST2ProductExp>
 
-  /*!
-   * Partial specialisation
-   */
+  //! \brief partial specialisation in 2D
   template <TensorConcept TensorResultType>
   struct Expr<TensorResultType, StensorT2toST2ProductExpr<2u>>
       : public TensorConceptBase<
@@ -90,9 +88,9 @@ namespace tfel::math {
      * \param[in] a : first term of the product
      * \param[in] b : second term of the product
      */
-    template <StensorConcept StensorType, typename T2toST2Type>
-    TFEL_MATH_INLINE Expr(const StensorType& a, const T2toST2Type& b) {
-      static_assert(implementsT2toST2Concept<T2toST2Type>());
+    template <StensorConcept StensorType, T2toST2Concept T2toST2Type>
+    TFEL_HOST_DEVICE constexpr Expr(const StensorType& a,
+                                    const T2toST2Type& b) noexcept {
       static_assert(getSpaceDimension<T2toST2Type>() == 2u);
       static_assert(getSpaceDimension<StensorType>() == 2u);
       this->v[0] =
@@ -110,21 +108,19 @@ namespace tfel::math {
      * \brief access operator
      * \param[in] i : index
      */
-    TFEL_MATH_INLINE const value_type& operator()(
-        const unsigned short i) const {
+    TFEL_HOST_DEVICE constexpr const value_type& operator()(
+        const unsigned short i) const noexcept {
       return this->v[i];
     }  // end of operator()
     /*!
      * \return the runtime properties
      */
-    TFEL_MATH_INLINE RunTimeProperties getRunTimeProperties() const {
+    TFEL_HOST_DEVICE constexpr auto getRunTimeProperties() const noexcept {
       return RunTimeProperties();
     }
   };  // end of struct Expr<TensorResultType,T2SStensorT2toST2ProductExp>
 
-  /*!
-   * Partial specialisation
-   */
+  //! \brief partial specialisation in 3D
   template <TensorConcept TensorResultType>
   struct Expr<TensorResultType, StensorT2toST2ProductExpr<3u>>
       : public TensorConceptBase<
@@ -141,9 +137,9 @@ namespace tfel::math {
      * \param[in] a : first term of the product
      * \param[in] b : second term of the product
      */
-    template <StensorConcept StensorType, typename T2toST2Type>
-    TFEL_MATH_INLINE Expr(const StensorType& a, const T2toST2Type& b) {
-      static_assert(implementsT2toST2Concept<T2toST2Type>());
+    template <StensorConcept StensorType, T2toST2Concept T2toST2Type>
+    TFEL_HOST_DEVICE constexpr Expr(const StensorType& a,
+                                    const T2toST2Type& b) noexcept {
       static_assert(getSpaceDimension<T2toST2Type>() == 3u);
       static_assert(getSpaceDimension<StensorType>() == 3u);
       this->v[0] = b(5, 0) * a[5] + b(4, 0) * a[4] + b(3, 0) * a[3] +
@@ -169,14 +165,14 @@ namespace tfel::math {
      * \brief access operator
      * \param[in] i : index
      */
-    TFEL_MATH_INLINE const value_type& operator()(
-        const unsigned short i) const {
+    TFEL_HOST_DEVICE constexpr const value_type& operator()(
+        const unsigned short i) const noexcept {
       return this->v[i];
     }  // end of operator()
     /*!
      * \return the runtime properties
      */
-    TFEL_MATH_INLINE RunTimeProperties getRunTimeProperties() const {
+    TFEL_HOST_DEVICE constexpr auto getRunTimeProperties() const noexcept {
       return RunTimeProperties();
     }
   };  // end of struct Expr<TensorResultType,T2SStensorT2toST2ProductExp>
diff --git a/include/TFEL/Math/T2toST2/T2toST2Concept.hxx b/include/TFEL/Math/T2toST2/T2toST2Concept.hxx
index f83a1e4e4..022aee39e 100644
--- a/include/TFEL/Math/T2toST2/T2toST2Concept.hxx
+++ b/include/TFEL/Math/T2toST2/T2toST2Concept.hxx
@@ -32,41 +32,40 @@ namespace tfel::math {
    * \brief Helper class to characterise t2tost2.
    */
   struct T2toST2Tag {};
-
+  /*!
+   * \brief an helper class that simply exposes publically a member named
+   * ConceptTag as an alias to T2toST2Tag.
+   *
+   * The main reason for this alias is to properly implement the `ConceptRebind`
+   * metafunction.
+   */
   template <class T>
-  struct T2toST2Concept {
-    typedef T2toST2Tag ConceptTag;
-
-   protected:
-    T2toST2Concept() = default;
-    T2toST2Concept(T2toST2Concept&&) = default;
-    T2toST2Concept(const T2toST2Concept&) = default;
-    T2toST2Concept& operator=(const T2toST2Concept&) = default;
-    ~T2toST2Concept() = default;
+  struct T2toST2ConceptBase {
+    using ConceptTag = T2toST2Tag;
   };
-
   /*!
-   * \brief an helper function which returns if the given type implements the
-   * `T2toST2Concept`.
-   * \tparam T2toST2Type: type tested
+   * \brief definition of the T2toST2Concept
+   * a class matching the stensor concept must expose the `T2toST2Tag` and have
+   * access operators.
    */
-  template <typename T2toST2Type>
-  TFEL_HOST_DEVICE constexpr bool implementsT2toST2Concept() {
-    return tfel::meta::implements<T2toST2Type, T2toST2Concept>();
-  }  // end of implementsT2toST2Concept
-
-  //! paratial specialisation for T2toST2
+  template <typename T>
+  concept T2toST2Concept =
+      (std::is_same_v<typename std::decay_t<T>::ConceptTag, T2toST2Tag>)&&  //
+      (requires(const T t, const unsigned short i, const unsigned short j) {
+        t(i, j);
+      });
+  //! \brief partial specialisation for T2toST2
   template <typename Type>
   struct ConceptRebind<T2toST2Tag, Type> {
-    using type = T2toST2Concept<Type>;
+    using type = T2toST2ConceptBase<Type>;
   };
-
-  template <typename T2toST2Type>
-  typename std::enable_if<
-      implementsT2toST2Concept<T2toST2Type>(),
-      typename tfel::typetraits::AbsType<numeric_type<T2toST2Type>>::type>::type
-  abs(const T2toST2Type&);
-
+  /*!
+   * \return the sum of the absolute values of all components of an
+   * linear application transforming a symmetric tensor in an unsymmetric tensor
+   * \param[in] t: linear application transforming a symmetric tensor in an
+   * unsymmetric tensor
+   */
+  TFEL_HOST_DEVICE constexpr void abs(const T2toST2Concept auto&) noexcept;
   /*!
    * \brief compute de derivative of the push-forward of a symmetric
    * second order tensor with respect to the deformation gradient, i.e.
@@ -77,17 +76,16 @@ namespace tfel::math {
    * \param[in]  : orginal tensor
    * \param[in]  : deformation gradient
    */
-  template <typename T2toST2ResultType,
+  template <T2toST2Concept T2toST2ResultType,
             StensorConcept StensorType,
             TensorConcept TensorType>
-  typename std::enable_if<implementsT2toST2Concept<T2toST2ResultType>() &&
-                              tfel::typetraits::IsFundamentalNumericType<
-                                  numeric_type<TensorType>>::cond &&
-                              isAssignableTo<numeric_type<StensorType>,
-                                             numeric_type<T2toST2ResultType>>(),
-                          void>::type
+  TFEL_HOST_DEVICE constexpr void
   computePushForwardDerivativeWithRespectToDeformationGradient(
-      T2toST2ResultType&, const StensorType&, const TensorType&);
+      T2toST2ResultType&, const StensorType&, const TensorType&) noexcept  //
+      requires(tfel::typetraits::IsFundamentalNumericType<
+               numeric_type<TensorType>>::cond&&
+                   isAssignableTo<numeric_type<StensorType>,
+                                  numeric_type<T2toST2ResultType>>());
   /*!
    * \brief compute the Cauchy stress derivative from the Kirchhoff
    * stress derivative with respect to the deformation gradient
@@ -96,13 +94,17 @@ namespace tfel::math {
    * \param[in]  s    : Cauchy stress
    * \param[in]  F    : deformation gradient
    */
-  template <typename T2toST2ResultType,
-            typename T2toST2Type,
+  template <T2toST2Concept T2toST2ResultType,
+            T2toST2Concept T2toST2Type,
             StensorConcept StensorType,
             TensorConcept TensorType>
-  typename std::enable_if<
-      implementsT2toST2Concept<T2toST2ResultType>() &&
-          implementsT2toST2Concept<T2toST2Type>() &&
+  TFEL_HOST_DEVICE constexpr void
+  computeCauchyStressDerivativeFromKirchhoffStressDerivative(
+      T2toST2ResultType&,
+      const T2toST2Type&,
+      const StensorType&,
+      const TensorType&) noexcept  //
+      requires(
           getSpaceDimension<T2toST2ResultType>() ==
               getSpaceDimension<T2toST2Type>() &&
           getSpaceDimension<T2toST2ResultType>() ==
@@ -114,12 +116,7 @@ namespace tfel::math {
           isAssignableTo<typename ComputeBinaryResult<numeric_type<T2toST2Type>,
                                                       numeric_type<StensorType>,
                                                       OpPlus>::Result,
-                         numeric_type<T2toST2ResultType>>(),
-      void>::type
-  computeCauchyStressDerivativeFromKirchhoffStressDerivative(T2toST2ResultType&,
-                                                             const T2toST2Type&,
-                                                             const StensorType&,
-                                                             const TensorType&);
+                         numeric_type<T2toST2ResultType>>());
   /*!
    * \brief compute the Kirchhoff stress derivative from the Cauchy
    * stress derivative with respect to the deformation gradient
@@ -128,13 +125,17 @@ namespace tfel::math {
    * \param[in]  s    : Cauchy stress
    * \param[in]  F    : deformation gradient
    */
-  template <typename T2toST2ResultType,
-            typename T2toST2Type,
+  template <T2toST2Concept T2toST2ResultType,
+            T2toST2Concept T2toST2Type,
             StensorConcept StensorType,
             TensorConcept TensorType>
-  typename std::enable_if<
-      implementsT2toST2Concept<T2toST2ResultType>() &&
-          implementsT2toST2Concept<T2toST2Type>() &&
+  TFEL_HOST_DEVICE constexpr void
+  computeKirchhoffStressDerivativeFromCauchyStressDerivative(
+      T2toST2ResultType&,
+      const T2toST2Type&,
+      const StensorType&,
+      const TensorType&) noexcept  //
+      requires(
           getSpaceDimension<T2toST2ResultType>() ==
               getSpaceDimension<T2toST2Type>() &&
           getSpaceDimension<T2toST2ResultType>() ==
@@ -146,12 +147,16 @@ namespace tfel::math {
           isAssignableTo<typename ComputeBinaryResult<numeric_type<T2toST2Type>,
                                                       numeric_type<StensorType>,
                                                       OpPlus>::Result,
-                         numeric_type<T2toST2ResultType>>(),
-      void>::type
-  computeKirchhoffStressDerivativeFromCauchyStressDerivative(T2toST2ResultType&,
-                                                             const T2toST2Type&,
-                                                             const StensorType&,
-                                                             const TensorType&);
+                         numeric_type<T2toST2ResultType>>());
+  /*!
+   * \brief an helper function which returns if the given type implements the
+   * `T2toST2Concept`.
+   * \tparam T2toST2Type: type tested
+   */
+  template <T2toST2Concept T2toST2Type>
+  [[deprecated]] TFEL_HOST_DEVICE constexpr bool implementsT2toST2Concept() {
+    return T2toST2Concept<T2toST2Type>;
+  }  // end of implementsT2toST2Concept
 
 }  // end of namespace tfel::math
 
diff --git a/include/TFEL/Math/T2toST2/T2toST2Concept.ixx b/include/TFEL/Math/T2toST2/T2toST2Concept.ixx
index 8df41c86e..ff7094f24 100644
--- a/include/TFEL/Math/T2toST2/T2toST2Concept.ixx
+++ b/include/TFEL/Math/T2toST2/T2toST2Concept.ixx
@@ -20,37 +20,35 @@
 
 namespace tfel::math {
 
-  template <typename T2toST2Type>
-  typename std::enable_if<
-      implementsT2toST2Concept<T2toST2Type>(),
-      typename tfel::typetraits::AbsType<numeric_type<T2toST2Type>>::type>::type
-  abs(const T2toST2Type& v) {
+  TFEL_HOST_DEVICE constexpr auto abs(const T2toST2Concept auto& t) noexcept {
+    using T2toST2Type = decltype(t);
     using NumType = numeric_type<T2toST2Type>;
     using AbsNumType = typename tfel::typetraits::AbsType<NumType>::type;
-    AbsNumType a(0);
+    auto a = AbsNumType{};
     constexpr auto ssize =
         StensorDimeToSize<getSpaceDimension<T2toST2Type>()>::value;
     constexpr auto tsize =
         TensorDimeToSize<getSpaceDimension<T2toST2Type>()>::value;
     for (unsigned short i = 0; i < ssize; ++i) {
       for (unsigned short j = 0; j < tsize; ++j) {
-        a += abs(v(i, j));
+        a += abs(t(i, j));
       }
     }
     return a;
-  }
+  }  // end of abs
 
-  template <typename T2toST2ResultType,
+  template <T2toST2Concept T2toST2ResultType,
             StensorConcept StensorType,
             TensorConcept TensorType>
-  typename std::enable_if<implementsT2toST2Concept<T2toST2ResultType>() &&
-                              tfel::typetraits::IsFundamentalNumericType<
-                                  numeric_type<TensorType>>::cond &&
-                              isAssignableTo<numeric_type<StensorType>,
-                                             numeric_type<T2toST2ResultType>>(),
-                          void>::type
+  TFEL_HOST_DEVICE constexpr void
   computePushForwardDerivativeWithRespectToDeformationGradient(
-      T2toST2ResultType& dTdF, const StensorType& S, const TensorType& F) {
+      T2toST2ResultType& dTdF,
+      const StensorType& S,
+      const TensorType& F) noexcept  //
+      requires(tfel::typetraits::IsFundamentalNumericType<
+               numeric_type<TensorType>>::cond&&
+                   isAssignableTo<numeric_type<StensorType>,
+                                  numeric_type<T2toST2ResultType>>()) {
     using value_type = numeric_type<T2toST2ResultType>;
     constexpr auto N = getSpaceDimension<T2toST2ResultType>();
     static_assert(getSpaceDimension<StensorType>() == N);
@@ -148,13 +146,17 @@ namespace tfel::math {
     }
   }  // end of computePushForwardDerivativeWithRespectToDeformationGradient
 
-  template <typename T2toST2ResultType,
-            typename T2toST2Type,
+  template <T2toST2Concept T2toST2ResultType,
+            T2toST2Concept T2toST2Type,
             StensorConcept StensorType,
             TensorConcept TensorType>
-  typename std::enable_if<
-      implementsT2toST2Concept<T2toST2ResultType>() &&
-          implementsT2toST2Concept<T2toST2Type>() &&
+  TFEL_HOST_DEVICE constexpr void
+  computeCauchyStressDerivativeFromKirchhoffStressDerivative(
+      T2toST2ResultType& ds,
+      const T2toST2Type& dt_K,
+      const StensorType& s,
+      const TensorType& F) noexcept  //
+      requires(
           getSpaceDimension<T2toST2ResultType>() ==
               getSpaceDimension<T2toST2Type>() &&
           getSpaceDimension<T2toST2ResultType>() ==
@@ -166,25 +168,23 @@ namespace tfel::math {
           isAssignableTo<typename ComputeBinaryResult<numeric_type<T2toST2Type>,
                                                       numeric_type<StensorType>,
                                                       OpPlus>::Result,
-                         numeric_type<T2toST2ResultType>>(),
-      void>::type
-  computeCauchyStressDerivativeFromKirchhoffStressDerivative(
-      T2toST2ResultType& ds,
-      const T2toST2Type& dt_K,
-      const StensorType& s,
-      const TensorType& F) {
+                         numeric_type<T2toST2ResultType>>()) {
     const auto iJ = 1 / det(F);
     const auto dJ = computeDeterminantDerivative(F);
     ds = iJ * (dt_K - (s ^ dJ));
   }
 
-  template <typename T2toST2ResultType,
-            typename T2toST2Type,
+  template <T2toST2Concept T2toST2ResultType,
+            T2toST2Concept T2toST2Type,
             StensorConcept StensorType,
             TensorConcept TensorType>
-  typename std::enable_if<
-      implementsT2toST2Concept<T2toST2ResultType>() &&
-          implementsT2toST2Concept<T2toST2Type>() &&
+  TFEL_HOST_DEVICE constexpr void
+  computeKirchhoffStressDerivativeFromCauchyStressDerivative(
+      T2toST2ResultType& dt_K,
+      const T2toST2Type& ds,
+      const StensorType& s,
+      const TensorType& F) noexcept  //
+      requires(
           getSpaceDimension<T2toST2ResultType>() ==
               getSpaceDimension<T2toST2Type>() &&
           getSpaceDimension<T2toST2ResultType>() ==
@@ -196,13 +196,7 @@ namespace tfel::math {
           isAssignableTo<typename ComputeBinaryResult<numeric_type<T2toST2Type>,
                                                       numeric_type<StensorType>,
                                                       OpPlus>::Result,
-                         numeric_type<T2toST2ResultType>>(),
-      void>::type
-  computeKirchhoffStressDerivativeFromCauchyStressDerivative(
-      T2toST2ResultType& dt_K,
-      const T2toST2Type& ds,
-      const StensorType& s,
-      const TensorType& F) {
+                         numeric_type<T2toST2ResultType>>()) {
     const auto J = det(F);
     const auto dJ = computeDeterminantDerivative(F);
     dt_K = J * ds + (s ^ dJ);
diff --git a/include/TFEL/Math/T2toST2/T2toST2ConceptIO.hxx b/include/TFEL/Math/T2toST2/T2toST2ConceptIO.hxx
index f30134a0e..83024df61 100644
--- a/include/TFEL/Math/T2toST2/T2toST2ConceptIO.hxx
+++ b/include/TFEL/Math/T2toST2/T2toST2ConceptIO.hxx
@@ -19,27 +19,27 @@
 
 namespace tfel::math {
 
-  // Serialisation operator
-  template <typename T2toST2Type>
-  std::enable_if_t<implementsT2toST2Concept<T2toST2Type>(), std::ostream&>
-  operator<<(std::ostream& os, const T2toST2Type& s) {
+  //! \brief serialisation operator
+  TFEL_HOST std::ostream& operator<<(std::ostream& os,
+                                     const T2toST2Concept auto& t) {
+    using T2toST2Type = decltype(t);
     constexpr auto N = getSpaceDimension<T2toST2Type>();
-    constexpr auto stensor_size = StensorDimeToSize<N>::value;
-    constexpr auto tensor_size = TensorDimeToSize<N>::value;
+    constexpr auto ssize = StensorDimeToSize<N>::value;
+    constexpr auto tsize = TensorDimeToSize<N>::value;
     os << "[";
-    for (unsigned short i = 0; i < stensor_size;) {
+    for (unsigned short i = 0; i < ssize;) {
       if (i != 0) {
         os << " [";
       } else {
         os << "[";
       }
-      for (unsigned short j = 0; j < tensor_size;) {
-        os << s(i, j);
-        if (++j != tensor_size) {
+      for (unsigned short j = 0; j < tsize;) {
+        os << t(i, j);
+        if (++j != tsize) {
           os << ",";
         }
       }
-      if (++i != stensor_size) {
+      if (++i != ssize) {
         os << "]\n";
       } else {
         os << "]";
diff --git a/include/TFEL/Math/T2toST2/T2toST2ConceptOperations.hxx b/include/TFEL/Math/T2toST2/T2toST2ConceptOperations.hxx
index 5a0149032..cfbb759af 100644
--- a/include/TFEL/Math/T2toST2/T2toST2ConceptOperations.hxx
+++ b/include/TFEL/Math/T2toST2/T2toST2ConceptOperations.hxx
@@ -144,14 +144,14 @@ namespace tfel::math {
                                       Expr<Result, UnaryOperation<A, OpNeg>>>;
   };
 
-  template <typename T2>
-  TFEL_HOST_DEVICE constexpr auto
-  operator|(const StensorConcept auto& a, const T2& b) noexcept requires(
-      implementsT2toST2Concept<T2>() &&
-      !isInvalid<BinaryOperationResult<decltype(a), T2, OpMult>>()) {
-    typedef BinaryOperationHandler<decltype(a), T2, OpMult> Handle;
+  TFEL_HOST_DEVICE constexpr auto operator|(
+      const StensorConcept auto& a,
+      const T2toST2Concept auto& b) noexcept  //
+      requires(!isInvalid<
+               BinaryOperationResult<decltype(a), decltype(b), OpMult>>()) {
+    typedef BinaryOperationHandler<decltype(a), decltype(b), OpMult> Handle;
     return Handle(a, b);
-  }
+  }  // end of operator |
 
 }  // end of namespace tfel::math
 
diff --git a/include/TFEL/Math/T2toST2/T2toST2T2toT2ProductExpr.hxx b/include/TFEL/Math/T2toST2/T2toST2T2toT2ProductExpr.hxx
index 17d2da7de..88c1e222a 100644
--- a/include/TFEL/Math/T2toST2/T2toST2T2toT2ProductExpr.hxx
+++ b/include/TFEL/Math/T2toST2/T2toST2T2toT2ProductExpr.hxx
@@ -31,10 +31,10 @@ namespace tfel::math {
   };  // end of struct T2toST2T2toT2ProductExpr
 
   //! \brief partial specialisation in 1D
-  template <typename T2toST2ResultType>
+  template <T2toST2Concept T2toST2ResultType>
   struct TFEL_VISIBILITY_LOCAL
       Expr<T2toST2ResultType, T2toST2T2toT2ProductExpr<1u>>
-      : public T2toST2Concept<
+      : public T2toST2ConceptBase<
             Expr<T2toST2ResultType, T2toST2T2toT2ProductExpr<1u>>>,
         public array_holder<
             StensorDimeToSize<getSpaceDimension<T2toST2ResultType>()>::value *
@@ -49,10 +49,9 @@ namespace tfel::math {
      * \param[in] a : first term of the product
      * \param[in] b : second term of the product
      */
-    template <typename T2toST2Type, T2toT2Concept T2toT2Type2>
+    template <T2toST2Concept T2toST2Type, T2toT2Concept T2toT2Type2>
     TFEL_HOST_DEVICE constexpr Expr(const T2toST2Type& a,
                                     const T2toT2Type2& b) noexcept {
-      static_assert(implementsT2toST2Concept<T2toST2Type>());
       static_assert(getSpaceDimension<T2toST2Type>() == 1u);
       static_assert(getSpaceDimension<T2toT2Type2>() == 1u);
       this->v[0] = a(0, 0) * b(0, 0) + a(0, 1) * b(1, 0) + a(0, 2) * b(2, 0);
@@ -81,10 +80,10 @@ namespace tfel::math {
   };  // end of struct Expr<T2toST2ResultType,T2toST2T2toT2ProductExpr>
 
   //! \brief partial specialisation in 2D
-  template <typename T2toST2ResultType>
+  template <T2toST2Concept T2toST2ResultType>
   struct TFEL_VISIBILITY_LOCAL
       Expr<T2toST2ResultType, T2toST2T2toT2ProductExpr<2u>>
-      : public T2toST2Concept<
+      : public T2toST2ConceptBase<
             Expr<T2toST2ResultType, T2toST2T2toT2ProductExpr<2u>>>,
         public array_holder<
             StensorDimeToSize<getSpaceDimension<T2toST2ResultType>()>::value *
@@ -99,10 +98,9 @@ namespace tfel::math {
      * \param[in] a : first term of the product
      * \param[in] b : second term of the product
      */
-    template <typename T2toST2Type, T2toT2Concept T2toT2Type2>
+    template <T2toST2Concept T2toST2Type, T2toT2Concept T2toT2Type2>
     TFEL_HOST_DEVICE constexpr Expr(const T2toST2Type& a,
                                     const T2toT2Type2& b) noexcept {
-      static_assert(implementsT2toST2Concept<T2toST2Type>());
       static_assert(getSpaceDimension<T2toST2Type>() == 2u);
       static_assert(getSpaceDimension<T2toT2Type2>() == 2u);
       this->v[0] = a(0, 0) * b(0, 0) + a(0, 1) * b(1, 0) + a(0, 2) * b(2, 0) +
@@ -162,10 +160,10 @@ namespace tfel::math {
   };  // end of struct Expr<T2toST2ResultType,T2toST2T2toT2ProductExpr>
 
   //! \brief partial specialisation in 3D
-  template <typename T2toST2ResultType>
+  template <T2toST2Concept T2toST2ResultType>
   struct TFEL_VISIBILITY_LOCAL
       Expr<T2toST2ResultType, T2toST2T2toT2ProductExpr<3u>>
-      : public T2toST2Concept<
+      : public T2toST2ConceptBase<
             Expr<T2toST2ResultType, T2toST2T2toT2ProductExpr<3u>>>,
         public array_holder<
             StensorDimeToSize<getSpaceDimension<T2toST2ResultType>()>::value *
@@ -180,10 +178,9 @@ namespace tfel::math {
      * \param[in] a : first term of the product
      * \param[in] b : second term of the product
      */
-    template <typename T2toST2Type, T2toT2Concept T2toT2Type2>
+    template <T2toST2Concept T2toST2Type, T2toT2Concept T2toT2Type2>
     TFEL_HOST_DEVICE constexpr Expr(const T2toST2Type& a,
                                     const T2toT2Type2& b) noexcept {
-      static_assert(implementsT2toST2Concept<T2toST2Type>());
       static_assert(getSpaceDimension<T2toST2Type>() == 3u);
       static_assert(getSpaceDimension<T2toT2Type2>() == 3u);
       this->v[0] = a(0, 0) * b(0, 0) + a(0, 1) * b(1, 0) + a(0, 2) * b(2, 0) +
diff --git a/include/TFEL/Math/T2toST2/T2toST2TensorProductExpr.hxx b/include/TFEL/Math/T2toST2/T2toST2TensorProductExpr.hxx
index 047fd4c72..8a75ddcb4 100644
--- a/include/TFEL/Math/T2toST2/T2toST2TensorProductExpr.hxx
+++ b/include/TFEL/Math/T2toST2/T2toST2TensorProductExpr.hxx
@@ -22,14 +22,14 @@
 
 namespace tfel::math {
 
-  //! Empty structure used for partial specialisation of the
-  //! Expr class
+  /*!
+   * \brief empty structure used for partial specialisation of the
+   * Expr class
+   */
   template <unsigned short N>
   struct T2toST2TensorProductExpr {};  // end of struct T2toST2TensorProductExpr
 
-  /*!
-   * Partial specialisation
-   */
+  //! \brief partial specialisation in 1D
   template <StensorConcept StensorResultType>
   struct Expr<StensorResultType, T2toST2TensorProductExpr<1u>>
       : public StensorConceptBase<
@@ -44,10 +44,9 @@ namespace tfel::math {
      * \param[in] a : first term of the product
      * \param[in] b : second term of the product
      */
-    template <typename T2toST2Type, TensorConcept TensorType>
+    template <T2toST2Concept T2toST2Type, TensorConcept TensorType>
     TFEL_HOST_DEVICE constexpr Expr(const T2toST2Type& a,
                                     const TensorType& b) noexcept {
-      static_assert(implementsT2toST2Concept<T2toST2Type>());
       static_assert(getSpaceDimension<T2toST2Type>() == 1u);
       static_assert(getSpaceDimension<TensorType>() == 1u);
       this->v[0] = a(0, 0) * b(0) + a(0, 1) * b(1) + a(0, 2) * b(2);
@@ -56,6 +55,14 @@ namespace tfel::math {
     }  // end of Expr
     //! using operator []
     using array_holder<3u, value_type>::operator[];
+    /*!
+     * \brief access operator
+     * \param[in] i : index
+     */
+    TFEL_HOST_DEVICE constexpr const value_type& operator[](
+        const unsigned short i) const noexcept {
+      return this->operator()(i);
+    }  // end of operator()
     /*!
      * \brief access operator
      * \param[in] i : index
@@ -65,15 +72,12 @@ namespace tfel::math {
       return this->v[i];
     }  // end of operator()
     //! \return the runtime properties
-    TFEL_HOST_DEVICE constexpr RunTimeProperties getRunTimeProperties() const
-        noexcept {
+    TFEL_HOST_DEVICE constexpr auto getRunTimeProperties() const noexcept {
       return RunTimeProperties();
     }
   };  // end of struct Expr<StensorResultType,T2ST2toST2TensorProductExp>
 
-  /*!
-   * Partial specialisation
-   */
+  //! \brief partial specialisation in 2D
   template <StensorConcept StensorResultType>
   struct Expr<StensorResultType, T2toST2TensorProductExpr<2u>>
       : public StensorConceptBase<
@@ -88,10 +92,9 @@ namespace tfel::math {
      * \param[in] a : first term of the product
      * \param[in] b : second term of the product
      */
-    template <typename T2toST2Type, TensorConcept TensorType>
+    template <T2toST2Concept T2toST2Type, TensorConcept TensorType>
     TFEL_HOST_DEVICE constexpr Expr(const T2toST2Type& a,
                                     const TensorType& b) noexcept {
-      static_assert(implementsT2toST2Concept<T2toST2Type>());
       static_assert(getSpaceDimension<T2toST2Type>() == 2u);
       static_assert(getSpaceDimension<TensorType>() == 2u);
       this->v[0] = a(0, 4) * b[4] + a(0, 3) * b[3] + a(0, 2) * b[2] +
@@ -105,6 +108,14 @@ namespace tfel::math {
     }  // end of Expr
     //! using operator []
     using array_holder<4u, value_type>::operator[];
+    /*!
+     * \brief access operator
+     * \param[in] i : index
+     */
+    TFEL_HOST_DEVICE constexpr const value_type& operator[](
+        const unsigned short i) const noexcept {
+      return this->operator()(i);
+    }  // end of operator()
     /*!
      * \brief access operator
      * \param[in] i : index
@@ -114,15 +125,12 @@ namespace tfel::math {
       return this->v[i];
     }  // end of operator()
     //! \return the runtime properties
-    TFEL_HOST_DEVICE constexpr RunTimeProperties getRunTimeProperties() const
-        noexcept {
+    TFEL_HOST_DEVICE constexpr auto getRunTimeProperties() const noexcept {
       return RunTimeProperties();
     }
   };  // end of struct Expr<StensorResultType,T2ST2toST2TensorProductExp>
 
-  /*!
-   * Partial specialisation
-   */
+  //! \brief partial specialisation in 2D
   template <StensorConcept StensorResultType>
   struct Expr<StensorResultType, T2toST2TensorProductExpr<3u>>
       : public StensorConceptBase<
@@ -137,10 +145,9 @@ namespace tfel::math {
      * \param[in] a : first term of the product
      * \param[in] b : second term of the product
      */
-    template <typename T2toST2Type, TensorConcept TensorType>
+    template <T2toST2Concept T2toST2Type, TensorConcept TensorType>
     TFEL_HOST_DEVICE constexpr Expr(const T2toST2Type& a,
                                     const TensorType& b) noexcept {
-      static_assert(implementsT2toST2Concept<T2toST2Type>());
       static_assert(getSpaceDimension<T2toST2Type>() == 3u);
       static_assert(getSpaceDimension<TensorType>() == 3u);
       this->v[0] = a(0, 8) * b[8] + a(0, 7) * b[7] + a(0, 6) * b[6] +
@@ -164,6 +171,14 @@ namespace tfel::math {
     }  // end of Expr
     //! using operator []
     using array_holder<6u, value_type>::operator[];
+    /*!
+     * \brief access operator
+     * \param[in] i : index
+     */
+    TFEL_HOST_DEVICE constexpr const value_type& operator[](
+        const unsigned short i) const noexcept {
+      return this->operator()(i);
+    }  // end of operator()
     /*!
      * \brief access operator
      * \param[in] i : index
@@ -173,8 +188,7 @@ namespace tfel::math {
       return this->v[i];
     }  // end of operator()
     //! \return the runtime properties
-    TFEL_HOST_DEVICE constexpr RunTimeProperties getRunTimeProperties() const
-        noexcept {
+    TFEL_HOST_DEVICE constexpr auto getRunTimeProperties() const noexcept {
       return RunTimeProperties();
     }
   };  // end of struct Expr<StensorResultType,T2ST2toST2TensorProductExp>
diff --git a/include/TFEL/Math/T2toST2/t2tost2.ixx b/include/TFEL/Math/T2toST2/t2tost2.ixx
index 8e7be7f25..23a3f4ebb 100644
--- a/include/TFEL/Math/T2toST2/t2tost2.ixx
+++ b/include/TFEL/Math/T2toST2/t2tost2.ixx
@@ -28,30 +28,37 @@ namespace tfel::math {
 
   template <unsigned short N, typename T>
   template <TensorConcept TensorType>
-  TFEL_HOST_DEVICE std::enable_if_t<
+  TFEL_HOST_DEVICE constexpr auto
+  t2tost2<N, T>::dCdF(const TensorType& F) noexcept requires(
       getSpaceDimension<TensorType>() == N &&
-          isAssignableTo<numeric_type<TensorType>, T>(),
-      Expr<t2tost2<N, T>, RightCauchyGreenTensorDerivativeExpr<N>>>
-  t2tost2<N, T>::dCdF(const TensorType& F) {
+      isAssignableTo<numeric_type<TensorType>, T>()) {
     return Expr<t2tost2<N, T>, RightCauchyGreenTensorDerivativeExpr<N>>(F);
   }  // end of t2tost2::dCdF
 
   template <unsigned short N, typename T>
   template <TensorConcept TensorType>
-  TFEL_HOST_DEVICE std::enable_if_t<
+  TFEL_HOST_DEVICE constexpr auto
+  t2tost2<N, T>::dBdF(const TensorType& F) noexcept requires(
       getSpaceDimension<TensorType>() == N &&
-          isAssignableTo<numeric_type<TensorType>, T>(),
-      Expr<t2tost2<N, T>, LeftCauchyGreenTensorDerivativeExpr<N>>>
-  t2tost2<N, T>::dBdF(const TensorType& F) {
+      isAssignableTo<numeric_type<TensorType>, T>()) {
     return Expr<t2tost2<N, T>, LeftCauchyGreenTensorDerivativeExpr<N>>(F);
   }  // end of t2tost2::dBdF
 
   template <unsigned short N, typename T>
-  template <typename InputIterator>
-  TFEL_HOST_DEVICE void t2tost2<N, T>::copy(const InputIterator src) {
+  TFEL_HOST_DEVICE constexpr void t2tost2<N, T>::import(
+      const base_type<T>* const src) noexcept {
+    tfel::fsalgo::transform<StensorDimeToSize<N>::value *
+                            TensorDimeToSize<N>::value>::exe(src, this->begin(),
+                                                             [](const auto& v) {
+                                                               return T(v);
+                                                             });
+  }  // end of import
+
+  template <unsigned short N, typename T>
+  TFEL_HOST_DEVICE constexpr void t2tost2<N, T>::copy(const auto p) noexcept {
     tfel::fsalgo::copy<StensorDimeToSize<N>::value *
-                       TensorDimeToSize<N>::value>::exe(src, *this);
-  }
+                       TensorDimeToSize<N>::value>::exe(p, *this);
+  }  // end of copy
 
   TFEL_HOST_DEVICE constexpr auto convertToT2toST2(
       const T2toT2Concept auto& t) noexcept {
@@ -195,22 +202,20 @@ namespace tfel::math {
     }
   }
 
-  template <typename T2toST2Type,
+  template <T2toST2Concept T2toST2Type,
             StensorConcept StensorType,
             TensorConcept TensorType>
-  TFEL_HOST_DEVICE std::enable_if_t<
-      implementsT2toST2Concept<T2toST2Type>() &&
-          getSpaceDimension<T2toST2Type>() ==
-              getSpaceDimension<StensorType>() &&
-          getSpaceDimension<T2toST2Type>() == getSpaceDimension<TensorType>() &&
-          tfel::typetraits::IsFundamentalNumericType<
-              numeric_type<TensorType>>::cond,
-      t2tost2<getSpaceDimension<T2toST2Type>(),
-              typename ComputeBinaryResult<numeric_type<T2toST2Type>,
-                                           numeric_type<StensorType>,
-                                           OpPlus>::Result>>
+  TFEL_HOST_DEVICE constexpr auto
   computeCauchyStressDerivativeFromKirchhoffStressDerivative(
-      const T2toST2Type& dt, const StensorType& s, const TensorType& F) {
+      const T2toST2Type& dt,
+      const StensorType& s,
+      const TensorType& F) noexcept  //
+      requires(getSpaceDimension<T2toST2Type>() ==
+                   getSpaceDimension<StensorType>() &&
+               getSpaceDimension<T2toST2Type>() ==
+                   getSpaceDimension<TensorType>() &&
+               tfel::typetraits::IsFundamentalNumericType<
+                   numeric_type<TensorType>>::cond) {
     using stress =
         typename ComputeBinaryResult<numeric_type<T2toST2Type>,
                                      numeric_type<StensorType>, OpPlus>::Result;
@@ -219,22 +224,20 @@ namespace tfel::math {
     return r;
   }
 
-  template <typename T2toST2Type,
+  template <T2toST2Concept T2toST2Type,
             StensorConcept StensorType,
             TensorConcept TensorType>
-  TFEL_HOST_DEVICE std::enable_if_t<
-      implementsT2toST2Concept<T2toST2Type>() &&
-          getSpaceDimension<T2toST2Type>() ==
-              getSpaceDimension<StensorType>() &&
-          getSpaceDimension<T2toST2Type>() == getSpaceDimension<TensorType>() &&
-          tfel::typetraits::IsFundamentalNumericType<
-              numeric_type<TensorType>>::cond,
-      t2tost2<getSpaceDimension<T2toST2Type>(),
-              typename ComputeBinaryResult<numeric_type<T2toST2Type>,
-                                           numeric_type<StensorType>,
-                                           OpPlus>::Result>>
+  TFEL_HOST_DEVICE constexpr auto
   computeKirchhoffStressDerivativeFromCauchyStressDerivative(
-      const T2toST2Type& ds, const StensorType& s, const TensorType& F) {
+      const T2toST2Type& ds,
+      const StensorType& s,
+      const TensorType& F) noexcept  //
+      requires(getSpaceDimension<T2toST2Type>() ==
+                   getSpaceDimension<StensorType>() &&
+               getSpaceDimension<T2toST2Type>() ==
+                   getSpaceDimension<TensorType>() &&
+               tfel::typetraits::IsFundamentalNumericType<
+                   numeric_type<TensorType>>::cond) {
     using stress =
         typename ComputeBinaryResult<numeric_type<T2toST2Type>,
                                      numeric_type<StensorType>, OpPlus>::Result;
@@ -243,24 +246,43 @@ namespace tfel::math {
     return r;
   }
 
-  template <typename T2toST2ResultType,
-            typename T2toST2Type,
+  template <T2toST2Concept T2toST2Type,
+            StensorConcept StensorType,
+            TensorConcept TensorType>
+  TFEL_HOST_DEVICE constexpr auto computePushForwardDerivative(
+      const T2toST2Type& K,
+      const StensorType& S,
+      const TensorType& F) noexcept  //
+      requires(getSpaceDimension<StensorType>() ==
+                   getSpaceDimension<T2toST2Type>() &&
+               getSpaceDimension<TensorType>() ==
+                   getSpaceDimension<T2toST2Type>() &&
+               tfel::typetraits::IsFundamentalNumericType<
+                   numeric_type<TensorType>>::cond) {
+    using stress =
+        typename ComputeBinaryResult<numeric_type<T2toST2Type>,
+                                     numeric_type<StensorType>, OpPlus>::Result;
+    t2tost2<getSpaceDimension<T2toST2Type>(), stress> r;
+    computePushForwardDerivative(r, K, S, F);
+    return r;
+  }
+
+  template <T2toST2Concept T2toST2ResultType,
+            T2toST2Concept T2toST2Type,
             StensorConcept StensorType,
             TensorConcept TensorType>
-  TFEL_HOST_DEVICE typename std::enable_if<
-      implementsT2toST2Concept<T2toST2ResultType>() &&
-          implementsT2toST2Concept<T2toST2Type>() &&
-          tfel::typetraits::IsFundamentalNumericType<
-              numeric_type<TensorType>>::cond &&
-          isAssignableTo<typename ComputeBinaryResult<numeric_type<T2toST2Type>,
-                                                      numeric_type<StensorType>,
-                                                      OpPlus>::Result,
-                         numeric_type<T2toST2ResultType>>(),
-      void>::type
-  computePushForwardDerivative(T2toST2ResultType& dT_dF,
-                               const T2toST2Type& dS_dF,
-                               const StensorType& S,
-                               const TensorType& F) {
+  TFEL_HOST_DEVICE constexpr void computePushForwardDerivative(
+      T2toST2ResultType& dT_dF,
+      const T2toST2Type& dS_dF,
+      const StensorType& S,
+      const TensorType& F) noexcept  //
+      requires(tfel::typetraits::IsFundamentalNumericType<
+               numeric_type<TensorType>>::cond&&
+                   isAssignableTo<
+                       typename ComputeBinaryResult<numeric_type<T2toST2Type>,
+                                                    numeric_type<StensorType>,
+                                                    OpPlus>::Result,
+                       numeric_type<T2toST2ResultType>>()) {
     constexpr auto N = getSpaceDimension<T2toST2ResultType>();
     static_assert(getSpaceDimension<T2toST2Type>() == N);
     static_assert(getSpaceDimension<StensorType>() == N);
@@ -274,55 +296,21 @@ namespace tfel::math {
     dT_dF = d1 + d2 * dS_dF;
   }  // end of computePushForwardDerivative
 
-  template <typename T2toST2Type,
-            StensorConcept StensorType,
-            TensorConcept TensorType>
-  TFEL_HOST_DEVICE std::enable_if_t<
-      implementsT2toST2Concept<T2toST2Type>() &&
-          getSpaceDimension<StensorType>() ==
-              getSpaceDimension<T2toST2Type>() &&
-          getSpaceDimension<TensorType>() == getSpaceDimension<T2toST2Type>() &&
-          tfel::typetraits::IsFundamentalNumericType<
-              numeric_type<TensorType>>::cond,
-      t2tost2<getSpaceDimension<T2toST2Type>(),
-              typename ComputeBinaryResult<numeric_type<T2toST2Type>,
-                                           numeric_type<StensorType>,
-                                           OpPlus>::Result>>
-  computePushForwardDerivative(const T2toST2Type& K,
-                               const StensorType& S,
-                               const TensorType& F) {
-    using stress =
-        typename ComputeBinaryResult<numeric_type<T2toST2Type>,
-                                     numeric_type<StensorType>, OpPlus>::Result;
-    t2tost2<getSpaceDimension<T2toST2Type>(), stress> r;
-    computePushForwardDerivative(r, K, S, F);
-    return r;
-  }
-
-  template <typename T2toST2Type>
-  TFEL_HOST_DEVICE
-      std::enable_if_t<((implementsT2toST2Concept<T2toST2Type>()) &&
-                        (getSpaceDimension<T2toST2Type>() == 1u)),
-                       t2tost2<1u, numeric_type<T2toST2Type>>>
-      change_basis(const T2toST2Type& s,
-                   const rotation_matrix<numeric_type<T2toST2Type>>&) {
-    return s;
-  }  // end of change_basis
-
-  template <typename T2toST2Type>
-  TFEL_HOST_DEVICE std::enable_if_t<
-      ((implementsT2toST2Concept<T2toST2Type>()) &&
-       (getSpaceDimension<T2toST2Type>() != 1u)),
-      t2tost2<getSpaceDimension<T2toST2Type>(), numeric_type<T2toST2Type>>>
-  change_basis(const T2toST2Type& s,
-               const rotation_matrix<numeric_type<T2toST2Type>>& r) {
+  template <T2toST2Concept T2toST2Type>
+  TFEL_HOST_DEVICE constexpr auto change_basis(
+      const T2toST2Type& s,
+      const rotation_matrix<numeric_type<T2toST2Type>>& r) noexcept {
     constexpr auto N = getSpaceDimension<T2toST2Type>();
-    using t2tot2 = tfel::math::t2tot2<N, numeric_type<T2toST2Type>>;
-    using st2tost2 = tfel::math::st2tost2<N, numeric_type<T2toST2Type>>;
-    const auto sr = st2tost2::fromRotationMatrix(r);
-    const auto sir = t2tot2::fromRotationMatrix(transpose(r));
-    return sr * s * sir;
-    return s;
+    using Result = t2tost2<N, numeric_type<T2toST2Type>>;
+    if constexpr (N == 1) {
+      return Result{s};
+    } else {
+      using t2tot2 = tfel::math::t2tot2<N, numeric_type<T2toST2Type>>;
+      using st2tost2 = tfel::math::st2tost2<N, numeric_type<T2toST2Type>>;
+      const auto sr = st2tost2::fromRotationMatrix(r);
+      const auto sir = t2tot2::fromRotationMatrix(transpose(r));
+      return Result{sr * s * sir};
+    }
   }  // end of change_basis
 
 }  // end of namespace tfel::math
diff --git a/include/TFEL/Math/T2toT2/T2toT2Concept.hxx b/include/TFEL/Math/T2toT2/T2toT2Concept.hxx
index 6bf517bf9..2d78dc0bf 100644
--- a/include/TFEL/Math/T2toT2/T2toT2Concept.hxx
+++ b/include/TFEL/Math/T2toT2/T2toT2Concept.hxx
@@ -32,7 +32,7 @@ namespace tfel::math {
   struct T2toT2Tag {};
   /*!
    * \brief an helper class that simply exposes publically a member named
-   * ConceptTag as an alias to ST2toST2Tag.
+   * ConceptTag as an alias to T2toT2Tag.
    *
    * The main reason for this alias is to properly implement the `ConceptRebind`
    * metafunction.
@@ -42,9 +42,10 @@ namespace tfel::math {
     using ConceptTag = T2toT2Tag;
   };
   /*!
-   * \brief definition of the ST2toST2Concept
-   * a class matching the stensor concept must expose the `ST2toST2Tag` and have
-   * access operators.
+   * \brief definition of the T2toT2Concept
+   *
+   * A class matching this concept must expose the `T2toT2Tag` and have
+   * appropriate access operators.
    */
   template <typename T>
   concept T2toT2Concept =
diff --git a/include/TFEL/Math/T2toT2/t2tot2.ixx b/include/TFEL/Math/T2toT2/t2tot2.ixx
index 84871847f..f6833f9c3 100644
--- a/include/TFEL/Math/T2toT2/t2tot2.ixx
+++ b/include/TFEL/Math/T2toT2/t2tot2.ixx
@@ -361,15 +361,22 @@ namespace tfel::math {
   }  // end of K
 
   template <unsigned short N, typename T>
-  template <
-      typename T2toST2Type,
-      std::enable_if_t<((implementsT2toST2Concept<T2toST2Type>()) &&
-                        (isAssignableTo<numeric_type<T2toST2Type>, T>()) &&
-                        (getSpaceDimension<T2toST2Type>() == N)),
-                       bool>>
-  constexpr t2tot2<N, T>::t2tot2(const T2toST2Type& s) {
+  template <T2toST2Concept T2toST2Type>
+  constexpr t2tot2<N, T>::t2tot2(const T2toST2Type& s) noexcept requires(
+      (isAssignableTo<numeric_type<T2toST2Type>, T>()) &&
+      (getSpaceDimension<T2toST2Type>() == N)) {
     convert(*this, s);
-  }  // end of t2tot2<N, T>::t2tot2
+  }  // end of t2tot2
+
+  template <unsigned short N, typename T>
+  TFEL_HOST_DEVICE constexpr void t2tot2<N, T>::import(
+      const base_type<T>* const src) noexcept {
+    tfel::fsalgo::transform<TensorDimeToSize<N>::value *
+                            TensorDimeToSize<N>::value>::exe(src, this->begin(),
+                                                             [](const auto& v) {
+                                                               return T(v);
+                                                             });
+  }  // end of import
 
   template <unsigned short N, typename T>
   TFEL_HOST_DEVICE constexpr void t2tot2<N, T>::copy(const auto p) noexcept {
@@ -447,13 +454,11 @@ namespace tfel::math {
     }
   }  // end of computeDeterminantSecondDerivative
 
-  template <typename T, typename T2toST2Type>
-  TFEL_HOST_DEVICE
-      std::enable_if_t<((implementsT2toST2Concept<T2toST2Type>()) &&
-                        (isAssignableTo<numeric_type<T2toST2Type>, T>())),
-                       void>
-      convert(t2tot2<getSpaceDimension<T2toST2Type>(), T>& d,
-              const T2toST2Type& s) {
+  template <typename T, T2toST2Concept T2toST2Type>
+  TFEL_HOST_DEVICE constexpr void convert(
+      t2tot2<getSpaceDimension<T2toST2Type>(), T>& d,
+      const T2toST2Type& s) noexcept  //
+      requires(isAssignableTo<numeric_type<T2toST2Type>, T>()) {
     constexpr auto N = getSpaceDimension<T2toST2Type>();
     static_assert((N == 1) || (N == 2) || (N == 3));
     if constexpr (N == 1) {
diff --git a/include/TFEL/Math/Tensor/tensor.ixx b/include/TFEL/Math/Tensor/tensor.ixx
index 7c4b80a63..31b9ae272 100644
--- a/include/TFEL/Math/Tensor/tensor.ixx
+++ b/include/TFEL/Math/Tensor/tensor.ixx
@@ -137,9 +137,9 @@ namespace tfel::math {
     tfel::fsalgo::copy<TensorDimeToSize<N>::value>::exe(p, *this);
   }  // end of copy
 
-  template <unsigned short N, typename T, typename OutputIterator>
-  TFEL_HOST_DEVICE TFEL_MATH_INLINE2 std::enable_if_t<isScalar<T>(), void>
-  exportToBaseTypeArray(const tensor<N, T>& t, OutputIterator p) {
+  template <unsigned short N, typename T>
+  TFEL_HOST_DEVICE constexpr void exportToBaseTypeArray(
+      const tensor<N, T>& t, const auto p) noexcept requires(isScalar<T>()) {
     tfel::fsalgo::transform<TensorDimeToSize<N>::value>::exe(
         t.cbegin(), p, [](const auto& v) { return base_type_cast(v); });
   }
diff --git a/include/TFEL/Math/st2tost2.hxx b/include/TFEL/Math/st2tost2.hxx
index 8932c8397..4a7e42c83 100644
--- a/include/TFEL/Math/st2tost2.hxx
+++ b/include/TFEL/Math/st2tost2.hxx
@@ -147,19 +147,17 @@ namespace tfel::math {
      * convert a T2toST2 to a ST2toST2
      * \param[in] src : T2toST2 to be converted
      */
-    template <typename T2toST2Type>
-    static TFEL_MATH_INLINE std::enable_if_t<
-        implementsT2toST2Concept<T2toST2Type>() &&
-            getSpaceDimension<T2toST2Type>() == N &&
-            isAssignableTo<numeric_type<T2toST2Type>, ValueType>(),
-        Expr<st2tost2<N, ValueType>, ConvertT2toST2ToST2toST2Expr<N>>>
-    convert(const T2toST2Type&);
+    template <T2toST2Concept T2toST2Type>
+    TFEL_HOST_DEVICE static constexpr auto
+    convert(const T2toST2Type&) noexcept requires(
+        getSpaceDimension<T2toST2Type>() == N &&
+        isAssignableTo<numeric_type<T2toST2Type>, ValueType>());
     /*!
      * build the equivalent st2tost2 from a rotation matrix
      * \param[in] r : rotation matrix
      */
-    static tfel::math::st2tost2<N, base_type<ValueType>> fromRotationMatrix(
-        const rotation_matrix<ValueType>&);
+    TFEL_HOST_DEVICE static constexpr auto fromRotationMatrix(
+        const rotation_matrix<ValueType>&) noexcept;
     /*!
      * \brief compute the derivative of the symmetric tensor product:
      * \f[
@@ -168,11 +166,10 @@ namespace tfel::math {
      * \param[in] s: second tensor of the product
      */
     template <StensorConcept StensorType>
-    static TFEL_MATH_INLINE std::enable_if_t<
+    TFEL_HOST_DEVICE static constexpr auto
+    stpd(const StensorType&) noexcept requires(
         getSpaceDimension<StensorType>() == N &&
-            isAssignableTo<numeric_type<StensorType>, ValueType>(),
-        tfel::math::st2tost2<N, ValueType>>
-    stpd(const StensorType&);
+        isAssignableTo<numeric_type<StensorType>, ValueType>());
     //
     static constexpr st2tost2 Id() noexcept;
     static constexpr st2tost2 IxI() noexcept;
@@ -216,16 +213,14 @@ namespace tfel::math {
    * \param[in] r : rotation matrix
    */
   template <ST2toST2Concept ST2toST2Type>
-  TFEL_HOST_DEVICE constexpr st2tost2<getSpaceDimension<ST2toST2Type>(),
-                                      numeric_type<ST2toST2Type>>
-  change_basis(const ST2toST2Type&,
-               const rotation_matrix<numeric_type<ST2toST2Type>>&) noexcept;
+  TFEL_HOST_DEVICE constexpr auto change_basis(
+      const ST2toST2Type&,
+      const rotation_matrix<numeric_type<ST2toST2Type>>&) noexcept;
   /*!
    * \return the invert of a st2tost2
    * \param[in] s : st2tost2 to be inverted
    */
   TFEL_HOST constexpr auto invert(const ST2toST2Concept auto&);
-
   /*!
    * \return the push-forward of a fourth order tensor:
    * \f[
diff --git a/include/TFEL/Math/st2tot2.hxx b/include/TFEL/Math/st2tot2.hxx
index d0deb0937..0209475ca 100644
--- a/include/TFEL/Math/st2tot2.hxx
+++ b/include/TFEL/Math/st2tot2.hxx
@@ -39,13 +39,11 @@ namespace tfel::math {
    * \brief partial specialisation of the `DerivativeTypeDispatcher`
    * metafunction.
    */
-  template <typename ST2toT2Type, typename ScalarType>
+  template <ST2toT2Concept ST2toT2Type, typename ScalarType>
   struct DerivativeTypeDispatcher<ST2toT2Tag,
                                   ScalarTag,
                                   ST2toT2Type,
                                   ScalarType> {
-    static_assert(implementsST2toT2Concept<ST2toT2Type>(),
-                  "template argument ST2toT2Type is not a st2tot2");
     static_assert(isScalar<ScalarType>(),
                   "template argument ScalarType is not a scalar");
     static_assert(isScalar<numeric_type<ST2toT2Type>>(),
@@ -59,13 +57,11 @@ namespace tfel::math {
    * \brief partial specialisation of the `DerivativeTypeDispatcher`
    * metafunction.
    */
-  template <typename ScalarType, typename ST2toT2Type>
+  template <typename ScalarType, ST2toT2Concept ST2toT2Type>
   struct DerivativeTypeDispatcher<ScalarTag,
                                   ST2toT2Tag,
                                   ScalarType,
                                   ST2toT2Type> {
-    static_assert(implementsST2toT2Concept<ST2toT2Type>(),
-                  "template argument ST2toT2Type is not a st2tot2");
     static_assert(isScalar<ScalarType>(),
                   "template argument ScalarType is not a scalar");
     static_assert(isScalar<numeric_type<ST2toT2Type>>(),
@@ -100,7 +96,7 @@ namespace tfel::math {
 
   template <unsigned short N, typename ValueType>
   struct st2tot2
-      : ST2toT2Concept<st2tot2<N, ValueType>>,
+      : ST2toT2ConceptBase<st2tot2<N, ValueType>>,
         GenericFixedSizeArray<
             st2tot2<N, ValueType>,
             FixedSizeRowMajorMatrixPolicy<TensorDimeToSize<N>::value,
@@ -118,51 +114,47 @@ namespace tfel::math {
      * \return the left part of the derivative of a tensor product
      */
     template <StensorConcept StensorType>
-    static TFEL_MATH_INLINE std::enable_if_t<
+    static TFEL_HOST_DEVICE constexpr auto
+    tpld(const StensorType&) noexcept requires(
         getSpaceDimension<StensorType>() == N &&
-            isAssignableTo<numeric_type<StensorType>, ValueType>(),
-        Expr<st2tot2<N, ValueType>, StensorProductLeftDerivativeExpr<N>>>
-    tpld(const StensorType&);
+        isAssignableTo<numeric_type<StensorType>, ValueType>());
     /*!
      * \param[in] B : second tensor of the product
      * \param[in] C : derivative of the first tensor
      * \return the left part of the derivative of a tensor product
      */
     template <StensorConcept StensorType, ST2toST2Concept ST2toST2Type>
-    static TFEL_MATH_INLINE std::enable_if_t<
+    static TFEL_HOST_DEVICE constexpr auto
+    tpld(const StensorType&, const ST2toST2Type&) noexcept requires(
         getSpaceDimension<StensorType>() == N &&
-            getSpaceDimension<ST2toST2Type>() == N &&
-            isAssignableTo<BinaryOperationResult<numeric_type<StensorType>,
-                                                 numeric_type<ST2toST2Type>,
-                                                 OpMult>,
-                           ValueType>(),
-        Expr<st2tot2<N, ValueType>, StensorProductLeftDerivativeExpr<N>>>
-    tpld(const StensorType&, const ST2toST2Type&);
+        getSpaceDimension<ST2toST2Type>() == N &&
+        isAssignableTo<BinaryOperationResult<numeric_type<StensorType>,
+                                             numeric_type<ST2toST2Type>,
+                                             OpMult>,
+                       ValueType>());
     /*!
      * \param[in] A : first tensor of the product
      * \return the right part of the derivative of a tensor product
      */
     template <StensorConcept StensorType>
-    TFEL_HOST_DEVICE static constexpr std::enable_if_t<
+    TFEL_HOST_DEVICE static constexpr auto
+    tprd(const StensorType&) noexcept requires(
         getSpaceDimension<StensorType>() == N &&
-            isAssignableTo<numeric_type<StensorType>, ValueType>(),
-        Expr<st2tot2<N, ValueType>, StensorProductRightDerivativeExpr<N>>>
-    tprd(const StensorType&) noexcept;
+        isAssignableTo<numeric_type<StensorType>, ValueType>());
     /*!
      * \param[in] A : first tensor of the product
      * \param[in] C : derivative of the first tensor
      * \return the right part of the derivative of a tensor product
      */
     template <StensorConcept StensorType, ST2toST2Concept ST2toST2Type>
-    TFEL_HOST_DEVICE static constexpr std::enable_if_t<
+    TFEL_HOST_DEVICE static constexpr auto
+    tprd(const StensorType&, const ST2toST2Type&) noexcept requires(
         getSpaceDimension<StensorType>() == N &&
-            getSpaceDimension<ST2toST2Type>() == N &&
-            isAssignableTo<BinaryOperationResult<numeric_type<StensorType>,
-                                                 numeric_type<ST2toST2Type>,
-                                                 OpMult>,
-                           ValueType>(),
-        Expr<st2tot2<N, ValueType>, StensorProductRightDerivativeExpr<N>>>
-    tprd(const StensorType&, const ST2toST2Type&) noexcept;
+        getSpaceDimension<ST2toST2Type>() == N &&
+        isAssignableTo<BinaryOperationResult<numeric_type<StensorType>,
+                                             numeric_type<ST2toST2Type>,
+                                             OpMult>,
+                       ValueType>());
     //
     TFEL_MATH_FIXED_SIZE_ARRAY_DEFAULT_METHODS(st2tot2,
                                                GenericFixedSizeArrayBase);
@@ -170,8 +162,9 @@ namespace tfel::math {
     using GenericFixedSizeArrayBase::operator[];
     using GenericFixedSizeArrayBase::operator();
     //! \brief import values from an external memory location
-    void import(const base_type<ValueType>* const);
-    //
+    TFEL_HOST_DEVICE constexpr void import(
+        const base_type<ValueType>* const) noexcept;
+    //! \brief copy from a range
     TFEL_HOST_DEVICE constexpr void copy(const auto) noexcept;
   };
 
diff --git a/include/TFEL/Math/stensor.hxx b/include/TFEL/Math/stensor.hxx
index 01010f0e5..83e294ec1 100644
--- a/include/TFEL/Math/stensor.hxx
+++ b/include/TFEL/Math/stensor.hxx
@@ -767,10 +767,9 @@ namespace tfel::math {
    * \return the rotated symmetric tensor
    */
   template <StensorConcept StensorType>
-  TFEL_HOST_DEVICE TFEL_MATH_INLINE2
-      stensor<getSpaceDimension<StensorType>(), numeric_type<StensorType>>
-      change_basis(const StensorType&,
-                   const rotation_matrix<numeric_type<StensorType>>&);
+  TFEL_HOST_DEVICE constexpr auto change_basis(
+      const StensorType&,
+      const rotation_matrix<numeric_type<StensorType>>&) noexcept;
   /*!
    * \return the inverse of a symmetric tensor
    * \param[in] s: symmetric tensor
diff --git a/include/TFEL/Math/t2tost2.hxx b/include/TFEL/Math/t2tost2.hxx
index c12120ebb..c6b226fdb 100644
--- a/include/TFEL/Math/t2tost2.hxx
+++ b/include/TFEL/Math/t2tost2.hxx
@@ -46,13 +46,11 @@ namespace tfel::math {
    * \brief partial specialisation of the `DerivativeTypeDispatcher`
    * metafunction.
    */
-  template <typename T2toST2Type, typename ScalarType>
+  template <T2toST2Concept T2toST2Type, typename ScalarType>
   struct DerivativeTypeDispatcher<T2toST2Tag,
                                   ScalarTag,
                                   T2toST2Type,
                                   ScalarType> {
-    static_assert(implementsT2toST2Concept<T2toST2Type>(),
-                  "template argument T2toST2Type is not a t2tost2");
     static_assert(isScalar<ScalarType>(),
                   "template argument ScalarType is not a scalar");
     static_assert(isScalar<numeric_type<T2toST2Type>>(),
@@ -66,13 +64,11 @@ namespace tfel::math {
    * \brief partial specialisation of the `DerivativeTypeDispatcher`
    * metafunction.
    */
-  template <typename ScalarType, typename T2toST2Type>
+  template <typename ScalarType, T2toST2Concept T2toST2Type>
   struct DerivativeTypeDispatcher<ScalarTag,
                                   T2toST2Tag,
                                   ScalarType,
                                   T2toST2Type> {
-    static_assert(implementsT2toST2Concept<T2toST2Type>(),
-                  "template argument T2toST2Type is not a t2tost2");
     static_assert(isScalar<ScalarType>(),
                   "template argument ScalarType is not a scalar");
     static_assert(isScalar<numeric_type<T2toST2Type>>(),
@@ -107,7 +103,7 @@ namespace tfel::math {
   //
   template <unsigned short N, typename ValueType>
   struct t2tost2
-      : T2toST2Concept<t2tost2<N, ValueType>>,
+      : T2toST2ConceptBase<t2tost2<N, ValueType>>,
         GenericFixedSizeArray<
             t2tost2<N, ValueType>,
             FixedSizeRowMajorMatrixPolicy<StensorDimeToSize<N>::value,
@@ -126,22 +122,20 @@ namespace tfel::math {
      * with respect to the deformation gradient
      */
     template <TensorConcept TensorType>
-    TFEL_HOST_DEVICE static TFEL_MATH_INLINE std::enable_if_t<
+    TFEL_HOST_DEVICE static TFEL_HOST_DEVICE constexpr auto
+    dCdF(const TensorType&) noexcept requires(
         getSpaceDimension<TensorType>() == N &&
-            isAssignableTo<numeric_type<TensorType>, ValueType>(),
-        Expr<t2tost2<N, ValueType>, RightCauchyGreenTensorDerivativeExpr<N>>>
-    dCdF(const TensorType&);
+        isAssignableTo<numeric_type<TensorType>, ValueType>());
     /*!
      * \param[in] F : deformation gradient
      * \return the derivative of the Cauchy left symmetric tensor
      * with respect to the deformation gradient
      */
     template <TensorConcept TensorType>
-    TFEL_HOST_DEVICE static TFEL_MATH_INLINE std::enable_if_t<
+    TFEL_HOST_DEVICE static TFEL_HOST_DEVICE constexpr auto
+    dBdF(const TensorType&) noexcept requires(
         getSpaceDimension<TensorType>() == N &&
-            isAssignableTo<numeric_type<TensorType>, ValueType>(),
-        Expr<t2tost2<N, ValueType>, LeftCauchyGreenTensorDerivativeExpr<N>>>
-    dBdF(const TensorType&);
+        isAssignableTo<numeric_type<TensorType>, ValueType>());
     //
     TFEL_MATH_FIXED_SIZE_ARRAY_DEFAULT_METHODS(t2tost2,
                                                GenericFixedSizeArrayBase);
@@ -149,13 +143,10 @@ namespace tfel::math {
     using GenericFixedSizeArrayBase::operator[];
     using GenericFixedSizeArrayBase::operator();
     //! \brief import values
-    template <typename T2>
-    TFEL_HOST_DEVICE
-        std::enable_if_t<std::is_convertible_v<T2, base_type<ValueType>>, void>
-        import(const T2* const);
-
-    template <typename InputIterator>
-    TFEL_HOST_DEVICE TFEL_MATH_INLINE2 void copy(const InputIterator src);
+    TFEL_HOST_DEVICE constexpr void import(
+        const base_type<ValueType>* const) noexcept;
+    //! \brief copy from a range
+    TFEL_HOST_DEVICE constexpr void copy(const auto) noexcept;
   };
 
   /*!
@@ -174,30 +165,14 @@ namespace tfel::math {
   using ConstT2toST2View = ConstView<t2tost2<N, T>>;
 
   /*!
-   * \return change the basis of a 1D t2tost2
-   * \param[in] s: t2tost2
-   * \param[in] r: rotation matrix
-   */
-  template <typename T2toST2Type>
-  TFEL_HOST_DEVICE
-      std::enable_if_t<((implementsT2toST2Concept<T2toST2Type>()) &&
-                        (getSpaceDimension<T2toST2Type>() == 1u)),
-                       t2tost2<1u, numeric_type<T2toST2Type>>>
-      change_basis(const T2toST2Type&,
-                   const rotation_matrix<numeric_type<T2toST2Type>>&);
-  /*!
-   * \return change the basis of a 2D or 3D t2tost2
+   * \return change the basis of a t2tost2
    * \param[in] s: t2tost2
    * \param[in] r: rotation matrix
    */
-  template <typename T2toST2Type>
-  TFEL_HOST_DEVICE std::enable_if_t<
-      ((implementsT2toST2Concept<T2toST2Type>()) &&
-       (getSpaceDimension<T2toST2Type>() != 1u)),
-      t2tost2<getSpaceDimension<T2toST2Type>(), numeric_type<T2toST2Type>>>
-  change_basis(const T2toST2Type& s,
-               const rotation_matrix<numeric_type<T2toST2Type>>&);
-
+  template <T2toST2Concept T2toST2Type>
+  TFEL_HOST_DEVICE constexpr auto change_basis(
+      const T2toST2Type&,
+      const rotation_matrix<numeric_type<T2toST2Type>>&) noexcept;
   /*!
    * \brief convert a t2tot2 tensor to a t2tost2
    * \param[in] t: fourth order tensor to be converted
@@ -219,23 +194,20 @@ namespace tfel::math {
    * \param[in]  s:  Cauchy stress
    * \param[in]  F:  deformation gradient
    */
-  template <typename T2toST2Type,
+  template <T2toST2Concept T2toST2Type,
             StensorConcept StensorType,
             TensorConcept TensorType>
-  TFEL_HOST_DEVICE std::enable_if_t<
-      implementsT2toST2Concept<T2toST2Type>() &&
-          getSpaceDimension<T2toST2Type>() ==
-              getSpaceDimension<StensorType>() &&
-          getSpaceDimension<T2toST2Type>() == getSpaceDimension<TensorType>() &&
-          tfel::typetraits::IsFundamentalNumericType<
-              numeric_type<TensorType>>::cond,
-      t2tost2<getSpaceDimension<T2toST2Type>(),
-              typename ComputeBinaryResult<numeric_type<T2toST2Type>,
-                                           numeric_type<StensorType>,
-                                           OpPlus>::Result>>
-  computeCauchyStressDerivativeFromKirchhoffStressDerivative(const T2toST2Type&,
-                                                             const StensorType&,
-                                                             const TensorType&);
+  TFEL_HOST_DEVICE constexpr auto
+  computeCauchyStressDerivativeFromKirchhoffStressDerivative(
+      const T2toST2Type&,
+      const StensorType&,
+      const TensorType&) noexcept  //
+      requires(getSpaceDimension<T2toST2Type>() ==
+                   getSpaceDimension<StensorType>() &&
+               getSpaceDimension<T2toST2Type>() ==
+                   getSpaceDimension<TensorType>() &&
+               tfel::typetraits::IsFundamentalNumericType<
+                   numeric_type<TensorType>>::cond);
   /*!
    * \brief compute the Cauchy stress derivative from the Kirchhoff stress
    * derivative
@@ -243,23 +215,20 @@ namespace tfel::math {
    * \param[in]  s:  Cauchy stress
    * \param[in]  F:  deformation gradient
    */
-  template <typename T2toST2Type,
+  template <T2toST2Concept T2toST2Type,
             StensorConcept StensorType,
             TensorConcept TensorType>
-  TFEL_HOST_DEVICE std::enable_if_t<
-      implementsT2toST2Concept<T2toST2Type>() &&
-          getSpaceDimension<T2toST2Type>() ==
-              getSpaceDimension<StensorType>() &&
-          getSpaceDimension<T2toST2Type>() == getSpaceDimension<TensorType>() &&
-          tfel::typetraits::IsFundamentalNumericType<
-              numeric_type<TensorType>>::cond,
-      t2tost2<getSpaceDimension<T2toST2Type>(),
-              typename ComputeBinaryResult<numeric_type<T2toST2Type>,
-                                           numeric_type<StensorType>,
-                                           OpPlus>::Result>>
-  computeKirchhoffStressDerivativeFromCauchyStressDerivative(const T2toST2Type&,
-                                                             const StensorType&,
-                                                             const TensorType&);
+  TFEL_HOST_DEVICE constexpr auto
+  computeKirchhoffStressDerivativeFromCauchyStressDerivative(
+      const T2toST2Type&,
+      const StensorType&,
+      const TensorType&) noexcept  //
+      requires(getSpaceDimension<T2toST2Type>() ==
+                   getSpaceDimension<StensorType>() &&
+               getSpaceDimension<T2toST2Type>() ==
+                   getSpaceDimension<TensorType>() &&
+               tfel::typetraits::IsFundamentalNumericType<
+                   numeric_type<TensorType>>::cond);
   /*!
    * \return the derivative of the push-forward of a symmetric
    * second order tensor with respect to the deformation gradient
@@ -274,23 +243,19 @@ namespace tfel::math {
    * \param[in]  : orginal tensor (second Piola-Kirschoff stress)
    * \param[in]  : deformation gradient
    */
-  template <typename T2toST2Type,
+  template <T2toST2Concept T2toST2Type,
             StensorConcept StensorType,
             TensorConcept TensorType>
-  TFEL_HOST_DEVICE std::enable_if_t<
-      implementsT2toST2Concept<T2toST2Type>() &&
-          getSpaceDimension<StensorType>() ==
-              getSpaceDimension<T2toST2Type>() &&
-          getSpaceDimension<TensorType>() == getSpaceDimension<T2toST2Type>() &&
-          tfel::typetraits::IsFundamentalNumericType<
-              numeric_type<TensorType>>::cond,
-      t2tost2<getSpaceDimension<T2toST2Type>(),
-              typename ComputeBinaryResult<numeric_type<T2toST2Type>,
-                                           numeric_type<StensorType>,
-                                           OpPlus>::Result>>
-  computePushForwardDerivative(const T2toST2Type&,
-                               const StensorType&,
-                               const TensorType&);
+  TFEL_HOST_DEVICE constexpr auto computePushForwardDerivative(
+      const T2toST2Type&,
+      const StensorType&,
+      const TensorType&) noexcept  //
+      requires(getSpaceDimension<StensorType>() ==
+                   getSpaceDimension<T2toST2Type>() &&
+               getSpaceDimension<TensorType>() ==
+                   getSpaceDimension<T2toST2Type>() &&
+               tfel::typetraits::IsFundamentalNumericType<
+                   numeric_type<TensorType>>::cond);
   /*!
    * \brief compute de derivative of the push-forward of a symmetric
    * second order tensor with respect to the deformation gradient
@@ -306,24 +271,22 @@ namespace tfel::math {
    * \param[in]  : orginal tensor (second Piola-Kirschoff stress)
    * \param[in]  : deformation gradient
    */
-  template <typename T2toST2ResultType,
-            typename T2toST2Type,
+  template <T2toST2Concept T2toST2ResultType,
+            T2toST2Concept T2toST2Type,
             StensorConcept StensorType,
             TensorConcept TensorType>
-  TFEL_HOST_DEVICE typename std::enable_if<
-      implementsT2toST2Concept<T2toST2ResultType>() &&
-          implementsT2toST2Concept<T2toST2Type>() &&
-          tfel::typetraits::IsFundamentalNumericType<
-              numeric_type<TensorType>>::cond &&
-          isAssignableTo<typename ComputeBinaryResult<numeric_type<T2toST2Type>,
-                                                      numeric_type<StensorType>,
-                                                      OpPlus>::Result,
-                         numeric_type<T2toST2ResultType>>(),
-      void>::type
-  computePushForwardDerivative(T2toST2ResultType&,
-                               const T2toST2Type&,
-                               const StensorType&,
-                               const TensorType&);
+  TFEL_HOST_DEVICE constexpr void computePushForwardDerivative(
+      T2toST2ResultType&,
+      const T2toST2Type&,
+      const StensorType&,
+      const TensorType&) noexcept  //
+      requires(tfel::typetraits::IsFundamentalNumericType<
+               numeric_type<TensorType>>::cond&&
+                   isAssignableTo<
+                       typename ComputeBinaryResult<numeric_type<T2toST2Type>,
+                                                    numeric_type<StensorType>,
+                                                    OpPlus>::Result,
+                       numeric_type<T2toST2ResultType>>());
 
 }  // end of namespace tfel::math
 
diff --git a/include/TFEL/Math/t2tot2.hxx b/include/TFEL/Math/t2tot2.hxx
index a060dfefa..514a05875 100644
--- a/include/TFEL/Math/t2tot2.hxx
+++ b/include/TFEL/Math/t2tot2.hxx
@@ -184,20 +184,20 @@ namespace tfel::math {
      * \brief constructor from a t2tost2
      * \param[in] v : values
      */
-    template <typename T2toST2Type,
-              std::enable_if_t<
-                  ((implementsT2toST2Concept<T2toST2Type>()) &&
-                   (isAssignableTo<numeric_type<T2toST2Type>, ValueType>()) &&
-                   (getSpaceDimension<T2toST2Type>() == N)),
-                  bool> = true>
-    TFEL_HOST_DEVICE constexpr t2tot2(const T2toST2Type&);
+    template <T2toST2Concept T2toST2Type>
+    TFEL_HOST_DEVICE constexpr t2tot2(const T2toST2Type&) noexcept requires(
+        (isAssignableTo<numeric_type<T2toST2Type>, ValueType>()) &&
+        (getSpaceDimension<T2toST2Type>() == N));
     //
     TFEL_MATH_FIXED_SIZE_ARRAY_DEFAULT_METHODS(t2tot2,
                                                GenericFixedSizeArrayBase);
     //! exposing GenericFixedSizeArray access operator
     using GenericFixedSizeArrayBase::operator[];
     using GenericFixedSizeArrayBase::operator();
-
+    //! \brief import values from an external memory location
+    TFEL_HOST_DEVICE constexpr void import(
+        const base_type<ValueType>* const) noexcept;
+    //! \brief copy from a range
     TFEL_HOST_DEVICE constexpr void copy(const auto) noexcept;
   };
 
@@ -225,7 +225,6 @@ namespace tfel::math {
   TFEL_HOST_DEVICE constexpr auto change_basis(
       const T2toT2Type&,
       const rotation_matrix<numeric_type<T2toT2Type>>&) noexcept;
-
   /*!
    * \return compute the derivative of the velocity gradient
    * \param[in] F: deformation gradient
@@ -248,12 +247,11 @@ namespace tfel::math {
   computeDeterminantSecondDerivative(const TensorType&) noexcept requires(
       isScalar<numeric_type<TensorType>>());
 
-  template <typename T, typename T2toST2Type>
-  TFEL_HOST_DEVICE
-      std::enable_if_t<((implementsT2toST2Concept<T2toST2Type>()) &&
-                        (isAssignableTo<numeric_type<T2toST2Type>, T>())),
-                       void>
-      convert(t2tot2<getSpaceDimension<T2toST2Type>(), T>&, const T2toST2Type&);
+  template <typename T, T2toST2Concept T2toST2Type>
+  TFEL_HOST_DEVICE constexpr void convert(
+      t2tot2<getSpaceDimension<T2toST2Type>(), T>&,
+      const T2toST2Type&) noexcept  //
+      requires(isAssignableTo<numeric_type<T2toST2Type>, T>());
 
 }  // end of namespace tfel::math
 
diff --git a/include/TFEL/Math/tensor.hxx b/include/TFEL/Math/tensor.hxx
index 25f0fc30a..912185579 100644
--- a/include/TFEL/Math/tensor.hxx
+++ b/include/TFEL/Math/tensor.hxx
@@ -158,9 +158,9 @@ namespace tfel::math {
   template <unsigned short N, typename T>
   using ConstTensorView = ConstView<tensor<N, T>>;
 
-  template <unsigned short N, typename T, typename OutputIterator>
-  TFEL_HOST_DEVICE TFEL_MATH_INLINE2 std::enable_if_t<isScalar<T>(), void>
-  exportToBaseTypeArray(const tensor<N, T>&, OutputIterator);
+  template <unsigned short N, typename T>
+  TFEL_HOST_DEVICE constexpr void exportToBaseTypeArray(
+      const tensor<N, T>&, const auto) noexcept requires(isScalar<T>());
   /*!
    * \return the invert of a tensor
    * \param[in] t : tensor to be inverted