// 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 _GENTENSOR3_H_
#define _GENTENSOR3_H_

#include "genTensorBase.h"
#include "genTensor0.h"
#include "genTensor1.h"
#include "genTensor2.h"
#include <iostream>
#include <iomanip>

template<class scalar=double,int N=3> class genTensor3 : public genTensorBase<N>
{
 protected:
  scalar _val[N*N*N];
 public:
  // operator on the data list
  inline scalar & operator[](int i) {return _val[i];}
  inline scalar operator[](int i) const {return _val[i];}
  // operator on the tensor
  inline int getIndex(int i, int j, int k) const {return genTensorBase<N>::getIndex(i,j,k);}
  inline scalar & operator()(int i, int j, int k) {return _val[getIndex(i,j,k)];}
  inline scalar operator()(int i, int j, int k) const {return _val[getIndex(i,j,k)];}

  // default constructor, null tensor
  genTensor3(const scalar v=scalar())
  {
    for (int i=0; i<N; ++i)
      for (int j=0; j<N; ++j)
        for (int k=0; k<N; ++k)
          _val[getIndex(i,j,k)] = v;
  }
  genTensor3(const genTensor3<scalar,N> &other)
  {
    for (int i=0; i<N*N*N; ++i) _val[i] = other._val[i];
  }
  genTensor3(const scalar* array)
  {
    for (int i=0; i<N*N*N; ++i) _val[i]=array[i];
  }
  genTensor3(const std::vector<scalar> &array)
  {
    for (int i=0; i<N*N*N; ++i)
      _val[i] = array[i];
  }
  template <int N2> genTensor3(const genTensor3<scalar,N2> &in)
  {
    int i=0;
    const int Nmin = (N<N2)?N:N2;
    for (; i<Nmin; ++i)
    {
      int j=0;
      for (; j<Nmin; ++j)
      {
        int k=0;
        for (; k<Nmin; ++k)
          _val[getIndex(i,j,k)] = in(i,j,k);
        for (; k<N; ++k)
          _val[getIndex(i,j,k)] = scalar();
      }
      for (; j<N; ++j)
        for (int k=0; k<N; ++k) 
          _val[getIndex(i,j,k)] = scalar();
    }
    for (; i<N; ++i)
      for (int j=0; j<N; ++j)
        for (int k=0; k<N; ++k)
          _val[getIndex(i,j,k)] = scalar();
  }

  scalar normalize()
  {
    scalar nrm = norm(*this);
    if(nrm)
      for (int i=0; i<N*N*N; ++i)
        _val[i] /= nrm;
    return nrm;
  }

  genTensor3<scalar,N> transpose (int n, int m) const
  {
    genTensor3<scalar,N> ithis;
    if ((n==0 && m==1) || (n==1 && m==0))
    {
      for (int i=0; i<N; ++i)
        for (int j=0; j<N; ++j)
          for (int k=0; k<N; ++k)
            ithis(i,j,k) = (*this)(j,i,k);
      return ithis;
    }
    if ((n==0 && m==2) || (n==2 && m==0))
    {
      for (int i=0; i<N; ++i)
        for (int j=0; j<N; ++j)
          for (int k=0; k<N; ++k)
            ithis(i,j,k) = (*this)(k,j,i);
      return ithis;
    }
    if ((n==1 && m==2) || (n==2 && m==1))
    {
      for (int i=0; i<N; ++i)
        for (int j=0; j<N; ++j)
          for (int k=0; k<N; ++k)
            ithis(i,j,k) = (*this)(i,k,j);
      return ithis;
    }
    return ithis += (*this);
  }

  // operator on the data list
  genTensor3<scalar,N> operator+(const genTensor3<scalar,N> &other) const
  {
    genTensor3<scalar> res(*this);
    for (int i=0; i<N*N*N; ++i) res._val[i] += other._val[i];
    return res;
  }
  genTensor3<scalar,N> operator-(const genTensor3<scalar,N> &other) const
  {
    genTensor3<scalar> res(*this);
    for (int i=0; i<N*N*N; ++i) res._val[i] -= other._val[i];
    return res;
  }
  genTensor3<scalar,N> & operator+=(const genTensor3<scalar,N> &other)
  {
    for (int i=0; i<N*N*N; ++i) _val[i] += other._val[i];
    return *this;
  }
  genTensor3<scalar,N> & operator-=(const genTensor3<scalar,N> &other)
  {
    for (int i=0; i<N*N*N; ++i) _val[i] -= other._val[i];
    return *this;
  }
//// misuse of the operator *
//   genTensor3<scalar,N> & operator*=(const genTensor3<scalar,N> &other)
//   {
//     for (int i=0; i<N*N*N; ++i) _val[i] *= other._val[i];
//     return *this;
//   }
  genTensor3<scalar,N> & operator*=(const scalar &s)
  {
    for (int i=0; i<N*N*N; ++i) _val[i] *= s;
    return *this;
  }

  void print(std::string name="") const
  {
    std::cout <<" tensor3 " << name << std::endl << *this;
  }
};

template <class scalar,int N> std::ostream & operator << (std::ostream & output,const genTensor3<scalar,N> & t)
{
  output.precision(5);
  output << std::setiosflags( std::ios::showpos );
  output << std::setiosflags( std::ios::scientific );

  for (int i=0; i<N; ++i){
    for (int j=0; j<N; ++j){
      for (int k=0; k<N; ++k)
        output << t(i,j,k) << " ";
      output << std::endl; }
    output << std::endl; }

  output << std::resetiosflags( std::ios::showpos );
  output << std::resetiosflags( std::ios::scientific );
  return output;
}


template <class scalar,int N> inline scalar dot(const genTensor3<scalar,N> &a, const genTensor3<scalar,N> &b)
{
  scalar prod = scalar();
  for (int i=0; i<N; ++i)
    for (int j=0; j<N; ++j)
      for (int k=0; k<N; ++k)
        prod += a(i,j,k)*b(i,j,k);
  return prod;
}

template<class scalar,int N1,int N2> inline scalar dot(const genTensor3<scalar,N1> &a, const genTensor3<scalar,N2> &b)
{
  scalar prod = scalar();
  const int Nmin=(N1>N2)?N2:N1;
  for (int i=0; i<Nmin; ++i)
    for (int j=0; j<Nmin; ++j)
      for (int k=0; k<Nmin; ++k)
        prod += a(i,j,k)*b(i,j,k);
  return prod;
}

template<class scalar,int N> inline scalar norm(const genTensor3<scalar,N> &t)
{
  return sqrt(dot(t,t));
}


// tensor product
template <class scalar,int N> inline void tensprod(const genTensor1<scalar,N> &a, const genTensor2<scalar,N> &b, genTensor3<scalar,N> &c)
{
  for (int i=0; i<N; ++i)
    for (int j=0; j<N; ++j)
      for (int k=0; k<N; ++k)
        c(i,j,k) = a(i)*b(j,k);
}
template <class scalar,int N> inline void tensprod(const genTensor2<scalar,N> &a, const genTensor1<scalar,N> &b, genTensor3<scalar,N> &c)
{
  for (int i=0; i<N; ++i)
    for (int j=0; j<N; ++j)
      for (int k=0; k<N; ++k)
        c(i,j,k) = a(i,j)*b(k);
}


// full contracted product
template <class scalar,int N> inline genTensor3<scalar,N> operator*(const genTensor3<scalar,N> &t,const scalar s)
{
  genTensor3<scalar,N> val(t);
  val *= s;
  return val;
}

template <class scalar,int N> inline genTensor3<scalar,N> operator*(const scalar s,const genTensor3<scalar,N> &t)
{
  genTensor3<scalar,N> val(t);
  val *= s;
  return val;
}

template <class scalar,int N> inline genTensor3<scalar,N> operator*(const genTensor3<scalar,N> &t,const genTensor0<scalar,N> s)
{
  genTensor3<scalar,N> val(t);
  val *= s;
  return val;
}

template <class scalar,int N> inline genTensor3<scalar,N> operator*(const genTensor0<scalar,N> s,const genTensor3<scalar,N> &t)
{
  genTensor3<scalar,N> val(t);
  val *= s;
  return val;
}

template <class scalar,int N> inline genTensor2<scalar,N> operator*(const genTensor3<scalar,N> &t, const genTensor1<scalar,N> &v)
{
  genTensor2<scalar,N> val;
  for (int i=0; i<N; ++i)
    for (int j=0; j<N; ++j)
      for (int k=0; k<N; ++k)
        val(i,j) += t(i,j,k)*v(k);
  return val;
}

template <class scalar,int N> inline genTensor2<scalar,N> operator*( const genTensor1<scalar,N> &v , const genTensor3<scalar,N> &t)
{
  genTensor2<scalar,N> val;
  for (int i=0; i<N; ++i)
    for (int j=0; j<N; ++j)
      for (int k=0; k<N; ++k)
        val(j,k) += v(i)*t(i,j,k);
  return val;
}

template <class scalar,int N> inline genTensor1<scalar,N> operator*(const genTensor3<scalar,N> &t, const genTensor2<scalar,N> &m)
{
  genTensor1<scalar,N> val;
  for (int i=0; i<N; ++i)
    for (int j=0; j<N; ++j)
      for (int k=0; k<N; ++k)
        val(i) += t(i,j,k)*m(k,j);
  return val;
}

template <class scalar,int N> inline genTensor1<scalar,N> operator*( const genTensor2<scalar,N> &m , const genTensor3<scalar,N> &t)
{
  genTensor1<scalar,N> val;
  for (int i=0; i<N; ++i)
    for (int j=0; j<N; ++j)
      for (int k=0; k<N; ++k)
        val(k) += m(j,i)*t(i,j,k);
  return val;
}

template <class scalar,int N> inline scalar operator*(const genTensor3<scalar,N> &m, const genTensor3<scalar,N> &t)
{
  scalar val = scalar();
  for (int i=0; i<N; ++i)
    for (int j=0; j<N; ++j)
      for (int k=0; k<N; ++k)
        val += m(i,j,k)*t(k,j,i);
  return val;
}

#endif // _GENTENSOR3_H_