//#####################################################################
// Copyright 2013.
// This file is part of PhysBAM whose distribution is governed by the license contained in the accompanying file PHYSBAM_COPYRIGHT.txt.
//#####################################################################
// Class VEC_ID_TENSOR
//##################################################################### 
#ifndef __VEC_ID_TENSOR__
#define __VEC_ID_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=v_i delta_jk
// uu = 1: T_ijk=v_j delta_ik
// uu = 2: T_ijk=v_k delta_ij
template<class T,int uu,int mm,int nn>
class VEC_ID_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<T,mm> v;
    explicit VEC_ID_TENSOR(VECTOR<T,mm> v=(VECTOR<T,mm>())): v(v) {}

    VEC_ID_TENSOR operator-() const
    {return VEC_ID_TENSOR(-v);}

    VEC_ID_TENSOR operator*(T a) const
    {return VEC_ID_TENSOR(v*a);}

    VEC_ID_TENSOR operator/(T a) const
    {return VEC_ID_TENSOR(v/a);}

    template<class RW> void Read(UNTYPED_ISTREAM<RW>& input)
    {Read_Binary(input,v);}
    template<class RW> void Write(UNTYPED_OSTREAM<RW>& output) const
    {Write_Binary(output,v);}
    void Fill_Uniform(RANDOM_NUMBERS<T>& rand,T a,T b)
    {rand.Fill_Uniform(v,a,b);}
};
template<int i,int j,int k,class T,int u,int m,int n> auto get(const VEC_ID_TENSOR<T,u,m,n>& a)
{
    if constexpr(u==0 && j==k) return a.v(i);
    else if constexpr(u==1 && i==k) return a.v(j);
    else if constexpr(u==2 && i==j) return a.v(k);
    else return T{};
}
}
#include <Core/Tensors/PRIMITIVE_TENSORS.h>
#endif
