// GenFem - A high-level finite element library
// Copyright (C) 2010-2026 Eric Bechet
//
// See the LICENSE file for license information and contributions.
// Please report all bugs and problems to <bechet@cadxfem.org>.

#ifndef _GEN_TENSOROPS__H_
#define _GEN_TENSOROPS__H_

#include "genTensorInverse.h"
#include "genTensorFullContrProd.h"
#include "genTensorContrProd.h"
#include "genTensorTensProd.h"
#include "genTensorTranspose.h"
#include "genTensorBuild.h"
#include "genTensorAdd.h"
#include <complex>


template<class T> class SqrtOp
{
public :
};

template<class scalar,int N> class SqrtOp<genTensor0<scalar,N> > // only for scalars.
{
public :
  typedef genTensor0<scalar,N> ValType;
  void operator() (const genTensor0<scalar,N> &a, ValType &r) const
  {
    r()=sqrt(a());
  }
};

// new genTensor
template<class scalar,int N> class SqrtOp<genTensor<0,scalar,N> > // only for scalars.
{
public :
  typedef genTensor<0,scalar,N> ValType;
  void operator() (const genTensor<0,scalar,N> &a, ValType &r) const
  {
    r()=sqrt(a());
  }
};


template<class T1,class T2> class CrossProdOp
{
public :
};

template<class scalar,int N1, int N2> class CrossProdOp<genTensor1<scalar,N1>,genTensor1<scalar,N2> > // only for genTensor1.
{
public :
  typedef genTensor1<scalar,3> ValType;
  void operator ()(const genTensor1<scalar,N1> &a, const genTensor1<scalar,N2> &b, ValType &r) const 
  {
    r=crossprod(a,b);
  }
};

// new genTensor
template<class scalar,int N1,int N2> class CrossProdOp<genTensor<1,scalar,N1>,genTensor<1,scalar,N2> > // only for genTensor1.
{
public :
  typedef genTensor<1,scalar,3> ValType;
  void operator ()(const genTensor<1,scalar,N1> &a, const genTensor<1,scalar,N2> &b, ValType &r) const 
  {
    r=crossprod(a,b);
  }
};

template<class T> class TraceOp
{
public :
};

template<class scalar,int N> class TraceOp<genTensor2<scalar,N> > // only for genTensor2.
{
public :
  typedef genTensor0<scalar,N> ValType;
  void operator() (const genTensor2<scalar,N> &a, ValType &r) const
  {
    r()=trace(a);
  }
};

// new genTensor
template<class scalar,int N> class TraceOp<genTensor<2,scalar,N> > // only for genTensor2.
{
public :
  typedef genTensor<0,scalar,N> ValType;
  void operator() (const genTensor<2,scalar,N> &a, ValType &r) const
  {
    r()=trace(a);
  }
};

template<class T> class ConjOpScal // default : scalars are considered as real numbers
{
public :
  void operator() (const T &a, T &r) const
  {
    r=a;
  }
};

template<class T> class ConjOpScal<std::complex<T> > // special case if complex numbers
{
public :
  void operator() (const std::complex<T> &a, std::complex<T> &r) const
  {
    r=std::conj(a);
  }
};


template<class T> class ConjOp
{
public :
};

// here we can take the conjugate of either complex variables or real variables


// new genTensor
template<int order,class scalar,int N> class ConjOp<genTensor<order,scalar,N> >
{
public :
  typedef genTensor<order,scalar,N> ValType;
  void operator() (const genTensor<order,scalar,N> &a, ValType &r) const
  {
    for (int i = 0; i < ValType::Nele; ++i)
      ConjOpScal<scalar>()(a[i],r[i]);
  }
};

template<class scalar,int N> class ConjOp<genTensor0<scalar,N> >
{
public :
  typedef genTensor0<scalar,N> ValType;
  void operator() (const genTensor0<scalar,N> &a, ValType &r) const
  {
    ConjOpScal<scalar>()(a(),r());
  }
};

template<class scalar,int N> class ConjOp<genTensor1<scalar,N> >
{
public :
  typedef genTensor1<scalar,N> ValType;
  void operator() (const genTensor1<scalar,N> &a, ValType &r) const
  {
    for (int i = 0; i < N; ++i)
      ConjOpScal<scalar>()(a(i),r(i));
  }
};

template<class scalar,int N> class ConjOp<genTensor2<scalar,N> >
{
public :
  typedef genTensor2<scalar,N> ValType;
  void operator() (const genTensor2<scalar,N> &a, ValType &r) const
  {
    for (int i = 0; i < N; ++i)
      for (int j = 0; j < N; ++j)
        ConjOpScal<scalar>()(a(i,j),r(i,j));
  }
};

template<class scalar,int N> class ConjOp<genTensor3<scalar,N> >
{
public :
  typedef genTensor3<scalar,N> ValType;
  void operator() (const genTensor3<scalar,N> &a, ValType &r) const
  {
    for (int i = 0; i < N; ++i)
      for (int j = 0; j < N; ++j)
        for (int k = 0; k < N; ++k)
          ConjOpScal<scalar>()(a(i,j,k),r(i,j,k));
  }
};

template<class scalar,int N> class ConjOp<genTensor4<scalar,N> >
{
public :
  typedef genTensor4<scalar,N> ValType;
  void operator() (const genTensor4<scalar,N> &a, ValType &r) const
  {
    for (int i = 0; i < N; ++i)
      for (int j = 0; j < N; ++j)
        for (int k = 0; k < N; ++k)
          for (int l = 0; l < N; ++l)
            ConjOpScal<scalar>()(a(i,j,k,l),r(i,j,k,l));
  }
};


template<class scalar> class ReOpScal // real part
{
public :
  void operator() (const scalar &a, typename ScalarTraits<scalar>::RealType &r) const
  {
      r=std::real(a);
  }
};

template<class T> class ReOp
{
public :
};

// here we can take the real part of either complex variables or real variables


// new genTensor
template<int order,class scalar,int N> class ReOp<genTensor<order,scalar,N> >
{
public :
  typedef genTensor<order,typename ScalarTraits<scalar>::RealType,N> ValType;
  void operator() (const genTensor<order,scalar,N> &a, ValType &r) const
  {
    for (int i = 0; i < ValType::Nele; ++i)
      ReOpScal<scalar>()(a[i],r[i]);
  }
};

template<class scalar,int N> class ReOp<genTensor0<scalar,N> >
{
public :
  typedef genTensor0<typename ScalarTraits<scalar>::RealType,N> ValType;
  void operator() (const genTensor0<scalar,N> &a, ValType &r) const
  {
    ReOpScal<scalar>()(a(),r());
  }
};

template<class scalar,int N> class ReOp<genTensor1<scalar,N> >
{
public :
  typedef genTensor1<typename ScalarTraits<scalar>::RealType,N> ValType;
  void operator() (const genTensor1<scalar,N> &a, ValType &r) const
  {
    for (int i = 0; i < N; ++i)
      ReOpScal<scalar>()(a(i),r(i));
  }
};

template<class scalar,int N> class ReOp<genTensor2<scalar,N> >
{
public :
  typedef genTensor2<typename ScalarTraits<scalar>::RealType,N> ValType;
  void operator() (const genTensor2<scalar,N> &a, ValType &r) const
  {
    for (int i = 0; i < N; ++i)
      for (int j = 0; j < N; ++j)
        ReOpScal<scalar>()(a(i,j),r(i,j));
  }
};

template<class scalar,int N> class ReOp<genTensor3<scalar,N> >
{
public :
  typedef genTensor3<typename ScalarTraits<scalar>::RealType,N> ValType;
  void operator() (const genTensor3<scalar,N> &a, ValType &r) const
  {
    for (int i = 0; i < N; ++i)
      for (int j = 0; j < N; ++j)
        for (int k = 0; k < N; ++k)
          ReOpScal<scalar>()(a(i,j,k),r(i,j,k));
  }
};

template<class scalar,int N> class ReOp<genTensor4<scalar,N> >
{
public :
  typedef genTensor4<typename ScalarTraits<scalar>::RealType,N> ValType;
  void operator() (const genTensor4<scalar,N> &a, ValType &r) const
  {
    for (int i = 0; i < N; ++i)
      for (int j = 0; j < N; ++j)
        for (int k = 0; k < N; ++k)
          for (int l = 0; l < N; ++l)
            ReOpScal<scalar>()(a(i,j,k,l),r(i,j,k,l));
  }
};

template<class scalar> class ImOpScal // im part : return zero for real types, and the imaginary part if complex
{
public :
  void operator() (const scalar &a,typename ScalarTraits<scalar>::RealType &r) const
  {
      r=std::imag(a);
  }
};

template<class T> class ImOp
{
public :
};

// here we can take the imaginary part of either complex variables or real variables

// new genTensor
template<int order,class scalar,int N> class ImOp<genTensor<order,scalar,N> >
{
public :
  typedef genTensor<order,typename ScalarTraits<scalar>::RealType,N> ValType;
  void operator() (const genTensor<order,scalar,N> &a, ValType &r) const
  {
    for (int i = 0; i < ValType::Nele; ++i)
      ImOpScal<scalar>()(a[i],r[i]);
  }
};

template<class scalar,int N> class ImOp<genTensor0<scalar,N> >
{
public :
  typedef genTensor0<typename ScalarTraits<scalar>::RealType,N> ValType;
  void operator() (const genTensor0<scalar,N> &a, ValType &r) const
  {
    ImOpScal<scalar>()(a(),r());
  }
};

template<class scalar,int N> class ImOp<genTensor1<scalar,N> >
{
public :
  typedef genTensor1<typename ScalarTraits<scalar>::RealType,N> ValType;
  void operator() (const genTensor1<scalar,N> &a, ValType &r) const
  {
    for (int i = 0; i < N; ++i)
      ImOpScal<scalar>()(a(i),r(i));
  }
};

template<class scalar,int N> class ImOp<genTensor2<scalar,N> >
{
public :
  typedef genTensor2<typename ScalarTraits<scalar>::RealType,N> ValType;
  void operator() (const genTensor2<scalar,N> &a, ValType &r) const
  {
    for (int i = 0; i < N; ++i)
      for (int j = 0; j < N; ++j)
        ImOpScal<scalar>()(a(i,j),r(i,j));
  }
};

template<class scalar,int N> class ImOp<genTensor3<scalar,N> >
{
public :
  typedef genTensor3<typename ScalarTraits<scalar>::RealType,N> ValType;
  void operator() (const genTensor3<scalar,N> &a, ValType &r) const
  {
    for (int i = 0; i < N; ++i)
      for (int j = 0; j < N; ++j)
        for (int k = 0; k < N; ++k)
          ImOpScal<scalar>()(a(i,j,k),r(i,j,k));
  }
};

template<class scalar,int N> class ImOp<genTensor4<scalar,N> >
{
public :
  typedef genTensor4<typename ScalarTraits<scalar>::RealType,N> ValType;
  void operator() (const genTensor4<scalar,N> &a, ValType &r) const
  {
    for (int i = 0; i < N; ++i)
      for (int j = 0; j < N; ++j)
        for (int k = 0; k < N; ++k)
          for (int l = 0; l < N; ++l)
            ImOpScal<scalar>()(a(i,j,k,l),r(i,j,k,l));
  }
};


template<class scalar> class AbsOpScal // modulus (or absolute value)
{
public :
  void operator() (const scalar &a, typename ScalarTraits<scalar>::RealType &r) const
  {
      r=std::abs(a);
  }
};

template<class T> class AbsOp
{
public :
};

// here we can take the modulus of either complex variables or real variables

// new genTensor
template<int order,class scalar,int N> class AbsOp<genTensor<order,scalar,N> >
{
public :
  typedef genTensor<order,typename ScalarTraits<scalar>::RealType,N> ValType;
  void operator() (const genTensor<order,scalar,N> &a, ValType &r) const
  {
    for (int i = 0; i < ValType::Nele; ++i)
      AbsOpScal<scalar>()(a[i],r[i]);
  }
};

template<class scalar,int N> class AbsOp<genTensor0<scalar,N> >
{
public :
  typedef genTensor0<typename ScalarTraits<scalar>::RealType,N> ValType;
  void operator() (const genTensor0<scalar,N> &a, ValType &r) const
  {
    AbsOpScal<scalar>()(a(),r());
  }
};

template<class scalar,int N> class AbsOp<genTensor1<scalar,N> >
{
public :
  typedef genTensor1<typename ScalarTraits<scalar>::RealType,N> ValType;
  void operator() (const genTensor1<scalar,N> &a, ValType &r) const
  {
    for (int i = 0; i < N; ++i)
      AbsOpScal<scalar>()(a(i),r(i));
  }
};

template<class scalar,int N> class AbsOp<genTensor2<scalar,N> >
{
public :
  typedef genTensor2<typename ScalarTraits<scalar>::RealType,N> ValType;
  void operator() (const genTensor2<scalar,N> &a, ValType &r) const
  {
    for (int i = 0; i < N; ++i)
      for (int j = 0; j < N; ++j)
        AbsOpScal<scalar>()(a(i,j),r(i,j));
  }
};

template<class scalar,int N> class AbsOp<genTensor3<scalar,N> >
{
public :
  typedef genTensor3<typename ScalarTraits<scalar>::RealType,N> ValType;
  void operator() (const genTensor3<scalar,N> &a, ValType &r) const
  {
    for (int i = 0; i < N; ++i)
      for (int j = 0; j < N; ++j)
        for (int k = 0; k < N; ++k)
          AbsOpScal<scalar>()(a(i,j,k),r(i,j,k));
  }
};

template<class scalar,int N> class AbsOp<genTensor4<scalar,N> >
{
public :
  typedef genTensor4<typename ScalarTraits<scalar>::RealType,N> ValType;
  void operator() (const genTensor4<scalar,N> &a, ValType &r) const
  {
    for (int i = 0; i < N; ++i)
      for (int j = 0; j < N; ++j)
        for (int k = 0; k < N; ++k)
          for (int l = 0; l < N; ++l)
            AbsOpScal<scalar>()(a(i,j,k,l),r(i,j,k,l));
  }
};


template<class scalar> class ArgOpScal // argument :  for real numbers = 0
{
public :
  void operator() (const scalar &a, typename ScalarTraits<scalar>::RealType &r) const
  {
      r=std::arg(a);
  }
};

template<class T> class ArgOp
{
public :
};

// here we can take the modulus of either complex variables or real variables

// new genTensor
template<int order,class scalar,int N> class ArgOp<genTensor<order,scalar,N> >
{
public :
  typedef genTensor<order,typename ScalarTraits<scalar>::RealType,N> ValType;
  void operator() (const genTensor<order,scalar,N> &a, ValType &r) const
  {
    for (int i = 0; i < ValType::Nele; ++i)
      ArgOpScal<scalar>()(a[i],r[i]);
  }
};

template<class scalar,int N> class ArgOp<genTensor0<scalar,N> >
{
public :
  typedef genTensor0<typename ScalarTraits<scalar>::RealType,N> ValType;
  void operator() (const genTensor0<scalar,N> &a, ValType &r) const
  {
    ArgOpScal<scalar>()(a(),r());
  }
};

template<class scalar,int N> class ArgOp<genTensor1<scalar,N> >
{
public :
  typedef genTensor1<typename ScalarTraits<scalar>::RealType,N> ValType;
  void operator() (const genTensor1<scalar,N> &a, ValType &r) const
  {
    for (int i = 0; i < N; ++i)
      ArgOpScal<scalar>()(a(i),r(i));
  }
};

template<class scalar,int N> class ArgOp<genTensor2<scalar,N> >
{
public :
  typedef genTensor2<typename ScalarTraits<scalar>::RealType,N> ValType;
  void operator() (const genTensor2<scalar,N> &a, ValType &r) const
  {
    for (int i = 0; i < N; ++i)
      for (int j = 0; j < N; ++j)
        ArgOpScal<scalar>()(a(i,j),r(i,j));
  }
};

template<class scalar,int N> class ArgOp<genTensor3<scalar,N> >
{
public :
  typedef genTensor3<typename ScalarTraits<scalar>::RealType,N> ValType;
  void operator() (const genTensor3<scalar,N> &a, ValType &r) const
  {
    for (int i = 0; i < N; ++i)
      for (int j = 0; j < N; ++j)
        for (int k = 0; k < N; ++k)
          ArgOpScal<scalar>()(a(i,j,k),r(i,j,k));
  }
};

template<class scalar,int N> class ArgOp<genTensor4<scalar,N> >
{
public :
  typedef genTensor4<typename ScalarTraits<scalar>::RealType,N> ValType;
  void operator() (const genTensor4<scalar,N> &a, ValType &r) const
  {
    for (int i = 0; i < N; ++i)
      for (int j = 0; j < N; ++j)
        for (int k = 0; k < N; ++k)
          for (int l = 0; l < N; ++l)
            ArgOpScal<scalar>()(a(i,j,k,l),r(i,j,k,l));
  }
};

#endif // _GEN_TENSOROPS__H_