// 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>.
//
// Initial design: Eric Bechet (rev.999)


#ifndef _GEN_TERM_APPLY__H_
#define _GEN_TERM_APPLY__H_

#include "genTerm.h"
#include "genTraits.h"
#include "genMatrix.h"


template<class T,int n> class genTerm;
template<class T,int n> class diffTerm;

////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// applies an operator (binary) to two genTerms stemming from the same function space, yelding a new genTerm. //
////////////////////////////////////////////////////////////////////////////////////////////////////////////////

template<class T1,class T2,int n1,int n2, class OP> class ApplyTermBinary  : public genTerm<typename OP::ValType,n1>
{
public:
  typedef typename OP::ValType T;
  typedef typename TensorialTraits<T>::ValType ValType;
  typedef typename ContainerTraits<ValType,n1>::Container ContainerValType;
  static const int Nsp=n1;
protected:
  typename genTerm<T1,n1>::ConstHandle gt1;
  typename genTerm<T2,n2>::ConstHandle gt2;
  OP op;

public:
  ApplyTermBinary(const typename genTerm<T1,n1>::ConstHandle &gt1_, const typename genTerm<T2,n2>::ConstHandle &gt2_, const OP &op_) : gt1( gt1_),gt2( gt2_),op(op_)
  {
    if (n2==0) return; // special case
    if (n1==n2)
    {
      for (int i=0;i<n1;++i)
      {
        bool samespace=(gt1->getIncidentSpaceTag(i)==gt2->getIncidentSpaceTag(i));
        if (!samespace)
        {
          std::cerr << "ApplyTerm2 only allowed for genTerms stemming from identical functions spaces" << std::endl;
          throw "ApplyTerm2 only allowed for genTerms stemming from identical functions spaces";
        }
      }
    }
    else 
    {
      std::cerr << "ApplyTerm2 only allowed for genTerms havin the same matricial type" << std::endl;
      throw "ApplyTerm2 only allowed for genTerms havin the same matricial type";
    }
  }
  virtual int getNumKeys(MElement* ele, int k=0) const
  {
    if (n1==0) return 0;
    return gt1->getNumKeys(ele,k);
  }
  virtual void getKeys(MElement* ele, std::vector<Dof> &keys, int k=0) const
  {
    if (n1==0) return;
    gt1->getKeys(ele,keys,k);
  }
  virtual int getIncidentSpaceTag(int k=0) const
  { 
    if (n1==0) return 0;
    return gt1->getIncidentSpaceTag(k);
  }

  virtual void get(MElement* ele, int npts, IntPt* GP, std::vector<ContainerValType> &vvals) const;

  virtual ApplyTermBinary<T1,T2,n1,n2,OP>* clone () const {return new ApplyTermBinary<T1,T2,n1,n2,OP>(gt1,gt2,op);}
};



// version with two template parameter (no user defined template parameter left)
template< template<class TA,class TB> class OP,
          class GT1,
          class GT2 > 
std::shared_ptr<genTerm<typename OP<typename GT1::ValType,typename GT2::ValType>::ValType,GT1::Nsp> >
  Apply(const std::shared_ptr<GT1> &gt1, const std::shared_ptr<GT2> &gt2, OP<typename GT1::ValType,typename GT2::ValType> op=OP<typename GT1::ValType,typename GT2::ValType>())
  {
    return std::shared_ptr<genTerm<typename OP<typename GT1::ValType,typename GT2::ValType>::ValType,GT1::Nsp> >
    (new ApplyTermBinary<typename GT1::ValType,typename GT2::ValType,GT1::Nsp,GT2::Nsp,OP<typename GT1::ValType,typename GT2::ValType> >(gt1,gt2,op));
  }

// version with three template parameter (remains one 'user defined' : T)
template< template<class TA,class TB, class TC> class OP,
          class T,
          class GT1,
          class GT2 > 
std::shared_ptr<genTerm<typename OP<typename GT1::ValType,typename GT2::ValType, T>::ValType,GT1::Nsp> >
  Apply(const std::shared_ptr<GT1> &gt1, const std::shared_ptr<GT2> &gt2, OP<typename GT1::ValType,typename GT2::ValType, T> op=OP<typename GT1::ValType,typename GT2::ValType, T>())
  {
    return std::shared_ptr<genTerm<typename OP<typename GT1::ValType,typename GT2::ValType, T>::ValType,GT1::Nsp> >
    (new ApplyTermBinary<typename GT1::ValType,typename GT2::ValType,GT1::Nsp,GT2::Nsp,OP<typename GT1::ValType,typename GT2::ValType, T> >(gt1,gt2,op));
  }

  

template < class OP,
           class T1,
           class T2> 
void Apply(const genScalar<T1> &cont1, const genScalar<T2> &cont2, genScalar<typename OP::ValType> &res, OP op=OP())
{
  op(cont1(),cont2(),res());
}

template < class OP,
           class T1,
           class T2> 
void Apply(const genVector<T1> &cont1, const genVector<T2> &cont2, genVector<typename OP::ValType> &res, OP op=OP())
{
  res.resize(cont1.size(),false);
  for (int i=0;i<cont1.size();++i)
    op(cont1(i),cont2(i),res(i));
}

template < class OP,
           class T1,
           class T2> 
void Apply(const genVector<T1> &cont1, const genScalar<T2> &cont2, genVector<typename OP::ValType> &res, OP op=OP())
{
  res.resize(cont1.size(),false);
  for (int i=0;i<cont1.size();++i)
    op(cont1(i),cont2(),res(i));
}

template < class OP,
           class T1,
           class T2> 
void Apply(const genMatrix<T1> &cont1, const genMatrix<T2> &cont2, genMatrix<typename OP::ValType> &res, OP op=OP())
{
  res.resize(cont1.size1(), cont1.size2(),false);
  for (int i=0;i<cont1.size1();++i)
    for (int j=0;j<cont1.size2();++j)
      op(cont1(i,j),cont2(i,j),res(i,j));
}

template < class OP,
           class T1,
           class T2> 
void Apply(const genMatrix<T1> &cont1, const genScalar<T2> &cont2, genMatrix<typename OP::ValType> &res, OP op=OP())
{
  res.resize(cont1.size1(), cont1.size2(),false);
  for (int i=0;i<cont1.size1();++i)
    for (int j=0;j<cont1.size2();++j)
      op(cont1(i,j),cont2(),res(i,j));
}


template<class T1,class T2,int n1, int n2,class OP>
void ApplyTermBinary<T1,T2,n1,n2,OP>::get(MElement* ele, int npts, IntPt* GP, std::vector<ContainerValType> &vvals) const
{
  std::vector< typename ContainerTraits<T1,n1>::Container > vvals1(npts);
  std::vector< typename ContainerTraits<T2,n2>::Container > vvals2(npts);
  gt1->get(ele,npts,GP,vvals1);
  gt2->get(ele,npts,GP,vvals2);

  for (int i=0;i<npts;++i)
    Apply(vvals1[i],vvals2[i],vvals[i],op);
}





////////////////////////////////////////////////////////////////////////
// applies an operator (unary) to one genTerm, yelding a new genTerm. //
////////////////////////////////////////////////////////////////////////

template<class T1,int n, class OP> class ApplyTermUnary : public genTerm<typename OP::ValType,n>
{
public:
  typedef typename OP::ValType T;
  typedef typename TensorialTraits<T>::ValType ValType;
  typedef typename ContainerTraits<ValType,n>::Container ContainerValType;
  static const int Nsp=n;
protected:
  typename genTerm<T1,n>::ConstHandle gt1;
  OP op;

public:
  ApplyTermUnary(const typename genTerm<T1,n>::ConstHandle &gt1_,const OP &op_) : gt1( gt1_),op(op_) {}
  virtual int getNumKeys(MElement* ele,int k=0) const
  {
    if (n==0) return 0;
    return gt1->getNumKeys(ele,k);
  }
  virtual void getKeys(MElement* ele, std::vector<Dof> &keys, int k=0) const
  {
    if (n==0) return;
    gt1->getKeys(ele,keys,k);
  }
  virtual int getIncidentSpaceTag(int k=0) const
  { 
    if (n==0) return 0;
    return gt1->getIncidentSpaceTag(k);
  }

  virtual void get(MElement* ele, int npts, IntPt* GP, std::vector<ContainerValType> &vvals) const;

  virtual ApplyTermUnary<T1,n,OP>* clone () const {return new ApplyTermUnary<T1,n,OP>(gt1,op);}
};
/*
template<class T1,class Tp1,int n, class OP> class ApplyTermUnary1  : public genTerm<typename OP::ValType,n>
{
public:
  typedef typename OP::ValType T;
  typedef typename TensorialTraits<T>::ValType ValType;
  typedef typename ContainerTraits<ValType,n>::Container ContainerValType;
  static const int Nsp=n;
protected:
  typename genTerm<T1,n>::ConstHandle gt1;
  typename genTerm<Tp1,0>::ConstHandle p1;
  OP op;

public:
  ApplyTermUnary1(const typename genTerm<T1,n>::ConstHandle &gt1_, const typename genTerm<Tp1,0>::ConstHandle &p1_, const OP &op_) : gt1( gt1_),p1(p1_),op(op_) {}
  virtual int getNumKeys(MElement* ele,int k=0) const
  {
    if (n==0) return 0;
    return gt1->getNumKeys(ele,k);
  }
  virtual void getKeys(MElement* ele, std::vector<Dof> &keys, int k=0) const
  {
    if (n==0) return;
    gt1->getKeys(ele,keys,k);
  }
  virtual int getIncidentSpaceTag(int k=0) const
  { 
    if (n==0) return 0;
    return gt1->getIncidentSpaceTag(k);
  }

  virtual void get(MElement* ele, int npts, IntPt* GP, std::vector<ContainerValType> &vvals) const;

  virtual ApplyTermUnary1<T1,Tp1,n,OP>* clone () const {return new ApplyTermUnary1<T1,Tp1,n,OP>(gt1,p1,op);}
};

*/


// version with only one template parameter (no user defined template parameter left)
template< template<class TA> class OP,
          class GT1> 
std::shared_ptr<genTerm<typename OP<typename GT1::ValType>::ValType,GT1::Nsp> >
  Apply(const std::shared_ptr<GT1> &gt1, OP<typename GT1::ValType> op=OP<typename GT1::ValType>())
  {
    return std::shared_ptr<genTerm<typename OP<typename GT1::ValType>::ValType,GT1::Nsp> >
    (new ApplyTermUnary<typename GT1::ValType,GT1::Nsp,OP<typename GT1::ValType> >(gt1,op));
  }


// version with two template parameter (remains one 'user defined' : T)
template< template<class TA, class TB> class OP,
          class T,class GT1> 
std::shared_ptr<genTerm<typename OP<typename GT1::ValType, T>::ValType,GT1::Nsp> >
  Apply(const std::shared_ptr<GT1> &gt1, OP<typename GT1::ValType, T> op=OP<typename GT1::ValType, T>())
  {
    return std::shared_ptr<genTerm<typename OP<typename GT1::ValType, T>::ValType,GT1::Nsp> >
    (new ApplyTermUnary<typename GT1::ValType,GT1::Nsp,OP<typename GT1::ValType, T> >(gt1,op));
  }

/*
template< template<class TA, class TB> class OP,
          class GT1,class GTp1> 
std::shared_ptr<genTerm<typename OP<typename GT1::ValType, typename GTp1::ValType>::ValType,GT1::Nsp> >
  Apply(const std::shared_ptr<GT1> &gt1,const std::shared_ptr<GTp1> &p1, OP<typename GT1::ValType, typename GTp1::ValType> op=OP<typename GT1::ValType,typename GTp1::ValType>())
  {
    return std::shared_ptr<genTerm<typename OP<typename GT1::ValType, typename GTp1::ValType>::ValType,GT1::Nsp> >
    (new ApplyTermUnary1<typename GT1::ValType,typename GTp1::ValType,GT1::Nsp,OP<typename GT1::ValType, typename GTp1::ValType> >(gt1,p1,op));
  }

  */


// version with three template parameters (remains two 'user defined' : T1, T2)
template< template<class TA, class TB, class TC> class OP,
          class T1, class T2, class GT1> 
std::shared_ptr<genTerm<typename OP<typename GT1::ValType, T1, T2>::ValType,GT1::Nsp> >
  Apply(const std::shared_ptr<GT1> &gt1, OP<typename GT1::ValType, T1, T2> op=OP<typename GT1::ValType, T1, T2>())
  {
    return std::shared_ptr<genTerm<typename OP<typename GT1::ValType, T1, T2>::ValType,GT1::Nsp> >
    (new ApplyTermUnary<typename GT1::ValType,GT1::Nsp,OP<typename GT1::ValType, T1, T2> >(gt1,op));
  }

// no user parameter

template < class OP,
           class T >
void Apply(const genScalar<T> &cont, genScalar<typename OP::ValType> &res, OP op=OP())
{
  op(cont(),res());
}

template < class OP,
           class T > 
void Apply(const genVector<T> &cont, genVector<typename OP::ValType> &res, OP op=OP())
{
  res.resize(cont.size(),false);
  for (int i=0;i<cont.size();++i)
    op(cont(i),res(i));
}

template < class OP,
           class T > 
void Apply(const genMatrix<T> &cont, genMatrix<typename OP::ValType> &res, OP op=OP())
{
  res.resize(cont.size1(), cont.size2(),false);
  for (int i=0;i<cont.size1();++i)
    for (int j=0;j<cont.size2();++j)
      op(cont(i,j),res(i,j));
}

// one user parameter

template < class OP,
           class T , class Tp1> 
void Apply(const genScalar<T> &cont, genScalar<Tp1> &p1, genScalar<typename OP::ValType> &res, OP op=OP())
{
  op(cont(),p1(),res());
}


template < class OP,
           class T , class Tp1> 
void Apply(const genVector<T> &cont, genScalar<Tp1> &p1, genVector<typename OP::ValType> &res, OP op=OP())
{
  res.resize(cont.size(),false);
  for (int i=0;i<cont.size();++i)
    op(cont(i),p1(),res(i));
}



template < class OP,
           class T , class Tp1>
void Apply(const genMatrix<T> &cont, genScalar<Tp1> &p1, genMatrix<typename OP::ValType> &res, OP op=OP())
{
  res.resize(cont.size1(), cont.size2(),false);
  for (int i=0;i<cont.size1();++i)
    for (int j=0;j<cont.size2();++j)
      op(cont(i,j),p1(),res(i,j));
}

// two user parameters

template < class OP,
           class T , class Tp1,class Tp2> 
void Apply(const genScalar<T> &cont, genScalar<Tp1> &p1, genScalar<Tp2> &p2, genScalar<typename OP::ValType> &res, OP op=OP())
{
  op(cont(),p1(),p2(),p2(),res());
}


template < class OP,
           class T , class Tp1,class Tp2> 
void Apply(const genVector<T> &cont, genScalar<Tp1> &p1, genScalar<Tp2> &p2, genVector<typename OP::ValType> &res, OP op=OP())
{
  res.resize(cont.size(),false);
  for (int i=0;i<cont.size();++i)
    op(cont(i),p1(),p2(),res(i));
}


template < class OP,
           class T , class Tp1,class Tp2>
void Apply(const genMatrix<T> &cont, genScalar<Tp1> &p1, genScalar<Tp2> &p2, genMatrix<typename OP::ValType> &res, OP op=OP())
{
  res.resize(cont.size1(), cont.size2(),false);
  for (int i=0;i<cont.size1();++i)
    for (int j=0;j<cont.size2();++j)
      op(cont(i,j),p1(),p2(),res(i,j));
}




template<class T1,int n,class OP> 
void ApplyTermUnary<T1,n,OP>::get(MElement* ele, int npts, IntPt* GP, std::vector<ContainerValType> &vvals) const
{
  std::vector< typename ContainerTraits<T1,n>::Container > vvals1(npts);
  gt1->get(ele,npts,GP,vvals1);
  for (int i=0;i<npts;++i)
    Apply(vvals1[i],vvals[i],op);
}
/*
template<class T1,class Tp1,int n,class OP> 
void ApplyTermUnary1<T1,Tp1,n,OP>::get(MElement* ele, int npts, IntPt* GP, std::vector<ContainerValType> &vvals) const
{
  std::vector< typename ContainerTraits<T1,n>::Container > vvals1(npts);
  std::vector< typename ContainerTraits<Tp1,0>::Container > vvalsp1(npts);
  gt1->get(ele,npts,GP,vvals1);
  p1->get(ele,npts,GP,vvalsp1);
  for (int i=0;i<npts;++i)
    Apply(vvals1[i],vvalsp1[i], vvals[i],op);
}
*/

#endif //GEN_TERM_APPLY