TriangleTransformation.hpp

#ifndef TRIANGLE_TRANSFORMATION_HPP
#define TRIANGLE_TRANSFORMATION_HPP

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

template <size_t GivenDimension>
class TriangleTransformation;

template <>
class TriangleTransformation<2>
{
 public:
  constexpr static size_t Dimension = 2;

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

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

  double
  jacobianDeterminant() const
  {
    return m_jacobian_determinant;
  }

  PUGS_INLINE
  TriangleTransformation(const TinyVector<Dimension>& a, const TinyVector<Dimension>& b, const TinyVector<Dimension>& c)
  {
    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_shift = a;

    m_jacobian_determinant = det(m_jacobian);
  }

  ~TriangleTransformation() = default;
};

template <size_t GivenDimension>
class TriangleTransformation
{
 public:
  constexpr static size_t Dimension = GivenDimension;
  static_assert(Dimension == 3, "Triangle transformation is defined in dimension 2 or 3");

 private:
  TinyMatrix<Dimension, 2> m_jacobian;
  TinyVector<Dimension> m_shift;
  double m_area_variation_norm;

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

  double
  areaVariationNorm() const
  {
    return m_area_variation_norm;
  }

  PUGS_INLINE
  TriangleTransformation(const TinyVector<Dimension>& a, const TinyVector<Dimension>& b, const TinyVector<Dimension>& c)
  {
    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_shift = a;

    m_area_variation_norm = l2Norm(crossProduct(b - a, c - a));
  }

  ~TriangleTransformation() = default;
};

#endif   // TRIANGLE_TRANSFORMATION_HPP