//#####################################################################
// Copyright 2013.
// This file is part of PhysBAM whose distribution is governed by the license contained in the accompanying file PHYSBAM_COPYRIGHT.txt.
//#####################################################################
// Class SYMMETRIC_TENSOR
//##################################################################### 
#ifndef __SYMMETRIC_TENSOR__
#define __SYMMETRIC_TENSOR__

#include <Core/Math_Tools/cube.h>
#include <Core/Matrices/MATRIX.h>
#include <Core/Matrices/PRIMITIVE_MATRICES.h>
#include <Core/Matrices/SYMMETRIC_MATRIX.h>
#include <Core/Tensors/TENSORS_FORWARD.h>
#include <Core/Vectors/VECTOR.h>
#include <cmath>
namespace PhysBAM{

// uu = unique index (0..2)
// uu=0: T_ijk=T_ikj
// uu=1: T_ijk=T_kji
// uu=2: T_ijk=T_jik
template<class T,int uu,int mm,int nn>
class SYMMETRIC_TENSOR
{
public:
    static const bool is_tensor=true;
    static_assert(uu>=0 && uu<=2);
    typedef T SCALAR;
    enum {u=uu,m=uu==0?mm:nn,n=uu==1?mm:nn,p=uu==2?mm:nn,um=mm,un=nn};
    VECTOR<SYMMETRIC_MATRIX<T,nn>,mm> x;

    SYMMETRIC_TENSOR operator-() const
    {SYMMETRIC_TENSOR t;
    for(int i=0;i<mm;i++) t.x(i)=-x(i);
    return t;}

    SYMMETRIC_TENSOR operator+(const SYMMETRIC_TENSOR& A) const
    {SYMMETRIC_TENSOR s;
    for(int i=0;i<mm;i++) s.x(i)=x(i)+A.x(i);
    return s;}

    SYMMETRIC_TENSOR operator-(const SYMMETRIC_TENSOR& A) const
    {SYMMETRIC_TENSOR s;
    for(int i=0;i<mm;i++) s.x(i)=x(i)-A.x(i);
    return s;}

    SYMMETRIC_TENSOR& operator+=(const SYMMETRIC_TENSOR& A)
    {for(int i=0;i<mm;i++) x(i)+=A.x(i);
    return *this;}

    SYMMETRIC_TENSOR& operator-=(const SYMMETRIC_TENSOR& A)
    {for(int i=0;i<mm;i++) x(i)-=A.x(i);
    return *this;}

    SYMMETRIC_TENSOR operator*(T a) const
    {SYMMETRIC_TENSOR t;
    for(int i=0;i<mm;i++) t.x(i)=x(i)*a;
    return t;}

    SYMMETRIC_TENSOR operator/(T a) const
    {SYMMETRIC_TENSOR t;
    for(int i=0;i<mm;i++) t.x(i)=x(i)/a;
    return t;}

    const T& operator()(int i,int j,int k) const
    {
        if constexpr(u==0) return x(i)(j,k);
        else if constexpr(u==1) return x(j)(i,k);
        else return x(k)(i,j);
    }

    T& operator()(int i,int j,int k)
    {
        if constexpr(u==0) return x(i)(j,k);
        else if constexpr(u==1) return x(j)(i,k);
        else return x(k)(i,j);
    }

    template<class RW> void Read(UNTYPED_ISTREAM<RW>& input)
    {Read_Binary(input,x);}
    template<class RW> void Write(UNTYPED_OSTREAM<RW>& output) const
    {Write_Binary(output,x);}
    void Fill_Uniform(RANDOM_NUMBERS<T>& rand,T a,T b)
    {rand.Fill_Uniform(x,a,b);}
};

template<int i,int j,int k,class T,int u,int m,int n> auto get(const SYMMETRIC_TENSOR<T,u,m,n>& a)
{
    if constexpr(u==0) return get<j,k>(a.x(i));
    else if constexpr(u==1) return get<i,k>(a.x(j));
    else return get<i,j>(a.x(k));
}

template<class T,int u,int m,int n>
SYMMETRIC_TENSOR<T,u,m,n> operator*(T a,const SYMMETRIC_TENSOR<T,u,m,n>& s)
{return s*a;}

template<class T,int m,int n,int p> inline std::ostream&
operator<<(std::ostream& o,const SYMMETRIC_TENSOR<T,m,n,p>& A)
{o<<"("<<A.x<<")";return o;}
}
#include <Core/Tensors/PRIMITIVE_TENSORS.h>
#endif
