#ifndef TETRAHEDRON_TRANSFORMATION_HPP
#define TETRAHEDRON_TRANSFORMATION_HPP

#include <algebra/TinyMatrix.hpp>
#include <algebra/TinyVector.hpp>

class TetrahedronTransformation
{
 public:
  constexpr static size_t Dimension = 3;

 private:
  constexpr static size_t NumberOfPoints = Dimension + 1;

  TinyMatrix<Dimension> m_jacobian;
  TinyVector<Dimension> m_shift;
  double m_jacobian_determinant;

 public:
  PUGS_INLINE
  TinyVector<Dimension>
  operator()(const TinyVector<Dimension>& x) const
  {
    return m_jacobian * x + m_shift;
  }

  PUGS_INLINE
  double
  jacobianDeterminant() const
  {
    return m_jacobian_determinant;
  }

  PUGS_INLINE
  double
  jacobianDeterminant(const TinyVector<3>&) const
  {
    return m_jacobian_determinant;
  }

  PUGS_INLINE
  TetrahedronTransformation(const TinyVector<Dimension>& a,
                            const TinyVector<Dimension>& b,
                            const TinyVector<Dimension>& c,
                            const TinyVector<Dimension>& d)
  {
    for (size_t i = 0; i < Dimension; ++i) {
      m_jacobian(i, 0) = b[i] - a[i];
      m_jacobian(i, 1) = c[i] - a[i];
      m_jacobian(i, 2) = d[i] - a[i];
    }

    m_shift = a;

    m_jacobian_determinant = det(m_jacobian);
  }

  ~TetrahedronTransformation() = default;
};

#endif   // TETRAHEDRON_TRANSFORMATION_HPP