diff --git a/include/expolit/differentiation.hh b/include/expolit/differentiation.hh
index 94e8e4b..a1a496c 100644
--- a/include/expolit/differentiation.hh
+++ b/include/expolit/differentiation.hh
@@ -1,143 +1,143 @@
 /**
  * @file
  *
  * @author Lucas Frérot <lucas.frerot@epfl.ch>
  *
  * @section LICENSE
  *
- * Copyright (©)  2018-2020 EPFL  (Ecole Polytechnique  Fédérale de
+ * Copyright (©)  2018-2021 EPFL  (Ecole Polytechnique  Fédérale de
  * Lausanne)  Laboratory (LSMS  -  Laboratoire de  Simulation  en Mécanique  des
  * Solides)
  *
  * Expolit is free software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as  published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Expolit is  distributed in the  hope that it  will be useful, but  WITHOUT
  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
  * FOR A PARTICULAR PURPOSE. See  the GNU  Lesser General  Public License for
  * more details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Expolit. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 
 #ifndef DIFFERENTIATION_HH
 #define DIFFERENTIATION_HH
 
 #include "exponential.hh"
 #include "litteral.hh"
 #include "operations.hh"
 #include "polynomial.hh"
 #include "types.hh"
 #include <algorithm>
 
 namespace expolit {
 
 namespace detail {
 template <std::size_t i>
 struct for_diff {
   template <typename CT, std::size_t n>
   constexpr static void diff(std::array<CT, n - 1>& res,
                              const std::array<CT, n>& a) {
     std::get<i>(res) = a[i + 1] * static_cast<CT>(i + 1);
     for_diff<i - 1>::diff(res, a);
   }
 };
 
 template <>
 struct for_diff<0> {
   template <typename CT, std::size_t n>
   constexpr static void diff(std::array<CT, n - 1>& res,
                              const std::array<CT, n>& a) {
     std::get<0>(res) = a[1] * static_cast<CT>(1);
   }
 };
 }  // namespace detail
 
 /// Differentiate polynomial
 template <typename CT, UInt order>
 constexpr auto differentiate(const Polynomial<CT, order>& p) {
   Polynomial<CT, order - 1> dp;
   detail::for_diff<dp.terms - 1>::diff(dp.coeffs, p.coeffs);
   return dp;
 }
 
 /// Differentiate constant
 template <typename CT>
 constexpr auto differentiate(const Constant<CT>& /*p*/) {
   return Constant<CT>({0});
 }
 
 /// Differentiate litteral
 template <typename Tag>
 constexpr auto differentiate(const Litteral<Tag>& /*p*/) {
   return 0;
 }
 
 /// Differentiate exponential
 template <typename T>
 constexpr auto differentiate(const Exponential<T>& e) {
   return differentiate(e.expression) * e;
 }
 
 /// Differentiate generic product
 template <typename T, typename U>
 constexpr auto differentiate(const Product<T, U>& p) {
   auto&& o = p.operands;
   return differentiate(o.first) * o.second + o.first * differentiate(o.second);
 }
 
 /// Differentiate constant + expression
 template <typename CT, typename U>
 constexpr auto differentiate(const Sum<Constant<CT>, U>& p) {
   return differentiate(p.operands.second);
 }
 
 /// Differentiate expression + constant
 template <typename T, typename U, typename CT>
 constexpr auto differentiate(const Sum<U, Constant<CT>>& p) {
   return differentiate(p.commute());
 }
 
 /// Differentiate litteral + expression
 template <typename Tag, typename U>
 constexpr auto differentiate(const Sum<Litteral<Tag>, U>& p) {
   return differentiate(p.operands.second);
 }
 
 /// Differentiate expression + litteral
 template <typename Tag, typename U>
 constexpr auto differentiate(const Sum<U, Litteral<Tag>>& p) {
   return differentiate(p.commute());
 }
 
 /// Differentiate constant * expression
 template <typename CT, typename U>
 constexpr auto differentiate(const Product<Constant<CT>, U>& p) {
   return p.operands.first * differentiate(p.operands.second);
 }
 
 /// Differentiate expression * constant
 template <typename CT, typename U>
 constexpr auto differentiate(const Product<U, Constant<CT>>& p) {
   return differentiate(p.commute());
 }
 
 /// Differentiate litteral * expression
 template <typename Tag, typename U>
 constexpr auto differentiate(const Product<Litteral<Tag>, U>& p) {
   return p.operands.first * differentiate(p.operands.second);
 }
 
 /// Differentiate expression * litteral
 template <typename Tag, typename U>
 constexpr auto differentiate(const Product<U, Litteral<Tag>>& p) {
   return differentiate(p.commute());
 }
 
 }  // namespace expolit
 
 #endif
diff --git a/include/expolit/expolit b/include/expolit/expolit
index 1e7cb15..d9511bd 100644
--- a/include/expolit/expolit
+++ b/include/expolit/expolit
@@ -1,38 +1,38 @@
 /**
  * @file
  *
  * @author Lucas Frérot <lucas.frerot@epfl.ch>
  *
  * @section LICENSE
  *
- * Copyright (©)  2018-2020 EPFL  (Ecole Polytechnique  Fédérale de
+ * Copyright (©)  2018-2021 EPFL  (Ecole Polytechnique  Fédérale de
  * Lausanne)  Laboratory (LSMS  -  Laboratoire de  Simulation  en Mécanique  des
  * Solides)
  *
  * Expolit is free software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as  published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Expolit is  distributed in the  hope that it  will be useful, but  WITHOUT
  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
  * FOR A PARTICULAR PURPOSE. See  the GNU  Lesser General  Public License for
  * more details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Expolit. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 
 #ifndef EXPOLIT_HEADER
 #define EXPOLIT_HEADER
 
 #include "types.hh"
 #include "operations.hh"
 #include "polynomial.hh"
 #include "exponential.hh"
 #include "litteral.hh"
 #include "differentiation.hh"
 #include "integration.hh"
 
 #endif
diff --git a/include/expolit/exponential.hh b/include/expolit/exponential.hh
index 3047a72..b600ac5 100644
--- a/include/expolit/exponential.hh
+++ b/include/expolit/exponential.hh
@@ -1,74 +1,74 @@
 /**
  * @file
  *
  * @author Lucas Frérot <lucas.frerot@epfl.ch>
  *
  * @section LICENSE
  *
- * Copyright (©)  2018-2020 EPFL  (Ecole Polytechnique  Fédérale de
+ * Copyright (©)  2018-2021 EPFL  (Ecole Polytechnique  Fédérale de
  * Lausanne)  Laboratory (LSMS  -  Laboratoire de  Simulation  en Mécanique  des
  * Solides)
  *
  * Expolit is free software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as  published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Expolit is  distributed in the  hope that it  will be useful, but  WITHOUT
  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
  * FOR A PARTICULAR PURPOSE. See  the GNU  Lesser General  Public License for
  * more details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Expolit. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 
 #ifndef EXPONENTIAL_HH
 #define EXPONENTIAL_HH
 
 #include "types.hh"
 #include "operations.hh"
 #include <cmath>
 #include <ostream>
 #include <iostream>
 
 namespace expolit {
 
 namespace detail {
 // Necessary for proper dispatch
 using std::exp;
 
 struct ExponentialFunctor {
   template <typename Val>
   constexpr auto operator()(Val&& v) const {
     return exp(v);  // Should dispatch in case expression is passed
   }
 };
 }
 
 /**
  * @brief Exponential expression
  *
  * Handles exponential of arbitrary expression
  */
 template <typename T>
 using Exponential = UnaryOperation<T, detail::ExponentialFunctor>;
 
 template <typename Derived>
 constexpr auto exp(const Expression<Derived>& e) {
   return Exponential<std::remove_cv_t<Derived>>(e.downcast());
 }
 
 /* -------------------------------------------------------------------------- */
 /* Output */
 /* -------------------------------------------------------------------------- */
 template <typename T>
 std::ostream& operator<<(std::ostream& os, const Exponential<T>& e) {
   os << "exp(" << e.expression << ")";
   return os;
 }
 
 }  // namespace expolit
 
 #endif
diff --git a/include/expolit/integration.hh b/include/expolit/integration.hh
index 04ad777..ba61b4f 100644
--- a/include/expolit/integration.hh
+++ b/include/expolit/integration.hh
@@ -1,155 +1,155 @@
 /**
  * @file
  *
  * @author Lucas Frérot <lucas.frerot@epfl.ch>
  *
  * @section LICENSE
  *
- * Copyright (©)  2018-2020 EPFL  (Ecole Polytechnique  Fédérale de
+ * Copyright (©)  2018-2021 EPFL  (Ecole Polytechnique  Fédérale de
  * Lausanne)  Laboratory (LSMS  -  Laboratoire de  Simulation  en Mécanique  des
  * Solides)
  *
  * Expolit is free software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as  published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Expolit is  distributed in the  hope that it  will be useful, but  WITHOUT
  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
  * FOR A PARTICULAR PURPOSE. See  the GNU  Lesser General  Public License for
  * more details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Expolit. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 
 #ifndef INTEGRATION_HH
 #define INTEGRATION_HH
 
 #include "differentiation.hh"
 #include "exponential.hh"
 #include "litteral.hh"
 #include "polynomial.hh"
 #include "types.hh"
 
 namespace expolit {
 
 namespace detail {
 template <std::size_t i>
 struct for_integ {
   template <typename CT, std::size_t n>
   static constexpr void integ(std::array<CT, n+1>& res,
 			      const std::array<CT, n>& a) {
     std::get<i>(res) = a[i-1] / static_cast<CT>(i);
     for_integ<i-1>::integ(res, a);
   }
 };
 
 template <>
 struct for_integ<1> {
   template <typename CT, std::size_t n>
   static constexpr void integ(std::array<CT, n+1>& res,
 			      const std::array<CT, n>& a) {
     std::get<1>(res) = a[0];
   }
 };
 }  // namespace detail
 
 /// Integrate polynomial
 template <typename CT, UInt order>
 constexpr auto integrate(const Polynomial<CT, order>& p) {
   Polynomial<CT, order + 1> ip;
   detail::for_integ<ip.terms-1>::integ(ip.coeffs, p.coeffs);
   return ip;
 }
 
 /// Integrate litteral
 template <typename Tag>
 constexpr auto integrate(const Litteral<Tag>& q) {
   constexpr Polynomial<Int, 1> id({0, 1});
   return q * id;
 }
 
 /// Integrate constant * expression
 template <typename CT, typename Derived>
 constexpr auto integrate(const Product<Constant<CT>, Derived>& p) {
   return p.operands.first * integrate(p.operands.second);
 }
 
 /// Integrate expression * constant
 template <typename CT, typename Derived>
 constexpr auto integrate(const Product<Derived, Constant<CT>>& p) {
   return integrate(p.commute());
 }
 
 /// Integrate litteral * expression
 template <typename Tag, typename Derived>
 constexpr auto integrate(const Product<Litteral<Tag>, Derived>& p) {
   return p.operands.first * integrate(p.operands.second);
 }
 
 /// Integrate expression * litteral
 template <typename Tag, typename Derived>
 constexpr auto integrate(const Product<Derived, Litteral<Tag>>& p) {
   return integrate(p.commute());
 }
 
 /// Integrate litteral / expression
 template <typename Tag, typename Derived>
 constexpr auto integrate(const Division<Derived, Litteral<Tag>>& p) {
   return integrate(p.operands.first) / p.operands.second;
 }
 
 
 /// Integrate exponential of order 1 polynomial
 template <typename CT>
 constexpr auto integrate(const Exponential<Polynomial<CT, 1>>& e) {
   return Constant<CT>({1 / e.expression.coeffs.back()}) * e;
 }
 
 /// Integrate exponential of litteral * order 1 polynomial
 template <typename Tag, typename CT>
 constexpr auto
 integrate(const Exponential<Product<Litteral<Tag>, Polynomial<CT, 1>>>& e) {
   return e * Constant<CT>({1 / e.expression.operands.second.coeffs.back()}) /
          e.expression.operands.first;
 }
 
 /// Integrate exponential of order 1 polynomial * litteral
 template <typename Tag, typename CT>
 constexpr auto
 integrate(const Exponential<Product<Polynomial<CT, 1>, Litteral<Tag>>>& e) {
   return integrate(e.op(e.expression.commute()));
 }
 
 /// Integration by parts of expression * polynomial
 template <typename Derived, typename CT, UInt order>
 constexpr auto integrate(const Product<Derived, Polynomial<CT, order>>& p) {
   auto u = integrate(p.operands.first);
   return u * p.operands.second -
          integrate(u * differentiate(p.operands.second));
 }
 
 /// Integration by parts of polynomial * expression
 template <typename Derived, typename CT, UInt order>
 constexpr auto integrate(const Product<Polynomial<CT, order>, Derived>& p) {
   return integrate(p.commute());
 }
 
 /// Integration of sum
 template <typename Der1, typename Der2>
 constexpr auto integrate(const Sum<Der1, Der2>& s) {
   return integrate(s.operands.first) + integrate(s.operands.second);
 }
 
 /// Compute definite integral
 template <typename T, typename Derived>
 constexpr auto definite_integral(const std::pair<T, T>& bounds,
                                  const Expression<Derived>& exp) {
   auto integ = integrate(exp.downcast());
   return integ(bounds.second) - integ(bounds.first);
 }
 
 }  // namespace expolit
 #endif
diff --git a/include/expolit/litteral.hh b/include/expolit/litteral.hh
index 67073b0..096c02c 100644
--- a/include/expolit/litteral.hh
+++ b/include/expolit/litteral.hh
@@ -1,80 +1,80 @@
 /**
  * @file
  *
  * @author Lucas Frérot <lucas.frerot@epfl.ch>
  *
  * @section LICENSE
  *
- * Copyright (©)  2018-2020 EPFL  (Ecole Polytechnique  Fédérale de
+ * Copyright (©)  2018-2021 EPFL  (Ecole Polytechnique  Fédérale de
  * Lausanne)  Laboratory (LSMS  -  Laboratoire de  Simulation  en Mécanique  des
  * Solides)
  *
  * Expolit is free software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as  published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Expolit is  distributed in the  hope that it  will be useful, but  WITHOUT
  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
  * FOR A PARTICULAR PURPOSE. See  the GNU  Lesser General  Public License for
  * more details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Expolit. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 
 #ifndef LITTERAL_HH
 #define LITTERAL_HH
 
 #include "operations.hh"
 #include "polynomial.hh"
 #include "types.hh"
 
 namespace expolit {
 template <typename Tag>
 struct Litteral : Expression<Litteral<Tag>> {
   using tag = Tag;
   Litteral() = default;
 };
 
 /// Substitue litteral
 template <typename Tag, typename New>
 constexpr auto substitute(const Expression<New>& val,
 			  const Litteral<Tag>& /*t*/) {
   return val.downcast();
 }
 
 /// Substitute litteral in binary operation
 template <typename Tag, typename New, typename Der1, typename Der2,
           typename Op>
 constexpr auto substitute(const Expression<New>& val,
                           const BinaryOperation<Der1, Der2, Op>& op) {
   return op.op(substitute<Tag>(val, op.operands.first),
                substitute<Tag>(val, op.operands.second));
 }
 
 /// Substitute litteral in polynomial (aka do nothing)
 template <typename Tag, typename New, typename CT, UInt order>
 constexpr auto substitute(const Expression<New>& /*val*/,
 			  const Polynomial<CT, order>& p) {
   return p;
 }
 
 /// Substitute litteral in polynomial (aka do nothing)
 template <typename Tag, typename New, typename T, typename Op>
 constexpr auto substitute(const Expression<New>& val,
 			  const UnaryOperation<T, Op>& p) {
   return p.op(substitute<Tag>(val, p.expression));
 }
 
 /// Substitute generic expression by downcast
 // template <typename Tag, typename New, typename Derived>
 // constexpr auto substitute(const Expression<New>& val,
 // 			  const Expression<Derived>& e) {
 //   return substitute<Tag>(val, e.downcast());
 // }
 
 }  // namespace expolit
 
 #endif
diff --git a/include/expolit/operations.hh b/include/expolit/operations.hh
index 38954d3..78d005e 100644
--- a/include/expolit/operations.hh
+++ b/include/expolit/operations.hh
@@ -1,189 +1,189 @@
 /**
  * @file
  *
  * @author Lucas Frérot <lucas.frerot@epfl.ch>
  *
  * @section LICENSE
  *
- * Copyright (©)  2018-2020 EPFL  (Ecole Polytechnique  Fédérale de
+ * Copyright (©)  2018-2021 EPFL  (Ecole Polytechnique  Fédérale de
  * Lausanne)  Laboratory (LSMS  -  Laboratoire de  Simulation  en Mécanique  des
  * Solides)
  *
  * Expolit is free software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as  published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Expolit is  distributed in the  hope that it  will be useful, but  WITHOUT
  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
  * FOR A PARTICULAR PURPOSE. See  the GNU  Lesser General  Public License for
  * more details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Expolit. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 
 #ifndef OPERATIONS_HH
 #define OPERATIONS_HH
 
 #include "types.hh"
 #include <functional>
 #include <ostream>
 #include <utility>
 
 namespace expolit {
 /**
  * @brief Type for unary operation
  */
 template <typename T, typename Op>
 struct UnaryOperation : Expression<UnaryOperation<T, Op>> {
   /// Construct from an expression
   constexpr UnaryOperation(const Expression<T>& t) : expression(t.downcast()) {}
   /// Evaluate operation
   template <typename Val>
   constexpr auto operator()(Val&& v) const {
     return op(expression(v));
   }
 
   const T expression;
   const Op op = Op();
 };
 
 /**
  * @brief Type for binary operation
  *
  * Handles binary operations between two arbitrary expressions
  */
 template <typename T, typename U, typename Op>
 struct BinaryOperation : Expression<BinaryOperation<T, U, Op>> {
   /// Construct from two expressions
   constexpr BinaryOperation(const Expression<T>& t, const Expression<U>& u)
       : operands(t.downcast(), u.downcast()) {}
 
   /// Evaluate operation
   template <typename Val>
   constexpr auto operator()(Val&& v) const {
     return op(operands.first(v), operands.second(v));
   }
 
   constexpr auto commute() const {
     return BinaryOperation<U, T, Op>(operands.second, operands.first);
   }
 
   /// Product operarands
   const std::pair<T, U> operands;
   const Op op = Op();
 };
 
 /// Symmetric of T, U
 template <typename T, typename U, typename Op>
 using SBinaryOperation = BinaryOperation<U, T, Op>;
 
 /// Product expression
 template <typename T, typename U>
 using Product = BinaryOperation<T, U, std::multiplies<>>;
 
 /// Sum expression
 template <typename T, typename U>
 using Sum = BinaryOperation<T, U, std::plus<>>;
 
 /// Difference expression
 template <typename T, typename U>
 using Difference = BinaryOperation<T, U, std::minus<>>;
 
 /// Division expression
 template <typename T, typename U>
 using Division = BinaryOperation<T, U, std::divides<>>;
 
 /* -------------------------------------------------------------------------- */
 /* Algebraic rules */
 /* -------------------------------------------------------------------------- */
 /// Product operator
 template <typename Der1, typename Der2>
 constexpr auto operator*(const Expression<Der1>& e1,
                          const Expression<Der2>& e2) {
   return Product<Der1, Der2>(e1, e2);
 }
 
 /// Sum operator
 template <typename Der1, typename Der2>
 constexpr auto operator+(const Expression<Der1>& e1,
                          const Expression<Der2>& e2) {
   return Sum<Der1, Der2>(e1, e2);
 }
 
 /// Difference operator
 template <typename Der1, typename Der2>
 constexpr auto operator-(const Expression<Der1>& e1,
                          const Expression<Der2>& e2) {
   return Difference<Der1, Der2>(e1, e2);
 }
 
 /// Division operator
 template <typename Der1, typename Der2>
 constexpr auto operator/(const Expression<Der1>& e1,
                          const Expression<Der2>& e2) {
   return Division<Der1, Der2>(e1, e2);
 }
 
 namespace detail {
 template <UInt n>
 struct power_acc {
   template <typename DerAcc, typename Derived>
   static constexpr auto accumulate(const Expression<DerAcc>& acc,
                                    const Expression<Derived>& e) {
     return power_acc<n - 1>::accumulate(acc.downcast() * e.downcast(), e);
   }
 };
 
 template <>
 struct power_acc<1> {
   template <typename DerAcc, typename Derived>
   static constexpr auto accumulate(const Expression<DerAcc>& acc,
                                    const Expression<Derived>&) {
     return acc.downcast();
   }
 };
 
 }  // namespace detail
 
 template <UInt n, typename Derived>
 constexpr auto pow(const Expression<Derived>& e) {
   return detail::power_acc<n>::accumulate(e, e);
 }
 
 /* -------------------------------------------------------------------------- */
 /* Output */
 /* -------------------------------------------------------------------------- */
 
 template <typename T, typename U>
 std::ostream& operator<<(std::ostream& os, const Product<T, U>& p) {
   os << "(" << p.operands.first << ")"
      << " * "
      << "(" << p.operands.second << ")";
   return os;
 }
 
 template <typename T, typename U>
 std::ostream& operator<<(std::ostream& os, const Division<T, U>& p) {
   os << "(" << p.operands.first << ")"
      << " / "
      << "(" << p.operands.second << ")";
   return os;
 }
 
 template <typename T, typename U>
 std::ostream& operator<<(std::ostream& os, const Difference<T, U>& p) {
   os << p.operands.first << " - " << p.operands.second;
   return os;
 }
 
 template <typename T, typename U>
 std::ostream& operator<<(std::ostream& os, const Sum<T, U>& p) {
   os << p.operands.first << " + " << p.operands.second;
   return os;
 }
 
 }  // namespace expolit
 #endif
diff --git a/include/expolit/polynomial.hh b/include/expolit/polynomial.hh
index c42525c..fd002b9 100644
--- a/include/expolit/polynomial.hh
+++ b/include/expolit/polynomial.hh
@@ -1,327 +1,327 @@
 /**
  * @file
  *
  * @author Lucas Frérot <lucas.frerot@epfl.ch>
  *
  * @section LICENSE
  *
- * Copyright (©)  2018-2020 EPFL  (Ecole Polytechnique  Fédérale de
+ * Copyright (©)  2018-2021 EPFL  (Ecole Polytechnique  Fédérale de
  * Lausanne)  Laboratory (LSMS  -  Laboratoire de  Simulation  en Mécanique  des
  * Solides)
  *
  * Expolit is free software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as  published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Expolit is  distributed in the  hope that it  will be useful, but  WITHOUT
  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
  * FOR A PARTICULAR PURPOSE. See  the GNU  Lesser General  Public License for
  * more details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Expolit. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 
 #ifndef POLYNOMIALS_HH
 #define POLYNOMIALS_HH
 
 #include "operations.hh"
 #include "types.hh"
 #include <algorithm>
 #include <array>
 #include <numeric>
 #include <ostream>
 
 namespace expolit {
 
 template <typename From, typename To>
 using compatible = std::enable_if_t<std::is_convertible<From, To>::value>;
 
 /**
  * @brief Polynomial type
  *
  * Type for polynomial expressions of arbitrary order, variable type and
  * coefficient type
  */
 template <typename CoeffType, UInt order>
 struct Polynomial : Expression<Polynomial<CoeffType, order>> {
   /// Number of terms
   constexpr static std::size_t terms = order + 1;
 
   constexpr Polynomial() = default;
   constexpr explicit Polynomial(const std::array<CoeffType, terms>& coeffs)
       : coeffs(coeffs) {}
 
   /// Evaluating polynomial
   template <typename Val, typename = compatible<Val, CoeffType>>
   constexpr auto operator()(Val z) const {
     using Res = decltype(std::declval<CoeffType>() * z);
     Res res(coeffs.front());
     Val power(z);
     for (UInt i = 1; i < terms; ++i) {
       res += coeffs[i] * power;
       power *= z;
     }
     return res;
   }
 
   /// Storing coefficients
   std::array<CoeffType, terms> coeffs = {0};
 };
 
 template <typename CoeffType>
 using Constant = Polynomial<CoeffType, 0>;
 
 namespace detail {
 template <typename Poly>
 struct get_terms
     : std::integral_constant<
           std::size_t, std::remove_reference_t<std::remove_cv_t<Poly>>::terms> {
 };
 }  // namespace detail
 
 /* -------------------------------------------------------------------------- */
 /* Algebraic rules */
 /* -------------------------------------------------------------------------- */
 
 /// Sum with scalar
 template <typename CT, UInt order, typename DT, typename = compatible<DT, CT>>
 constexpr auto operator+(DT x, const Polynomial<CT, order>& p) {
   Polynomial<decltype(std::declval<CT>() + x), order> res(p);
   std::get<0>(res.coeffs) += x;
   return res;
 }
 
 /// Sum with scalar
 template <typename CT, UInt order, typename DT, typename = compatible<DT, CT>>
 constexpr auto operator+(const Polynomial<CT, order>& p, DT x) {
   return x + p;
 }
 
 namespace detail {
 template <std::size_t i>
 struct for_scalar {
   template <typename CT, std::size_t n, typename DT>
   constexpr static void prod(std::array<CT, n>& res, DT x) {
     std::get<i>(res) *= x;
     for_scalar<i - 1>::prod(res, x);
   }
 
   template <typename CT, std::size_t n, typename DT>
   constexpr static void div(std::array<CT, n>& res, DT x) {
     std::get<i>(res) /= x;
     for_scalar<i - 1>::prod(res, x);
   }
 };
 
 template <>
 struct for_scalar<0> {
   template <typename CT, std::size_t n, typename DT>
   constexpr static void prod(std::array<CT, n>& res, DT x) {
     std::get<0>(res) *= x;
   }
 
   template <typename CT, std::size_t n, typename DT>
   constexpr static void div(std::array<CT, n>& res, DT x) {
     std::get<0>(res) /= x;
   }
 };
 }  // namespace detail
 
 /// Product with scalar
 template <typename CT, UInt order, typename DT, typename = compatible<DT, CT>>
 constexpr auto operator*(DT x, const Polynomial<CT, order>& p) {
   Polynomial<decltype(std::declval<CT>() * x), order> res(p);
   detail::for_scalar<res.terms - 1>::prod(res.coeffs, x);
   return res;
 }
 
 /// Product with scalar
 template <typename CT, UInt order, typename DT, typename = compatible<DT, CT>>
 constexpr auto operator*(const Polynomial<CT, order>& p, DT x) {
   return x * p;
 }
 
 /// Division with scalar
 template <typename CT, UInt order, typename DT, typename = compatible<DT, CT>>
 constexpr auto operator/(const Polynomial<CT, order>& p, DT x) {
   Polynomial<decltype(std::declval<CT>() / x), order> res(p);
   detail::for_scalar<res.terms - 1>::div(res.coeffs, x);
   return res;
 }
 
 /// Division with constant
 template <typename CT, UInt order, typename CT2, typename = compatible<CT2, CT>>
 constexpr auto operator/(const Polynomial<CT, order>& p,
                          const Constant<CT2>& x) {
   Polynomial<decltype(std::declval<CT>() / x.coeffs.front()), order> res(p);
   detail::for_scalar<res.terms - 1>::div(res.coeffs, std::get<0>(x.coeffs));
   return res;
 }
 
 namespace detail {
 template <std::size_t i, std::size_t j>
 struct coeff_product {
   template <typename CT1, typename CT2, typename CT3, std::size_t n,
             std::size_t m, std::size_t l>
   static constexpr void prod(std::array<CT1, n>& res,
                              const std::array<CT2, m>& a,
                              const std::array<CT3, l>& b) {
     std::get<i + j>(res) += a[i] * b[j];
     coeff_product<i, j - 1>::prod(res, a, b);
   }
 };
 
 template <std::size_t i>
 struct coeff_product<i, 0> {
   template <typename CT1, typename CT2, typename CT3, std::size_t n,
             std::size_t m, std::size_t l>
   static constexpr void prod(std::array<CT1, n>& res,
                              const std::array<CT2, m>& a,
                              const std::array<CT3, l>& b) {
     std::get<i>(res) += a[i] * b[0];
     coeff_product<i - 1, l - 1>::prod(res, a, b);
   }
 };
 
 template <>
 struct coeff_product<0, 0> {
   template <typename CT1, typename CT2, typename CT3, std::size_t n,
             std::size_t m, std::size_t l>
   static constexpr void prod(std::array<CT1, n>& res,
                              const std::array<CT2, m>& a,
                              const std::array<CT3, l>& b) {
     std::get<0>(res) += a[0] * b[0];
   }
 };
 }  // namespace detail
 
 /// Product of polynomials
 template <typename CT1, typename CT2, UInt o1, UInt o2>
 constexpr auto operator*(const Polynomial<CT1, o1>& p1,
                          const Polynomial<CT2, o2>& p2) {
   Polynomial<decltype(std::declval<CT1>() * std::declval<CT2>()), o1 + o2> res;
   detail::coeff_product<detail::get_terms<decltype(p1)>::value - 1,
                         detail::get_terms<decltype(p2)>::value -
                             1>::prod(res.coeffs, p1.coeffs, p2.coeffs);
   return res;
 }
 
 namespace detail {
 template <std::size_t i>
 struct coeff_sum {
   template <typename CT1, typename CT2, typename CT3, std::size_t n,
             std::size_t m, std::size_t l>
   static constexpr void sum(std::array<CT1, n>& res,
                             const std::array<CT2, m>& a,
                             const std::array<CT3, l>& b) {
     std::get<i>(res) += ((i < m) ? a[i] : 0) + ((i < l) ? b[i] : 0);
     coeff_sum<i - 1>::sum(res, a, b);
   }
 };
 
 template <>
 struct coeff_sum<0> {
   template <typename CT1, typename CT2, typename CT3, std::size_t n,
             std::size_t m, std::size_t l>
   static constexpr void sum(std::array<CT1, n>& res,
                             const std::array<CT2, m>& a,
                             const std::array<CT3, l>& b) {
     std::get<0>(res) += a[0] + b[0];
   }
 };
 }  // namespace detail
 
 /// Sum of polynomials
 template <typename CT1, typename CT2, UInt o1, UInt o2>
 constexpr auto operator+(const Polynomial<CT1, o1>& p1,
                          const Polynomial<CT2, o2>& p2) {
   Polynomial<decltype(std::declval<CT1>() + std::declval<CT2>()),
              std::max(o1, o2)>
       res;
   detail::coeff_sum<res.terms - 1>::sum(res.coeffs, p1.coeffs, p2.coeffs);
   return res;
 }
 
 /* -------------------------------------------------------------------------- */
 /* Simplification rules */
 /* -------------------------------------------------------------------------- */
 
 #define POLY_SIMPLIFICATION(opsign, std_op)                                    \
   template <typename CT1, typename CT2, UInt o1, UInt o2, typename Derived>    \
   constexpr auto operator opsign(                                              \
       const Polynomial<CT1, o1>& p1,                                           \
       const BinaryOperation<Polynomial<CT2, o2>, Derived, std_op>& p2) {       \
     return p2.op(p2.op(p1, p2.operands.first), p2.operands.second);            \
   }                                                                            \
                                                                                \
   template <typename CT1, typename CT2, UInt o1, UInt o2, typename Derived>    \
   constexpr auto operator opsign(                                              \
       const BinaryOperation<Polynomial<CT2, o2>, Derived, std_op>& p2,         \
       const Polynomial<CT1, o1>& p1) {                                         \
     return p2.op(p1, p2);                                                      \
   }                                                                            \
                                                                                \
   template <typename CT1, typename CT2, UInt o1, UInt o2, typename Derived>    \
   constexpr auto operator opsign(                                              \
       const Polynomial<CT1, o1>& p1,                                           \
       const SBinaryOperation<Polynomial<CT2, o2>, Derived, std_op>& p2) {      \
     return p2.op(p1, p2.commute());                                            \
   }                                                                            \
                                                                                \
   template <typename CT1, typename CT2, UInt o1, UInt o2, typename Derived>    \
   constexpr auto operator opsign(                                              \
       const SBinaryOperation<Polynomial<CT2, o2>, Derived, std_op>& p2,        \
       const Polynomial<CT1, o1>& p1) {                                         \
     return p2.op(p1, p2.commute());                                            \
   }
 
 POLY_SIMPLIFICATION(+, std::plus<>)
 POLY_SIMPLIFICATION(*, std::multiplies<>)
 
 #undef POLY_SIMPLIFICATION
 
 template <typename CT1, UInt o, typename CT2, typename Derived>
 constexpr auto operator/(const Product<Polynomial<CT1, o>, Derived>& p,
                          const Constant<CT2>& c) {
   return p.operands.second * (p.operands.first / c);
 }
 
 template <typename CT1, UInt o, typename CT2, typename Derived>
 constexpr auto operator/(const Product<Derived, Polynomial<CT1, o>>& p,
                          const Constant<CT2>& c) {
   return p.operands.second * (p.operands.first / c);
 }
 
 template <typename CT1, UInt o, typename CT2, typename Derived>
 constexpr auto operator*(const Polynomial<CT1, o>& p,
                          const Division<Derived, Constant<CT2>>& d) {
   return (p / d.operands.second) * d.operands.first;
 }
 
 template <typename CT1, UInt o, typename CT2, typename Derived>
 constexpr auto operator*(const Division<Derived, Constant<CT2>>& d,
                          const Polynomial<CT1, o>& p) {
   return (p / d.operands.second) * d.operands.first;
 }
 
 /* -------------------------------------------------------------------------- */
 /* Output */
 /* -------------------------------------------------------------------------- */
 
 template <typename CT, UInt order>
 std::ostream& operator<<(std::ostream& os, const Polynomial<CT, order>& p) {
   os << p.coeffs[0];
 
   if (p.terms >= 2) {
     os << " + " << p.coeffs[1] << "*z";
     for (UInt i = 2; i < p.terms; ++i)
       os << " + " << p.coeffs[i] << "*z^" << i;
   }
   return os;
 }
 
 }  // namespace expolit
 #endif
diff --git a/include/expolit/types.hh b/include/expolit/types.hh
index 68b998a..729d0f3 100644
--- a/include/expolit/types.hh
+++ b/include/expolit/types.hh
@@ -1,52 +1,52 @@
 /**
  * @file
  *
  * @author Lucas Frérot <lucas.frerot@epfl.ch>
  *
  * @section LICENSE
  *
- * Copyright (©)  2018-2020 EPFL  (Ecole Polytechnique  Fédérale de
+ * Copyright (©)  2018-2021 EPFL  (Ecole Polytechnique  Fédérale de
  * Lausanne)  Laboratory (LSMS  -  Laboratoire de  Simulation  en Mécanique  des
  * Solides)
  *
  * Expolit is free software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as  published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Expolit is  distributed in the  hope that it  will be useful, but  WITHOUT
  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
  * FOR A PARTICULAR PURPOSE. See  the GNU  Lesser General  Public License for
  * more details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Expolit. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 
 #ifndef TYPES_HH
 #define TYPES_HH
 
 #include <cstddef>
 #include <ostream>
 #include <type_traits>
 #include <utility>
 
 namespace expolit {
 using Int = int;
 using UInt = unsigned int;
 
 /// Expression class for expression templates
 template <typename Derived>
 struct Expression {
   constexpr const Derived& downcast() const { return static_cast<const Derived&>(*this); }
   constexpr Derived& downcast() { return static_cast<Derived&>(*this); }
 };
 
 template <typename Derived>
 std::ostream& operator<<(std::ostream& os, const Expression<Derived>& e) {
   os << e.downcast();
   return os;
 }
 }  // namespace expolit
 #endif
diff --git a/tests/tests_gtest.cc b/tests/tests_gtest.cc
index b7eb93b..a603731 100644
--- a/tests/tests_gtest.cc
+++ b/tests/tests_gtest.cc
@@ -1,175 +1,175 @@
 /**
  * @file
  *
  * @author Lucas Frérot <lucas.frerot@epfl.ch>
  *
  * @section LICENSE
  *
- * Copyright (©)  2018-2020 EPFL  (Ecole Polytechnique  Fédérale de
+ * Copyright (©)  2018-2021 EPFL  (Ecole Polytechnique  Fédérale de
  * Lausanne)  Laboratory (LSMS  -  Laboratoire de  Simulation  en Mécanique  des
  * Solides)
  *
  * Expolit is free software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as  published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Expolit is  distributed in the  hope that it  will be useful, but  WITHOUT
  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
  * FOR A PARTICULAR PURPOSE. See  the GNU  Lesser General  Public License for
  * more details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Expolit. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 
 #include <expolit/expolit>
 #include <complex>
 #include <gtest/gtest.h>
 
 using namespace expolit;
 using Real = double;
 
 static constexpr Polynomial<Int, 1> Zi({0, 1});
 static constexpr Polynomial<std::complex<Real>, 1> Zc({0, 1});
 static constexpr Polynomial<Real, 1> Zr({0, 1});
 
 using lin = std::remove_cv_t<decltype(Zr)>;
 
 template <typename CT, UInt order>
 constexpr void compare(const Polynomial<CT, order>& p1,
                        const Polynomial<CT, order>& p2) {
   for (UInt i = 0; i < p1.terms; ++i)
     ASSERT_NEAR(std::abs(p1.coeffs[i] - p2.coeffs[i]), 0, 1e-14)
         << "Polynomials are not equal (i = " << i << ")";
 }
 
 template <typename T, UInt order>
 constexpr void compare(const Polynomial<Int, order>& p1,
                        const Polynomial<Int, order>& p2) {
   for (UInt i = 0; i < p1.terms; ++i)
     ASSERT_EQ(p1.coeffs[i], p2.coeffs[i])
         << "Polynomials are not equal (i = " << i << ")";
 }
 
 TEST(Polynomial, initialization) {
   constexpr auto p = 1 + 2 * pow<2>(Zi) + pow<3>(Zi);
   constexpr Int z = 3;
   constexpr Int sol = 1 + 2 * z * z + z * z * z;
   constexpr Int res = p(z);
   static_assert(sol == res, "Bad coefficients");
 
   // This guy cannot be consexpr because std::complex is stupid
   auto p2 = 1. + Zc * Zc;
   std::complex<Real> z2(0, 1);
   ASSERT_NEAR(std::abs(p2(z2)), 0, 1e-14) << "Bad coefficents (complex)";
 }
 
 TEST(Polynomial, sum) {
   constexpr auto p = 1 + 2 * pow<2>(Zi) + pow<3>(Zi);
   constexpr auto p2 = Zi + 2 * Zi * Zi;
   constexpr auto res = 1 + Zi + 4 * Zi * Zi + Zi * Zi * Zi;
   compare(p + p2, res);
 }
 
 TEST(Polynomial, product) {
   constexpr auto p = 1 + 2 * pow<2>(Zi) + pow<3>(Zi);
   constexpr auto p2 = Zi + 2 * pow<2>(Zi);
   constexpr auto res =
       Zi + 2 * pow<2>(Zi) + 2 * pow<3>(Zi) + 5 * pow<4>(Zi) + 2 * pow<5>(Zi);
   compare(p * p2, res);
 }
 
 TEST(Polynomial, derivative) {
   constexpr auto p = 1 + 2 * pow<2>(Zi) + pow<3>(Zi);
   constexpr auto res = 4 * Zi + 3 * Zi * Zi;
   compare(differentiate(p), res);
 }
 
 TEST(Polynomial, integration) {
   constexpr auto p = 1. + 2. * pow<2>(Zr) + pow<3>(Zr);
   constexpr auto res = Zr + 2 / 3. * pow<3>(Zr) + 1 / 4. * pow<4>(Zr);
   compare(integrate(p), res);
 }
 
 TEST(Exponential, evaluation) {
   constexpr auto e = exp(2. * Zr);
   constexpr Real z = 3.4;
   const Real res = e(z);  // cannot be constexpr because of exponential
   ASSERT_NEAR(std::exp(2 * z), res, 1e-14) << "Exponential evaluation fail";
 }
 
 TEST(Exponential, derivative) {
   constexpr auto e = exp(2. * Zr);
   constexpr auto diff = differentiate(e);
   constexpr auto res = Constant<Real>({2}) * e;
   compare(res.operands.first, diff.operands.first);
   compare(e.expression, diff.operands.second.expression);
 }
 
 TEST(Exponential, integration) {
   constexpr auto e = exp(2. * Zr);
   constexpr auto integ = integrate(e);
   constexpr auto res = Constant<Real>({1 / 2.}) * e;
   compare(res.operands.first, integ.operands.first);
   compare(e.expression, integ.operands.second.expression);
 }
 
 TEST(Integration, parts) {
   constexpr auto e = exp(2. * Zr);
   constexpr auto expr = e * (2. * Zr);
   constexpr auto res = Zr * e - Constant<Real>({1 / 2.}) * e;
   constexpr auto integ = integrate(expr);
   compare(integ.operands.first.operands.first,
           res.operands.first.operands.first);
   compare(integ.operands.second.operands.first,
           res.operands.second.operands.first);
 }
 
 TEST(Integration, definite_integral) {
   constexpr auto p = Zr;
   constexpr std::pair<Real, Real> bounds(0, 1);
   constexpr auto res = definite_integral(bounds, p);
   ASSERT_NEAR(0.5, res, 1e-14) << "Definite integral fail";
 
   constexpr auto pi = 2 * Zi;
   static_assert(definite_integral(std::make_pair(0, 1), pi) == 1,
                 "Definite integral fail");
 }
 
 using Q = Litteral<struct Q_>;
 
 TEST(Litteral, substitution) {
   constexpr Q q;
   constexpr auto p = q * Zi;
   constexpr auto c = Constant<Int>({2});
   constexpr auto p_sub = substitute<Q::tag>(c, p);
   compare(2 * Zi, p_sub);
 
   constexpr auto expo = exp(q * Zr);
   constexpr auto x_sub = substitute<Q::tag>(c, expo);
   compare(2. * Zr, x_sub.expression);
 }
 
 TEST(Litteral, differentiation) {
   constexpr Q q;
   constexpr auto c = Constant<Real>({2});
   constexpr auto expo = exp(q * Zr);
   constexpr auto dexpo = differentiate(expo);
   constexpr auto res = substitute<Q::tag>(c, dexpo);
   compare(c, res.operands.first);
   compare(2. * Zr, res.operands.second.expression);
 }
 
 TEST(Litteral, integration) {
   constexpr Q q;
   constexpr auto c = Constant<Real>({2});
   constexpr auto res = substitute<Q::tag>(c, integrate(2 * Zr * exp(q * Zr)));
   constexpr auto sol =
       Zr * exp(2. * Zr) - Constant<Real>({1 / 2.}) * exp(2. * Zr);
   compare(sol.operands.first.operands.first,
           res.operands.first.operands.first);
   compare(sol.operands.second.operands.first,
           res.operands.second.operands.first);
 }