diff --git a/src/core/static_types.hh b/src/core/static_types.hh
index e4a2ff8..b53d0a7 100644
--- a/src/core/static_types.hh
+++ b/src/core/static_types.hh
@@ -1,810 +1,810 @@
 /*
  *  SPDX-License-Indentifier: AGPL-3.0-or-later
  *
  *  Copyright (©) 2016-2024 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *  Copyright (©) 2020-2024 Lucas Frérot
  *
  *  This program is free software: you can redistribute it and/or modify
  *  it under the terms of the GNU Affero General Public License as published
  *  by the Free Software Foundation, either version 3 of the License, or
  *  (at your option) any later version.
  *
  *  This program 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 Affero General Public License for more details.
  *
  *  You should have received a copy of the GNU Affero General Public License
  *  along with this program.  If not, see <https://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #ifndef STATIC_TYPES_HH
 #define STATIC_TYPES_HH
 /* -------------------------------------------------------------------------- */
 #include "tamaas.hh"
-#include <thrust/detail/config/host_device.h>
+#include <thrust/detail/config/config.h>
 #include <thrust/sort.h>
 #include <type_traits>
 
 namespace tamaas {
 
 /* -------------------------------------------------------------------------- */
 
 namespace detail {
 /// \cond HIDDEN_SYMBOLS
 template <UInt acc, UInt n, UInt... ns>
 struct product_tail_rec : product_tail_rec<acc * n, ns...> {};
 /// \endcond
 
 template <UInt acc, UInt n>
 struct product_tail_rec<acc, n> : std::integral_constant<UInt, acc * n> {};
 
 /// \cond HIDDEN_SYMBOLS
 template <UInt N, UInt n, UInt... ns>
 struct get_rec : get_rec<N - 1, ns...> {};
 /// \endcond
 
 template <UInt n, UInt... ns>
 struct get_rec<0, n, ns...> : std::integral_constant<UInt, n> {};
 }  // namespace detail
 
 template <UInt... ns>
 struct product : detail::product_tail_rec<1, ns...> {};
 
 template <UInt N, UInt... ns>
 struct get : detail::get_rec<N, ns...> {};
 
 template <typename T>
 struct is_arithmetic : std::is_arithmetic<T> {};
 
 template <typename T>
 struct is_arithmetic<complex<T>> : std::true_type {};
 
 template <UInt dim>
 struct voigt_size;
 
 template <>
 struct voigt_size<3> : std::integral_constant<UInt, 6> {};
 
 template <>
 struct voigt_size<2> : std::integral_constant<UInt, 3> {};
 
 template <>
 struct voigt_size<1> : std::integral_constant<UInt, 1> {};
 
 /* -------------------------------------------------------------------------- */
 
 /**
  * @brief Static Array
  *
  * This class is meant to be a small and fast object for intermediate
  * calculations, possibly on wrapped memory belonging to a grid. Support type
  * show be either a pointer or a C array. It should not contain any virtual
  * method.
  */
 template <typename DataType, typename SupportType, UInt _size>
 class StaticArray {
   static_assert(std::is_array<SupportType>::value or
                     std::is_pointer<SupportType>::value,
                 "the support type of StaticArray should be either a pointer or "
                 "a C-array");
 
   using T = DataType;
   using T_bare = typename std::remove_cv_t<T>;
 
 public:
   using value_type = T;
   static constexpr UInt size = _size;
 
 public:
   StaticArray() = default;
   ~StaticArray() = default;
   StaticArray(const StaticArray&) = delete;
   StaticArray(StaticArray&&) = delete;
   StaticArray& operator=(StaticArray&&) = delete;
 
 public:
   /// Access operator
   __device__ __host__ auto operator()(UInt i) -> T& { return _mem[i]; }
 
   /// Access operator
   __device__ __host__ auto operator()(UInt i) const -> const T& {
     return _mem[i];
   }
 
   /// Scalar product
   template <typename DT, typename ST>
   __device__ __host__ auto dot(const StaticArray<DT, ST, size>& o) const
       -> T_bare {
     decltype(T_bare(0) * DT(0)) res = 0;
     for (UInt i = 0; i < size; ++i)
       res += (*this)(i)*o(i);
     return res;
   }
 
   /// L2 norm squared
   __device__ __host__ T_bare l2squared() const { return this->dot(*this); }
 
   /// L2 norm
   __device__ __host__ T_bare l2norm() const { return std::sqrt(l2squared()); }
 
   /// Sum of all elements
   __device__ __host__ T_bare sum() const {
     T_bare res = 0;
     for (UInt i = 0; i < size; ++i)
       res += (*this)(i);
     return res;
   }
 
 #define VECTOR_OP(op)                                                          \
   template <typename DT, typename ST>                                          \
   __device__ __host__ void operator op(const StaticArray<DT, ST, size>& o) {   \
     for (UInt i = 0; i < size; ++i)                                            \
       (*this)(i) op o(i);                                                      \
   }
 
   VECTOR_OP(+=)
   VECTOR_OP(-=)
   VECTOR_OP(*=)
   VECTOR_OP(/=)
 #undef VECTOR_OP
 
 #define SCALAR_OP(op)                                                          \
   template <typename T1>                                                       \
   __device__ __host__ std::enable_if_t<is_arithmetic<T1>::value, StaticArray&> \
   operator op(const T1& x) {                                                   \
     for (UInt i = 0; i < size; ++i)                                            \
       (*this)(i) op x;                                                         \
     return *this;                                                              \
   }
 
   SCALAR_OP(+=)
   SCALAR_OP(-=)
   SCALAR_OP(*=)
   SCALAR_OP(/=)
   SCALAR_OP(=)
 #undef SCALAR_OP
 
   /// Overriding the implicit copy operator
   __device__ __host__ StaticArray& operator=(const StaticArray& o) {
     copy(o);
     return *this;
   }
 
   template <typename DT, typename ST>
   __device__ __host__ void operator=(const StaticArray<DT, ST, size>& o) {
     copy(o);
   }
 
   template <typename DT, typename ST>
   __device__ __host__ StaticArray& copy(const StaticArray<DT, ST, size>& o) {
     for (UInt i = 0; i < size; ++i)
       (*this)(i) = o(i);
     return *this;
   }
 
   T* begin() { return _mem; }
   const T* begin() const { return _mem; }
 
   T* end() { return _mem + size; }
   const T* end() const { return _mem + size; }
 
 private:
   template <typename U>
   using valid_size_t = std::enable_if_t<(size > 0), U>;
 
 public:
   __host__ __device__ valid_size_t<T&> front() { return *_mem; }
 
   __host__ __device__ valid_size_t<const T&> front() const { return *_mem; }
 
   __host__ __device__ valid_size_t<T&> back() { return _mem[size - 1]; }
 
   __host__ __device__ valid_size_t<const T&> back() const {
     return _mem[size - 1];
   }
 
 protected:
   SupportType _mem;
 };
 
 /**
  * @brief Static Tensor
  *
  * This class implements a multi-dimensional tensor behavior.
  */
 template <typename DataType, typename SupportType = DataType*, UInt... dims>
 class StaticTensor
     : public StaticArray<DataType, SupportType, product<dims...>::value> {
   using parent = StaticArray<DataType, SupportType, product<dims...>::value>;
   using T = DataType;
 
 public:
   static constexpr UInt dim = sizeof...(dims);
   using parent::operator=;
 
 private:
   template <typename... Idx>
   __device__ __host__ static UInt unpackOffset(UInt offset, UInt index,
                                                Idx... rest) {
     constexpr UInt size = sizeof...(rest);
     offset += index;
     offset *= get<dim - size, dims...>::value;
     return unpackOffset(offset, rest...);
   }
 
   template <typename... Idx>
   __device__ __host__ static UInt unpackOffset(UInt offset, UInt index) {
     return offset + index;
   }
 
 public:
   template <typename... Idx>
   __device__ __host__ const T& operator()(Idx... idx) const {
     return parent::operator()(unpackOffset(0, idx...));
   }
 
   template <typename... Idx>
   __device__ __host__ T& operator()(Idx... idx) {
     return parent::operator()(unpackOffset(0, idx...));
   }
 };
 
 /* -------------------------------------------------------------------------- */
 /* Common Static Types */
 /* -------------------------------------------------------------------------- */
 // Forward declaration
 template <typename DataType, typename SupportType, UInt n>
 class StaticVector;
 template <typename DataType, typename SupportType, UInt n>
 class StaticSymMatrix;
 
 /* -------------------------------------------------------------------------- */
 
 template <typename DataType, typename SupportType, UInt n, UInt m>
 class StaticMatrix : public StaticTensor<DataType, SupportType, n, m> {
   using T = DataType;
   using T_bare = typename std::remove_cv_t<T>;
 
 public:
   using StaticTensor<DataType, SupportType, n, m>::operator=;
 
   // /// Initialize from a symmetric matrix
   template <typename DT, typename ST>
   __device__ __host__ std::enable_if_t<n == m>
   fromSymmetric(const StaticSymMatrix<DT, ST, n>& o);
   /// Outer product of two vectors
   template <typename DT1, typename ST1, typename DT2, typename ST2>
   __device__ __host__ void outer(const StaticVector<DT1, ST1, n>& a,
                                  const StaticVector<DT2, ST2, m>& b);
 
   template <typename DT1, typename ST1, typename DT2, typename ST2, UInt l>
   __device__ __host__ void mul(const StaticMatrix<DT1, ST1, n, l>& a,
                                const StaticMatrix<DT2, ST2, l, m>& b) {
     (*this) = T(0);
     for (UInt i = 0; i < n; ++i)
       for (UInt j = 0; j < m; ++j)
         for (UInt k = 0; k < l; ++k)
           (*this)(i, j) += a(i, k) * b(k, j);
   }
 
   __device__ __host__ std::enable_if_t<n == m, T_bare> trace() const {
     T_bare res{0};
     for (UInt i = 0; i < n; ++i)
       res += (*this)(i, i);
     return res;
   }
 
   template <typename DT1, typename ST1>
   __device__ __host__ std::enable_if_t<n == m>
   deviatoric(const StaticMatrix<DT1, ST1, n, m>& mat, Real factor = n) {
     auto norm_trace = mat.trace() / factor;
     for (UInt i = 0; i < n; ++i)
       for (UInt j = 0; j < m; ++j)
         (*this)(i, j) = mat(i, j) - (i == j) * norm_trace;
   }
 };
 
 /* -------------------------------------------------------------------------- */
 
 /// Vector class with size determined at compile-time
 template <typename DataType, typename SupportType, UInt n>
 class StaticVector : public StaticTensor<DataType, SupportType, n> {
   using T = std::remove_cv_t<DataType>;
 
 public:
   using StaticTensor<DataType, SupportType, n>::operator=;
   /// Matrix-vector product
   template <bool transpose, typename DT1, typename ST1, typename DT2,
             typename ST2, UInt m>
   __device__ __host__ void mul(const StaticMatrix<DT1, ST1, n, m>& mat,
                                const StaticVector<DT2, ST2, m>& vec) {
     *this = T(0);
     for (UInt i = 0; i < n; ++i)
       for (UInt j = 0; j < m; ++j)
         (*this)(i) +=
             ((transpose) ? mat(j, i) : mat(i, j)) * vec(j);  // can be optimized
   }
 };
 
 /* -------------------------------------------------------------------------- */
 
 /// Symmetric matrix in Voigt notation
 template <typename DataType, typename SupportType, UInt n>
 class StaticSymMatrix
     : public StaticVector<DataType, SupportType, voigt_size<n>::value> {
   using parent = StaticVector<DataType, SupportType, voigt_size<n>::value>;
   using T = std::remove_cv_t<DataType>;
 
 private:
   template <typename DT, typename ST, typename BinOp>
   __device__ __host__ void sym_binary(const StaticMatrix<DT, ST, n, n>& m,
                                       BinOp&& op) {
     for (UInt i = 0; i < n; ++i)
       op((*this)(i), m(i, i));
 
     const auto a = 0.5 * std::sqrt(2);
     for (UInt i = 0, b = n; i < n; ++i)
       for (UInt j = i + 1; j < n; ++j)
         op((*this)(b++), a * (m(i, j) + m(j, i)));
 
     // Old Mendel-Voigt loops
     // for (UInt j = n - 1, b = n; j > 0; --j)
     //   for (Int i = j - 1; i >= 0; --i)
     //     op((*this)(b++), a * (m(i, j) + m(j, i)));
   }
 
 public:
   /// Copy values from matrix and symmetrize
   template <typename DT, typename ST>
   __device__ __host__ void symmetrize(const StaticMatrix<DT, ST, n, n>& m) {
     sym_binary(m, [](auto&& v, auto&& w) { v = w; });
   }
 
   /// Add values from symmetrized matrix
   template <typename DT, typename ST>
   __device__ __host__ void operator+=(const StaticMatrix<DT, ST, n, n>& m) {
     sym_binary(m, [](auto&& v, auto&& w) { v += w; });
   }
 
   __device__ __host__ auto trace() const {
     std::remove_cv_t<DataType> res = 0;
     for (UInt i = 0; i < n; ++i)
       res += (*this)(i);
     return res;
   }
 
   template <typename DT, typename ST>
   __device__ __host__ void deviatoric(const StaticSymMatrix<DT, ST, n>& m,
                                       Real factor = n) {
     auto tr = m.trace() / factor;
     for (UInt i = 0; i < n; ++i)
       (*this)(i) = m(i) - tr;
     for (UInt i = n; i < voigt_size<n>::value; ++i)
       (*this)(i) = m(i);
   }
 
   using parent::operator+=;
   using parent::operator=;
 };
 
 /* -------------------------------------------------------------------------- */
 
 // Implementation of constructor from symmetric matrix
 template <typename DataType, typename SupportType, UInt n, UInt m>
 template <typename DT, typename ST>
 __device__ __host__ std::enable_if_t<n == m>
 StaticMatrix<DataType, SupportType, n, m>::fromSymmetric(
     const StaticSymMatrix<DT, ST, n>& o) {
   for (UInt i = 0; i < n; ++i)
     (*this)(i, i) = o(i);
 
   // We use Mendel-Voigt notation for the vector representation
   // Specifically from
   // https://thelfer.github.io/tfel/web/tensors.html#vector-notations-for-symmetric-tensors
   const auto a = 1. / std::sqrt(2);
   for (UInt i = 0, b = n; i < n; ++i)
     for (UInt j = i + 1; j < n; ++j)
       (*this)(i, j) = (*this)(j, i) = a * o(b++);
 
   // Old Mendel-Voigt
   // for (UInt j = n - 1, b = n; j > 0; --j)
   //   for (Int i = j - 1; i >= 0; --i)
   //     (*this)(i, j) = (*this)(j, i) = a * o(b++);
 }
 
 // Implementation of outer product
 template <typename DataType, typename SupportType, UInt n, UInt m>
 template <typename DT1, typename ST1, typename DT2, typename ST2>
 __device__ __host__ void StaticMatrix<DataType, SupportType, n, m>::outer(
     const StaticVector<DT1, ST1, n>& a, const StaticVector<DT2, ST2, m>& b) {
   for (UInt i = 0; i < n; ++i)
     for (UInt j = 0; j < m; ++j)
       (*this)(i, j) = a(i) * b(j);
 }
 
 /* -------------------------------------------------------------------------- */
 /* On the stack static types */
 /* -------------------------------------------------------------------------- */
 template <template <typename, typename, UInt...> class StaticParent,
           UInt... dims>
 struct static_size_helper : product<dims...> {};
 
 template <UInt n>
 struct static_size_helper<StaticSymMatrix, n> : voigt_size<n> {};
 
 template <template <typename, typename, UInt...> class StaticParent, typename T,
           UInt... dims>
 class Tensor
     : public StaticParent<
           T, T[static_size_helper<StaticParent, dims...>::value], dims...> {
   static constexpr UInt size = static_size_helper<StaticParent, dims...>::value;
   using parent = StaticParent<T, T[size], dims...>;
 
 public:
   using parent::operator=;
   using parent::copy;
 
   /// Default constructor
   Tensor() = default;
   /// Construct with default value
   __device__ __host__ Tensor(T val) { *this = val; }
   /// Construct from array
   __device__ __host__ Tensor(const std::array<T, size>& arr) {
     // we use size to ensure static loop unrolling
     for (UInt i = 0; i < size; ++i)
       this->_mem[i] = arr[i];
   }
 
   /// Copy from array
   __device__ __host__ Tensor& operator=(const std::array<T, size>& arr) {
     // we use size to ensure static loop unrolling
     for (UInt i = 0; i < size; ++i)
       (*this)(i) = arr[i];
   }
 
   /// Construct by copy from static tensor
   template <typename DT, typename ST>
   __device__ __host__ Tensor(const StaticParent<DT, ST, dims...>& o) {
     copy(o);
   }
 
   __device__ __host__ Tensor(const Tensor& o) : parent() { copy(o); }
   __device__ __host__ Tensor& operator=(const Tensor& o) {
     copy(o);
     return *this;
   }
 
   __device__ __host__ Tensor(Tensor&& o) noexcept { copy(o); }
 };
 
 template <typename T, UInt n, UInt m>
 using Matrix = Tensor<StaticMatrix, T, n, m>;
 
 template <typename T, UInt n>
 using SymMatrix = Tensor<StaticSymMatrix, T, n>;
 
 template <typename T, UInt n>
 using Vector = Tensor<StaticVector, T, n>;
 
 /* -------------------------------------------------------------------------- */
 /* Proxy Static Types */
 /* -------------------------------------------------------------------------- */
 
 /// Proxy type for tensor
 template <template <typename, typename, UInt...> class StaticParent, typename T,
           UInt... dims>
 class TensorProxy : public StaticParent<T, T*, dims...> {
   using parent = StaticParent<T, T*, dims...>;
 
 public:
   /// Explicit construction from data location
   __device__ __host__ explicit TensorProxy(T* spot) { this->_mem = spot; }
   /// Explicit construction from lvalue-reference
   __device__ __host__ explicit TensorProxy(T& spot) : TensorProxy(&spot) {}
   /// Construction from static tensor
   template <typename DataType, typename SupportType>
   __device__ __host__
   TensorProxy(StaticParent<DataType, SupportType, dims...>& o)
       : TensorProxy(o.begin()) {}
 
   using parent::operator=;
 
   __device__ __host__ TensorProxy(const TensorProxy& o) { this->_mem = o._mem; }
   __device__ __host__ TensorProxy& operator=(const TensorProxy& o) {
     this->copy(o);
     return *this;
   }
   __device__ __host__ TensorProxy(TensorProxy&& o) noexcept
       : TensorProxy(exchange(o._mem, nullptr)) {}
 
 public:
   using stack_type = Tensor<StaticParent, T, dims...>;
 };
 
 template <typename T, UInt n, UInt m>
 using MatrixProxy = TensorProxy<StaticMatrix, T, n, m>;
 
 template <typename T, UInt n>
 using SymMatrixProxy = TensorProxy<StaticSymMatrix, T, n>;
 
 template <typename T, UInt n>
 using VectorProxy = TensorProxy<StaticVector, T, n>;
 
 /* -------------------------------------------------------------------------- */
 /* -------------------------------------------------------------------------- */
 /* Arithmetic operators creating temporaries */
 /* -------------------------------------------------------------------------- */
 /* -------------------------------------------------------------------------- */
 
 /* -------------------------------------------------------------------------- */
 /* Simple operators */
 /* -------------------------------------------------------------------------- */
 
 template <typename DT1, typename ST1, typename DT2, typename ST2, UInt dim>
 __device__ __host__ Vector<decltype(DT1(0) + DT2(0)), dim>
 operator+(const StaticVector<DT1, ST1, dim>& a,
           const StaticVector<DT2, ST2, dim>& b) {
   Vector<decltype(DT1(0) + DT2(0)), dim> res(a);
   res += b;
   return res;
 }
 
 template <typename DT1, typename ST1, typename DT2, typename ST2, UInt dim>
 __device__ __host__ Vector<decltype(DT1(0) - DT2(0)), dim>
 operator-(const StaticVector<DT1, ST1, dim>& a,
           const StaticVector<DT2, ST2, dim>& b) {
   Vector<decltype(DT1(0) - DT2(0)), dim> res(a);
   res -= b;
   return res;
 }
 
 template <typename DT1, typename ST1, UInt dim>
 __device__ __host__ Vector<decltype(DT1(0)), dim>
 operator-(const StaticVector<DT1, ST1, dim>& a) {
   Vector<decltype(DT1(0)), dim> res(a);
   res *= -1;
   return res;
 }
 
 template <typename DT1, typename ST1, typename DT2, typename ST2, UInt n,
           UInt m>
 __device__ __host__ Matrix<decltype(DT1(0) + DT2(0)), n, m>
 operator+(const StaticMatrix<DT1, ST1, n, m>& a,
           const StaticMatrix<DT2, ST2, n, m>& b) {
   Matrix<decltype(DT1(0) + DT2(0)), n, m> res(a);
   res += b;
   return res;
 }
 
 template <typename DT1, typename ST1, typename DT2, typename ST2, UInt n,
           UInt m>
 __device__ __host__ Matrix<decltype(DT1(0) - DT2(0)), n, m>
 operator-(const StaticMatrix<DT1, ST1, n, m>& a,
           const StaticMatrix<DT2, ST2, n, m>& b) {
   Matrix<decltype(DT1(0) - DT2(0)), n, m> res(a);
   res -= b;
   return res;
 }
 
 template <typename DT1, typename ST1, UInt n, UInt m>
 __device__ __host__ Matrix<decltype(DT1(0)), n, m>
 operator-(const StaticMatrix<DT1, ST1, n, m>& a) {
   Matrix<decltype(DT1(0)), n, m> res(a);
   res *= -1;
   return res;
 }
 
 template <typename DT1, typename ST1, typename DT2, typename ST2, UInt n>
 __device__ __host__ SymMatrix<decltype(DT1(0) + DT2(0)), n>
 operator+(const StaticSymMatrix<DT1, ST1, n>& a,
           const StaticSymMatrix<DT2, ST2, n>& b) {
   SymMatrix<decltype(DT1(0) + DT2(0)), n> res(a);
   res += b;
   return res;
 }
 
 template <typename DT1, typename ST1, typename DT2, typename ST2, UInt n>
 __device__ __host__ SymMatrix<decltype(DT1(0) - DT2(0)), n>
 operator-(const StaticSymMatrix<DT1, ST1, n>& a,
           const StaticSymMatrix<DT2, ST2, n>& b) {
   SymMatrix<decltype(DT1(0) - DT2(0)), n> res(a);
   res -= b;
   return res;
 }
 
 template <typename DT1, typename ST1, UInt dim>
 __device__ __host__ SymMatrix<decltype(DT1(0)), dim>
 operator-(const StaticSymMatrix<DT1, ST1, dim>& a) {
   SymMatrix<decltype(DT1(0)), dim> res(a);
   res *= -1;
   return res;
 }
 
 template <typename DT, typename ST, typename T, UInt n,
           typename = std::enable_if_t<is_arithmetic<T>::value>>
 __device__ __host__ Vector<decltype(DT(0) * T(0)), n>
 operator*(const StaticVector<DT, ST, n>& a, const T& b) {
   Vector<decltype(DT(0) * T(0)), n> res{a};
   res *= b;
   return res;
 }
 
 // symmetry
 template <typename DT, typename ST, typename T, UInt n,
           typename = std::enable_if_t<is_arithmetic<T>::value>>
 __device__ __host__ Vector<decltype(DT(0) * T(0)), n>
 operator*(const T& b, const StaticVector<DT, ST, n>& a) {
   return a * b;
 }
 
 template <typename DT, typename ST, typename T, UInt n, UInt m,
           typename = std::enable_if_t<is_arithmetic<T>::value>>
 __device__ __host__ Matrix<decltype(DT(0) * T(0)), n, m>
 operator*(const StaticMatrix<DT, ST, n, m>& a, const T& b) {
   Matrix<decltype(DT(0) * T(0)), n, m> res{a};
   res *= b;
   return res;
 }
 
 // symmetry
 template <typename DT, typename ST, typename T, UInt n, UInt m,
           typename = std::enable_if_t<is_arithmetic<T>::value, void>>
 __device__ __host__ Matrix<decltype(DT(0) * T(0)), n, m>
 operator*(const T& b, const StaticMatrix<DT, ST, n, m>& a) {
   return a * b;
 }
 
 template <typename DT, typename ST, typename T, UInt n,
           typename = std::enable_if_t<is_arithmetic<T>::value>>
 __device__ __host__ SymMatrix<decltype(DT(0) * T(0)), n>
 operator*(const StaticSymMatrix<DT, ST, n>& a, const T& b) {
   SymMatrix<decltype(DT(0) * T(0)), n> res{a};
   res *= b;
   return res;
 }
 
 // symmetry
 template <typename DT, typename ST, typename T, UInt n,
           typename = std::enable_if_t<is_arithmetic<T>::value>>
 __device__ __host__ SymMatrix<decltype(DT(0) * T(0)), n>
 operator*(const T& b, const StaticSymMatrix<DT, ST, n>& a) {
   return a * b;
 }
 
 /* -------------------------------------------------------------------------- */
 /* Linear algebra operators */
 /* -------------------------------------------------------------------------- */
 
 /// Matrix-vector multiplication
 template <typename DT1, typename ST1, typename DT2, typename ST2, UInt n,
           UInt m>
 __device__ __host__ Vector<decltype(DT1(0) * DT2(0)), n>
 operator*(const StaticMatrix<DT1, ST1, n, m>& a,
           const StaticVector<DT2, ST2, m>& b) {
   Vector<decltype(DT1(0) * DT2(0)), n> res;
   res.template mul<false>(a, b);
   return res;
 }
 
 /// Matrix-matrix multiplication
 template <typename DT1, typename ST1, typename DT2, typename ST2, UInt n,
           UInt m, UInt l>
 __device__ __host__ Matrix<decltype(DT1(0) * DT2(0)), n, m>
 operator*(const StaticMatrix<DT1, ST1, n, l>& a,
           const StaticMatrix<DT2, ST2, l, m>& b) {
   Matrix<decltype(DT1(0) * DT2(0)), n, m> res;
   res.mul(a, b);
   return res;
 }
 
 template <typename DT1, typename ST1, typename DT2, typename ST2, UInt n,
           UInt m>
 __device__ __host__ Matrix<decltype(DT1(0) * DT2(0)), n, m>
 outer(const StaticVector<DT1, ST1, n>& a, const StaticVector<DT2, ST2, m>& b) {
   Matrix<decltype(DT1(0) * DT2(0)), n, m> res;
   res.outer(a, b);
   return res;
 }
 
 /* -------------------------------------------------------------------------- */
 /* Dense/Sparse */
 /* -------------------------------------------------------------------------- */
 template <typename DT, typename ST, UInt n>
 __device__ __host__ Matrix<std::remove_cv_t<DT>, n, n>
 dense(const StaticSymMatrix<DT, ST, n>& m) {
   Matrix<std::remove_cv_t<DT>, n, n> res;
   res.fromSymmetric(m);
   return res;
 }
 
 template <typename DT, typename ST, UInt n>
 __device__ __host__ auto dense(const StaticVector<DT, ST, n>& v)
     -> Vector<DT, n> {
   return v;
 }
 
 template <typename DT, typename ST, UInt n>
 __device__ __host__ SymMatrix<std::remove_cv_t<DT>, n>
 symmetrize(const StaticMatrix<DT, ST, n, n>& m) {
   SymMatrix<std::remove_cv_t<DT>, n> res;
   res.symmetrize(m);
   return res;
 }
 
 /* -------------------------------------------------------------------------- */
 template <typename DT, typename ST>
 __device__ __host__ Vector<std::remove_cv_t<DT>, 3>
 invariants(const StaticSymMatrix<DT, ST, 3>& m) {
   return {
       {// I1 = tr(A)
        m.trace(),
        // I2 = 1/2 * (tr(A)^2 - tr(A^2))
        m(0) * m(1) + m(1) * m(2) + m(0) * m(2) -
            (m(3) * m(3) + m(4) * m(4) + m(5) * m(5)) / 2,
        // I3 = det(A)
        m(0) * m(1) * m(2) + m(5) * m(3) * m(4) / std::sqrt(Real(2)) -
            (m(4) * m(4) * m(1) + m(5) * m(5) * m(0) + m(3) * m(3) * m(2)) / 2}};
 }
 
 template <typename DT, typename ST>
 __device__ __host__ Vector<std::remove_cv_t<DT>, 3>
 eigenvalues(const StaticSymMatrix<DT, ST, 3>& m) {
   constexpr UInt n = 3;
   Vector<std::remove_cv_t<DT>, n> eigenv;
   auto inv = invariants(m);
   Real a = 1, b = -inv(0), c = inv(1), d = -inv(2);
   auto p = (3 * a * c - b * b) / (3 * a * a);
   auto q = (2 * b * b * b - 9 * a * b * c + 27 * a * a * d) / (27 * a * a * a);
 
   for (UInt k = 0; k < n; ++k)
     eigenv(k) =
         2. * std::sqrt(-p / 3.) *
             std::cos(1. / 3. *
                          std::acos(3. * q / (2. * p) * std::sqrt(-3. / p)) -
                      2. * M_PI * k / 3.) -
         b / (3. * a);
 
   // Sorting eigenvalues
   std::remove_cv_t<DT> max = eigenv(0);
   UInt max_k = 0;
   for (UInt k = 1; k < n; ++k) {
     if (eigenv(k) > max) {
       max = eigenv(k);
       max_k = k;
     }
   }
 
   thrust::swap(eigenv.back(), eigenv(max_k));
 
   if (eigenv(0) > eigenv(1))
     thrust::swap(eigenv(0), eigenv(1));
 
   return eigenv;
 }
 
 /* -------------------------------------------------------------------------- */
 /* Type traits */
 /* -------------------------------------------------------------------------- */
 
 template <class Type>
 struct is_proxy : std::false_type {};
 
 template <template <typename, typename, UInt...> class StaticParent, typename T,
           UInt... dims>
 struct is_proxy<TensorProxy<StaticParent, T, dims...>> : std::true_type {};
 
 template <typename T, UInt n, UInt m>
 struct is_proxy<MatrixProxy<T, n, m>> : std::true_type {};
 
 template <typename T, UInt n>
 struct is_proxy<SymMatrixProxy<T, n>> : std::true_type {};
 
 template <typename T, UInt n>
 struct is_proxy<VectorProxy<T, n>> : std::true_type {};
 
 }  // namespace tamaas
 
 #endif  // STATIC_TYPES_HH