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

#include "genTensorBase.h"
#include "genTensor0.h"
#include <iostream>
#include <iomanip>
#include <vector>

template<class scalar=double,int N=3>
class genTensor1 : public genTensorBase<N>
{
 protected:
  scalar _val[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) const {return genTensorBase<N>::getIndex(i);}
  inline scalar & operator()(int i) {return _val[i]; }
  inline scalar operator()(int i) const {return _val[i];}

  // default constructor, null tensor
  genTensor1(const scalar v=scalar())
  {
    for (int i=0; i<N; ++i) _val[i] = v;
  }
  genTensor1(const genTensor1<scalar,N> &other)
  {
    for (int i=0; i<N; ++i) _val[i] = other._val[i];
  }
  genTensor1(const scalar* array)
  {
    for (int i=0; i<N; ++i)
      _val[i] = array[i];
  }
  genTensor1(const std::vector<scalar> &array)
  {
    for (int i=0; i<N; ++i)
      _val[i] = array[i];
  }
  template<int N2> genTensor1(const genTensor1<scalar,N2> &in)
  {
    int i=0;
    const int Nmin = (N<N2)?N:N2;
    for (; i<Nmin; ++i)
      _val[i] = in[i];
    for (; i<N; ++i)
      _val[i] = scalar();
  }
  genTensor1(const scalar x, const scalar y, const scalar z)
  {
    _val[0] = x;
    if (N>1) _val[1] = y;
    if (N>2) _val[2] = z;
  }
  genTensor1(const scalar x, const scalar y)
  {
    _val[0] = x;
    if (N>1) _val[1] = y;
    if (N>2) _val[2] = 0;
  }

  inline scalar x() const { return _val[0]; }
  inline scalar y() const { if (N>1) return _val[1]; else return scalar();}
  inline scalar z() const { if (N>2) return _val[2]; else return scalar();}

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

  void negate()
  {
    for (int i=0;i<N ;++i)
      _val[i] = -_val[i];
  }

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

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

template <class scalar,int N> std::ostream & operator << (std::ostream &output,const genTensor1<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)
    output << t(i) << " "<< 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 genTensor1<scalar,N> &a, const genTensor1<scalar,N> &b)
{
  scalar prod = scalar();
  for (int i=0;i<N;++i)
    prod += a(i)*b(i);
  return prod;
}

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

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


template<class scalar,int N1,int N2> inline genTensor1<scalar,3> crossprod(const genTensor1<scalar,N1> &a, const genTensor1<scalar,N2> &b)
{ return genTensor1<scalar,3>(a.y() * b.z() - a.z() * b.y(),
                              a.z() * b.x() - a.x() * b.z(),
                              a.x() * b.y() - a.y() * b.x()); }

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


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

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

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

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

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

// misuse of the operator *
// template<class scalar,int N> inline genTensor1<scalar,N> operator*(const genTensor1<scalar,N> &v1, const genTensor1<scalar,N> &v2)
// { return genTensor1<scalar,N>(v1) *= v2;}

// template<class scalar,int N> inline genTensor1<scalar> operator+(const genTensor1<scalar,N> &a,const genTensor1<scalar,N> &b)
// { return genTensor1<scalar,N>(a) += b; }

// template<class scalar,int N> inline genTensor1<scalar,N> operator-(const genTensor1<scalar,N> &a,const genTensor1<scalar,N> &b)
// { return genTensor1<scalar,N>(a) -= b; }

#endif // _GENTENSOR1_H_