//#####################################################################
// 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;
    static constexpr int u=uu,m=uu==0?mm:nn,n=uu==1?mm:nn,p=uu==2?mm:nn,um=mm,un=nn;
    static constexpr int type=ten_sym+u;
    VECTOR<SYMMETRIC_MATRIX<T,nn>,mm> x;

    void From(const TENSOR<T,m,n,p>& t)
    {
        if constexpr(u==0)
        {
            for(int i=0;i<m;i++)
                x(i).From_Lower(t.x(i));
        }
        else if constexpr(u==1)
        {
            for(int j=0;j<n;j++)
                for(int i=0;i<m;i++)
                    for(int k=0;k<=i;k++)
                        x(j).Element_Lower(i,k)=t(i,j,k);
        }
        else
        {
            for(int k=0;k<p;k++)
                for(int i=0;i<m;i++)
                    for(int j=0;j<=i;j++)
                        x(k).Element_Lower(i,j)=t(i,j,k);
        }
    }
    
    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);
    }

    T* begin() // for stl
    {return x.begin()->begin();}

    const T* begin() const // for stl
    {return x.begin()->begin();}

    T* end() // for stl
    {return begin()+nn*(nn+1)/2*mm;}

    const T* end() const // for stl
    {return begin()+nn*(nn+1)/2*mm;}
    
    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<int i,int j,int k,class T,int u,int m,int n> void set(SYMMETRIC_TENSOR<T,u,m,n>& a, T x)
{
    if constexpr(u==0) return set<j,k>(a.x(i),x);
    else if constexpr(u==1) return set<i,k>(a.x(j),x);
    else return set<i,j>(a.x(k),x);
}

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
