// 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.6912 in gmsh)


#ifndef _GENMISCTERMS_H_
#define _GENMISCTERMS_H_

#include <vector>
#include <iterator>

#include "genTerm.h"
#include "genFSpace.h"
#include "genTensors.h"
#include "genTraits.h"
#include "genSimpleFunction.h"

template<class T2> class ScalarTermBase;
template<class T2> class LinearTermBase;
template<class T2> class PlusTerm;
template<class T2> class BilinearTermBase;

inline double dot(const double &a, const double &b)
{
  return a * b;
}


inline int delta(int i,int j) {if (i==j) return 1; else return 0;}


//--------------------------------------------------------------------------
// ScalarTermBase
//--------------------------------------------------------------------------

template<class S=double> class  ScalarTermBase : public genTerm<S,0>
{
public :
  typedef typename TensorialTraits<S>::ValType ValType;
  typedef typename ContainerTraits<ValType,0>::Container ContainerValType;
  static const int Nsp=0;

public :
  virtual ~ScalarTermBase() {}
//  virtual void get(MElement* ele, int npts, IntPt* GP, std::vector<ContainerValType> &vvals) const {};
  virtual void get(MElement* ele, int npts, IntPt* GP, ContainerValType &vals) const { genTerm<S,0>::get(ele,npts,GP,vals);}
  virtual ScalarTermBase<S>* clone () const =0;
};

template<class S=double> class ScalarTerm : public ScalarTermBase<S>
{
public :
  typedef typename TensorialTraits<S>::ValType ValType;
  typedef typename ContainerTraits<ValType,0>::Container ContainerValType;
  static const int Nsp=0;
  virtual ~ScalarTerm() {}
};

template<class S=double> class ScalarTermConstant : public ScalarTerm<S>
{
public :
  typedef typename TensorialTraits<S>::ValType ValType;
  typedef typename ContainerTraits<ValType,0>::Container ContainerValType;
  static const int Nsp=0;
protected :
  S cst;
  
public :
  ScalarTermConstant(S val_ = S()) : cst(val_) {}
  virtual ~ScalarTermConstant() {}
  virtual int getNumKeys(MElement* ele,int k=0) const {return 0;}
  virtual void getKeys(MElement* ele, std::vector<Dof> &keys, int k=0) const {}
  virtual int getIncidentSpaceTag(int k=0) const { return 0;}
  virtual void get(MElement* ele, int npts, IntPt* GP, std::vector<ContainerValType> &vvals) const;
  virtual void get(MElement* ele, int npts, IntPt* GP, ContainerValType &vals) const;
  virtual void get(MVertex* ver, ContainerValType &vals) const;
  virtual ScalarTermConstant<S>* clone () const {return new ScalarTermConstant<S>(cst);}
};

//--------------------------------------------------------------------------
// LinearTermBase
//--------------------------------------------------------------------------

template<class S=double> class  LinearTermBase : public genTerm<S,1>
{
public :
  typedef typename TensorialTraits<S>::ValType ValType;
  typedef typename ContainerTraits<ValType,1>::Container ContainerValType;
  static const int Nsp=1;
  
public :
  virtual ~LinearTermBase() {}
  virtual void get(MElement* ele, int npts, IntPt* GP, std::vector<ContainerValType > &vvals) const =0;
  virtual void get(MElement* ele, int npts, IntPt* GP, ContainerValType &vals) const { genTerm<S,1>::get(ele,npts,GP,vals);}
//  virtual LinearTermBase<S>* clone () const = 0;
};

template<class T,class S=double> class LinearTerm : public LinearTermBase<S>
{
public :
  typedef typename TensorialTraits<S>::ValType ValType;
  typedef typename ContainerTraits<ValType,1>::Container ContainerValType;
  static const int Nsp=1;
protected :
  typename genFSpace<T>::ConstHandle space1;
  
public :
  LinearTerm(typename genFSpace<T>::ConstHandle space1_) : space1(space1_) {}
  virtual ~LinearTerm() {}
  virtual int getNumKeys(MElement* ele,int k=0) const {return space1->getNumKeys(ele,k);}
  virtual void getKeys(MElement* ele, std::vector<Dof> &keys, int k=0) const { space1->getKeys(ele,keys,k); }
  virtual int getIncidentSpaceTag(int k=0) const { return space1->getIncidentSpaceTag(k);}
};

template<class T> class LoadTerm : public LinearTerm<T>
{
public :
  typedef typename TensorialTraits<T>::ValType ValType;
  typedef typename ContainerTraits<ValType,1>::Container ContainerValType;
  static const int Nsp=1;
protected :
  double _eqfac;
  genSimpleFunction<typename TensorialTraits<T>::ValType> &Load;

public :
  LoadTerm(typename genFSpace<T>::ConstHandle space1_, genSimpleFunction<typename TensorialTraits<T>::ValType> &Load_, double eqfac = 1.0) :
      LinearTerm<T>(space1_), _eqfac(eqfac), Load(Load_) {}
  virtual ~LoadTerm() {}
  virtual void get(MElement* ele, int npts, IntPt* GP, std::vector<ContainerValType> &vvals) const;
  virtual void get(MElement* ele, int npts, IntPt* GP, ContainerValType &vals) const;
  virtual LoadTerm<T>* clone () const { return new LoadTerm<T>(LinearTerm<T>::space1,Load,_eqfac);}
};

template<class T> class LoadTermOnBorder : public LinearTerm<T>
{
public :
  typedef typename TensorialTraits<T>::ValType ValType;
  typedef typename ContainerTraits<ValType,1>::Container ContainerValType;
  static const int Nsp=1;
protected :
  double _eqfac;
  genSimpleFunction<typename TensorialTraits<T>::ValType> &Load;

public :
  LoadTermOnBorder(typename genFSpace<T>::ConstHandle space1_, genSimpleFunction<typename TensorialTraits<T>::ValType> &Load_, double eqfac = 1.0) :
      LinearTerm<T>(space1_), _eqfac(eqfac), Load(Load_) {}
  virtual ~LoadTermOnBorder() {}
  virtual void get(MElement* ele, int npts, IntPt* GP, std::vector<ContainerValType> &vvals) const;
  virtual void get(MElement* ele, int npts, IntPt* GP, ContainerValType &vals) const;
  virtual LoadTermOnBorder<T>* clone () const { return new LoadTermOnBorder<T>(LinearTerm<T>::space1,Load,_eqfac);}
};


//--------------------------------------------------------------------------
// BilinearTermBase
//--------------------------------------------------------------------------

template<class S=double> class  BilinearTermBase : public genTerm<S,2>
{
public :
  typedef typename TensorialTraits<S>::ValType ValType;
  typedef typename ContainerTraits<ValType,2>::Container ContainerValType;
  static const int Nsp=2;
  
public :
  virtual ~BilinearTermBase() {}
  virtual void get(MElement* ele, int npts, IntPt* GP, std::vector<ContainerValType> &vvals) const = 0;
  virtual void get(MElement* ele, int npts, IntPt* GP, ContainerValType &vals) const { genTerm<S,2>::get(ele,npts,GP,vals);}
  virtual BilinearTermBase<S>* clone () const = 0;
};

template<class T1,class T2> class BilinearTerm : public BilinearTermBase<genTensor0<double> >
{
protected :
  typename genFSpace<T1>::ConstHandle space1;
  typename genFSpace<T2>::ConstHandle space2;
  static const int Nsp=2;
  
public :
  BilinearTerm(typename genFSpace<T1>::ConstHandle space1_, typename genFSpace<T2>::ConstHandle space2_) : space1(space1_), space2(space2_) {}
  virtual ~BilinearTerm() {}
  virtual int getNumKeys(MElement* ele,int k=0) const { if (k==0) return space1->getNumKeys(ele);else return space2->getNumKeys(ele);}
  virtual void getKeys(MElement* ele, std::vector<Dof> &keys, int k=0) const { if (k==0) return space1->getKeys(ele,keys); else return space2->getKeys(ele,keys); }
  virtual int getIncidentSpaceTag(int k=0) const { if (k==0) return space1->getIncidentSpaceTag(); else return space2->getIncidentSpaceTag();}
};

template<class T1,class T2> class LaplaceTerm : public BilinearTerm<T1,T2>
{
public :
  typedef typename TensorialTraitsBinary<T1,T2>::FullContrProdType ValType;
  typedef typename ContainerTraits<ValType,2>::Container ContainerValType;
  static const int Nsp=2;
  
public :
  LaplaceTerm(typename genFSpace<T1>::ConstHandle space1_, typename genFSpace<T2>::ConstHandle space2_) : BilinearTerm<T1,T2>(space1_, space2_)
  {}
  virtual ~LaplaceTerm() {}
  virtual void get(MElement* ele, int npts, IntPt* GP, std::vector<ContainerValType> &vvals) const
  {
    Msg::Error("LaplaceTerm<S1, S2> w/ S1 != S2 not implemented");
  }
  virtual void get(MElement* ele, int npts, IntPt* GP, ContainerValType &vals) const
  {
    Msg::Error("LaplaceTerm<S1, S2> w/ S1 != S2 not implemented");
  }
  virtual void get(MVertex* ver, ContainerValType &vals)
  {
    Msg::Error("LaplaceTerm<S1, S2> w/ S1 != S2 not implemented");
  }
  virtual LaplaceTerm<T1,T2>* clone () const {return new LaplaceTerm<T1,T2>(BilinearTerm<T1,T2>::space1,BilinearTerm<T1,T2>::space2);}
}; // class

template<class T> class LaplaceTerm<T,T> : public BilinearTerm<T,T> // symmetric
{
public :
  typedef typename TensorialTraitsBinary<T,T>::FullContrProdType ValType;
  typedef typename ContainerTraits<ValType,2>::Container ContainerValType;
  static const int Nsp=2;
protected:
  double diffusivity;
  
public :
  LaplaceTerm(typename genFSpace<T>::ConstHandle space1_, double diff = 1) : BilinearTerm<T,T>(space1_, space1_), diffusivity(diff) {}
  virtual ~LaplaceTerm() {}
  virtual void get(MElement* ele, int npts, IntPt* GP, std::vector<ContainerValType> &vvals) const{}
  virtual void get(MElement* ele, int npts, IntPt* GP, ContainerValType &vals) const;
  virtual LaplaceTerm<T,T>* clone () const {return new LaplaceTerm<T,T>(BilinearTerm<T,T>::space1,diffusivity);}
}; // class

class StabilizedLagrangeMultiplierTerm : public BilinearTerm<genTensor1<double>, genTensor0<double> >
{
public :
  typedef TensorialTraits<genTensor0<double> >::ValType ValType;
  typedef ContainerTraits<ValType,2>::Container ContainerValType;
  static const int Nsp=2;
  genTensor1<double>  _d;
  
public :
  StabilizedLagrangeMultiplierTerm(genFSpace<genTensor1<double> >::ConstHandle space1_, genFSpace<genTensor0<double> >::ConstHandle space2_ ,const genTensor1<double> &d) :
      BilinearTerm<genTensor1<double>, genTensor0<double> >(space1_, space2_) { for (int i = 0; i < 3; ++i) _d(i) = d(i); }
  virtual ~StabilizedLagrangeMultiplierTerm() {}
  virtual void get(MElement* ele, int npts, IntPt* GP, std::vector<ContainerValType> &vvals) const;
  virtual void get(MElement* ele, int npts, IntPt* GP, ContainerValType &vals) const;
  virtual StabilizedLagrangeMultiplierTerm* clone () const {return new StabilizedLagrangeMultiplierTerm(BilinearTerm<genTensor1<double>, genTensor0<double> >::space1,BilinearTerm<genTensor1<double>, genTensor0<double> >::space2,_d);}
};

class DispLagrangeMultiplierTerm: public BilinearTerm<genTensor1<double>, genTensor1<double> >
{
public :
  typedef TensorialTraits<genTensor0<double> >::ValType ValType;
  typedef ContainerTraits<ValType,2>::Container ContainerValType;
  static const int Nsp=2;
protected :
  double _eqfac;
  
public :
  DispLagrangeMultiplierTerm(genFSpace<genTensor1<double> >::ConstHandle space1_, genFSpace<genTensor1<double> >::ConstHandle space2_, double eqfac = 1.0) :
      BilinearTerm<genTensor1<double>, genTensor1<double> >(space1_, space2_), _eqfac(eqfac) {}
  virtual ~DispLagrangeMultiplierTerm() {}
  virtual void get(MElement* ele, int npts, IntPt* GP, std::vector<ContainerValType> &vvals) const;
  virtual void get(MElement* ele, int npts, IntPt* GP, ContainerValType &vals) const;
  virtual DispLagrangeMultiplierTerm* clone () const {return new DispLagrangeMultiplierTerm(BilinearTerm<genTensor1<double>,genTensor1<double> >::space1,BilinearTerm<genTensor1<double>, genTensor1<double> >::space2,_eqfac);}
};


template<int n1,int n2> class ElasticTerm : public genTerm<genTensor0<double>,n1+n2>
{
public :
  typedef typename TensorialTraits<genTensor0<double> >::ValType ValType;
  typedef typename ContainerTraits<ValType,n1+n2>::Container ContainerValType;
  static const int Nsp=n1+n2;
  genMatrix<double> H;
protected :
  typename diffTerm<genTensor1<double>,n1>::ConstHandle space1;
  typename diffTerm<genTensor1<double>,n2>::ConstHandle space2;
  bool sym;

public :
  ElasticTerm(typename diffTerm<genTensor1<double>,n1>::ConstHandle space1_, typename diffTerm<genTensor1<double>,n2>::ConstHandle space2_);
  ElasticTerm(typename diffTerm<genTensor1<double>,n1>::ConstHandle space1_);
  virtual ~ElasticTerm() {}
  virtual int getNumKeys(MElement* ele,int k=0) const
  {
    if ((n1==0)&&(n2==0)) return 0;
    if ((n1!=0)&&(n2==0)) return space1->getNumKeys(ele,k);
    if ((n1==0)&&(n2!=0)) return space2->getNumKeys(ele,k);
    if ((n1==1)&&(n2==1))
    {
      if (k==0) return space1->getNumKeys(ele);
      if (k==1) return space2->getNumKeys(ele);
    }
    return 0;
  }
  virtual void getKeys(MElement* ele, std::vector<Dof> &keys, int k=0) const
  {
    if ((n1==0)&&(n2==0)) return;
    if ((n1!=0)&&(n2==0)) space1->getKeys(ele,keys,k);
    if ((n1==0)&&(n2!=0)) space2->getKeys(ele,keys,k);
    if ((n1==1)&&(n2==1))
    {
      if (k==0) return space1->getKeys(ele,keys);
      if (k==1) return space2->getKeys(ele,keys);
    }
    return;
  }
  virtual int getIncidentSpaceTag(int k=0) const
  { 
    if ((n1==0)&&(n2==0)) return 0;
    if ((n1!=0)&&(n2==0)) return space1->getIncidentSpaceTag(k);
    if ((n1==0)&&(n2!=0)) return space2->getIncidentSpaceTag(k);
    if ((n1==1)&&(n2==1))
    {
      if (k==0) return space1->getIncidentSpaceTag();
      if (k==1) return space2->getIncidentSpaceTag();
    }
    return 0;
  }
  virtual void get(MElement* ele, int npts, IntPt* GP, std::vector<ContainerValType> &vvals) const;
  virtual void get(MElement* ele, int npts, IntPt* GP, ContainerValType &vals) const {genTerm<genTensor0<double>,n1+n2>::get(ele,npts,GP,vals);}
  virtual ElasticTerm<n1,n2>* clone () const {return new ElasticTerm(space1,space2 );}
  
  double get_H(int i, int j) const{ return H(i,j); }
  void print_H() const{H.print();}
};


#include "genMiscTerms.hpp"

#endif// _GENMISCTERMS_H_