Templates.hpp

Go to the documentation of this file.

/** \file Templates.hpp
 * \brief Templates declarations
 */
 
#ifndef __TEMPLATES_HPP__
#define __TEMPLATES_HPP__
 
namespace MoFEM {
 
/** @name Data layout helpers
 *  Traits describing tensor storage layout conventions.
 *  @{
 */
enum class DataLayout { GaussByCoeffs, CoeffsByGauss };
template <DataLayout DL> struct DataLayoutTraits;
template <> struct DataLayoutTraits<DataLayout::GaussByCoeffs> {
  static constexpr bool isGaussByCoeffs = true;
  static constexpr bool isCoeffsByGauss = false;
};
template <> struct DataLayoutTraits<DataLayout::CoeffsByGauss> {
  static constexpr bool isGaussByCoeffs = false;
  static constexpr bool isCoeffsByGauss = true;
};
/** @} */
 
/** @name Adaptor functions
 *  Views that wrap raw storage as vector or matrix adaptors.
 *  @{
 */
/**
 * @brief Get Vector adaptor
 *
 * \code
 *
 * double *a;
 * CHKERR VecGetArray(v,&a);
 *
 * for(int n = 0; n != nodes; ++n) {
 *
 *   auto a = getVectorAdaptor(&a[3*n], 3);
 *   double dot = inner_prod(a, a);
 *
 * }
 *
 * CHKERR VecRetsoreArray(v,&a);
 * \endcode
 *
 */
template <typename T1> inline auto getVectorAdaptor(T1 ptr, const size_t n) {
  typedef typename std::remove_pointer<T1>::type T;
  return VectorShallowArrayAdaptor<T>(n,
                                      ublas::shallow_array_adaptor<T>(n, ptr));
};
 
/**
 * @brief Get Matrix adaptor
 *
 * \code
 *
 * double *a;
 * CHKERR VecGetArray(v,&a);
 *
 * for(int n = 0; n != nodes; ++n) {
 *
 *   auto F = getMatrixAdaptor(&a[3*3*n], 3, 3);
 *   MatrixDouble C = prod(F, trans(F));
 *
 * }
 *
 * CHKERR VecRetsoreArray(v,&a);
 * \endcode
 *
 */
template <typename T1>
inline auto getMatrixAdaptor(T1 ptr, const size_t n, const size_t m) {
  typedef typename std::remove_pointer<T1>::type T;
  return MatrixShallowArrayAdaptor<T>(
      n, m, ublas::shallow_array_adaptor<T>(n * m, ptr));
};
/** @} */
 
/** @name Generic utility helpers
 *  Shared template utilities used across the header.
 *  @{
 */
/**
 * This small utility that cascades two key extractors will be
 * used throughout the boost example
 * <a
 * href=http://www.boost.org/doc/libs/1_53_0/libs/multi_index/example/complex_structs.cpp>
 * http://www.boost.org/doc/libs/1_53_0/libs/multi_index/example/complex_structs.cpp
 * </a>
 */
template <class KeyExtractor1, class KeyExtractor2> struct KeyFromKey {
public:
  typedef typename KeyExtractor1::result_type result_type;
 
  KeyFromKey(const KeyExtractor1 &key1_ = KeyExtractor1(),
             const KeyExtractor2 &key2_ = KeyExtractor2())
      : key1(key1_), key2(key2_) {}
 
  template <typename Arg> result_type operator()(Arg &arg) const {
    return key1(key2(arg));
  }
 
private:
  KeyExtractor1 key1;
  KeyExtractor2 key2;
};
 
template <typename id_type> struct LtBit {
  inline bool operator()(const id_type &valueA, const id_type &valueB) const {
    return valueA.to_ulong() < valueB.to_ulong();
  }
};
 
template <typename id_type> struct EqBit {
  inline bool operator()(const id_type &valueA, const id_type &valueB) const {
    return valueA.to_ulong() == valueB.to_ulong();
  }
};
 
template <typename id_type> struct HashBit {
  inline unsigned int operator()(const id_type &value) const {
    return value.to_ulong();
  }
};
 
template <class X> inline std::string toString(X x) {
  std::ostringstream buffer;
  buffer << x;
  return buffer.str();
}
/** @} */
 
/** @name Tensor0 functions
 *  Scalar tensor access helpers.
 *  @{
 */
template <int S, class V> struct GetFTensor0FromVecImpl {
  static inline auto get(V &data) {
    using T = std::remove_cv_t<std::remove_reference_t<decltype(data(0))>>;
    return FTensor::Tensor0<FTensor::PackPtr<T *, S>>(data.data().data());
  }
};
 
/**
* \brief Get tensor rank 0 (scalar) form data vector
 
Example how to use it.
\code
VectorDouble vec;
vec.resize(nb_gauss_pts,false);
vec.clear();
auto t0 = getFTensor0FromVec<1>(vec);
for(int gg = 0;gg!=nb_gauss_pts;gg++) {
 
  ++t0;
}
\endcode
 
*/
template <int S = 1, class V> static inline auto getFTensor0FromVec(V &data) {
  return GetFTensor0FromVecImpl<S, V>::get(data);
}
 
template <int S, class M> struct GetFTensor0FromMatImpl {
  static inline auto get(M &data) {
#ifndef NDEBUG
    // We can do rows, or columns, does not matter for rank 0 tensor, but we need to be sure that one of them is 1
    if (data.size1() != 1 && data.size2() != 1)
      THROW_MESSAGE("getFTensor0FromMat: wrong size of data matrix, number of "
                    "rows should be 1 but is " +
                    boost::lexical_cast<std::string>(data.size1()));
#endif
    using T = std::remove_cv_t<std::remove_reference_t<decltype(data(0, 0))>>;
    return FTensor::Tensor0<FTensor::PackPtr<T *, S>>(data.data().data());
  }
};
 
/**
* \brief Get tensor rank 0 (scalar) form data vector
 
Example how to use it.
\code
MatrixDouble vec;
vec.resize(1, nb_gauss_pts, false);
vec.clear();
auto t0 = getFTensor0FromMat<1>(vec);
for(int gg = 0;gg!=nb_gauss_pts;gg++) {
 
  ++t0;
}
\endcode
 
*/
template <int S = 1, class M> inline auto getFTensor0FromMat(M &data) {
  return GetFTensor0FromMatImpl<S, M>::get(data);
}
 
template <int S = 1, class T, class A>
inline auto getFTensor0FromMat(boost::weak_ptr<MatrixDouble> data_ptr) {
#ifndef NDEBUG
  if (data_ptr.expired())
    THROW_MESSAGE("getFTensor0FromMat: data pointer expired");
#endif
  return GetFTensor0FromMatImpl<S, MatrixDouble>::get(*data_ptr.lock());
}
 
template <int S = 1>
inline auto getFTensor0FromMat(boost::shared_ptr<MatrixDouble> data_ptr) {
#ifndef NDEBUG
  if (!data_ptr)
    THROW_MESSAGE("getFTensor0FromMat: data pointer is not available");
#endif
  return GetFTensor0FromMatImpl<S, MatrixDouble>::get(*data_ptr);
}
 
/**
 * \brief Get tensor rank 0 (scalar) from pointer
 */
template <int S = 1, class T> static inline auto getFTensor0FromPtr(T *ptr) {
  return FTensor::Tensor0<FTensor::PackPtr<T *, S>>(ptr);
}
/** @} */
 
/** @name Tensor1 functions
 *  Rank-1 tensor construction helpers.
 *  @{
 */
template <int Tensor_Dim, int S, class DL, class M>
struct GetFTensor1FromMatImpl;
 
template <int Tensor_Dim, int S, class M>
struct GetFTensor1FromMatImpl<Tensor_Dim, S,
                              DataLayoutTraits<DataLayout::GaussByCoeffs>, M> {
  static inline auto get(M &data, int rr, int cc) {
    // Preserve const matrix elements in the non-owning tensor view.
    using T = std::remove_reference_t<decltype(data(0, 0))>;
    constexpr int stride = S == -1 ? Tensor_Dim : S;
    static_assert(
        stride % Tensor_Dim == 0,
        "getFTensor1FromMat: stride should be a multiple of Tensor_Dim");
 
#ifndef NDEBUG
    if (data.size2() != Tensor_Dim)
      THROW_MESSAGE(
          "getFTensor1FromMat<" + boost::lexical_cast<std::string>(Tensor_Dim) +
          ">: wrong size of data matrix, number of columns should be " +
          boost::lexical_cast<std::string>(Tensor_Dim) + " but is " +
          boost::lexical_cast<std::string>(data.size2()));
    auto *first = &data(rr + 0, cc + 0);
    auto *last = &data(rr + 0, cc + 0 + Tensor_Dim - 1);
    if (last - first != Tensor_Dim - 1)
      THROW_MESSAGE("getFTensor1FromMat<" +
                    boost::lexical_cast<std::string>(Tensor_Dim) +
                    ">: row slice is not contiguous");
 
#endif
    return FTensor::Tensor1<FTensor::CursorPtr<T *, stride>, Tensor_Dim>(
        &data(rr + 0, cc + 0));
  }
};
 
template <int S, class M>
struct GetFTensor1FromMatImpl<3, S, DataLayoutTraits<DataLayout::CoeffsByGauss>,
                              M> {
  static inline auto get(M &data, int rr, int cc) {
    using T = std::remove_cv_t<std::remove_reference_t<decltype(data(0, 0))>>;
    constexpr int stride = S == -1 ? 1 : S;
 
#ifndef NDEBUG
    if (data.size1() % 3 != 0)
      THROW_MESSAGE(
          "getFTensor1FromMat<3>: wrong size of data matrix, number of "
          "rows modulo of 3 but is " +
          boost::lexical_cast<std::string>(data.size1()));
#endif
    return FTensor::Tensor1<FTensor::PackPtr<T *, stride>, 3>(
        &data(rr + 0, cc + 0), &data(rr + 1, cc + 0), &data(rr + 2, cc + 0));
  }
};
 
template <int S, class M>
struct GetFTensor1FromMatImpl<4, S, DataLayoutTraits<DataLayout::CoeffsByGauss>,
                              M> {
  static inline auto get(M &data, int rr, int cc) {
    using T = std::remove_cv_t<std::remove_reference_t<decltype(data(0, 0))>>;
    constexpr int stride = S == -1 ? 1 : S;
#ifndef NDEBUG
    if (data.size1() % 4 != 0)
      THROW_MESSAGE(
          "getFTensor1FromMat<4>: wrong size of data matrix, number of "
          "rows should be 4 but is " +
          boost::lexical_cast<std::string>(data.size1()));
#endif
    return FTensor::Tensor1<FTensor::PackPtr<T *, stride>, 4>(
        &data(rr + 0, cc + 0), &data(rr + 1, cc + 0), &data(rr + 2, cc + 0),
        &data(rr + 3, cc + 0));
  }
};
 
template <int S, class M>
struct GetFTensor1FromMatImpl<6, S, DataLayoutTraits<DataLayout::CoeffsByGauss>,
                              M> {
  static inline auto get(M &data, int rr, int cc) {
    using T = std::remove_cv_t<std::remove_reference_t<decltype(data(0, 0))>>;
    constexpr int stride = S == -1 ? 1 : S;
#ifndef NDEBUG
    if (data.size1() % 6 != 0)
      THROW_MESSAGE(
          "getFTensor1FromMat<6>: wrong size of data matrix, number of "
          "rows should be 6 but is " +
          boost::lexical_cast<std::string>(data.size1()));
#endif
    return FTensor::Tensor1<FTensor::PackPtr<T *, stride>, 6>(
        &data(rr + 0, cc + 0), &data(rr + 1, cc + 0), &data(rr + 2, cc + 0),
        &data(rr + 3, cc + 0), &data(rr + 4, cc + 0), &data(rr + 5, cc + 0));
  }
};
 
template <int S, class M>
struct GetFTensor1FromMatImpl<9, S, DataLayoutTraits<DataLayout::CoeffsByGauss>,
                              M> {
  static inline auto get(M &data, int rr, int cc) {
    using T = std::remove_cv_t<std::remove_reference_t<decltype(data(0, 0))>>;
    constexpr int stride = S == -1 ? 1 : S;
#ifndef NDEBUG
    if (data.size1() % 9 != 0)
      THROW_MESSAGE(
          "getFTensor1FromMat<9>: wrong size of data matrix, number of "
          "rows should be 9 but is " +
          boost::lexical_cast<std::string>(data.size1()));
#endif
    return FTensor::Tensor1<FTensor::PackPtr<T *, stride>, 9>(
        &data(rr + 0, cc + 0), &data(rr + 1, cc + 0), &data(rr + 2, cc + 0),
        &data(rr + 3, cc + 0), &data(rr + 4, cc + 0), &data(rr + 5, cc + 0),
        &data(rr + 6, cc + 0), &data(rr + 7, cc + 0), &data(rr + 8, cc + 0));
  }
};
 
template <int S, class M>
struct GetFTensor1FromMatImpl<2, S, DataLayoutTraits<DataLayout::CoeffsByGauss>,
                              M> {
  static inline auto get(M &data, int rr, int cc) {
    using T = std::remove_cv_t<std::remove_reference_t<decltype(data(0, 0))>>;
    constexpr int stride = S == -1 ? 1 : S;
#ifndef NDEBUG
    if (data.size1() % 2 != 0)
      THROW_MESSAGE(
          "getFTensor1FromMat<2>: wrong size of data matrix, number of "
          "row should be 2 but is " +
          boost::lexical_cast<std::string>(data.size1()));
#endif
    return FTensor::Tensor1<FTensor::PackPtr<T *, stride>, 2>(
        &data(rr + 0, cc + 0), &data(rr + 1, cc + 0));
  }
};
 
template <int S, class M>
struct GetFTensor1FromMatImpl<1, S, DataLayoutTraits<DataLayout::CoeffsByGauss>,
                              M> {
  static inline auto get(M &data, int rr, int cc) {
    using T = std::remove_cv_t<std::remove_reference_t<decltype(data(0, 0))>>;
    constexpr int stride = S == -1 ? 1 : S;
    return FTensor::Tensor1<FTensor::PackPtr<T *, stride>, 1>(
        &data(rr + 0, cc + 0));
  }
};
 
/**
 * \brief Get tensor rank 1 (vector) form data matrix
 */
template <int Tensor_Dim, int S = -1,
          class DL = DataLayoutTraits<DataLayout::GaussByCoeffs>,
          class M = MatrixDouble>
inline auto getFTensor1FromMat(M &data, int rr = 0, int cc = 0) {
  return GetFTensor1FromMatImpl<Tensor_Dim, S, DL, M>::get(data, rr, cc);
}
 
template <int Tensor_Dim, int S = -1,
          class DL = DataLayoutTraits<DataLayout::GaussByCoeffs>,
          class M = MatrixDouble>
using GetFTensor1FromMatType =
    decltype(GetFTensor1FromMatImpl<Tensor_Dim, S, DL, M>::get(
        std::declval<M &>(), 0, 0));
 
template <int Tensor_Dim, int S = -1,
          class DL = DataLayoutTraits<DataLayout::GaussByCoeffs>>
inline auto getFTensor1FromMat(boost::weak_ptr<MatrixDouble> data, int rr = 0,
                               int cc = 0) {
#ifndef NDEBUG
  if (!data.lock()) {
    CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL,
                      "Data matrix is not available");
  }
#endif
  return GetFTensor1FromMatImpl<Tensor_Dim, S, DL, MatrixDouble>::get(
      *data.lock(), rr, cc);
}
 
template <int Tensor_Dim, int S = -1,
          class DL = DataLayoutTraits<DataLayout::GaussByCoeffs>>
inline auto getFTensor1FromMat(boost::shared_ptr<MatrixDouble> data, int rr = 0,
                               int cc = 0) {
#ifndef NDEBUG
  if (!data) {
    CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL,
                      "Data matrix is not available");
  }
#endif
  return GetFTensor1FromMatImpl<Tensor_Dim, S, DL, MatrixDouble>::get(*data, rr,
                                                                      cc);
}
 
template <int Tensor_Dim, int S = Tensor_Dim, typename T = double>
inline auto getFTensor1FromPtr(T *ptr) {
  return FTensor::Tensor1<FTensor::CursorPtr<T *, S>, Tensor_Dim>(ptr);
}
 
template <int Tensor_Dim, int S = Tensor_Dim>
inline auto getFTensor1FromPtr(std::complex<double> *ptr) {
  return getFTensor1FromPtr<Tensor_Dim, S, std::complex<double>>(ptr);
}
 
#ifdef WITH_ADOL_C
template <int Tensor_Dim, int S = Tensor_Dim>
inline auto getFTensor1FromPtr(adouble *ptr) {
  return getFTensor1FromPtr<Tensor_Dim, S, adouble>(ptr);
}
#endif
 
template <int Tensor_Dim, int S = Tensor_Dim, typename T = VectorDouble>
inline auto getFTensor1FromArray(T &data) {
  static_assert(S % Tensor_Dim == 0,
                "getFTensor1FromArray(VectorDouble&) requires S to be a "
                "multiple of Tensor_Dim");
 
#ifndef NDEBUG
  if (data.size() % S != 0) {
    CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL,
                      "getFTensor1FromArray: data size should be divisible by "
                      "pack stride dimension");
  }
#endif
  return getFTensor1FromPtr<Tensor_Dim, S>(data.data().data());
}
 
/**
 * @brief Get FTensor1 from array
 *
 * \todo Generalize for different arrays and data types
 *
 * @tparam DIM
 * @param data
 * @param rr
 * @return FTensor::Tensor1<FTensor::PackPtr<double *, DIM>, DIM>
 */
template <int DIM, int S>
inline auto getFTensor1FromArrayDiag(MatrixDouble &data, const size_t rr);
 
template <>
inline auto getFTensor1FromArrayDiag<2, 2>(MatrixDouble &data,
                                           const size_t rr) {
  return FTensor::Tensor1<FTensor::PackPtr<double *, 2>, 2>{&data(rr + 0, 0),
                                                            &data(rr + 1, 1)};
}
 
template <>
inline auto getFTensor1FromArrayDiag<3, 3>(MatrixDouble &data,
                                           const size_t rr) {
  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>{
      &data(rr + 0, 0), &data(rr + 1, 1), &data(rr + 2, 2)};
}
/** @} */
 
/** @name Tensor2 functions
 *  Rank-2 tensor construction and array helper routines.
 *  @{
 */
template <int Tensor_Dim1, int Tensor_Dim2, int S, class DL, class M>
struct GetFTensor2FromMatImpl {
  static inline auto get(M &data, int rr = 0, int cc = 0) {
    // Preserve const matrix elements in the non-owning tensor view.
    using T = std::remove_reference_t<decltype(data(0, 0))>;
    constexpr int stride =
        S == -1 ? (DL::isGaussByCoeffs ? Tensor_Dim1 * Tensor_Dim2 : 1) : S;
    if constexpr (DL::isGaussByCoeffs) {
      static_assert(stride % (Tensor_Dim1 * Tensor_Dim2) == 0,
                    "getFTensor2FromMat: stride should be a multiple of "
                    "Tensor_Dim1 * Tensor_Dim2");
 
#ifndef NDEBUG
      if (data.size2() != Tensor_Dim1 * Tensor_Dim2) {
        THROW_MESSAGE(
            "getFTensor2FromMat<" +
            boost::lexical_cast<std::string>(Tensor_Dim1) + "," +
            boost::lexical_cast<std::string>(Tensor_Dim2) +
            ">: wrong size of columns in data matrix, should be " +
            boost::lexical_cast<std::string>(Tensor_Dim1 * Tensor_Dim2) +
            " but is " + boost::lexical_cast<std::string>(data.size2()));
      }
      auto *first = &data(rr + 0, cc + 0);
      auto *last = &data(rr + 0, cc + 0 + Tensor_Dim1 * Tensor_Dim2 - 1);
      if (last - first != Tensor_Dim1 * Tensor_Dim2 - 1) {
        THROW_MESSAGE("getFTensor2FromMat<" +
                      boost::lexical_cast<std::string>(Tensor_Dim1) + "," +
                      boost::lexical_cast<std::string>(Tensor_Dim2) +
                      ">: row slice is not contiguous");
      }
#endif
      return FTensor::Tensor2<FTensor::CursorPtr<T *, stride>, Tensor_Dim1,
                              Tensor_Dim2>(&data(rr + 0, cc + 0));
    } else if constexpr (DL::isCoeffsByGauss) {
      std::array<T *, Tensor_Dim1 * Tensor_Dim2> ptrs;
#ifndef NDEBUG
      if (data.size1() != Tensor_Dim1 * Tensor_Dim2) {
        THROW_MESSAGE(
            "getFTensor2FromMat<" +
            boost::lexical_cast<std::string>(Tensor_Dim1) + "," +
            boost::lexical_cast<std::string>(Tensor_Dim2) +
            ">: wrong size of rows in data matrix, should be " +
            boost::lexical_cast<std::string>(Tensor_Dim1 * Tensor_Dim2) +
            " but is " + boost::lexical_cast<std::string>(data.size1()));
      }
#endif
      for (auto i = 0; i != Tensor_Dim1 * Tensor_Dim2; ++i)
        ptrs[i] = &data(rr + i, cc + 0);
      return FTensor::Tensor2<FTensor::PackPtr<T *, stride>, Tensor_Dim1,
                              Tensor_Dim2>(ptrs);
    } else {
      static_assert(!std::is_same<M, M>::value,
                    "Unsupported data layout for getFTensor2FromMatImpl");
    }
  }
};
 
/**
 * \brief Get tensor rank 2 (matrix) form data matrix
 */
template <int Tensor_Dim1, int Tensor_Dim2, int S = -1,
          class DL = DataLayoutTraits<DataLayout::GaussByCoeffs>,
          class M = MatrixDouble>
inline auto getFTensor2FromMat(M &data) {
  return GetFTensor2FromMatImpl<Tensor_Dim1, Tensor_Dim2, S, DL, M>::get(data);
}
 
template <int Tensor_Dim0, int Tensor_Dim1, int S = -1,
          class DL = DataLayoutTraits<DataLayout::GaussByCoeffs>,
          class M = MatrixDouble>
using GetFTensor2FromMatType =
    decltype(GetFTensor2FromMatImpl<Tensor_Dim0, Tensor_Dim1, S, DL, M>::get(
        std::declval<M &>(), 0, 0));
 
/**
 * \brief Get tensor rank 2 (matrix) form data matrix
 */
template <int Tensor_Dim1, int Tensor_Dim2, int S = -1,
          class DL = DataLayoutTraits<DataLayout::GaussByCoeffs>>
inline auto getFTensor2FromMat(boost::weak_ptr<MatrixDouble> data_ptr) {
#ifndef NDEBUG
  if (!data_ptr.lock()) {
    CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL,
                      "Data matrix is not available");
  }
#endif
  return GetFTensor2FromMatImpl<Tensor_Dim1, Tensor_Dim2, S, DL,
                                MatrixDouble>::get(*(data_ptr.lock()));
}
 
/**
 * \brief Get tensor rank 2 (matrix) form data matrix
 */
template <int Tensor_Dim1, int Tensor_Dim2, int S = -1,
          class DL = DataLayoutTraits<DataLayout::GaussByCoeffs>>
inline auto getFTensor2FromMat(boost::shared_ptr<MatrixDouble> data_ptr) {
#ifndef NDEBUG
  if (!data_ptr) {
    CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL,
                      "Data matrix is not available");
  }
#endif
  return GetFTensor2FromMatImpl<Tensor_Dim1, Tensor_Dim2, S, DL,
                                MatrixDouble>::get(*data_ptr);
}
 
template <int Tensor_Dim1, int Tensor_Dim2, int S = -1>
inline auto getFTensor2FromVec(VectorDouble &data) {
  auto matrix_data =
      getMatrixAdaptor(&data[0], data.size() / Tensor_Dim1 * Tensor_Dim2,
                       Tensor_Dim1 * Tensor_Dim2);
  return GetFTensor2FromMatImpl<Tensor_Dim1, Tensor_Dim2, S,
                                DataLayoutTraits<DataLayout::GaussByCoeffs>,
                                decltype(matrix_data)>::get(matrix_data);
}
 
template <int Tensor_Dim1, int Tensor_Dim2, int S = -1,
          class DL = DataLayoutTraits<DataLayout::GaussByCoeffs>>
inline auto getFTensor2FromPtr(double *ptr) {
  if constexpr (DL::isGaussByCoeffs) {
    auto matrix_data = getMatrixAdaptor(ptr, 1, Tensor_Dim1 * Tensor_Dim2);
    return GetFTensor2FromMatImpl<Tensor_Dim1, Tensor_Dim2, S, DL,
                                  decltype(matrix_data)>::get(matrix_data);
  } else if constexpr (DL::isCoeffsByGauss) {
    auto matrix_data = getMatrixAdaptor(ptr, Tensor_Dim1 * Tensor_Dim2, 1);
    return GetFTensor2FromMatImpl<Tensor_Dim1, Tensor_Dim2, S, DL,
                                  decltype(matrix_data)>::get(matrix_data);
  }
}
 
template <int Tensor_Dim1, int Tensor_Dim2, int S = -1,
          class DL = DataLayoutTraits<DataLayout::GaussByCoeffs>>
inline auto getFTensor2FromPtr(std::complex<double> *ptr) {
  if constexpr (DL::isGaussByCoeffs) {
    auto matrix_data = getMatrixAdaptor(ptr, 1, Tensor_Dim1 * Tensor_Dim2);
    return GetFTensor2FromMatImpl<Tensor_Dim1, Tensor_Dim2, S, DL,
                                  decltype(matrix_data)>::get(matrix_data);
  } else if constexpr (DL::isCoeffsByGauss) {
    auto matrix_data = getMatrixAdaptor(ptr, Tensor_Dim1 * Tensor_Dim2, 1);
    return GetFTensor2FromMatImpl<Tensor_Dim1, Tensor_Dim2, S, DL,
                                  decltype(matrix_data)>::get(matrix_data);
  }
}
 
/**
 * @brief Make Tensor2 for HVec base from pointer
 *
 * @tparam DIM
 * @param ptr
 * @return FTensor::Tensor2<FTensor::PackPtr<double *, DIM1 * DIM2>, DIM1, DIM2>
 */
template <int DIM1, int DIM2> inline auto getFTensor2HVecFromPtr(double *ptr);
 
template <> inline auto getFTensor2HVecFromPtr<3, 2>(double *ptr) {
  return FTensor::Tensor2<FTensor::PackPtr<double *, 6>, 3, 2>(
      ptr + HVEC0_0, ptr + HVEC0_1,
 
      ptr + HVEC1_0, ptr + HVEC1_1,
 
      ptr + HVEC2_0, ptr + HVEC2_1);
}
 
template <> inline auto getFTensor2HVecFromPtr<3, 3>(double *ptr) {
  return FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3>(
      ptr + HVEC0_0, ptr + HVEC0_1, ptr + HVEC0_2,
 
      ptr + HVEC1_0, ptr + HVEC1_1, ptr + HVEC1_2,
 
      ptr + HVEC2_0, ptr + HVEC2_1, ptr + HVEC2_2);
}
/**
 * @brief Get FTensor2 from array
 *
 * \note Generalize for other data types
 *
 * @tparam DIM1
 * @tparam DIM2
 * @tparam S
 * @param data
 * @return FTensor::Tensor2<FTensor::PackPtr<T *, S>, DIM1, DIM2>
 */
template <int DIM1, int DIM2, int S, class M> struct GetFTensor2FromArrayImpl {
  GetFTensor2FromArrayImpl() = delete;
  inline static auto get(M &data, const size_t rr, const size_t cc = 0) {
    return make(data, rr, cc, std::make_index_sequence<DIM1 * DIM2>{});
  }
 
private:
  template <size_t... Is>
  inline static auto make(M &data, const size_t rr, const size_t cc,
                          std::index_sequence<Is...>) {
    using T = std::remove_cv_t<std::remove_reference_t<decltype(data(0, 0))>>;
    return FTensor::Tensor2<FTensor::PackPtr<T *, S>, DIM1, DIM2>{
        &data(rr + Is / DIM2, cc + Is % DIM2)...};
  }
};
 
/**
 * @brief Get FTensor2 from array
 *
 * \note Generalize for other data types
 *
 * @tparam DIM1
 * @tparam DIM2
 * @tparam S
 * @param data
 * @return FTensor::Tensor2<T *, DIM1, DIM2>
 */
template <int DIM1, int DIM2, class M> struct GetFTensor2FromArrayRawPtrImpl;
 
template <class M> struct GetFTensor2FromArrayRawPtrImpl<2, 2, M> {
  GetFTensor2FromArrayRawPtrImpl() = delete;
  inline static auto get(M &data, const size_t rr, const size_t cc,
                         const int ss = 0) {
    using T = std::remove_cv_t<std::remove_reference_t<decltype(data(0, 0))>>;
    return FTensor::Tensor2<T *, 2, 2>(
        &data(rr + 0, cc + 0), &data(rr + 0, cc + 1),
 
        &data(rr + 1, cc + 0), &data(rr + 1, cc + 1), ss);
  }
};
 
template <class M> struct GetFTensor2FromArrayRawPtrImpl<3, 3, M> {
  GetFTensor2FromArrayRawPtrImpl() = delete;
  inline static auto get(M &data, const size_t rr, const size_t cc,
                         const int ss = 0) {
    using T = std::remove_cv_t<std::remove_reference_t<decltype(data(0, 0))>>;
    return FTensor::Tensor2<T *, 3, 3>(
        &data(rr + 0, cc + 0), &data(rr + 0, cc + 1), &data(rr + 0, cc + 2),
        &data(rr + 1, cc + 0), &data(rr + 1, cc + 1), &data(rr + 1, cc + 2),
        &data(rr + 2, cc + 0), &data(rr + 2, cc + 1), &data(rr + 2, cc + 2),
        ss);
  }
};
 
template <int DIM1, int DIM2, int S>
inline auto getFTensor2FromArray(MatrixDouble &data, const size_t rr,
                                 const size_t cc = 0) {
  return GetFTensor2FromArrayImpl<DIM1, DIM2, S, MatrixDouble>::get(data, rr,
                                                                    cc);
}
 
template <int DIM1, int DIM2>
inline auto getFTensor2FromArray(MatrixDouble &data, const size_t rr,
                                 const size_t cc, const int ss) {
  return GetFTensor2FromArrayRawPtrImpl<DIM1, DIM2, MatrixDouble>::get(data, rr,
                                                                       cc, ss);
}
 
#ifdef WITH_ADOL_C
template <int DIM1, int DIM2, int S>
inline auto getFTensor2FromArray(MatrixADouble &data, const size_t rr) {
  return GetFTensor2FromArrayImpl<DIM1, DIM2, S, MatrixADouble>::get(data, rr);
}
#endif
/** @} */
 
/** @name Tensor2_symmetric functions
 *  Symmetric rank-2 tensor construction and conversion helpers.
 *  @{
 */
template <int Tensor_Dim, int S, class DL, class M>
struct GetFTensor2SymmetricFromMatImpl {};
 
template <int Tensor_Dim, int S, class M>
struct GetFTensor2SymmetricFromMatImpl<
    Tensor_Dim, S, DataLayoutTraits<DataLayout::GaussByCoeffs>, M> {
  static inline auto get(M &data, int rr = 0, int cc = 0) {
    // Preserve const matrix elements in the non-owning tensor view.
    using T = std::remove_reference_t<decltype(data(0, 0))>;
    constexpr int symm_size = (Tensor_Dim * (Tensor_Dim + 1)) / 2;
    constexpr int stride = S == -1 ? symm_size : S;
    static_assert(stride % symm_size == 0,
                  "getFTensor2SymmetricFromMat: stride should be a multiple of "
                  "symm_size");
 
#ifndef NDEBUG
    if (data.size2() != symm_size)
      THROW_MESSAGE(
          "getFTensor2SymmetricFromMat<" +
          boost::lexical_cast<std::string>(Tensor_Dim) +
          ">: wrong size of data matrix, number of columns should be " +
          boost::lexical_cast<std::string>(symm_size) + " but is " +
          boost::lexical_cast<std::string>(data.size2()));
    auto *first = &data(rr + 0, cc + 0);
    auto *last = &data(rr + 0, cc + 0 + symm_size - 1);
    if (last - first != symm_size - 1)
      THROW_MESSAGE("getFTensor2SymmetricFromMat<" +
                    boost::lexical_cast<std::string>(Tensor_Dim) +
                    ">: row slice is not contiguous");
#endif
    return FTensor::Tensor2_symmetric<FTensor::CursorPtr<T *, stride>,
                                      Tensor_Dim>(&data(rr + 0, cc + 0));
  }
};
 
template <int S, class M>
struct GetFTensor2SymmetricFromMatImpl<
    3, S, DataLayoutTraits<DataLayout::CoeffsByGauss>, M> {
  static inline auto get(M &data, int rr = 0, int cc = 0) {
    using T = std::remove_cv_t<std::remove_reference_t<decltype(data(0, 0))>>;
    constexpr int stride = S == -1 ? 1 : S;
#ifndef NDEBUG
    constexpr int symm_size = 6;
    if (data.size1() != symm_size)
      THROW_MESSAGE(
          "getFTensor2SymmetricFromMat<3>: wrong size of data matrix, "
          "number of rows should be 6 but is " +
          boost::lexical_cast<std::string>(data.size1()));
#endif
    return FTensor::Tensor2_symmetric<FTensor::PackPtr<T *, stride>, 3>(
        &data(rr + 0, cc + 0), &data(rr + 1, cc + 0), &data(rr + 2, cc + 0),
        &data(rr + 3, cc + 0), &data(rr + 4, cc + 0), &data(rr + 5, cc + 0));
  }
};
 
template <int S, class M>
struct GetFTensor2SymmetricFromMatImpl<
    2, S, DataLayoutTraits<DataLayout::CoeffsByGauss>, M> {
  static inline auto get(M &data, int rr = 0, int cc = 0) {
    using T = std::remove_cv_t<std::remove_reference_t<decltype(data(0, 0))>>;
    constexpr int stride = S == -1 ? 1 : S;
#ifndef NDEBUG
    constexpr int symm_size = 3;
    if (data.size1() != symm_size)
      THROW_MESSAGE(
          "getFTensor2SymmetricFromMat<2>: wrong size of data matrix, "
          "number of rows should be 3 but is " +
          boost::lexical_cast<std::string>(data.size1()));
#endif
    return FTensor::Tensor2_symmetric<FTensor::PackPtr<T *, stride>, 2>(
        &data(rr + 0, cc + 0), &data(rr + 1, cc + 0), &data(rr + 2, cc + 0));
  }
};
 
/**
 * \brief Get symmetric tensor rank 2 (matrix) form data matrix
 */
template <int Tensor_Dim, int S = -1,
          class DL = DataLayoutTraits<DataLayout::GaussByCoeffs>,
          class M = MatrixDouble>
static inline auto getFTensor2SymmetricFromMat(M &data) {
  return GetFTensor2SymmetricFromMatImpl<Tensor_Dim, S, DL, M>::get(data);
}
 
template <int Tensor_Dim, int S = -1,
          class DL = DataLayoutTraits<DataLayout::GaussByCoeffs>,
          class M = MatrixDouble>
using GetFTensor2SymmetricFromMatType =
    decltype(GetFTensor2SymmetricFromMatImpl<Tensor_Dim, S, DL, M>::get(
        std::declval<M &>(), 0, 0));
 
template <int Tensor_Dim, int S = -1,
          class DL = DataLayoutTraits<DataLayout::GaussByCoeffs>>
static inline auto
getFTensor2SymmetricFromMat(boost::weak_ptr<MatrixDouble> data_ptr) {
#ifndef NDEBUG
  if (data_ptr.expired()) {
    CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL,
                      "Data matrix is not available");
  }
#endif
  return GetFTensor2SymmetricFromMatImpl<Tensor_Dim, S, DL, MatrixDouble>::get(
      *data_ptr.lock());
}
 
template <int Tensor_Dim, int S = -1,
          class DL = DataLayoutTraits<DataLayout::GaussByCoeffs>>
static inline auto
getFTensor2SymmetricFromMat(boost::shared_ptr<MatrixDouble> data_ptr) {
#ifndef NDEBUG
  if (!data_ptr) {
    CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL,
                      "Data matrix is not available");
  }
#endif
  return GetFTensor2SymmetricFromMatImpl<Tensor_Dim, S, DL, MatrixDouble>::get(
      *data_ptr);
}
 
template <int Tensor_Dim> inline auto getFTensor2SymmetricFromPtr(double *ptr) {
  auto matrix_data =
      getMatrixAdaptor(ptr, 1, Tensor_Dim * (Tensor_Dim + 1) / 2);
  return GetFTensor2SymmetricFromMatImpl<
      Tensor_Dim, -1, DataLayoutTraits<DataLayout::GaussByCoeffs>,
      decltype(matrix_data)>::get(matrix_data, 0, 0);
}
 
#ifdef WITH_ADOL_C
template <int Tensor_Dim>
inline auto getFTensor2SymmetricFromPtr(adouble *ptr) {
  auto matrix_data =
      getMatrixAdaptor(ptr, 1, Tensor_Dim * (Tensor_Dim + 1) / 2);
  return GetFTensor2SymmetricFromMatImpl<
      Tensor_Dim, -1, DataLayoutTraits<DataLayout::GaussByCoeffs>,
      decltype(matrix_data)>::get(matrix_data, 0, 0);
}
#endif
/**
 * @brief Make symmetric Tensor2 from pointer, taking lower triangle of matrix
 *
 * @tparam DIM
 * @param ptr
 * @return FTensor::Tensor2<FTensor::PackPtr<double *, DIM1 * DIM2>, DIM1, DIM2>
 */
template <int DIM> inline auto getFTensor2SymmetricLowerFromPtr(double *ptr);
 
template <> inline auto getFTensor2SymmetricLowerFromPtr<3>(double *ptr) {
  return FTensor::Tensor2_symmetric<FTensor::PackPtr<double *, 9>, 3>(
      ptr + HVEC0_0, ptr + HVEC0_1, ptr + HVEC0_2,
 
      ptr + HVEC1_0, ptr + HVEC1_1,
 
      ptr + HVEC2_2);
};
 
template <> inline auto getFTensor2SymmetricLowerFromPtr<2>(double *ptr) {
  return FTensor::Tensor2_symmetric<FTensor::PackPtr<double *, 4>, 2>(
      ptr + 0, ptr + 1, ptr + 3);
};
// Template function for conversion of symmetric tensor to non-symmetric
 
template <int DIM, class T>
inline auto to_non_symm(const FTensor::Tensor2_symmetric<T, DIM> &symm) {
  FTensor::Tensor2<T, DIM, DIM> full;
  FTensor::Index<'i', DIM> i;
  FTensor::Index<'j', DIM> j;
  full(i, j) = symm(i, j);
  return full;
}
/** @} */
 
/** @name Tensor3 functions
 *  Rank-3 tensor construction helpers.
 *  @{
 */
template <int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2, int S, class DL,
          class M>
struct GetFTensor3FromMatImpl;
 
template <int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2, int S, class M>
struct GetFTensor3FromMatImpl<Tensor_Dim0, Tensor_Dim1, Tensor_Dim2, S,
                              DataLayoutTraits<DataLayout::GaussByCoeffs>, M> {
  static inline auto get(M &data, int rr = 0, int cc = 0) {
    using T = std::remove_cv_t<std::remove_reference_t<decltype(data(0, 0))>>;
    constexpr int tensor_size = Tensor_Dim0 * Tensor_Dim1 * Tensor_Dim2;
    constexpr int stride = S == -1 ? tensor_size : S;
    static_assert(
        stride % tensor_size == 0,
        "getFTensor3FromMat: stride should be a multiple of tensor_size");
 
#ifndef NDEBUG
    if (data.size2() != tensor_size) {
      THROW_MESSAGE("getFTensor3FromMat<" +
                    boost::lexical_cast<std::string>(Tensor_Dim0) + "," +
                    boost::lexical_cast<std::string>(Tensor_Dim1) + "," +
                    boost::lexical_cast<std::string>(Tensor_Dim2) +
                    ">: wrong size of columns in data matrix, should be " +
                    boost::lexical_cast<std::string>(tensor_size) + " but is " +
                    boost::lexical_cast<std::string>(data.size2()));
    }
    auto *first = &data(rr + 0, cc + 0);
    auto *last = &data(rr + 0, cc + 0 + tensor_size - 1);
    if (last - first != tensor_size - 1)
      THROW_MESSAGE("getFTensor3FromMat<" +
                    boost::lexical_cast<std::string>(Tensor_Dim0) + "," +
                    boost::lexical_cast<std::string>(Tensor_Dim1) + "," +
                    boost::lexical_cast<std::string>(Tensor_Dim2) +
                    ">: row slice is not contiguous");
#endif
    return FTensor::Tensor3<FTensor::CursorPtr<T *, stride>, Tensor_Dim0,
                            Tensor_Dim1, Tensor_Dim2>{&data(rr + 0, cc + 0)};
  }
};
 
template <int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2, int S, class M>
struct GetFTensor3FromMatImpl<Tensor_Dim0, Tensor_Dim1, Tensor_Dim2, S,
                              DataLayoutTraits<DataLayout::CoeffsByGauss>, M> {
  static inline auto get(M &data, int rr = 0, int cc = 0) {
    using T = std::remove_cv_t<std::remove_reference_t<decltype(data(0, 0))>>;
    constexpr int tensor_size = Tensor_Dim0 * Tensor_Dim1 * Tensor_Dim2;
    constexpr int stride = S == -1 ? 1 : S;
    std::array<T *, tensor_size> ptrs;
#ifndef NDEBUG
    if (data.size1() != tensor_size) {
      THROW_MESSAGE("getFTensor3FromMat<" +
                    boost::lexical_cast<std::string>(Tensor_Dim0) + "," +
                    boost::lexical_cast<std::string>(Tensor_Dim1) + "," +
                    boost::lexical_cast<std::string>(Tensor_Dim2) +
                    ">: wrong size of rows in data matrix, should be " +
                    boost::lexical_cast<std::string>(tensor_size) + " but is " +
                    boost::lexical_cast<std::string>(data.size1()));
    }
#endif
    for (auto i = 0; i != tensor_size; ++i)
      ptrs[i] = &data(rr + i, cc + 0);
    return FTensor::Tensor3<FTensor::PackPtr<T *, stride>, Tensor_Dim0,
                            Tensor_Dim1, Tensor_Dim2>(ptrs);
  }
};
 
/**
 * @brief Get tensor rank 3 (non symmetries) form data matrix
 *
 * @tparam Tensor_Dim0 dimension of first index
 * @tparam Tensor_Dim1 dimension of second index
 * @tparam Tensor_Dim2 dimension of third index
 * @tparam S shift size
 * @tparam DL data layout
 * @tparam M matrix type
 * @param data data container
 * @return FTensor::Tensor3<FTensor::PackPtr<T *, S>, Tensor_Dim0,
                               Tensor_Dim1, Tensor_Dim2>
 */
template <int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2, int S = -1,
          class DL = DataLayoutTraits<DataLayout::GaussByCoeffs>,
          class M = MatrixDouble>
static inline auto getFTensor3FromMat(M &data, int rr = 0, int cc = 0) {
  return GetFTensor3FromMatImpl<Tensor_Dim0, Tensor_Dim1, Tensor_Dim2, S, DL,
                                M>::get(data, rr, cc);
}
 
template <int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2, int S = -1,
          class DL = DataLayoutTraits<DataLayout::GaussByCoeffs>,
          class M = MatrixDouble>
using GetFTensor3FromMatType =
    decltype(GetFTensor3FromMatImpl<Tensor_Dim0, Tensor_Dim1, Tensor_Dim2, S,
                                    DL, M>::get(std::declval<M &>(), 0, 0));
 
template <int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2, int S = -1,
          class DL = DataLayoutTraits<DataLayout::GaussByCoeffs>>
static inline auto getFTensor3FromMat(boost::weak_ptr<MatrixDouble> data_ptr,
                                      int rr = 0, int cc = 0) {
#ifndef NDEBUG
  if (!data_ptr.lock()) {
    CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL,
                      "Data matrix is not available");
  }
#endif
  return GetFTensor3FromMatImpl<Tensor_Dim0, Tensor_Dim1, Tensor_Dim2, S, DL,
                                MatrixDouble>::get(*data_ptr.lock(), rr, cc);
}
 
template <int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2, int S = -1,
          class DL = DataLayoutTraits<DataLayout::GaussByCoeffs>>
static inline auto getFTensor3FromMat(boost::shared_ptr<MatrixDouble> data_ptr,
                                      int rr = 0, int cc = 0) {
#ifndef NDEBUG
  if (!data_ptr) {
    CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL,
                      "Data matrix is not available");
  }
#endif
  return GetFTensor3FromMatImpl<Tensor_Dim0, Tensor_Dim1, Tensor_Dim2, S, DL,
                                MatrixDouble>::get(*data_ptr, rr, cc);
}
/*
 * @brief Make Tensor3 from pointer
 *
 * @tparam DIM
 * @param ptr
 * @return FTensor::Tensor3<FTensor::PackPtr<double *, DIM1 * DIM2* DIM3>, DIM1,
 * DIM2, DIM3>
 */
template <int DIM1, int DIM2, int DIM3>
inline FTensor::Tensor3<FTensor::PackPtr<double *, DIM1 * DIM2 * DIM3>, DIM1,
                        DIM2, DIM3>
getFTensor3FromPtr(double *ptr) {
  static_assert(DIM1 == DIM1 || DIM2 != DIM2 || DIM3 != DIM3,
                "Such getFTensor3FromPtr is not implemented");
}
 
template <>
inline FTensor::Tensor3<FTensor::PackPtr<double *, 12>, 3, 2, 2>
getFTensor3FromPtr<3, 2, 2>(double *ptr) {
  return FTensor::Tensor3<FTensor::PackPtr<double *, 12>, 3, 2, 2>(
      ptr + 0, ptr + 1, ptr + 2, ptr + 3, ptr + 4, ptr + 5, ptr + 6, ptr + 7,
      ptr + 8, ptr + 9, ptr + 10, ptr + 11);
}
 
template <>
inline FTensor::Tensor3<FTensor::PackPtr<double *, 27>, 3, 3, 3>
getFTensor3FromPtr<3, 3, 3>(double *ptr) {
  return FTensor::Tensor3<FTensor::PackPtr<double *, 27>, 3, 3, 3>(
      ptr + 0, ptr + 1, ptr + 2, ptr + 3, ptr + 4, ptr + 5, ptr + 6, ptr + 7,
      ptr + 8, ptr + 9, ptr + 10, ptr + 11, ptr + 12, ptr + 13, ptr + 14,
      ptr + 15, ptr + 16, ptr + 17, ptr + 18, ptr + 19, ptr + 20, ptr + 21,
      ptr + 22, ptr + 23, ptr + 24, ptr + 25, ptr + 26);
}
/** @} */
 
/** @name Dg functions
 *  Helpers for third-order tensors symmetric in the first two indices.
 *  @{
 */
template <int Tensor_Dim01, int Tensor_Dim2, int S, class DL, class M>
struct GetFTensor3DgFromMatImpl {};
 
template <int Tensor_Dim01, int Tensor_Dim2, int S, class M>
struct GetFTensor3DgFromMatImpl<Tensor_Dim01, Tensor_Dim2, S,
                                DataLayoutTraits<DataLayout::GaussByCoeffs>,
                                M> {
  static inline auto get(M &data, int rr = 0, int cc = 0) {
    using T = std::remove_cv_t<std::remove_reference_t<decltype(data(0, 0))>>;
    constexpr int symm_size_01 = (Tensor_Dim01 * (Tensor_Dim01 + 1)) / 2;
    constexpr int dg_size = symm_size_01 * Tensor_Dim2;
    constexpr int stride = S == -1 ? dg_size : S;
    static_assert(
        stride % dg_size == 0,
        "getFTensor3DgFromMat: stride should be a multiple of dg_size");
 
#ifndef NDEBUG
    if (data.size2() != dg_size) {
      THROW_MESSAGE(
          "getFTensor3DgFromMat<" +
          boost::lexical_cast<std::string>(Tensor_Dim01) + ", " +
          boost::lexical_cast<std::string>(Tensor_Dim2) +
          ">: wrong size of data matrix, number of columns should be " +
          boost::lexical_cast<std::string>(dg_size) + " but is " +
          boost::lexical_cast<std::string>(data.size2()));
    }
    auto *first = &data(rr + 0, cc + 0);
    auto *last = &data(rr + 0, cc + 0 + dg_size - 1);
    if (last - first != dg_size - 1)
      THROW_MESSAGE("getFTensor3DgFromMat<" +
                    boost::lexical_cast<std::string>(Tensor_Dim01) + ", " +
                    boost::lexical_cast<std::string>(Tensor_Dim2) +
                    ">: row slice is not contiguous");
#endif
    return FTensor::Dg<FTensor::CursorPtr<T *, stride>, Tensor_Dim01,
                       Tensor_Dim2>{&data(rr + 0, cc + 0)};
  }
};
 
template <int S, class M>
struct GetFTensor3DgFromMatImpl<
    1, 1, S, DataLayoutTraits<DataLayout::CoeffsByGauss>, M> {
  static inline auto get(M &data, int rr = 0, int cc = 0) {
    using T = std::remove_cv_t<std::remove_reference_t<decltype(data(0, 0))>>;
    constexpr int stride = S == -1 ? 1 : S;
 
#ifndef NDEBUG
    constexpr int dg_size = 1;
    if (data.size1() != dg_size)
      THROW_MESSAGE(
          "getFTensor3DgFromMat<1, 1>: wrong size of data matrix, number "
          "of rows should be 1 but is " +
          boost::lexical_cast<std::string>(data.size1()));
#endif
    return FTensor::Dg<FTensor::PackPtr<T *, stride>, 1, 1>{
        &data(rr + 0, cc + 0)};
  }
};
 
template <int S, class M>
struct GetFTensor3DgFromMatImpl<
    2, 2, S, DataLayoutTraits<DataLayout::CoeffsByGauss>, M> {
  static inline auto get(M &data, int rr = 0, int cc = 0) {
    using T = std::remove_cv_t<std::remove_reference_t<decltype(data(0, 0))>>;
    constexpr int stride = S == -1 ? 1 : S;
#ifndef NDEBUG
    constexpr int dg_size = 6;
    if (data.size1() != dg_size) {
      THROW_MESSAGE(
          "getFTensor3DgFromMat<2, 2>: wrong size of data matrix, number "
          "of rows should be 6 but is " +
          boost::lexical_cast<std::string>(data.size1()));
    }
#endif
    return FTensor::Dg<FTensor::PackPtr<T *, stride>, 2, 2>{
        &data(rr + 0, cc + 0), &data(rr + 1, cc + 0), &data(rr + 2, cc + 0),
        &data(rr + 3, cc + 0), &data(rr + 4, cc + 0), &data(rr + 5, cc + 0)};
  }
};
 
template <int S, class M>
struct GetFTensor3DgFromMatImpl<
    3, 3, S, DataLayoutTraits<DataLayout::CoeffsByGauss>, M> {
  static inline auto get(M &data, int rr = 0, int cc = 0) {
    using T = std::remove_cv_t<std::remove_reference_t<decltype(data(0, 0))>>;
    constexpr int stride = S == -1 ? 1 : S;
#ifndef NDEBUG
    constexpr int dg_size = 18;
    if (data.size1() != dg_size) {
      THROW_MESSAGE(
          "getFTensor3DgFromMat<3, 3>: wrong size of data matrix, number "
          "of rows should be 18 but is " +
          boost::lexical_cast<std::string>(data.size1()));
    }
#endif
    return FTensor::Dg<FTensor::PackPtr<T *, stride>, 3, 3>{
        &data(rr + 0, cc + 0),  &data(rr + 1, cc + 0),  &data(rr + 2, cc + 0),
        &data(rr + 3, cc + 0),  &data(rr + 4, cc + 0),  &data(rr + 5, cc + 0),
        &data(rr + 6, cc + 0),  &data(rr + 7, cc + 0),  &data(rr + 8, cc + 0),
        &data(rr + 9, cc + 0),  &data(rr + 10, cc + 0), &data(rr + 11, cc + 0),
        &data(rr + 12, cc + 0), &data(rr + 13, cc + 0), &data(rr + 14, cc + 0),
        &data(rr + 15, cc + 0), &data(rr + 16, cc + 0), &data(rr + 17, cc + 0)};
  }
};
 
/**
 * @brief Get symmetric tensor rank 3  on the first two indices from
 * form data matrix
 *
 * @tparam Tensor_Dim01 dimension of first two indicies
 * @tparam Tensor_Dim2 dimension of last index
 * @tparam S Memory shift
 * @tparam DL data layout
 * @tparam M matrix type
 * @param data data container
 * @return FTensor::Dg<FTensor::PackPtr<T *, S>, Tensor_Dim01, Tensor_Dim2>
 */
template <int Tensor_Dim01, int Tensor_Dim2, int S = -1,
          class DL = DataLayoutTraits<DataLayout::GaussByCoeffs>,
          class M = MatrixDouble>
static inline auto getFTensor3DgFromMat(M &data) {
  return GetFTensor3DgFromMatImpl<Tensor_Dim01, Tensor_Dim2, S, DL, M>::get(
      data);
}
 
template <int Tensor_Dim01, int Tensor_Dim2, int S = -1,
          class DL = DataLayoutTraits<DataLayout::GaussByCoeffs>,
          class M = MatrixDouble>
using GetFTensor3DgFromMatType =
    decltype(GetFTensor3DgFromMatImpl<Tensor_Dim01, Tensor_Dim2, S, DL, M>::get(
        std::declval<M &>(), 0, 0));
 
template <int Tensor_Dim01, int Tensor_Dim2, int S = -1,
          class DL = DataLayoutTraits<DataLayout::GaussByCoeffs>>
static inline auto
getFTensor3DgFromMat(boost::weak_ptr<MatrixDouble> data_ptr) {
#ifndef NDEBUG
  if (!data_ptr.lock()) {
    CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL,
                      "Data matrix is not available");
  }
#endif
  return GetFTensor3DgFromMatImpl<Tensor_Dim01, Tensor_Dim2, S, DL,
                                  MatrixDouble>::get(*(data_ptr.lock()));
}
 
template <int Tensor_Dim01, int Tensor_Dim2, int S = -1,
          class DL = DataLayoutTraits<DataLayout::GaussByCoeffs>>
static inline auto
getFTensor3DgFromMat(boost::shared_ptr<MatrixDouble> data_ptr) {
#ifndef NDEBUG
  if (!data_ptr) {
    CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL,
                      "Data matrix is not available");
  }
#endif
  return GetFTensor3DgFromMatImpl<Tensor_Dim01, Tensor_Dim2, S, DL,
                                  MatrixDouble>::get(*data_ptr);
}
 
template <int Tensor_Dim01, int Tensor_Dim2, int S, class T = double>
static inline auto getFTensor3DgFromPtr(T *ptr) {
  constexpr auto dg_size =
      ((Tensor_Dim01 * (Tensor_Dim01 + 1)) / 2) * Tensor_Dim2;
  auto matrix_data = getMatrixAdaptor(ptr, 1, dg_size);
  return GetFTensor3DgFromMatImpl<Tensor_Dim01, Tensor_Dim2, S,
                                  DataLayoutTraits<DataLayout::GaussByCoeffs>,
                                  decltype(matrix_data)>::get(matrix_data);
}
/** @} */
 
/** @name Tensor4 functions
 *  Rank-4 tensor construction helpers.
 *  @{
 */
template <class T, int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2,
          int Tensor_Dim3, int S, std::size_t... Is>
static inline auto makeFTensor4FromPtrArray(
    const std::array<T *, Tensor_Dim0 * Tensor_Dim1 * Tensor_Dim2 * Tensor_Dim3>
        &ptrs,
    std::index_sequence<Is...>) {
  return FTensor::Tensor4<FTensor::PackPtr<T *, S>, Tensor_Dim0, Tensor_Dim1,
                          Tensor_Dim2, Tensor_Dim3>(ptrs[Is]...);
}
 
template <int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2, int Tensor_Dim3,
          int S, class DL, class M>
struct GetFTensor4FromMatImpl {
  static inline auto get(M &data, int rr = 0, int cc = 0) {
    using T = std::remove_cv_t<std::remove_reference_t<decltype(data(0, 0))>>;
    constexpr int tensor_size =
        Tensor_Dim0 * Tensor_Dim1 * Tensor_Dim2 * Tensor_Dim3;
    constexpr int stride =
        S == -1 ? (DL::isGaussByCoeffs ? tensor_size : 1) : S;
    if constexpr (DL::isGaussByCoeffs) {
      static_assert(
          stride % tensor_size == 0,
          "getFTensor4FromMat: stride should be a multiple of tensor_size");
 
#ifndef NDEBUG
      if (data.size2() != tensor_size) {
        THROW_MESSAGE("getFTensor4FromMat<" +
                      boost::lexical_cast<std::string>(Tensor_Dim0) + "," +
                      boost::lexical_cast<std::string>(Tensor_Dim1) + "," +
                      boost::lexical_cast<std::string>(Tensor_Dim2) + "," +
                      boost::lexical_cast<std::string>(Tensor_Dim3) +
                      ">: wrong size of columns in data matrix, should be " +
                      boost::lexical_cast<std::string>(tensor_size) +
                      " but is " +
                      boost::lexical_cast<std::string>(data.size2()));
      }
      auto *first = &data(rr + 0, cc + 0);
      auto *last = &data(rr + 0, cc + 0 + tensor_size - 1);
      if (last - first != tensor_size - 1) {
        THROW_MESSAGE("getFTensor4FromMat<" +
                      boost::lexical_cast<std::string>(Tensor_Dim0) + "," +
                      boost::lexical_cast<std::string>(Tensor_Dim1) + "," +
                      boost::lexical_cast<std::string>(Tensor_Dim2) + "," +
                      boost::lexical_cast<std::string>(Tensor_Dim3) +
                      ">: row slice is not contiguous");
      }
#endif
      return FTensor::Tensor4<FTensor::CursorPtr<T *, stride>, Tensor_Dim0,
                              Tensor_Dim1, Tensor_Dim2, Tensor_Dim3>{
          &data(rr + 0, cc + 0)};
    } else if constexpr (DL::isCoeffsByGauss) {
      std::array<T *, tensor_size> ptrs;
#ifndef NDEBUG
      if (data.size1() != tensor_size) {
        THROW_MESSAGE("getFTensor4FromMat<" +
                      boost::lexical_cast<std::string>(Tensor_Dim0) + "," +
                      boost::lexical_cast<std::string>(Tensor_Dim1) + "," +
                      boost::lexical_cast<std::string>(Tensor_Dim2) + "," +
                      boost::lexical_cast<std::string>(Tensor_Dim3) +
                      ">: wrong size of rows in data matrix, should be " +
                      boost::lexical_cast<std::string>(tensor_size) +
                      " but is " +
                      boost::lexical_cast<std::string>(data.size1()));
      }
#endif
      for (auto i = 0; i != tensor_size; ++i)
        ptrs[i] = &data(rr + i, cc + 0);
      return makeFTensor4FromPtrArray<T, Tensor_Dim0, Tensor_Dim1, Tensor_Dim2,
                                      Tensor_Dim3, stride>(
          ptrs, std::make_index_sequence<tensor_size>{});
    } else {
      static_assert(!std::is_same<M, M>::value,
                    "Unsupported data layout for getFTensor4FromMatImpl");
    }
  }
};
 
/**
 * @brief Get tensor rank 4 (non symmetric) form data matrix
 *
 * @tparam Tensor_Dim0 dimension of first index
 * @tparam Tensor_Dim1 dimension of second index
 * @tparam Tensor_Dim2 dimension of third index
 * @tparam Tensor_Dim3 dimension of fourth index
 * @tparam S Memory shift
 * @tparam DL Data layout
 * @tparam M Matrix type
 * @param data data container
 * @return FTensor tensor rank 4 view matching the selected data layout
 */
template <int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2, int Tensor_Dim3,
          int S = -1, class DL = DataLayoutTraits<DataLayout::GaussByCoeffs>,
          class M = MatrixDouble>
static inline auto getFTensor4FromMat(M &data) {
  return GetFTensor4FromMatImpl<Tensor_Dim0, Tensor_Dim1, Tensor_Dim2,
                                Tensor_Dim3, S, DL, M>::get(data);
}
 
template <int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2, int Tensor_Dim3,
          int S = -1, class DL = DataLayoutTraits<DataLayout::GaussByCoeffs>,
          class M = MatrixDouble>
using GetFTensor4FromMatType =
    decltype(GetFTensor4FromMatImpl<Tensor_Dim0, Tensor_Dim1, Tensor_Dim2,
                                    Tensor_Dim3, S, DL,
                                    M>::get(std::declval<M &>(), 0, 0));
 
template <int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2, int Tensor_Dim3,
          int S = -1, class DL = DataLayoutTraits<DataLayout::GaussByCoeffs>>
static inline auto getFTensor4FromMat(boost::weak_ptr<MatrixDouble> data_ptr) {
#ifndef NDEBUG
  if (!data_ptr.lock()) {
    CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL,
                      "Data matrix is not available");
  }
#endif
  return GetFTensor4FromMatImpl<Tensor_Dim0, Tensor_Dim1, Tensor_Dim2,
                                Tensor_Dim3, S, DL,
                                MatrixDouble>::get(*data_ptr.lock());
}
 
template <int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2, int Tensor_Dim3,
          int S = -1, class DL = DataLayoutTraits<DataLayout::GaussByCoeffs>>
static inline auto
getFTensor4FromMat(boost::shared_ptr<MatrixDouble> data_ptr) {
#ifndef NDEBUG
  if (!data_ptr) {
    CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL,
                      "Data matrix is not available");
  }
#endif
  return GetFTensor4FromMatImpl<Tensor_Dim0, Tensor_Dim1, Tensor_Dim2,
                                Tensor_Dim3, S, DL,
                                MatrixDouble>::get(*data_ptr);
}
 
template <int DIM1, int DIM2, int DIM3, int DIM4, int S = -1, class T = double>
inline auto getFTensor4FromPtr(T *ptr) {
  auto matrix_data = getMatrixAdaptor(ptr, 1, DIM1 * DIM2 * DIM3 * DIM4);
  return GetFTensor4FromMatImpl<DIM1, DIM2, DIM3, DIM4, S,
                                DataLayoutTraits<DataLayout::GaussByCoeffs>,
                                decltype(matrix_data)>::get(matrix_data);
}
/** @} */
 
/** @name Ddg functions
 *  Helpers for fourth-order tensors symmetric in index pairs.
 *  @{
 */
template <int Tensor_Dim01, int Tensor_Dim23, int S, class DL, class M>
struct GetFTensor4DdgFromMatImpl {};
 
template <int Tensor_Dim01, int Tensor_Dim23, int S, class M>
struct GetFTensor4DdgFromMatImpl<Tensor_Dim01, Tensor_Dim23, S,
                                 DataLayoutTraits<DataLayout::GaussByCoeffs>,
                                 M> {
  static inline auto get(M &data, int rr = 0, int cc = 0) {
    using T = std::remove_cv_t<std::remove_reference_t<decltype(data(0, 0))>>;
    constexpr int symm_size_01 = (Tensor_Dim01 * (Tensor_Dim01 + 1)) / 2;
    constexpr int symm_size_23 = (Tensor_Dim23 * (Tensor_Dim23 + 1)) / 2;
    constexpr int ddg_size = symm_size_01 * symm_size_23;
    constexpr int stride = S == -1 ? ddg_size : S;
    static_assert(
        stride % ddg_size == 0,
        "Stride should be a multiple of the size of the symmetric tensor");
 
#ifndef NDEBUG
    if (data.size2() != ddg_size) {
      THROW_MESSAGE(
          "getFTensor4DdgFromMat<" +
          boost::lexical_cast<std::string>(Tensor_Dim01) + ", " +
          boost::lexical_cast<std::string>(Tensor_Dim23) +
          ">: wrong size of data matrix, number of columns should be " +
          boost::lexical_cast<std::string>(ddg_size) + " but is " +
          boost::lexical_cast<std::string>(data.size2()));
    }
    auto *first = &data(rr + 0, cc + 0);
    auto *last = &data(rr + 0, cc + 0 + ddg_size - 1);
    if (last - first != ddg_size - 1)
      THROW_MESSAGE("getFTensor4DdgFromMat<" +
                    boost::lexical_cast<std::string>(Tensor_Dim01) + ", " +
                    boost::lexical_cast<std::string>(Tensor_Dim23) +
                    ">: row slice is not contiguous");
#endif
    return FTensor::Ddg<FTensor::CursorPtr<T *, stride>, Tensor_Dim01,
                        Tensor_Dim23>{&data(rr + 0, cc + 0)};
  }
};
 
template <int S, class M>
struct GetFTensor4DdgFromMatImpl<
    1, 1, S, DataLayoutTraits<DataLayout::CoeffsByGauss>, M> {
  static inline auto get(M &data, int rr = 0, int cc = 0) {
    using T = std::remove_cv_t<std::remove_reference_t<decltype(data(0, 0))>>;
    constexpr int stride = S == -1 ? 1 : S;
#ifndef NDEBUG
    constexpr int ddg_size = 1;
    if (data.size1() != ddg_size)
      THROW_MESSAGE(
          "getFTensor4DdgFromMat<1, 1>: wrong size of data matrix, number "
          "of rows should be 1 but is " +
          boost::lexical_cast<std::string>(data.size1()));
#endif
    return FTensor::Ddg<FTensor::PackPtr<T *, stride>, 1, 1>{
        &data(rr + 0, cc + 0)};
  }
};
 
template <int S, class M>
struct GetFTensor4DdgFromMatImpl<
    2, 2, S, DataLayoutTraits<DataLayout::CoeffsByGauss>, M> {
  static inline auto get(M &data, int rr = 0, int cc = 0) {
    using T = std::remove_cv_t<std::remove_reference_t<decltype(data(0, 0))>>;
    constexpr int stride = S == -1 ? 1 : S;
#ifndef NDEBUG
    constexpr int ddg_size = 9;
    if (data.size1() != ddg_size) {
      THROW_MESSAGE(
          "getFTensor4DdgFromMat<2, 2>: wrong size of data matrix, number "
          "of rows should be 9 but is " +
          boost::lexical_cast<std::string>(data.size1()));
    }
#endif
    return FTensor::Ddg<FTensor::PackPtr<T *, stride>, 2, 2>{
        &data(rr + 0, cc + 0), &data(rr + 1, cc + 0), &data(rr + 2, cc + 0),
        &data(rr + 3, cc + 0), &data(rr + 4, cc + 0), &data(rr + 5, cc + 0),
        &data(rr + 6, cc + 0), &data(rr + 7, cc + 0), &data(rr + 8, cc + 0)};
  }
};
 
template <int S, class M>
struct GetFTensor4DdgFromMatImpl<
    3, 3, S, DataLayoutTraits<DataLayout::CoeffsByGauss>, M> {
  static inline auto get(M &data, int rr = 0, int cc = 0) {
    using T = std::remove_cv_t<std::remove_reference_t<decltype(data(0, 0))>>;
    constexpr int stride = S == -1 ? 1 : S;
#ifndef NDEBUG
    constexpr int ddg_size = 36;
    if (data.size1() != ddg_size) {
      THROW_MESSAGE(
          "getFTensor4DdgFromMat<3, 3>: wrong size of data matrix, number "
          "of rows should be 36 but is " +
          boost::lexical_cast<std::string>(data.size1()));
    }
#endif
    return FTensor::Ddg<FTensor::PackPtr<T *, stride>, 3, 3>{
        &data(rr + 0, cc + 0),  &data(rr + 1, cc + 0),  &data(rr + 2, cc + 0),
        &data(rr + 3, cc + 0),  &data(rr + 4, cc + 0),  &data(rr + 5, cc + 0),
        &data(rr + 6, cc + 0),  &data(rr + 7, cc + 0),  &data(rr + 8, cc + 0),
        &data(rr + 9, cc + 0),  &data(rr + 10, cc + 0), &data(rr + 11, cc + 0),
        &data(rr + 12, cc + 0), &data(rr + 13, cc + 0), &data(rr + 14, cc + 0),
        &data(rr + 15, cc + 0), &data(rr + 16, cc + 0), &data(rr + 17, cc + 0),
        &data(rr + 18, cc + 0), &data(rr + 19, cc + 0), &data(rr + 20, cc + 0),
        &data(rr + 21, cc + 0), &data(rr + 22, cc + 0), &data(rr + 23, cc + 0),
        &data(rr + 24, cc + 0), &data(rr + 25, cc + 0), &data(rr + 26, cc + 0),
        &data(rr + 27, cc + 0), &data(rr + 28, cc + 0), &data(rr + 29, cc + 0),
        &data(rr + 30, cc + 0), &data(rr + 31, cc + 0), &data(rr + 32, cc + 0),
        &data(rr + 33, cc + 0), &data(rr + 34, cc + 0), &data(rr + 35, cc + 0)};
  }
};
 
/**
 * @brief Get symmetric tensor rank 4  on first two and last indices from
 * form data matrix
 *
 * @tparam Tensor_Dim01 dimension of first two indicies
 * @tparam Tensor_Dim23 dimension of second two indicies
 * @tparam Memory shift
 * @tparam M matrix type
 * @param data data container
 * @return FTensor::Ddg<FTensor::PackPtr<T *, S>, Tensor_Dim01, TensorDim23>
 */
template <int Tensor_Dim01, int Tensor_Dim23, int S = -1,
          class DL = DataLayoutTraits<DataLayout::GaussByCoeffs>,
          class M = MatrixDouble>
static inline auto getFTensor4DdgFromMat(M &data) {
  return GetFTensor4DdgFromMatImpl<Tensor_Dim01, Tensor_Dim23, S, DL, M>::get(
      data);
}
 
template <int Tensor_Dim01, int Tensor_Dim23, int S = -1,
          class DL = DataLayoutTraits<DataLayout::GaussByCoeffs>,
          class M = MatrixDouble>
using GetFTensor4DdgFromMatType =
    decltype(GetFTensor4DdgFromMatImpl<Tensor_Dim01, Tensor_Dim23, S, DL,
                                       M>::get(std::declval<M &>(), 0, 0));
 
template <int Tensor_Dim01, int Tensor_Dim23, int S = -1,
          class DL = DataLayoutTraits<DataLayout::GaussByCoeffs>>
static inline auto getFTensor4DdgFromMat(boost::weak_ptr<MatrixDouble> data) {
#ifndef NDEBUG
  if (!data.lock()) {
    CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL,
                      "Data matrix is not available");
  }
#endif
  return GetFTensor4DdgFromMatImpl<Tensor_Dim01, Tensor_Dim23, S, DL,
                                   MatrixDouble>::get(*data.lock());
}
 
template <int Tensor_Dim01, int Tensor_Dim23, int S = -1,
          class DL = DataLayoutTraits<DataLayout::GaussByCoeffs>>
static inline auto
getFTensor4DdgFromMat(boost::shared_ptr<MatrixDouble> data_ptr) {
#ifndef NDEBUG
  if (!data_ptr) {
    CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL,
                      "Data matrix is not available");
  }
#endif
  return GetFTensor4DdgFromMatImpl<Tensor_Dim01, Tensor_Dim23, S, DL,
                                   MatrixDouble>::get(*data_ptr);
}
 
template <int Tensor_Dim01, int Tensor_Dim23, int S, class T = double>
static inline auto getFTensor4DdgFromPtr(T *ptr) {
  constexpr auto symm_size01 = (Tensor_Dim01 * (Tensor_Dim01 + 1)) / 2;
  constexpr auto symm_size23 = (Tensor_Dim23 * (Tensor_Dim23 + 1)) / 2;
  constexpr auto symm_size = symm_size01 * symm_size23;
  auto matrix_data = getMatrixAdaptor(ptr, 1, symm_size);
  return GetFTensor4DdgFromMatImpl<Tensor_Dim01, Tensor_Dim23, S,
                                   DataLayoutTraits<DataLayout::GaussByCoeffs>,
                                   decltype(matrix_data)>::get(matrix_data);
}
/** @} */
 
/** @name Matrix sizeing helpers
 *  Matrix size helpers
 *  @{
*/
 
template <class Tensor, class DL> struct MatrixSizeHelper;
 
template <class T, int S, int Tensor_Dim>
struct MatrixSizeHelper<
    FTensor::Tensor1<FTensor::CursorPtr<T *, S>, Tensor_Dim>,
    DataLayoutTraits<DataLayout::GaussByCoeffs>> {
  template <class M>
  static inline auto size(M &m, const size_t nb_integration_pts) {
    if (m.size1() != nb_integration_pts || m.size2() != Tensor_Dim) {
      m.resize(nb_integration_pts, Tensor_Dim, false);
    }
    return [&m]() {
      return GetFTensor1FromMatImpl<Tensor_Dim, S,
                                    DataLayoutTraits<DataLayout::GaussByCoeffs>,
                                    M>::get(m, 0, 0);
    };
  }
  template <class M>
  static inline auto get(M &m,
                         [[maybe_unused]] const size_t nb_integration_pts) {
#ifndef NDEBUG
    if (S == -1 || S == Tensor_Dim) {
      if (m.size1() != nb_integration_pts) {
        CHK_THROW_MESSAGE(
            MOFEM_DATA_INCONSISTENCY,
            "Matrix size is not compatible with the number of "
            "integration points: expected " +
                boost::lexical_cast<std::string>(nb_integration_pts) +
                " but got " + boost::lexical_cast<std::string>(m.size1()));
      }
    }
#endif
    return [&m]() {
      return GetFTensor1FromMatImpl<Tensor_Dim, S,
                                    DataLayoutTraits<DataLayout::GaussByCoeffs>,
                                    M>::get(m, 0, 0);
    };
  }
};
 
template <class T, int S, int Tensor_Dim0, int Tensor_Dim1>
struct MatrixSizeHelper<
    FTensor::Tensor2<FTensor::CursorPtr<T *, S>, Tensor_Dim0, Tensor_Dim1>,
    DataLayoutTraits<DataLayout::GaussByCoeffs>> {
  template <class M>
  static inline auto size(M &m, const size_t nb_integration_pts) {
    if (m.size1() != nb_integration_pts ||
        m.size2() != Tensor_Dim0 * Tensor_Dim1) {
      m.resize(nb_integration_pts, Tensor_Dim0 * Tensor_Dim1, false);
    }
    return [&m]() {
      return GetFTensor2FromMatImpl<Tensor_Dim0, Tensor_Dim1, S,
                                    DataLayoutTraits<DataLayout::GaussByCoeffs>,
                                    M>::get(m, 0, 0);
    };
  }
  template <class M>
  static inline auto get(M &m,
                         [[maybe_unused]] const size_t nb_integration_pts) {
#ifndef NDEBUG
    if (S == -1 || S == Tensor_Dim0 * Tensor_Dim1) {
      if (m.size1() != nb_integration_pts) {
        CHK_THROW_MESSAGE(
            MOFEM_DATA_INCONSISTENCY,
            "Matrix size is not compatible with the number of "
            "integration points: expected " +
                boost::lexical_cast<std::string>(nb_integration_pts) +
                " but got " + boost::lexical_cast<std::string>(m.size1()));
      }
    }
#endif
    return [&m]() {
      return GetFTensor2FromMatImpl<Tensor_Dim0, Tensor_Dim1, S,
                                    DataLayoutTraits<DataLayout::GaussByCoeffs>,
                                    M>::get(m, 0, 0);
    };
  }
};
 
template <class T, int S, int Tensor_Dim>
struct MatrixSizeHelper<
    FTensor::Tensor2_symmetric<FTensor::CursorPtr<T *, S>, Tensor_Dim>,
    DataLayoutTraits<DataLayout::GaussByCoeffs>> {
  template <class M>
  static inline auto size(M &m, const size_t nb_integration_pts) {
    constexpr int symm_size = (Tensor_Dim * (Tensor_Dim + 1)) / 2;
    if (m.size1() != nb_integration_pts || m.size2() != symm_size) {
      m.resize(nb_integration_pts, symm_size, false);
    }
    return [&m]() {
      return GetFTensor2SymmetricFromMatImpl<
          Tensor_Dim, S, DataLayoutTraits<DataLayout::GaussByCoeffs>,
          M>::get(m, 0, 0);
    };
  }
  template <class M>
  static inline auto get(M &m,
                         [[maybe_unused]] const size_t nb_integration_pts) {
#ifndef NDEBUG
    constexpr int symm_size = (Tensor_Dim * (Tensor_Dim + 1)) / 2;
    if (S == -1 || S == symm_size) {
      if (m.size1() != nb_integration_pts) {
        CHK_THROW_MESSAGE(
            MOFEM_DATA_INCONSISTENCY,
            "Matrix size is not compatible with the number of "
            "integration points: expected " +
                boost::lexical_cast<std::string>(nb_integration_pts) +
                " but got " + boost::lexical_cast<std::string>(m.size1()));
      }
    }
#endif
    return [&m]() {
      return GetFTensor2SymmetricFromMatImpl<
          Tensor_Dim, S, DataLayoutTraits<DataLayout::GaussByCoeffs>,
          M>::get(m, 0, 0);
    };
  }
};
 
template <class T, int S, int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2>
struct MatrixSizeHelper<FTensor::Tensor3<FTensor::CursorPtr<T *, S>,
                                         Tensor_Dim0, Tensor_Dim1, Tensor_Dim2>,
                        DataLayoutTraits<DataLayout::GaussByCoeffs>> {
  template <class M>
  static inline auto size(M &m, const size_t nb_integration_pts) {
    constexpr int tensor_size = Tensor_Dim0 * Tensor_Dim1 * Tensor_Dim2;
    if (m.size1() != nb_integration_pts || m.size2() != tensor_size) {
      m.resize(nb_integration_pts, tensor_size, false);
    }
    return [&m]() {
      return GetFTensor3FromMatImpl<Tensor_Dim0, Tensor_Dim1, Tensor_Dim2, S,
                                    DataLayoutTraits<DataLayout::GaussByCoeffs>,
                                    M>::get(m, 0, 0);
    };
  }
  template <class M>
  static inline auto get(M &m,
                         [[maybe_unused]] const size_t nb_integration_pts) {
#ifndef NDEBUG
    constexpr int tensor_size = Tensor_Dim0 * Tensor_Dim1 * Tensor_Dim2;
    if (S == -1 || S == tensor_size) {
      if (m.size1() != nb_integration_pts) {
        CHK_THROW_MESSAGE(
            MOFEM_DATA_INCONSISTENCY,
            "Matrix size is not compatible with the number of "
            "integration points: expected " +
                boost::lexical_cast<std::string>(nb_integration_pts) +
                " but got " + boost::lexical_cast<std::string>(m.size1()));
      }
    }
#endif
    return [&m]() {
      return GetFTensor3FromMatImpl<Tensor_Dim0, Tensor_Dim1, Tensor_Dim2, S,
                                    DataLayoutTraits<DataLayout::GaussByCoeffs>,
                                    M>::get(m, 0, 0);
    };
  }
};
 
template <class T, int S, int Tensor_Dim01, int Tensor_Dim2>
struct MatrixSizeHelper<
    FTensor::Dg<FTensor::CursorPtr<T *, S>, Tensor_Dim01, Tensor_Dim2>,
    DataLayoutTraits<DataLayout::GaussByCoeffs>> {
  template <class M>
  static inline auto size(M &m, const size_t nb_integration_pts) {
    constexpr int dg_size =
        ((Tensor_Dim01 * (Tensor_Dim01 + 1)) / 2) * Tensor_Dim2;
    if (m.size1() != nb_integration_pts || m.size2() != dg_size) {
      m.resize(nb_integration_pts, dg_size, false);
    }
    return [&m]() {
      return GetFTensor3DgFromMatImpl<
          Tensor_Dim01, Tensor_Dim2, S,
          DataLayoutTraits<DataLayout::GaussByCoeffs>, M>::get(m, 0, 0);
    };
  }
  template <class M>
  static inline auto get(M &m,
                         [[maybe_unused]] const size_t nb_integration_pts) {
#ifndef NDEBUG
    constexpr int dg_size =
        ((Tensor_Dim01 * (Tensor_Dim01 + 1)) / 2) * Tensor_Dim2;
    if (S == -1 || S == dg_size) {
      if (m.size1() != nb_integration_pts) {
        CHK_THROW_MESSAGE(
            MOFEM_DATA_INCONSISTENCY,
            "Matrix size is not compatible with the number of "
            "integration points: expected " +
                boost::lexical_cast<std::string>(nb_integration_pts) +
                " but got " + boost::lexical_cast<std::string>(m.size1()));
      }
    }
#endif
    return [&m]() {
      return GetFTensor3DgFromMatImpl<
          Tensor_Dim01, Tensor_Dim2, S,
          DataLayoutTraits<DataLayout::GaussByCoeffs>, M>::get(m, 0, 0);
    };
  }
};
 
template <class T, int S, int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2,
          int Tensor_Dim3>
struct MatrixSizeHelper<
    FTensor::Tensor4<FTensor::CursorPtr<T *, S>, Tensor_Dim0, Tensor_Dim1,
                     Tensor_Dim2, Tensor_Dim3>,
    DataLayoutTraits<DataLayout::GaussByCoeffs>> {
  template <class M>
  static inline auto size(M &m, const size_t nb_integration_pts) {
    constexpr int tensor_size =
        Tensor_Dim0 * Tensor_Dim1 * Tensor_Dim2 * Tensor_Dim3;
    if (m.size1() != nb_integration_pts || m.size2() != tensor_size) {
      m.resize(nb_integration_pts, tensor_size, false);
    }
    return [&m]() {
      return GetFTensor4FromMatImpl<
          Tensor_Dim0, Tensor_Dim1, Tensor_Dim2, Tensor_Dim3, S,
          DataLayoutTraits<DataLayout::GaussByCoeffs>, M>::get(m, 0, 0);
    };
  }
  template <class M>
  static inline auto get(M &m,
                         [[maybe_unused]] const size_t nb_integration_pts) {
#ifndef NDEBUG
    if (S == -1 || S == Tensor_Dim0 * Tensor_Dim1 * Tensor_Dim2 * Tensor_Dim3) {
      if (m.size1() != nb_integration_pts) {
        CHK_THROW_MESSAGE(
            MOFEM_DATA_INCONSISTENCY,
            "Matrix size is not compatible with the number of "
            "integration points: expected " +
                boost::lexical_cast<std::string>(nb_integration_pts) +
                " but got " + boost::lexical_cast<std::string>(m.size1()));
      }
    }
#endif
    return [&m]() {
      return GetFTensor4FromMatImpl<
          Tensor_Dim0, Tensor_Dim1, Tensor_Dim2, Tensor_Dim3, S,
          DataLayoutTraits<DataLayout::GaussByCoeffs>, M>::get(m, 0, 0);
    };
  }
};
 
template <class T, int S, int Tensor_Dim01, int Tensor_Dim23>
struct MatrixSizeHelper<
    FTensor::Ddg<FTensor::CursorPtr<T *, S>, Tensor_Dim01, Tensor_Dim23>,
    DataLayoutTraits<DataLayout::GaussByCoeffs>> {
  template <class M>
  static inline auto size(M &m, const size_t nb_integration_pts) {
    constexpr int ddg_size = ((Tensor_Dim01 * (Tensor_Dim01 + 1)) / 2) *
                             ((Tensor_Dim23 * (Tensor_Dim23 + 1)) / 2);
    if (m.size1() != nb_integration_pts || m.size2() != ddg_size) {
      m.resize(nb_integration_pts, ddg_size, false);
    }
    return [&m]() {
      return GetFTensor4DdgFromMatImpl<
          Tensor_Dim01, Tensor_Dim23, S,
          DataLayoutTraits<DataLayout::GaussByCoeffs>, M>::get(m, 0, 0);
    };
  }
  template <class M>
  static inline auto get(M &m,
                         [[maybe_unused]] const size_t nb_integration_pts) {
#ifndef NDEBUG
    constexpr int ddg_size = ((Tensor_Dim01 * (Tensor_Dim01 + 1)) / 2) *
                             ((Tensor_Dim23 * (Tensor_Dim23 + 1)) / 2);
    if (S == -1 || S == ddg_size) {
 
      if (m.size1() != nb_integration_pts) {
        CHK_THROW_MESSAGE(
            MOFEM_DATA_INCONSISTENCY,
            "Matrix size is not compatible with the number of "
            "integration points: expected " +
                boost::lexical_cast<std::string>(nb_integration_pts) +
                " but got " + boost::lexical_cast<std::string>(m.size1()));
      }
    }
#endif
    return [&m]() {
      return GetFTensor4DdgFromMatImpl<
          Tensor_Dim01, Tensor_Dim23, S,
          DataLayoutTraits<DataLayout::GaussByCoeffs>, M>::get(m, 0, 0);
    };
  }
};
 
/** @} */
 
/** @name Voigt notation helpers
 *  Conversions between tensor views and Voigt-like storage.
 *  @{
 */
// Template functions for conversion to Voigt notation for 2nd order tensors
 
template <int DIM, typename T> inline auto getVoigtVec(T &t_mat) {
  if constexpr (DIM == 3) {
    return std::array<double, 9>{t_mat(0, 0), t_mat(1, 1), t_mat(2, 2),
                                 t_mat(0, 1), t_mat(1, 0), t_mat(0, 2),
                                 t_mat(2, 0), t_mat(1, 2), t_mat(2, 1)};
  } else {
    return std::array<double, 5>{t_mat(0, 0), t_mat(1, 1), 1.0, t_mat(0, 1),
                                 t_mat(1, 0)};
  }
}
 
template <int DIM, typename T> inline auto getVoigtVecSymm(T &t_mat) {
  constexpr double sqr2 = boost::math::constants::root_two<double>();
 
  if constexpr (DIM == 3) {
    return std::array<double, 9>{t_mat(0, 0),
                                 t_mat(1, 1),
                                 t_mat(2, 2),
                                 sqr2 * t_mat(0, 1),
                                 sqr2 * t_mat(0, 2),
                                 sqr2 * t_mat(1, 2),
                                 0.0,
                                 0.0,
                                 0.0};
  } else {
    return std::array<double, 4>{t_mat(0, 0), t_mat(1, 1), 0.0,
                                 sqr2 * t_mat(0, 1)};
  }
}
 
template <typename T>
inline auto getVoigtVecAxisymm(T &t_mat, const double hoop_term) {
  array<double, 5> vec{t_mat(0, 0), t_mat(1, 1), 1. + hoop_term, t_mat(0, 1),
                       t_mat(1, 0)};
 
  return vec;
}
 
template <typename T>
inline auto getVoigtVecSymmAxisymm(T &t_mat, const double hoop_term) {
  constexpr double sqr2 = boost::math::constants::root_two<double>();
 
  array<double, 4> vec_sym{t_mat(0, 0), t_mat(1, 1), hoop_term,
                           sqr2 * t_mat(0, 1)};
  return vec_sym;
}
/** @} */
 
/** @name Linear algebra helpers
 *  LAPACK-backed matrix and eigenvalue utilities.
 *  @{
 */
// list of lapack wrappers
/**
 * @brief compute matrix inverse with lapack dgetri
 *
 * @param mat input square matrix / output inverse matrix
 * @return MoFEMErrorCode
 */
inline MoFEMErrorCode computeMatrixInverse(MatrixDouble &mat) {
  MoFEMFunctionBegin;
 
  const size_t M = mat.size1();
  const size_t N = mat.size2();
 
  if (M != N)
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "The input matrix for inverse computation is not square %zu != %zu",
            M, N);
 
  int *ipv = new int[N];
  int lwork = N * N;
  double *work = new double[lwork];
  int info;
  info = lapack_dgetrf(N, N, &*mat.data().begin(), N, ipv);
  if (info != 0)
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "lapack error info = %d",
            info);
  info = lapack_dgetri(N, &*mat.data().begin(), N, ipv, work, lwork);
  if (info != 0)
    SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
            "lapack error info = %d", info);
 
  delete[] ipv;
  delete[] work;
 
  MoFEMFunctionReturn(0);
}
/**
 * @brief solve linear system with lapack dgesv
 *
 * @param mat input lhs square matrix / output L and U from the factorization
 * @param f input rhs vector / output solution vector
 * @return MoFEMErrorCode
 */
inline MoFEMErrorCode solveLinearSystem(MatrixDouble &mat, VectorDouble &f) {
  MoFEMFunctionBegin;
 
  const size_t M = mat.size1();
  const size_t N = mat.size2();
 
  if (M == 0 || M != N)
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "The input matrix for inverse computation is not square %zu != %zu",
            M, N);
  if (f.size() != M)
    f.resize(M, false);
 
  const int nrhs = 1;
  int info;
  int *ipiv = new int[M];
  info = lapack_dgesv(M, nrhs, &*mat.data().begin(), M, ipiv,
                      &*f.data().begin(), M);
  if (info != 0) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
            "error lapack solve dgesv info = %d", info);
  }
 
  delete[] ipiv;
  MoFEMFunctionReturn(0);
}
 
/**
 * @brief Solve linear system of equations using Lapack
 *
 * @param mat
 * @param f
 * @return MoFEMErrorCode
 */
inline MoFEMErrorCode solveLinearSystem(const MatrixDouble &mat,
                                        VectorDouble &f) {
  MoFEMFunctionBegin;
  // copy matrix since on output lapack returns factorisation
  auto mat_copy = mat;
  CHKERR solveLinearSystem(mat_copy, f);
  MoFEMFunctionReturn(0);
}
 
/** check if two doubles are approximately equal */
template <typename T = double> inline bool isEq(const T a, const T b) {
  const double eps = 100 * std::numeric_limits<T>::epsilon();
  const auto abs_max =
      std::max(static_cast<T>(1), std::max(std::abs(a), std::abs(b)));
  return std::abs(a - b) <= eps * abs_max;
}
 
template <int DIM, typename T = double>
inline auto convertToShallow(FTensor::Tensor1<T, DIM> &vec) {
  return getVectorAdaptor(&vec(0), DIM);
}
 
template <int DIM1, int DIM2, typename T = double>
inline auto convertToShallow(FTensor::Tensor2<T, DIM1, DIM2> &mat) {
  return getMatrixAdaptor(&mat(0, 0), DIM1, DIM2);
}
 
/**
 * @brief Get the Uniq Nb objects
 * 
 * @tparam DIM 
 * @param ptr 
 * @return auto 
 */
template <int DIM, typename T = double>
inline auto getUniqNb(VectorShallowArrayAdaptor<T> vec, FTensor::Number<DIM>) {
  std::array<double, DIM> tmp;
  std::copy(vec.data().begin(), vec.data().begin() + DIM, tmp.begin());
  std::sort(tmp.begin(), tmp.end());
  std::size_t n_unique = 1;
  for (int i = 1; i < DIM; ++i) {
    if (!isEq(tmp[i], tmp[i - 1])) {
      ++n_unique;
    }
  }
  return n_unique;
}
 
template <int DIM, typename T = double>
inline auto sortEigenVals(VectorShallowArrayAdaptor<T> eig,
                          MatrixShallowArrayAdaptor<T> eigen_vec,
                          FTensor::Number<DIM>) {
  if constexpr (DIM == 2) {
    MoFEMFunctionReturnHot(0);
  }
  static_assert(DIM == 2 || DIM == 3, "That works only for DIM = 3");
  MoFEMFunctionBegin;
  int i = 0, j = 1, k = 2;
  if (isEq<T>(eig(0), eig(1))) {
    i = 0;
    j = 2;
    k = 1;
  } else if (isEq<T>(eig(0), eig(2))) {
    i = 0;
    j = 1;
    k = 2;
  } else if (isEq<T>(eig(1), eig(2))) {
    i = 1;
    j = 0;
    k = 2;
  }
  const std::array<T, 9> eigen_vec_c{
      eigen_vec(i, 0), eigen_vec(i, 1), eigen_vec(i, 2),
      eigen_vec(j, 0), eigen_vec(j, 1), eigen_vec(j, 2),
      eigen_vec(k, 0), eigen_vec(k, 1), eigen_vec(k, 2)};
  const std::array<T, 3> eig_c{eig(i), eig(j), eig(k)};
  std::copy(eigen_vec_c.begin(), eigen_vec_c.end(), eigen_vec.data().begin());
  std::copy(eig_c.begin(), eig_c.end(), eig.data().begin());
  MoFEMFunctionReturn(0);
}
 
template <int DIM, typename T = double>
inline auto getUniqNb(FTensor::Tensor1<T, DIM> vec) {
  return getUniqNb<DIM, T>(convertToShallow(vec), FTensor::Number<DIM>{});
}
 
template <int DIM, typename T = double>
inline auto sortEigenVals(FTensor::Tensor1<T, DIM> &eig,
                          FTensor::Tensor2<T, DIM, DIM> &eigen_vec) {
  return sortEigenVals<DIM, T>(convertToShallow(eig),
                               convertToShallow(eigen_vec),
                               FTensor::Number<DIM>{});
}
 
/**
 * @brief compute eigenvalues of a symmetric matrix using lapack dsyev
 *
 * LAPACK `dsyev` returns eigenvectors as columns of the column-major matrix
 * `Q`. Since MoFEM passes the raw buffer of `eigen_vec` directly to LAPACK
 * and then accesses the same buffer with C++ indexing, `eigen_vec(i, j)`
 * stores `Q(j, i)`. Therefore, in MoFEM convention, row `i` of `eigen_vec`
 * stores eigenvector `i`, and `Q = eigen_vec^T` in the usual linear algebra
 * notation with eigenvectors in columns.
 *
 * @param mat input symmetric matrix
 * @param eig output eigen values sorted
 * @param eigen_vec output matrix of eigen vectors in MoFEM convention
 * @return MoFEMErrorCode
 */
inline MoFEMErrorCode computeEigenValuesSymmetric(const MatrixDouble &mat,
                                                  VectorDouble &eig,
                                                  MatrixDouble &eigen_vec) {
  MoFEMFunctionBegin;
 
  const size_t M = mat.size1();
  const size_t N = mat.size2();
 
  if (M == 0 || M != N)
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "The input matrix for eigen value computation is not square %zu "
            "!= %zu",
            M, N);
  if (eig.size() != M)
    eig.resize(M, false);
 
  eigen_vec = mat;
  const int n = M;
  const int lda = M;
  const int size = (M + 2) * M;
  int lwork = size;
  double *work = new double[size];
 
  if (lapack_dsyev('V', 'U', n, &*eigen_vec.data().begin(), lda,
                   &*eig.data().begin(), work, lwork) > 0)
    SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
            "The algorithm failed to compute eigenvalues.");
 
  delete[] work;
  MoFEMFunctionReturn(0);
}
 
/**
 * @brief compute eigenvalues of a symmetric matrix using lapack dsyev
 *
 * LAPACK `dsyev` returns eigenvectors as columns of the column-major matrix
 * `Q`. Since MoFEM passes the raw buffer of `eigen_vec` directly to LAPACK
 * and then accesses the same buffer with C++ indexing, `eigen_vec(i, j)`
 * stores `Q(j, i)`. Therefore, in MoFEM convention, row `i` of `eigen_vec`
 * stores eigenvector `i`, and `Q = eigen_vec^T` in the usual linear algebra
 * notation with eigenvectors in columns.
 *
 * @tparam DIM
 * @param eigen_vec input / output DIM x DIM matrix of eigen vectors in MoFEM
 *        convention
 * @param eig output eigen values sorted
 * @return MoFEMErrorCode
 */
template <int DIM, typename T1, typename T2>
inline MoFEMErrorCode
computeEigenValuesSymmetric(FTensor::Tensor2<T1, DIM, DIM> &eigen_vec,
                            FTensor::Tensor1<T2, DIM> &eig) {
  MoFEMFunctionBegin;
 
  const int n = DIM;
  const int lda = DIM;
  const int lwork = (DIM + 2) * DIM;
  std::array<double, (DIM + 2) * DIM> work;
 
  if (lapack_dsyev('V', 'U', n, &eigen_vec(0, 0), lda, &eig(0), work.data(),
                   lwork) > 0)
    SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
            "The algorithm failed to compute eigenvalues.");
  MoFEMFunctionReturn(0);
}
 
/**
 * @brief compute eigenvalues of a symmetric tensor using lapack dsyev
 *
 * LAPACK `dsyev` returns eigenvectors as columns of the column-major matrix
 * `Q`. Since MoFEM passes the raw buffer of `eigen_vec` directly to LAPACK
 * and then accesses the same buffer with C++ indexing, `eigen_vec(i, j)`
 * stores `Q(j, i)`. Therefore, in MoFEM convention, row `i` of `eigen_vec`
 * stores eigenvector `i`, and `Q = eigen_vec^T` in the usual linear algebra
 * notation with eigenvectors in columns.
 *
 * @tparam DIM
 * @param mat input tensor pointer of size DIM x DIM
 * @param eig output eigen values sorted
 * @param eigen_vec output matrix of eigen vectors in MoFEM convention
 * @return MoFEMErrorCode
 */
template <int DIM, typename T1, typename T2, typename T3>
inline MoFEMErrorCode
computeEigenValuesSymmetric(const FTensor::Tensor2_symmetric<T1, DIM> &mat,
                            FTensor::Tensor1<T2, DIM> &eig,
                            FTensor::Tensor2<T3, DIM, DIM> &eigen_vec) {
  MoFEMFunctionBegin;
  for (int ii = 0; ii != DIM; ii++)
    for (int jj = 0; jj != DIM; jj++)
      eigen_vec(ii, jj) = mat(ii, jj);
 
  CHKERR computeEigenValuesSymmetric(eigen_vec, eig);
 
  MoFEMFunctionReturn(0);
}
 
/** @} */
 
/** @name Tensor algebra helpers
 *  Determinant and inverse routines for tensor-valued objects.
 *  @{
 */
/**
 * @brief Calculate the determinant of a 3x3 matrix or a tensor of rank 2
 *
 * @tparam T
 * @param t
 * @return double
 */
template <typename T> static inline auto determinantTensor3by3(T &t) {
  return t(0, 0) * t(1, 1) * t(2, 2) + t(1, 0) * t(2, 1) * t(0, 2) +
         t(2, 0) * t(0, 1) * t(1, 2) - t(0, 0) * t(2, 1) * t(1, 2) -
         t(2, 0) * t(1, 1) * t(0, 2) - t(1, 0) * t(0, 1) * t(2, 2);
}
 
/**
 * @brief Calculate the determinant of a 2x2 matrix or a tensor of rank 2
 *
 * @tparam T
 * @param t
 * @return double
 */
template <typename T> static inline auto determinantTensor2by2(T &t) {
  return t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0);
}
 
template <typename T, int DIM> struct DeterminantTensorImpl;
 
template <typename T> struct DeterminantTensorImpl<T, 3> {
  static inline auto get(T &t) { return determinantTensor3by3(t); }
};
 
template <typename T> struct DeterminantTensorImpl<T, 2> {
  static auto get(T &t) { return determinantTensor2by2(t); }
};
 
/**
 * @brief Calculate the determinant of a tensor of rank DIM
 */
template <typename T, int DIM>
static inline auto determinantTensor(FTensor::Tensor2<T, DIM, DIM> &t) {
  return DeterminantTensorImpl<FTensor::Tensor2<T, DIM, DIM>, DIM>::get(t);
}
 
/**
 * @brief Calculate the determinant of a tensor of rank DIM
 */
template <typename T, int DIM>
static inline auto determinantTensor(FTensor::Tensor2_symmetric<T, DIM> &t) {
  return DeterminantTensorImpl<FTensor::Tensor2_symmetric<T, DIM>, DIM>::get(t);
}
 
/**
 * \brief Calculate inverse of tensor rank 2 at integration points
 
 */
template <int Tensor_Dim, class T, class L, class A>
inline MoFEMErrorCode invertTensor3by3(ublas::matrix<T, L, A> &jac_data,
                                       ublas::vector<T, A> &det_data,
                                       ublas::matrix<T, L, A> &inv_jac_data) {
  MoFEMFunctionBeginHot;
  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
          "Specialization for this template not yet implemented");
  MoFEMFunctionReturnHot(0);
}
 
template <>
inline MoFEMErrorCode
invertTensor3by3<3, double, ublas::row_major, DoubleAllocator>(
    MatrixDouble &jac_data, VectorDouble &det_data, MatrixDouble &inv_jac_data);
 
/**
 * \brief Calculate determinant 3 by 3
 
 */
template <class T1, class T2>
inline MoFEMErrorCode determinantTensor3by3(T1 &t, T2 &det) {
  MoFEMFunctionBeginHot;
  det = determinantTensor3by3(t);
  MoFEMFunctionReturnHot(0);
}
 
/**
 * \brief Calculate determinant 2 by 2
 
 */
template <class T1, class T2>
inline MoFEMErrorCode determinantTensor2by2(T1 &t, T2 &det) {
  MoFEMFunctionBeginHot;
  det = determinantTensor2by2(t);
  MoFEMFunctionReturnHot(0);
}
 
/**
 * \brief Calculate matrix inverse 3 by 3
 
 */
template <class T1, class T2, class T3>
inline MoFEMErrorCode invertTensor3by3(T1 &t, T2 &det, T3 &inv_t) {
  MoFEMFunctionBeginHot;
  const auto a = t(0, 0), b = t(0, 1), c = t(0, 2);
  const auto d = t(1, 0), e = t(1, 1), f = t(1, 2);
  const auto g = t(2, 0), h = t(2, 1), i = t(2, 2);
  const auto inv_det = 1.0 / det;
  inv_t(0, 0) = (e * i - f * h) * inv_det;
  inv_t(0, 1) = (c * h - b * i) * inv_det;
  inv_t(0, 2) = (b * f - c * e) * inv_det;
  inv_t(1, 0) = (f * g - d * i) * inv_det;
  inv_t(1, 1) = (a * i - c * g) * inv_det;
  inv_t(1, 2) = (c * d - a * f) * inv_det;
  inv_t(2, 0) = (d * h - e * g) * inv_det;
  inv_t(2, 1) = (b * g - a * h) * inv_det;
  inv_t(2, 2) = (a * e - b * d) * inv_det;
  MoFEMFunctionReturnHot(0);
}
 
/**
 * \brief Calculate matrix inverse 2 by 2
 
 */
template <class T1, class T2, class T3>
inline MoFEMErrorCode invertTensor2by2(T1 &t, T2 &det, T3 &inv_t) {
  MoFEMFunctionBeginHot;
  const auto a = t(0, 0);
  const auto b = t(0, 1);
  const auto c = t(1, 0);
  const auto d = t(1, 1);
  inv_t(0, 0) = d / det;
  inv_t(0, 1) = -b / det;
  inv_t(1, 0) = -c / det;
  inv_t(1, 1) = a / det;
  MoFEMFunctionReturnHot(0);
}
 
/**
 * \brief Calculate matrix inverse, overload for symmetric tensor
 
 */
template <typename T1, typename T2, typename T3>
inline MoFEMErrorCode
invertTensor3by3(FTensor::Tensor2_symmetric<T1, 3> &t, T2 &det,
                 FTensor::Tensor2_symmetric<T3, 3> &inv_t) {
  MoFEMFunctionBeginHot;
  const auto a = t(0, 0);
  const auto b = t(0, 1);
  const auto c = t(0, 2);
  const auto d = t(1, 1);
  const auto e = t(1, 2);
  const auto f = t(2, 2);
  const auto inv_det = 1.0 / det;
  inv_t(0, 0) = (d * f - e * e) * inv_det;
  inv_t(0, 1) = (c * e - b * f) * inv_det;
  inv_t(0, 2) = (b * e - c * d) * inv_det;
  inv_t(1, 1) = (a * f - c * c) * inv_det;
  inv_t(1, 2) = (b * c - a * e) * inv_det;
  inv_t(2, 2) = (a * d - b * b) * inv_det;
  MoFEMFunctionReturnHot(0);
}
 
/**
 * \brief Calculate matrix inverse, overload for symmetric tensor
 
 */
template <typename T1, typename T2, typename T3>
inline MoFEMErrorCode
invertTensor2by2(FTensor::Tensor2_symmetric<T1, 2> &t, T2 &det,
                 FTensor::Tensor2_symmetric<T3, 2> &inv_t) {
  MoFEMFunctionBeginHot;
  const auto a = t(0, 0);
  const auto b = t(0, 1);
  const auto d = t(1, 1);
  const auto inv_det = 1.0 / det;
  inv_t(0, 0) = d * inv_det;
  inv_t(0, 1) = -b * inv_det;
  inv_t(1, 1) = a * inv_det;
  MoFEMFunctionReturnHot(0);
}
 
template <typename T1, typename T2, typename T3, int DIM>
struct InvertTensorImpl;
 
template <typename T1, typename T2, typename T3>
struct InvertTensorImpl<T1, T2, T3, 3> {
  inline static MoFEMErrorCode invert(T1 &t, T2 &det, T3 &inv_t) {
    return invertTensor3by3(t, det, inv_t);
  }
};
 
template <typename T1, typename T2, typename T3>
struct InvertTensorImpl<T1, T2, T3, 2> {
  inline static MoFEMErrorCode invert(T1 &t, T2 &det, T3 &inv_t) {
    return invertTensor2by2(t, det, inv_t);
  }
};
 
template <typename T1, typename T2, typename T3, int DIM>
static inline MoFEMErrorCode
invertTensor(FTensor::Tensor2<T1, DIM, DIM> &t, T2 &det,
             FTensor::Tensor2<T3, DIM, DIM> &inv_t) {
  return InvertTensorImpl<FTensor::Tensor2<T1, DIM, DIM>, T2,
                          FTensor::Tensor2<T3, DIM, DIM>, DIM>::invert(t, det,
                                                                       inv_t);
}
 
template <typename T1, typename T2, typename T3, int DIM>
static inline MoFEMErrorCode
invertTensor(FTensor::Tensor2_symmetric<T1, DIM> &t, T2 &det,
             FTensor::Tensor2_symmetric<T3, DIM> &inv_t) {
  return InvertTensorImpl<FTensor::Tensor2_symmetric<T1, DIM>, T2,
                          FTensor::Tensor2_symmetric<T3, DIM>,
                          DIM>::invert(t, det, inv_t);
}
/** @} */
 
/** @name Mesh and range helpers
 *  Utility types and helpers for mesh entities, ranges, and handles.
 *  @{
 */
/**
 * @brief Extract entity handle form multi-index container
 *
 */
struct RefEntExtractor {
  template <typename Iterator>
  static inline EntityHandle extract(const Iterator &it) {
    return (*it)->getEnt();
  }
};
 
/**
 * @brief Insert ordered mofem multi-index into range
 *
 * \note Inserted range has to be ordered.
 *
 * \code
 * auto hi_rit = refEntsPtr->upper_bound(start);
 * auto hi_rit = refEntsPtr->upper_bound(end);
 * Range to_erase;
 * insertOrdered(to_erase, RefEntExtractor(), rit, hi_rit);
 * \endcode
 *
 * @tparam Iterator
 * @param r
 * @param begin_iter
 * @param end_iter
 * @return moab::Range::iterator
 */
template <typename Extractor, typename Iterator>
moab::Range::iterator insertOrdered(Range &r, Extractor, Iterator begin_iter,
                                    Iterator end_iter) {
  moab::Range::iterator hint = r.begin();
  while (begin_iter != end_iter) {
    size_t j = 0;
    auto bi = Extractor::extract(begin_iter);
    Iterator pj = begin_iter;
    while (pj != end_iter && (bi + j) == Extractor::extract(pj)) {
      ++pj;
      ++j;
    }
    hint = r.insert(hint, bi, bi + (j - 1));
    begin_iter = pj;
  }
  return hint;
};
 
/**
 * @brief Do nothing, used to rebuild database
 *
 */
struct Modify_change_nothing {
  Modify_change_nothing() = default;
  template <typename T> inline void operator()(T &e) {}
};
 
/**
 * @brief Template used to reconstruct multi-index
 *
 * @tparam MI multi-index
 * @tparam Modifier
 * @param mi
 * @param mo
 * @return MoFEMErrorCode
 */
template <typename MI, typename MO = Modify_change_nothing>
inline MoFEMErrorCode reconstructMultiIndex(const MI &mi,
                                            MO &&mo = Modify_change_nothing()) {
  MoFEMFunctionBegin;
  for (auto it = mi.begin(); it != mi.end(); ++it) {
    if (!const_cast<MI &>(mi).modify(it, mo))
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "Houston we have a problem");
  }
  MoFEMFunctionReturn(0);
}
 
struct TempMeshset {
  TempMeshset(moab::Interface &moab) : moab(moab) {
    rval = moab.create_meshset(MESHSET_SET, meshset);
    MOAB_THROW(rval);
  }
  virtual ~TempMeshset() { delete_meshset(); }
  operator EntityHandle() const { return meshset; }
  auto get_ptr() { return &meshset; }
 
private:
  void delete_meshset() {
    rval = moab.delete_entities(&meshset, 1);
    MOAB_THROW(rval);
  }
  EntityHandle meshset;
  moab::Interface &moab;
};
 
/**
 * @brief  Create smart pointer to temporary meshset
 *
 */
inline auto get_temp_meshset_ptr(moab::Interface &moab) {
  return boost::make_shared<TempMeshset>(moab);
};
 
inline auto id_from_handle(const EntityHandle h) {
  return static_cast<EntityID>(h & MB_ID_MASK);
};
 
/**
 * @brief get type from entity handle
 *
 */
inline auto type_from_handle(const EntityHandle h) {
  return static_cast<EntityType>(h >> MB_ID_WIDTH);
};
 
/**
 * @brief get entity handle from type and id
 *
 */
inline auto ent_form_type_and_id(const EntityType type, const EntityID id) {
  return (static_cast<EntityHandle>(type) << MB_ID_WIDTH) | id;
};
 
/**
 * @brief get entity dimension form handle
 *
 */
inline auto dimension_from_handle(const EntityHandle h) {
  return moab::CN::Dimension(type_from_handle(h));
};
 
/**
 * @brief get entity type name from handle
 *
 */
inline auto type_name_from_handle(const EntityHandle h) {
  return moab::CN::EntityTypeName(type_from_handle(h));
};
 
/**
 * @brief get field bit id from bit number
 *
 */
inline auto field_bit_from_bit_number(const int bit_number) {
  return BitFieldId().set(bit_number - 1);
};
 
/**
 * @brief Insert ranges
 *
 * @tparam I
 * @param f
 * @param s
 * @param tester
 * @param inserter
 * @return auto
 */
template <typename I>
auto rangeInserter(const I f, const I s, boost::function<bool(I it)> tester,
                   boost::function<MoFEMErrorCode(I f, I s)> inserter) {
  MoFEMFunctionBegin;
 
  auto first = f;
  while (first != s)
    if (tester(first)) {
 
      auto second = first;
      ++second;
 
      while (second != s) {
        if (tester(second))
          ++second;
        else
          break;
      }
 
      CHKERR inserter(first, second);
 
      first = second;
      if (first != s)
        ++first;
 
    } else {
      ++first;
    }
 
  MoFEMFunctionReturn(0);
}
 
/**
 * @brief Create Array
 *
 * See:
 * <a
 * href="https://stackoverflow.com/questions/50942556/current-status-of-stdmake-array">See
 * stack overflow</a>
 *
 * @tparam Dest
 * @tparam Arg
 * @param arg
 * @return constexpr auto
 */
template <typename Dest = void, typename... Arg>
constexpr auto make_array(Arg &&...arg) {
  if constexpr (std::is_same<void, Dest>::value)
    return std::array<std::common_type_t<std::decay_t<Arg>...>, sizeof...(Arg)>{
        {std::forward<Arg>(arg)...}};
  else
    return std::array<Dest, sizeof...(Arg)>{{std::forward<Arg>(arg)...}};
}
/** @} */
 
/** @name FTensor index aliases
 *  Convenient aliases for commonly used FTensor indices.
 *  @{
 */
template <int DIM> using i_FTIndex = FTensor::Index<'i', DIM>;
template <int DIM> using j_FTIndex = FTensor::Index<'j', DIM>;
template <int DIM> using k_FTIndex = FTensor::Index<'k', DIM>;
template <int DIM> using l_FTIndex = FTensor::Index<'l', DIM>;
template <int DIM> using m_FTIndex = FTensor::Index<'m', DIM>;
template <int DIM> using n_FTIndex = FTensor::Index<'n', DIM>;
template <int DIM> using o_FTIndex = FTensor::Index<'o', DIM>;
template <int DIM> using p_FTIndex = FTensor::Index<'p', DIM>;
template <int DIM> using x_FTIndex = FTensor::Index<'x', DIM>;
template <int DIM> using y_FTIndex = FTensor::Index<'y', DIM>;
template <int DIM> using I_FTIndex = FTensor::Index<'I', DIM>;
template <int DIM> using J_FTIndex = FTensor::Index<'J', DIM>;
template <int DIM> using K_FTIndex = FTensor::Index<'K', DIM>;
template <int DIM> using L_FTIndex = FTensor::Index<'L', DIM>;
template <int DIM> using M_FTIndex = FTensor::Index<'M', DIM>;
template <int DIM> using N_FTIndex = FTensor::Index<'N', DIM>;
template <int DIM> using O_FTIndex = FTensor::Index<'O', DIM>;
template <int DIM> using Z_FTIndex = FTensor::Index<'Z', DIM>;
template <int DIM> using Y_FTIndex = FTensor::Index<'Y', DIM>;
 
#define FTENSOR_INDEX(DIM, I) I##_FTIndex<DIM> I;
 
#define FTENSOR_INDEX_DECL(r, DIM, I) FTENSOR_INDEX(DIM, I)
 
#define FTENSOR_INDEXES(DIM, ...)                                              \
  BOOST_PP_SEQ_FOR_EACH(FTENSOR_INDEX_DECL, DIM,                               \
                        BOOST_PP_VARIADIC_TO_SEQ(__VA_ARGS__))
 
template <typename T> struct parent_of;
/** @} */
 
} // namespace MoFEM
 
#endif //__TEMPLATES_HPP__