text: #ifndef DYNEARTHSOL3D_UTILS_HPP
  #define DYNEARTHSOL3D_UTILS_HPP
  
  #include <cmath>
  #include <iostream>
  #include <vector>
  
  static void print(std::ostream& os, const double& x)
  {
    os << x;
  }
  
  
  static void print(std::ostream& os, const int& x)
  {
    os << x;
  }
  
  
  static void print(std::ostream& os, const std::size_t& x)
  {
    os << x;
  }
  
  
  template <typename T1, typename T2>
  void print(std::ostream& os, const std::pair<T1,T2>& x)
  {
      os << x.first << ':' << x.second;
  }
  
  
  template <typename T>
  void print(std::ostream& os, const T& A, std::size_t size)
  {
    os << "[";
    for (std::size_t i = 0; i != size; ++i) {
      print(os, A[i]);
      if (i+1 != size)
        os << ", ";
    }
    os << "]";
  }
  
  
  template <typename Array>
  void print(std::ostream& os, const Array& A)
  {
    typename Array::const_iterator i;
    os << "[";
    for (i = A.begin(); i != A.end(); ++i) {
      print(os, *i);
      os << ", ";
    }
    os << "]";
  }
  
  #pragma acc routine seq
  static double trace(const double* s)
  {
  #ifdef THREED
      return s[0] + s[1] + s[2];
  #else
      return s[0] + s[1];
  #endif
  }
  
  
  static double second_invariant2(const double* t)
  {
  #ifdef THREED
      double a = (t[0] + t[1] + t[2]) / 3;
      return ( 0.5 * ((t[0]-a)*(t[0]-a) + (t[1]-a)*(t[1]-a) + (t[2]-a)*(t[2]-a))
               + t[3]*t[3] + t[4]*t[4] + t[5]*t[5] );
  #else
      return 0.25*(t[0]-t[1])*(t[0]-t[1]) + t[2]*t[2];
  #endif
  }
  
  
  static double second_invariant(const double* t)
  {
      /* second invariant of the deviatoric part of tensor t
       * defined as: td = deviatoric(t); sqrt( td(i,j) * td(i,j) / 2)
       */
      return std::sqrt(second_invariant2(t));
  }
  
  #endif
  