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


#ifndef _GENFUNCTION_SPACE_H_
#define _GENFUNCTION_SPACE_H_


#include "genDofIdManager.h"
#include "genTerm.h"
#include "genFilters.h"
#include "genSimpleFunction.h"

#include "SPoint3.h"
#include "Numeric.h"
#include "MElement.h"
#include "nodalBasis.h"
//double inv3x3(double [3][3], double [3][3]);


template<class T> class genFSpace : public diffTerm<T,1>
{
public:
  typedef typename TensorialTraits<T>::Scalar Scalar;
  typedef typename TensorialTraits<T>::ValType ValType;
  typedef typename TensorialTraits<T>::GradType GradType;
  typedef typename TensorialTraits<T>::HessType HessType;
  typedef typename ContainerTraits<ValType,1>::Container ContainerValType;
  typedef typename ContainerTraits<GradType,1>::Container ContainerGradType;
  typedef typename ContainerTraits<HessType,1>::Container ContainerHessType;
  typedef typename std::shared_ptr<genFSpace<T> > Handle;
  typedef typename std::shared_ptr<const genFSpace<T> > ConstHandle;
  static const int Nsp=1;

  virtual ~genFSpace() {}
  virtual int getNumKeys(MElement* ele, int k=0) const = 0; // if one needs the number of dofs
  virtual void getKeys(MElement* ele, std::vector<Dof> &keys, int k=0) const = 0;
  virtual int getIncidentSpaceTag(int k=0) const =0;
  virtual void setDegree(int i) {};
  virtual void f(MElement* ele, double u, double v, double w, ContainerValType &vals) const
  {
    IntPt Pt;
    Pt.pt[0]=u; Pt.pt[1]=v; Pt.pt[2]=w;
    Pt.weight=1.;

    std::vector<ContainerValType> vvals(1);

    get(ele,1,&Pt,vvals);
    vals=vvals[0];
  }
  virtual void gradf(MElement* ele, double u, double v, double w, ContainerGradType &grads) const
  {
    IntPt Pt;
    Pt.pt[0]=u; Pt.pt[1]=v; Pt.pt[2]=w;
    Pt.weight=1.;

    std::vector< ContainerGradType > vgrads(1);

    getgradf(ele,1,&Pt,vgrads);
    grads=vgrads[0];
  }
  virtual void gradfuvw(MElement* ele, double u, double v, double w, ContainerGradType &grads) const
  {
    IntPt Pt;
    Pt.pt[0]=u; Pt.pt[1]=v; Pt.pt[2]=w;
    Pt.weight=1.;

    std::vector< ContainerGradType > vgrads(1);

    getgradfuvw(ele,1,&Pt,vgrads);
    grads=vgrads[0];
  }
  virtual void hessf(MElement* ele, double u, double v, double w, ContainerHessType &hess) const
  {
    IntPt Pt;
    Pt.pt[0]=u; Pt.pt[1]=v; Pt.pt[2]=w;
    Pt.weight=1.;

    std::vector< ContainerHessType > vhesss(1);

    gethessfuvw(ele,1,&Pt,vhesss);
    hess=vhesss[0];
  }
  virtual void hessfuvw(MElement* ele, double u, double v, double w, ContainerHessType &hess) const
  {
    IntPt Pt;
    Pt.pt[0]=u; Pt.pt[1]=v; Pt.pt[2]=w;
    Pt.weight=1.;

    std::vector< ContainerHessType > vhesss(1);

    gethessfuvw(ele,1,&Pt,vhesss);
    hess=vhesss[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 getgradf(MElement* ele, int npts, IntPt* GP, std::vector< ContainerGradType > &vgrads) const {}
  virtual void getgradf(MElement* ele, int npts, IntPt* GP, ContainerGradType &grads) const {}
  virtual void getgradfuvw(MElement* ele, int npts, IntPt* GP, std::vector< ContainerGradType > &vgrads) const {}
  virtual void gethessf(MElement* ele, int npts, IntPt* GP, std::vector< ContainerHessType > &vhesss) const {}
  virtual void gethessf(MElement* ele, int npts, IntPt* GP, ContainerHessType &hesss) const {}
  virtual void gethessfuvw(MElement* ele, int npts, IntPt* GP, std::vector< ContainerHessType > &vhesss) const {}
  virtual genFSpace<T>* clone () const =0;
};


template<class T> class ScalarLagrangeFSpace : public genFSpace<T>
{
public:
  typedef typename TensorialTraits<T>::Scalar Scalar;
  typedef typename TensorialTraits<T>::ScalarType ScalarType;
  typedef typename TensorialTraits<T>::ValType ValType;
  typedef typename TensorialTraits<T>::GradType GradType;
  typedef typename TensorialTraits<T>::HessType HessType;
  typedef typename ContainerTraits<ValType,1>::Container ContainerValType;
  typedef typename ContainerTraits<GradType,1>::Container ContainerGradType;
  typedef typename ContainerTraits<HessType,1>::Container ContainerHessType;
  static const int Nsp=1;
protected:
  int iField; // field number (used to build dof keys)
  int degree;
private:
  virtual void getKeys(MVertex* ver, std::vector<Dof> &keys) const
  {
    keys.push_back(Dof(ver->getNum(), iField));
  }

public:
  ScalarLagrangeFSpace() : iField(DofIdManager::getInstance()->createField()), degree(-1) {} // degree=-1 -> geo element degree
  ScalarLagrangeFSpace(int id) : iField(id), degree(-1) {} // degree=-1 -> geo element degree
  ScalarLagrangeFSpace(const ScalarLagrangeFSpace<T> &other) : iField(other.iField), degree(other.degree) {}
  virtual ~ScalarLagrangeFSpace() {}
  virtual void setDegree(int i) {degree=i;}
  virtual int getId(void) const
  {
    return iField;
  }
  virtual int getNumKeys(MElement* ele, int k=0) const
  {
    const nodalBasis* fs = ele->getFunctionSpace(degree);
    return fs->getNumShapeFunctions();
  }
  virtual void getKeys(MElement* ele, std::vector<Dof> &keys, int k=0) const // appends ...
  {
    const nodalBasis* fs = ele->getFunctionSpace(degree);
    int ndofs = fs->getNumShapeFunctions();
//    int ndofs = ele->getNumShapeFunctions();

    keys.reserve(keys.size() + ndofs);
    for (int i=0; i < ndofs; ++i)
      getKeys(ele->getShapeFunctionNode(i), keys);
  }
  virtual int getIncidentSpaceTag(int k=0) const {return iField;}

  virtual void f(MElement* ele, double u, double v, double w, ContainerValType &vals) const
  {
    genFSpace<T>::f(ele, u, v, w, vals);
  }
  // Fonction renvoyant un vecteur contenant le grandient de chaque FF
  virtual void gradf(MElement* ele, double u, double v, double w, ContainerGradType &grads) const
  {
    genFSpace<T>::gradf(ele, u, v, w, grads);
  }
  virtual void gradfuvw(MElement* ele, double u, double v, double w, ContainerGradType &grads) const
  {
    genFSpace<T>::gradfuvw(ele, u, v, w, grads);
  }
  // Fonction renvoyant un vecteur contenant le hessien [][] de chaque FF dans l'espace ISOPARAMETRIQUE
  virtual void hessfuvw(MElement* ele, double u, double v, double w, ContainerHessType &hesss) const
  {
    genFSpace<T>::hessfuvw(ele, u, v, w, hesss);
  }

  // Fonction renvoyant un vecteur contenant, pour chaque point de Gauss, les vecteurs représentant la valeur de chaque FF
  virtual void get(MElement* ele, int npts, IntPt* GP, std::vector<ContainerValType> &vvals) const
  {
    const nodalBasis* fs = ele->getFunctionSpace(degree);
    int ndofs = fs->getNumShapeFunctions();
//    int ndofs = ele->getNumShapeFunctions();

    double u, v, w;
    int curpos = 0;
    double valsuvw[ndofs];
    for (int i=0; i<npts; ++i)
    {
      u = GP[i].pt[0]; v = GP[i].pt[1]; w = GP[i].pt[2];

      curpos = vvals[i].size();
      vvals[i].resize(curpos + ndofs, false);
      ele->getShapeFunctions(u, v, w, valsuvw, degree);
      for (int k=0; k<ndofs; ++k)
        vvals[i](curpos+k) = ValType(valsuvw[k]);
//      ele->getShapeFunctions(u, v, w, &(vvals[i](curpos)), degree);
    }
//use        void MElement::getShapeFunctions(double u, double v, double w, double s[], int o)
  }
  virtual void getgradf(MElement* ele, int npts, IntPt* GP, std::vector< ContainerGradType > &vgrads) const
  {
    const nodalBasis* fs = ele->getFunctionSpace(degree);
    int ndofs = fs->getNumShapeFunctions();

    double u, v, w;
    int curpos = 0;
    double gradsuvw[ndofs][3];
    double detJ;
    double jac[3][3];
    double invjac[3][3];
    for (int i=0; i<npts; ++i)
    {
      u = GP[i].pt[0]; v = GP[i].pt[1]; w = GP[i].pt[2];

      curpos = vgrads[i].size();
      vgrads[i].resize(curpos + ndofs,false);

      detJ = ele->getJacobian(u, v, w, jac); // redondant : on fait cet appel a l'exterieur
      inv3x3(jac, invjac);

      ele->getGradShapeFunctions(u, v, w, gradsuvw, degree);
      for (int j=0; j<ndofs; ++j)
        vgrads[i](curpos+j)=GradType(
        invjac[0][0] * gradsuvw[j][0] + invjac[0][1] * gradsuvw[j][1] + invjac[0][2] * gradsuvw[j][2],
        invjac[1][0] * gradsuvw[j][0] + invjac[1][1] * gradsuvw[j][1] + invjac[1][2] * gradsuvw[j][2],
        invjac[2][0] * gradsuvw[j][0] + invjac[2][1] * gradsuvw[j][1] + invjac[2][2] * gradsuvw[j][2]);
    }
  }
  virtual void getgradfuvw(MElement* ele, int npts, IntPt* GP, std::vector< ContainerGradType > &vgrads) const
  {
    const nodalBasis* fs = ele->getFunctionSpace(degree);
    int ndofs = fs->getNumShapeFunctions();

    double u, v, w;
    int curpos = 0;
    double gradsuvw[ndofs][3];
    for (int i=0; i<npts; ++i)
    {
      u = GP[i].pt[0]; v = GP[i].pt[1]; w = GP[i].pt[2];

//       if (elep!=ele)
//          ele->movePointFromParentSpaceToElementSpace(u, v, w);

      curpos = vgrads[i].size();
      vgrads[i].resize(curpos + ndofs, false);
      ele->getGradShapeFunctions(u, v, w, gradsuvw, degree);
      for (int j=0; j<ndofs; ++j)
        vgrads[i](curpos+j)=GradType(gradsuvw[j][0], gradsuvw[j][1], gradsuvw[j][2]);
    }
  }
  virtual void gethessfuvw(MElement* ele, int npts, IntPt* GP, std::vector< ContainerHessType > &vhesss) const
  {
    const nodalBasis* fs = ele->getFunctionSpace(degree);
    int ndofs = fs->getNumShapeFunctions();

    double u, v, w;
    int curpos = 0;
    double hessuvw[ndofs][3][3];
    for (int i=0; i<npts; ++i)
    {
      u = GP[i].pt[0]; v = GP[i].pt[1]; w = GP[i].pt[2];

//       if (elep!=ele)
//          ele->movePointFromParentSpaceToElementSpace(u, v, w);

      curpos = vhesss[i].size();
      vhesss[i].resize(curpos + ndofs, false);
      ele->getHessShapeFunctions(u, v, w, hessuvw, degree);
      HessType hesst;
      for (int j=0; j<ndofs; ++j)
      {
        hesst(0,0) = hessuvw[j][0][0]; hesst(0,1) = hessuvw[j][0][1]; hesst(0,2) = hessuvw[j][0][2];
        hesst(1,0) = hessuvw[j][1][0]; hesst(1,1) = hessuvw[j][1][1]; hesst(1,2) = hessuvw[j][1][2];
        hesst(2,0) = hessuvw[j][2][0]; hesst(2,1) = hessuvw[j][2][1]; hesst(2,2) = hessuvw[j][2][2];
        vhesss[i](curpos+j) = hesst;
      }
    }
  }

  virtual ScalarLagrangeFSpace<T>* clone () const {return new ScalarLagrangeFSpace<T>(*this);}
};


template<class T> class ScalarToAnyFSpace : public genFSpace<T>
{
public :
  typedef typename TensorialTraits<T>::Scalar Scalar;
  typedef typename TensorialTraits<T>::ScalarType ScalarType;
  typedef typename TensorialTraits<T>::ValType ValType;
  typedef typename TensorialTraits<T>::GradType GradType;
  typedef typename TensorialTraits<T>::HessType HessType;
  typedef typename ContainerTraits<ValType,1>::Container ContainerValType;
  typedef typename ContainerTraits<GradType,1>::Container ContainerGradType;
  typedef typename ContainerTraits<HessType,1>::Container ContainerHessType;
  static const int Nsp=1;
protected :
  typename genFSpace<ScalarType>::ConstHandle ScalarFS;
  std::vector<T> multipliers;
  std::vector<int> comp;

public :

  ScalarToAnyFSpace(const ScalarToAnyFSpace<T> &other) : ScalarFS(other.ScalarFS), multipliers(other.multipliers), comp(other.comp) {}

  ScalarToAnyFSpace(const typename genFSpace<ScalarType>::ConstHandle &SFS,
                    const T &mult1, int comp1_) : ScalarFS(SFS)
  {
    multipliers.push_back(mult1); comp.push_back(comp1_);
  }
  ScalarToAnyFSpace(const typename genFSpace<ScalarType>::ConstHandle &SFS,
                    const T &mult1, int comp1_, const T &mult2, int comp2_) : ScalarFS(SFS)
  {
    multipliers.push_back(mult1); comp.push_back(comp1_);
    multipliers.push_back(mult2); comp.push_back(comp2_);
  }
  ScalarToAnyFSpace(const typename genFSpace<ScalarType>::ConstHandle &SFS,
                    const T &mult1, int comp1_, const T &mult2, int comp2_, const T &mult3, int comp3_) : ScalarFS(SFS)
  {
    multipliers.push_back(mult1); comp.push_back(comp1_);
    multipliers.push_back(mult2); comp.push_back(comp2_);
    multipliers.push_back(mult3); comp.push_back(comp3_);
  }
//   virtual ~ScalarToAnyFSpace() {delete ScalarFS;}

  virtual int getNumKeys(MElement* ele, int k=0) const
  {
    return ScalarFS->getNumKeys(ele) * comp.size();
  }
  virtual void getKeys(MElement* ele, std::vector<Dof> &keys, int k=0) const
  {
    std::vector<Dof> bufk;
    ScalarFS->getKeys(ele, bufk);
    int nbdofs = bufk.size();
    int nbcomp = comp.size();
    int curpos = keys.size();
    keys.reserve(curpos + nbcomp * nbdofs);
    // Order given to keys : (numEntity,Component)
    //  comp. x |  comp. y |  comp. z
    // 0x 1x 2x | 0y 1y 2y | 0z 1z 2z
    for (int j=0; j<nbcomp; ++j)
    {
      for (int k=0; k<nbdofs; ++k)
      {
        int iField = bufk[k].getType();
        int vNum = bufk[k].getEntity();
        int newType = DofIdManager::getInstance()->createTypeWithTwoInts(iField, comp[j]);
        keys.push_back(Dof(vNum,newType));
      }
    }
  }
  virtual int getIncidentSpaceTag(int k=0) const {return ScalarFS->getIncidentSpaceTag(k);}

  virtual void f(MElement* ele, double u, double v, double w, ContainerValType &vals) const
  {
    genFSpace<T>::f(ele, u, v, w, vals);
  }
  virtual void gradf(MElement* ele, double u, double v, double w, ContainerGradType &grads) const
  {
    genFSpace<T>::gradf(ele, u, v, w, grads);
  }
  virtual void gradfuvw(MElement* ele, double u, double v, double w, ContainerGradType &grads) const
  {
    genFSpace<T>::gradfuvw(ele, u, v, w, grads);
  }
  virtual void hessfuvw(MElement* ele, double u, double v, double w, ContainerHessType &hesss) const
  {
    genFSpace<T>::hessfuvw(ele, u, v, w, hesss);
  }

  // Fonction renvoyant un vecteur contenant, pour chaque point de Gauss, les vecteurs représentant la valeur de chaque FF
  virtual void get(MElement* ele, int npts, IntPt* GP, std::vector<ContainerValType> &vvals) const
  {
    std::vector< typename ContainerTraits<ScalarType,1>::Container > vvalsd(npts);
    ScalarFS->get(ele, npts, GP, vvalsd);

    int nbcomp = comp.size();
    int nbdofs = ScalarFS->getNumKeys(ele);
    int curpos;
    for (int i=0; i<npts; ++i)
    {
      curpos = vvals[i].size();
      vvals[i].resize(curpos + nbcomp * nbdofs, false);
      int l=0;
      for (int j=0; j<nbcomp; ++j)
      {
        for (int k=0; k<nbdofs; ++k)
        {
          vvals[i](curpos+l) = multipliers[j] * vvalsd[i](k)();
          ++l;
        }
      }
    }
  }
  virtual void getgradf(MElement* ele, int npts, IntPt* GP, std::vector< ContainerGradType > &vgrads) const
  {
    std::vector< typename ContainerTraits<typename TensorialTraits<ScalarType>::GradType,1>::Container > vgradsd(npts);
    ScalarFS->getgradf(ele, npts, GP, vgradsd);

    int nbcomp = comp.size();
    int nbdofs = ScalarFS->getNumKeys(ele);
    int curpos;
    GradType grad;
    for (int i=0; i<npts; ++i)
    {
      curpos = vgrads[i].size();
      vgrads[i].resize(curpos + nbcomp * nbdofs, false);
      int l=0;
      for (int j=0; j<nbcomp; ++j)
      {
        for (int k=0; k<nbdofs; ++k)
        {
          tensprod(multipliers[j], vgradsd[i](k), grad);
          vgrads[i](curpos+l) = grad;
          ++l;
        }
      }
    }
  }
  virtual void getgradfuvw(MElement* ele, int npts, IntPt* GP, std::vector< ContainerGradType > &vgrads) const
  {
    std::vector< typename ContainerTraits<typename TensorialTraits<ScalarType>::GradType,1>::Container > vgradsd(npts);
    ScalarFS->getgradfuvw(ele, npts, GP, vgradsd);

    int nbcomp = comp.size();
    int nbdofs = ScalarFS->getNumKeys(ele);
    int curpos;
    GradType grad;
    for (int i=0; i<npts; ++i)
    {
      curpos = vgrads[i].size();
      vgrads[i].resize(curpos + nbcomp * nbdofs, false);
      int l=0;
      for (int j=0; j<nbcomp; ++j)
      {
        for (int k=0; k<nbdofs; ++k)
        {
          tensprod(multipliers[j], vgradsd[i](k), grad);
          vgrads[i](curpos+l) = grad;
          ++l;
        }
      }
    }
  }
  virtual void gethessfuvw(MElement* ele, int npts, IntPt* GP, std::vector< ContainerHessType > &vhesss) const {}

  virtual ScalarToAnyFSpace<T>* clone () const {return new ScalarToAnyFSpace(*this);}
};


template<class scalar, int N=3> class VectorLagrangeFSpace : public ScalarToAnyFSpace<genTensor1<scalar,N> >
{
public:
  typedef typename TensorialTraits<genTensor1<scalar,N> >::Scalar Scalar;
  typedef typename TensorialTraits<genTensor1<scalar,N> >::ScalarType ScalarType;
  typedef typename TensorialTraits<genTensor1<scalar,N> >::ValType ValType;
  typedef typename TensorialTraits<genTensor1<scalar,N> >::GradType GradType;
  typedef typename TensorialTraits<genTensor1<scalar,N> >::HessType HessType;
  typedef typename ContainerTraits<ValType,1>::Container ContainerValType;
  typedef typename ContainerTraits<GradType,1>::Container ContainerGradType;
  typedef typename ContainerTraits<HessType,1>::Container ContainerHessType;
  static const int Nsp=1;
  enum Along { VECTOR_X = 0, VECTOR_Y = 1, VECTOR_Z = 2 };
protected:
  static const ValType BasisVectors[3];

public:
  VectorLagrangeFSpace(const VectorLagrangeFSpace<scalar,N> &other) : ScalarToAnyFSpace<ValType>::ScalarToAnyFSpace(other) {}
  VectorLagrangeFSpace() :
    ScalarToAnyFSpace<ValType>::ScalarToAnyFSpace(
    typename genFSpace<ScalarType>::ConstHandle (new ScalarLagrangeFSpace<ScalarType>()),
    ValType(1.,0.,0.), VECTOR_X, ValType(0.,1.,0.), VECTOR_Y, ValType(0.,0.,1.), VECTOR_Z)
  {}
  VectorLagrangeFSpace(Along comp1) :
    ScalarToAnyFSpace<ValType>::ScalarToAnyFSpace(
    typename genFSpace<ScalarType>::ConstHandle (new ScalarLagrangeFSpace<ScalarType>()),
    BasisVectors[comp1], comp1)
  {}
  VectorLagrangeFSpace(Along comp1, Along comp2) :
    ScalarToAnyFSpace<ValType>::ScalarToAnyFSpace(
    typename genFSpace<ScalarType>::ConstHandle (new ScalarLagrangeFSpace<ScalarType>()),
    BasisVectors[comp1], comp1, BasisVectors[comp2], comp2)
  {}
  VectorLagrangeFSpace(Along comp1, Along comp2, Along comp3) :
    ScalarToAnyFSpace<ValType>::ScalarToAnyFSpace(
    typename genFSpace<ScalarType>::ConstHandle (new ScalarLagrangeFSpace<ScalarType>()),
    BasisVectors[comp1], comp1, BasisVectors[comp2], comp2, BasisVectors[comp3], comp3)
  {}
  VectorLagrangeFSpace(int id) :
    ScalarToAnyFSpace<ValType>::ScalarToAnyFSpace(
    typename genFSpace<ScalarType>::ConstHandle (new ScalarLagrangeFSpace<ScalarType>(id)),
    ValType(1.,0.,0.), VECTOR_X, ValType(0.,1.,0.), VECTOR_Y, ValType(0.,0.,1.), VECTOR_Z)
  {}
  VectorLagrangeFSpace(int id, Along comp1) :
    ScalarToAnyFSpace<ValType>::ScalarToAnyFSpace(
    typename genFSpace<ScalarType>::ConstHandle (new ScalarLagrangeFSpace<ScalarType>(id)),
    BasisVectors[comp1], comp1)
  {}
  VectorLagrangeFSpace(int id, Along comp1, Along comp2) :
    ScalarToAnyFSpace<ValType>::ScalarToAnyFSpace(
    typename genFSpace<ScalarType>::ConstHandle (new ScalarLagrangeFSpace<ScalarType>(id)),
    BasisVectors[comp1], comp1, BasisVectors[comp2], comp2)
  {}
  VectorLagrangeFSpace(int id, Along comp1, Along comp2, Along comp3) :
    ScalarToAnyFSpace<ValType>::ScalarToAnyFSpace(
    typename genFSpace<ScalarType>::ConstHandle (new ScalarLagrangeFSpace<ScalarType>(id)),
    BasisVectors[comp1], comp1, BasisVectors[comp2], comp2, BasisVectors[comp3], comp3)
  {}

  virtual ~VectorLagrangeFSpace() {}
  virtual VectorLagrangeFSpace<scalar,N>* clone () const {return new VectorLagrangeFSpace<scalar,N>(*this);}
};

template<class scalar, int N> const typename VectorLagrangeFSpace<scalar,N>::ValType VectorLagrangeFSpace<scalar,N>::BasisVectors[3] =
  {VectorLagrangeFSpace<scalar,N>::ValType(1, 0, 0), VectorLagrangeFSpace<scalar,N>::ValType(0, 1, 0), VectorLagrangeFSpace<scalar,N>::ValType(0, 0, 1)};


// new genTensor
static const double BasisVectorsinit[3][3]={{1,0,0},{0,1,0},{0,0,1}};
template<class scalar, int N=3> class VectorLagrangeFSpaceNG : public ScalarToAnyFSpace<genTensor<1,scalar,N> >
{
public:
  typedef typename TensorialTraits<genTensor<1,scalar,N> >::Scalar Scalar;
  typedef typename TensorialTraits<genTensor<1,scalar,N> >::ScalarType ScalarType;
  typedef typename TensorialTraits<genTensor<1,scalar,N> >::ValType ValType;
  typedef typename TensorialTraits<genTensor<1,scalar,N> >::GradType GradType;
  typedef typename TensorialTraits<genTensor<1,scalar,N> >::HessType HessType;
  typedef typename ContainerTraits<ValType,1>::Container ContainerValType;
  typedef typename ContainerTraits<GradType,1>::Container ContainerGradType;
  typedef typename ContainerTraits<HessType,1>::Container ContainerHessType;
  static const int Nsp=1;
  enum Along { VECTOR_X = 0, VECTOR_Y = 1, VECTOR_Z = 2 };
protected:
  static const ValType BasisVectors[3];
  
public:
  VectorLagrangeFSpaceNG(const VectorLagrangeFSpace<scalar,N> &other) : ScalarToAnyFSpace<ValType>::ScalarToAnyFSpace(other) {}
  VectorLagrangeFSpaceNG() :
    ScalarToAnyFSpace<ValType>::ScalarToAnyFSpace(
    typename genFSpace<ScalarType>::ConstHandle (new ScalarLagrangeFSpace<ScalarType>()),
    ValType(BasisVectorsinit[0]), VECTOR_X, ValType(BasisVectorsinit[1]), VECTOR_Y, ValType(BasisVectorsinit[2]), VECTOR_Z)
//    ValType(1.,0.,0.), VECTOR_X, ValType(0.,1.,0.), VECTOR_Y, ValType(0.,0.,1.), VECTOR_Z)    
  {}
  VectorLagrangeFSpaceNG(Along comp1) :
    ScalarToAnyFSpace<ValType>::ScalarToAnyFSpace(
    typename genFSpace<ScalarType>::ConstHandle (new ScalarLagrangeFSpace<ScalarType>()),
    BasisVectors[comp1], comp1)
  {}
  VectorLagrangeFSpaceNG(Along comp1, Along comp2) :
    ScalarToAnyFSpace<ValType>::ScalarToAnyFSpace(
    typename genFSpace<ScalarType>::ConstHandle (new ScalarLagrangeFSpace<ScalarType>()),
    BasisVectors[comp1], comp1, BasisVectors[comp2], comp2)
  {}
  VectorLagrangeFSpaceNG(Along comp1, Along comp2, Along comp3) :
    ScalarToAnyFSpace<ValType>::ScalarToAnyFSpace(
    typename genFSpace<ScalarType>::ConstHandle (new ScalarLagrangeFSpace<ScalarType>()),
    BasisVectors[comp1], comp1, BasisVectors[comp2], comp2, BasisVectors[comp3], comp3)
  {}
  VectorLagrangeFSpaceNG(int id) :
    ScalarToAnyFSpace<ValType>::ScalarToAnyFSpace(
    typename genFSpace<ScalarType>::ConstHandle (new ScalarLagrangeFSpace<ScalarType>(id)),
    ValType(BasisVectorsinit[0]), VECTOR_X, ValType(BasisVectorsinit[1]), VECTOR_Y, ValType(BasisVectorsinit[2]), VECTOR_Z)
//    ValType(1.,0.,0.), VECTOR_X, ValType(0.,1.,0.), VECTOR_Y, ValType(0.,0.,1.), VECTOR_Z)
  {}
  VectorLagrangeFSpaceNG(int id, Along comp1) :
    ScalarToAnyFSpace<ValType>::ScalarToAnyFSpace(
    typename genFSpace<ScalarType>::ConstHandle (new ScalarLagrangeFSpace<ScalarType>(id)),
    BasisVectors[comp1], comp1)
  {}
  VectorLagrangeFSpaceNG(int id, Along comp1, Along comp2) :
    ScalarToAnyFSpace<ValType>::ScalarToAnyFSpace(
    typename genFSpace<ScalarType>::ConstHandle (new ScalarLagrangeFSpace<ScalarType>(id)),
    BasisVectors[comp1], comp1, BasisVectors[comp2], comp2)
  {}
  VectorLagrangeFSpaceNG(int id, Along comp1, Along comp2, Along comp3) :
    ScalarToAnyFSpace<ValType>::ScalarToAnyFSpace(
    typename genFSpace<ScalarType>::ConstHandle (new ScalarLagrangeFSpace<ScalarType>(id)),
    BasisVectors[comp1], comp1, BasisVectors[comp2], comp2, BasisVectors[comp3], comp3)
  {}

  virtual ~VectorLagrangeFSpaceNG() {}
  virtual VectorLagrangeFSpaceNG<scalar,N>* clone () const {return new VectorLagrangeFSpaceNG<scalar,N>(*this);}
};

template<class scalar, int N> const typename VectorLagrangeFSpaceNG<scalar,N>::ValType VectorLagrangeFSpaceNG<scalar,N>::BasisVectors[3] =
  {VectorLagrangeFSpaceNG<scalar,N>::ValType(BasisVectorsinit[0]), VectorLagrangeFSpaceNG<scalar,N>::ValType(BasisVectorsinit[1]), VectorLagrangeFSpaceNG<scalar,N>::ValType(BasisVectorsinit[2])};
//  {VectorLagrangeFSpaceNG<scalar,N>::ValType(BasisVectorsinit[0]), VectorLagrangeFSpaceNG<scalar,N>::ValType(BasisVectorsinit[1]), VectorLagrangeFSpaceNG<scalar,N>::ValType(BasisVectorsinit[2])};


  
  
  
template<class T> class CompositeFSpace : public genFSpace<T>
{
public:
  typedef typename TensorialTraits<T>::ValType ValType;
  typedef typename TensorialTraits<T>::GradType GradType;
  typedef typename TensorialTraits<T>::HessType HessType;
  typedef typename ContainerTraits<ValType,1>::Container ContainerValType;
  typedef typename ContainerTraits<GradType,1>::Container ContainerGradType;
  typedef typename ContainerTraits<HessType,1>::Container ContainerHessType;
  typedef typename std::vector<typename genFSpace<T>::ConstHandle>::const_iterator iterFS;
  static const int Nsp=1;
protected:
  std::vector< typename genFSpace<T>::ConstHandle > FSs;
  int iField;

public:
  CompositeFSpace(const CompositeFSpace<T> &other) : FSs(other.FSs),iField(other.iField) {}
  CompositeFSpace() : iField(DofIdManager::getInstance()->createField()) {}
  CompositeFSpace(int id) : iField(id) {}
  CompositeFSpace(const typename genFSpace<T>::ConstHandle &FS) : iField(DofIdManager::getInstance()->createField()) { FSs.push_back(FS);}
  CompositeFSpace(const typename genFSpace<T>::ConstHandle &FS1, const typename genFSpace<T>::ConstHandle &FS2) : iField(DofIdManager::getInstance()->createField())
  {
    FSs.push_back(FS1);
    FSs.push_back(FS2);
  }
  CompositeFSpace(const typename genFSpace<T>::ConstHandle &FS1, const typename genFSpace<T>::ConstHandle &FS2, const typename genFSpace<T>::ConstHandle &FS3) : iField(DofIdManager::getInstance()->createField())
  {
    FSs.push_back(FS1);
    FSs.push_back(FS2);
    FSs.push_back(FS3);
  }
  CompositeFSpace(const typename genFSpace<T>::ConstHandle &FS1, const typename genFSpace<T>::ConstHandle &FS2, const typename genFSpace<T>::ConstHandle &FS3, const typename genFSpace<T>::ConstHandle &FS4) : iField(DofIdManager::getInstance()->createField())
  {
    FSs.push_back(FS1);
    FSs.push_back(FS2);
    FSs.push_back(FS3);
    FSs.push_back(FS4);
  }

  void insert(const typename genFSpace<T>::ConstHandle &FS)
  {
    FSs.push_back(FS);
  }

  virtual int getId(void) const
  {
    return iField;
  }
  virtual int getNumKeys(MElement* ele, int k=0) const
  {
    int ndofs = 0;
    for (iterFS it = FSs.begin(); it != FSs.end(); ++it)
      ndofs += (*it)->getNumKeys(ele);
    return ndofs;
  }
  virtual void getKeys(MElement* ele, std::vector<Dof> &keys, int k=0) const
  {
    for (iterFS it = FSs.begin(); it != FSs.end(); ++it)
      (*it)->getKeys(ele, keys);
  }
  virtual int getIncidentSpaceTag(int k=0) const {return iField;}

  virtual void f(MElement* ele, double u, double v, double w, ContainerValType &vals) const
  {
    genFSpace<T>::f(ele, u, v, w, vals);
  }
  virtual void gradf(MElement* ele, double u, double v, double w, ContainerGradType &grads) const
  {
    genFSpace<T>::gradf(ele, u, v, w, grads);
  }
  virtual void gradfuvw(MElement* ele, double u, double v, double w, ContainerGradType &grads) const
  {
    genFSpace<T>::gradfuvw(ele, u, v, w, grads);
  }
  virtual void hessfuvw(MElement* ele, double u, double v, double w, ContainerHessType &hess) const
  {
    genFSpace<T>::hessfuvw(ele, u, v, w, hess);
  }

  // Fonction renvoyant un vecteur contenant, pour chaque point de Gauss, les vecteurs représentant la valeur de chaque FF
  virtual void get(MElement* ele, int npts, IntPt* GP, std::vector<ContainerValType> &vvals) const
  {
    for (iterFS it = FSs.begin(); it != FSs.end(); ++it)
      (*it)->get(ele, npts, GP, vvals);
  }
  virtual void getgradf(MElement* ele, int npts, IntPt* GP, std::vector< ContainerGradType > &vgrads) const
  {
    for (iterFS it = FSs.begin(); it != FSs.end(); ++it)
      (*it)->getgradf(ele, npts, GP, vgrads);
  }
  virtual void getgradfuvw(MElement* ele, int npts, IntPt* GP, std::vector< ContainerGradType > &vgrads) const
  {
    for (iterFS it = FSs.begin(); it != FSs.end(); ++it)
      (*it)->getgradfuvw(ele, npts, GP, vgrads);
  }
  virtual void gethessfuvw(MElement* ele, int npts, IntPt* GP, std::vector< ContainerHessType > &vhesss) const
  {
    for (iterFS it = FSs.begin(); it != FSs.end(); ++it)
      (*it)->gethessfuvw(ele, npts, GP, vhesss);
  }

  virtual CompositeFSpace<T>* clone () const {return new CompositeFSpace<T>(*this);}
};


template<class T> class EnrichedFSpace : public genFSpace<T>
{
public:
  typedef typename TensorialTraits<T>::ScalarType ScalarType;
  typedef typename TensorialTraits<T>::ValType ValType;
  typedef typename TensorialTraits<T>::GradType GradType;
  typedef typename TensorialTraits<T>::HessType HessType;
  typedef typename ContainerTraits<ValType,1>::Container ContainerValType;
  typedef typename ContainerTraits<GradType,1>::Container ContainerGradType;
  typedef typename ContainerTraits<HessType,1>::Container ContainerHessType;
  static const int Nsp=1;
protected:
  typename genFSpace<T>::ConstHandle BaseFS;
  int iIndex; // Enrichment number (used to build dof keys)
  typename diffTerm<ScalarType,0>::ConstHandle enrichFunction;

public:
  EnrichedFSpace(const EnrichedFSpace<T> &other) : BaseFS(other.BaseFS), iIndex(other.iIndex), enrichFunction(other.enrichFunction) {}
  EnrichedFSpace(const typename genFSpace<T>::ConstHandle &BaseFS_, const typename diffTerm<ScalarType,0>::ConstHandle &enrichFunction_) :
    BaseFS(BaseFS_), iIndex(DofIdManager::getInstance()->createIndex(BaseFS->getIncidentSpaceTag())), enrichFunction(enrichFunction_){}
  virtual ~EnrichedFSpace() {}

  virtual int getNumKeys(MElement* ele, int k=0) const
  {
    int nbdofs = BaseFS->getNumKeys(ele);
    return nbdofs;
  }
  virtual void getKeys(MElement* ele, std::vector<Dof> &keys, int k=0) const
  {
    int normalk = BaseFS->getNumKeys(ele);

    std::vector<Dof> bufk;
    bufk.reserve(normalk);
    BaseFS->getKeys(ele, bufk);
    int normaldofs = bufk.size();
    int nbdofs = normaldofs;

    int curpos = keys.size();
    keys.reserve(curpos + nbdofs);

    // Order given to keys : (numEntity,Component,numIndex)
    //                numIndex 0               |                numIndex 1
    //   comp. x   |   comp. y   |   comp. z   |   comp. x   |   comp. y   |   comp. z
    // 0x0 1x0 2x0 | 0y0 1y0 2y0 | 0z0 1z0 2z0 | 0x1 1x1 2x1 | 0y1 1y1 2y1 | 0z1 1z1 2z1

    for (int k=0; k<nbdofs; ++k)
    {
      int vNum = bufk[k].getEntity();
      int iField, iComp;
      DofIdManager::getInstance()->getTwoIntsFromType(bufk[k].getType(), iField, iComp);
      int newType = DofIdManager::getInstance()->createTypeWithThreeInts(iField, iComp, iIndex);
      keys.push_back(Dof(vNum, newType));
    }
  }
  virtual int getIncidentSpaceTag(int k=0) const {return 0;}

  virtual void f(MElement* ele, double u, double v, double w, ContainerValType &vals) const
  {
    genFSpace<T>::f(ele, u, v, w, vals);
  }
  virtual void gradf(MElement* ele, double u, double v, double w, ContainerGradType &grads) const
  {
    genFSpace<T>::gradf(ele, u, v, w, grads);
  }
  virtual void gradfuvw(MElement* ele, double u, double v, double w, ContainerGradType &grads) const
  {
    genFSpace<T>::gradfuvw(ele, u, v, w, grads);
  }
  virtual void hessfuvw(MElement* ele, double u, double v, double w, ContainerHessType &hesss) const
  {
    genFSpace<T>::hessfuvw(ele, u, v, w, hesss);
  }

  virtual void get(MElement* ele, int npts, IntPt* GP, std::vector<ContainerValType> &vvals) const
  {
    // Get the spacebase valsd
    std::vector<ContainerValType> vvalsd(npts);
    BaseFS->get(ele, npts, GP, vvalsd);

    // Identification of enriched dofs
    int nbdofs = BaseFS->getNumKeys(ele);
    if (!nbdofs) return;

    // Enrichment function computation
    std::vector< genScalar<ScalarType> > vvals_enrich(npts);
    enrichFunction->get(ele, npts, GP, vvals_enrich);

    int curpos;
    for (int i = 0; i<npts; ++i)
    {
      curpos = vvals[i].size();           // if in composite function space
      vvals[i].resize(curpos + nbdofs, false);
      for (int j=0; j<nbdofs; ++j)
        vvals[i](curpos+j) = vvalsd[i](j) * vvals_enrich[i]();
    }
  }
  virtual void getgradf(MElement* ele, int npts, IntPt* GP, std::vector< ContainerGradType > &vgrads) const
  {
    // Get the spacebase gradsd
    std::vector< ContainerGradType > vgradsd(npts);
    BaseFS->getgradf(ele, npts, GP, vgradsd);

    // We need spacebase valsd to compute total gradient
    std::vector<ContainerValType> vvalsd(npts);
    BaseFS->get(ele, npts, GP, vvalsd);

    // Identification of enriched dofs
    int nbdofs = BaseFS->getNumKeys(ele);
    if (!nbdofs) return;

    // Enrichment function computation
    std::vector< genScalar<ScalarType> > vvals_enrich(npts);
    enrichFunction->get(ele, npts, GP, vvals_enrich);
    std::vector< genScalar< typename diffTerm<ScalarType,0>::GradType> > vgrads_enrich(npts);
    enrichFunction->getgradf(ele, npts, GP, vgrads_enrich);

    GradType GradFunc;
    int curpos;
    for (int i = 0; i<npts; ++i)
    {
      curpos = vgrads[i].size();           // if in composite function space
      vgrads[i].resize(curpos + nbdofs, false);
      for (int j=0; j<nbdofs; ++j)
      {
        tensprod(vvalsd[i](j), vgrads_enrich[i](), GradFunc);
        vgrads[i](curpos+j) = (vgradsd[i](j) * vvals_enrich[i]()) + GradFunc;
      }
    }
  }
  virtual void getgradfuvw(MElement* ele, int npts, IntPt* GP, std::vector< ContainerGradType > &vgrads) const {}
  virtual void gethessfuvw(MElement* ele, int npts, IntPt* GP, std::vector< ContainerHessType > &vhesss) const {}

  virtual EnrichedFSpace<T>* clone () const {return new EnrichedFSpace<T>(*this);}
};


/*template<class T> class xFemFSpace : public genFSpace<T>
{
public:
  typedef typename TensorialTraits<T>::ValType ValType;
  typedef typename TensorialTraits<T>::GradType GradType;
  typedef typename TensorialTraits<T>::HessType HessType;
  typedef typename ContainerTraits<ValType,1>::Container ContainerValType;
  typedef typename ContainerTraits<GradType,1>::Container ContainerGradType;
  typedef typename ContainerTraits<HessType,1>::Container ContainerHessType;
  static const int Nsp=1;
protected:
  typename genFSpace<T>::ConstHandle BaseFS;
  int iIndex; // Enrichment number (used to build dof keys)
  simpleFunctionOnElement<double>* enrichFunction;

public:
  xFemFSpace(const xFemFSpace<T>&other) : BaseFS(other.BaseFS), iIndex(other.iIndex), enrichFunction(other.enrichFunction) {}
  xFemFSpace(const typename genFSpace<T>::ConstHandle &BaseFS_, simpleFunctionOnElement<double>* enrichFunction_) :
    BaseFS(BaseFS_), iIndex(DofIdManager::getInstance()->createIndex(BaseFS->getIncidentSpaceTag())), enrichFunction(enrichFunction_){}
  virtual ~xFemFSpace() {}

  virtual int getNumKeys(MElement* ele) const
  {
    MElement* elep = ele->getParent();
    if (!elep) elep = ele;
    int nbdofs = BaseFS->getNumKeys(elep);
    return nbdofs;
  }
  virtual void getKeys(MElement* ele, std::vector<Dof> &keys) const
  {
    MElement* elep = ele->getParent();
    if (!elep) elep = ele;

    int normalk = BaseFS->getNumKeys(elep);

    std::vector<Dof> bufk;
    bufk.reserve(normalk);
    BaseFS->getKeys(elep, bufk);
    int normaldofs = bufk.size();
    int nbdofs = normaldofs;

    int curpos = keys.size();
    keys.reserve(curpos + nbdofs);

    // get keys so the order is ex:(a2x,a2y,a3x,a3y)
    // enriched dof tagged with ( i2 -> i2 + 1 )
    for (int k=0; k<nbdofs; ++k)
    {
      int iField, iComp;
      DofIdManager::getInstance()->getTwoIntsFromType(bufk[k].getType(), iField, iComp);
      int vNum = bufk[k].getEntity();
      int newType = DofIdManager::getInstance()->createTypeWithThreeInts(iField, iComp, iIndex);
      keys.push_back(Dof(vNum, newType));
    }
  }
  virtual int getIncidentSpaceTag(int k=0) const {return 0;}

  virtual void f(MElement* ele, double u, double v, double w, ContainerValType &vals) const
  {
    genFSpace<T>::f(ele, u, v, w, vals);
  }
  virtual void gradf(MElement* ele, double u, double v, double w, ContainerGradType &grads) const
  {
    genFSpace<T>::gradf(ele, u, v, w, grads);
  }
  virtual void gradfuvw(MElement* ele, double u, double v, double w, ContainerGradType &grads) const
  {
    genFSpace<T>::gradfuvw(ele, u, v, w, grads);
  }
  virtual void hessfuvw(MElement* ele, double u, double v, double w, ContainerHessType &hesss) const
  {
    genFSpace<T>::hessfuvw(ele, u, v, w, hesss);
  }

  virtual void get(MElement* ele, int npts, IntPt* GP, std::vector<ContainerValType> &vvals) const
  {
    // We need parent parameters
    MElement* elep = ele->getParent();
    if (!elep) elep = ele;

    // Get the spacebase valsd
    std::vector<ContainerValType> vvalsd(npts);
    BaseFS->get(elep, npts, GP, vvalsd);

    int nbdofs = 0;
    int curpos = 0;
    SPoint3 p;
    double u, v, w;
    double func;

    for (int i=0; i<npts; ++i)
    {
      nbdofs = vvalsd[i].size();
      curpos = vvals[i].size(); // if in composite function space
      vvals[i].resize(curpos + nbdofs, false);
      // then enriched dofs so the order is ex:(a2x,a2y,a3x,a3y)
      if (nbdofs > 0)
      { // if enriched
        // Enrichment function computation
        u = GP[i].pt[0]; v = GP[i].pt[1]; w = GP[i].pt[2];
        elep->pnt(u, v, w, p); // parametric to cartesian coordinates
        enrichFunction->setElement(elep);
        func = (*enrichFunction)(p.x(), p.y(), p.z());
        for (int j=0; j<nbdofs; ++j)
          vvals[i](curpos+j) = vvalsd[i](j) * func;
      }
    }
  }
  virtual void getgradf(MElement* ele, int npts, IntPt* GP, std::vector< ContainerGradType > &vgrads) const
  {
    // We need parent parameters
    MElement* elep = ele->getParent();
    if (!elep) elep = ele;

    // Get the spacebase gradsd
    std::vector< ContainerGradType > vgradsd(npts);
    BaseFS->getgradf(elep, npts, GP, vgradsd);

    // We need spacebase valsd to compute total gradient
    std::vector<ContainerValType> vvalsd(npts);
    BaseFS->get(elep, npts, GP, vvalsd);

    int nbdofs = 0;
    int curpos = 0;
    SPoint3 p;
    double u, v, w;
    double func;
    double df[3];
    GradType GradFunc;

    for (int i=0; i<npts; ++i)
    {
      nbdofs = vgradsd[i].size();
      curpos = vgrads[i].size(); // if in composite function space
      vgrads[i].resize(curpos + nbdofs, false);
      // then enriched dofs so the order is ex:(a2x,a2y,a3x,a3y)
      if (nbdofs > 0)
      { // if enriched
        u = GP[i].pt[0]; v = GP[i].pt[1]; w = GP[i].pt[2];
        elep->pnt(u, v, w, p);
        enrichFunction->setElement(elep);
        enrichFunction->gradient(p.x(), p.y(), p.z(), df[0], df[1], df[2]);
        ValType gradfuncenrich(df[0], df[1], df[2]);

        // Enrichment function computation
        enrichFunction->setElement(elep);
        func = (*enrichFunction)(p.x(), p.y(), p.z());

        for (int j=0; j<nbdofs; ++j)
        {
          tensprod(vvalsd[i](j), gradfuncenrich, GradFunc);
          vgrads[i](curpos+j) = vvalsd[i](j) * func + GradFunc;
        }
      }
    }
  }
  virtual void getgradfuvw(MElement* ele, int npts, IntPt* GP, std::vector< ContainerGradType > &vgrads) const {}
  virtual void gethessfuvw(MElement* ele, int npts, IntPt* GP, std::vector< ContainerHessType > &vhesss) const {}

  virtual xFemFSpace<T>* clone () const {return new xFemFSpace<T>(*this);}
};*/


template<class T> class FilteredFSpace : public genFSpace<T>
{
public :
  typedef typename TensorialTraits<T>::ValType ValType;
  typedef typename TensorialTraits<T>::GradType GradType;
  typedef typename TensorialTraits<T>::HessType HessType;
  typedef typename ContainerTraits<ValType,1>::Container ContainerValType;
  typedef typename ContainerTraits<GradType,1>::Container ContainerGradType;
  typedef typename ContainerTraits<HessType,1>::Container ContainerHessType;
  static const int Nsp=1;
protected:
  typename genFSpace<T>::ConstHandle BaseFS;
  genFilterDof* filterDof;
  genFilterElement* filterEle;

public:
  FilteredFSpace(const FilteredFSpace<T> &other) : BaseFS(other.BaseFS), filterDof(other.filterDof), filterEle(other.filterEle) {}
  FilteredFSpace<T>(const typename genFSpace<T>::ConstHandle &BaseFS_, genFilterDof* filterDof_) : BaseFS(BaseFS_), filterDof(filterDof_), filterEle(new genFilterElementTrivial) {}
  FilteredFSpace<T>(const typename genFSpace<T>::ConstHandle &BaseFS_, genFilterElement* filterEle_) : BaseFS(BaseFS_), filterDof(new genFilterDofTrivial), filterEle(filterEle_) {}
  FilteredFSpace<T>(const typename genFSpace<T>::ConstHandle &BaseFS_, genFilterDof* filterDof_, genFilterElement* filterEle_) : BaseFS(BaseFS_), filterDof(filterDof_), filterEle(filterEle_) {}
  virtual ~FilteredFSpace() {}
  virtual int getNumKeys(MElement* ele, int k=0) const
  {
    MElement* elep = ele->getParent();
    if (!elep) elep = ele;

    int nbfiltereddofs = 0;
    if ((*filterEle)(elep))
    {
      std::vector<Dof> bufk;
      BaseFS->getKeys(elep, bufk);
      int nbdofs = bufk.size();
      for (int i=0; i<nbdofs; ++i)
      {
        if ((*filterDof)(bufk[i]))
          ++nbfiltereddofs;
      }
    }
    return nbfiltereddofs;
  }
  virtual void getKeys(MElement* ele, std::vector<Dof> &keys, int k=0) const
  {
    MElement* elep = ele->getParent();
    if (!elep) elep = ele;

    if ((*filterEle)(elep))
    {
      std::vector<Dof> bufk;
      BaseFS->getKeys(elep, bufk);
      int nbdofs = bufk.size();
      for (int i=0; i<nbdofs; ++i)
      {
        if ((*filterDof)(bufk[i]))
          keys.push_back(bufk[i]);
      }
    }
  }
  virtual int getIncidentSpaceTag(int k=0) const {return BaseFS->getIncidentSpaceTag();}

  virtual void f(MElement* ele, double u, double v, double w, ContainerValType &vals) const
  {
    genFSpace<T>::f(ele, u, v, w, vals);
  }
  virtual void gradf(MElement* ele, double u, double v, double w, ContainerGradType &grads) const
  {
    genFSpace<T>::gradf(ele, u, v, w, grads);
  }
  virtual void gradfuvw(MElement* ele, double u, double v, double w, ContainerGradType &grads) const
  {
    genFSpace<T>::gradfuvw(ele, u, v, w, grads);
  }
  virtual void hessfuvw(MElement* ele, double u, double v, double w, ContainerHessType &hesss) const
  {
    genFSpace<T>::hessfuvw(ele, u, v, w, hesss);
  }

  virtual void get(MElement* ele, int npts, IntPt* GP, std::vector<ContainerValType> &vvals) const
  {
    // We need parent parameters
    MElement* elep = ele->getParent();
    if (!elep) elep = ele;

    if ((*filterEle)(elep))
    {
      // Get the BaseFS valsd
      std::vector<ContainerValType> vvalsd(npts);
      BaseFS->get(ele, npts, GP, vvalsd);

      // Get the BaseFS keys
      std::vector<Dof> bufk;
      BaseFS->getKeys(ele, bufk);
      int nbdofs = bufk.size();

      // Filter the BaseFS keys
      std::vector<int> selection;
      for (int k=0; k<nbdofs; ++k)
        if ((*filterDof)(bufk[k]))
          selection.push_back(k);
      int nbfiltereddofs = selection.size();

      int curpos = 0;
      for (int i=0; i<npts; ++i)
      {
        curpos = vvals[i].size(); // if in composite function space
        vvals[i].resize(curpos + nbfiltereddofs, false);
        for (int j=0; j<nbfiltereddofs; ++j)
          vvals[i](curpos+j) = vvalsd[i](selection[j]);
      }
    }
  }
  virtual void getgradf(MElement* ele, int npts, IntPt* GP, std::vector< ContainerGradType > &vgrads) const
  {
    // We need parent parameters
    MElement* elep = ele->getParent();
    if (!elep) elep = ele;

    if ((*filterEle)(elep))
    {
      // Get the spacebase gradsd
      std::vector< ContainerGradType > vgradsd(npts);
      BaseFS->getgradf(ele, npts, GP, vgradsd);

      // Get the BaseFS keys
      std::vector<Dof> bufk;
      BaseFS->getKeys(ele, bufk);
      int nbdofs = bufk.size();

      // Filter the BaseFS keys
      std::vector<int> selection;
      for (int k=0; k<nbdofs; ++k)
        if ((*filterDof)(bufk[k]))
          selection.push_back(k);
      int nbfiltereddofs = selection.size();

      int curpos = 0;
      for (int i=0; i<npts; ++i)
      {
        curpos = vgrads[i].size(); // if in composite function space
        vgrads[i].resize(curpos + nbfiltereddofs, false);
        for (int j=0; j<nbfiltereddofs; ++j)
          vgrads[i](curpos+j) = vgradsd[i](selection[j]);
      }
    }
  }
  virtual void getgradfuvw(MElement* ele, int npts, IntPt* GP, std::vector< ContainerGradType > &vgrads) const {}
  virtual void gethessfuvw(MElement* ele, int npts, IntPt* GP, std::vector< ContainerHessType > &vhesss) const {}

  virtual FilteredFSpace<T>* clone () const {return new FilteredFSpace<T>(*this);}
};




//-----------------------------------------------------------------------------------------------------------------------------------
// The ScalarToAnyFSpaceOriented class has a genTerm as template parameter, instead of a tensorial type.
// This enables to orient the local basis according to each element.
// Using a constantField as genTerm type us back to the case of ScalarToAnyFSpace.
//-----------------------------------------------------------------------------------------------------------------------------------


template<class GT> class ScalarToAnyFSpaceOriented : public genFSpace<typename GT::ValType>
{
public :
  typedef typename TensorialTraits<typename GT::ValType>::Scalar Scalar;
  typedef typename TensorialTraits<typename GT::ValType>::ScalarType ScalarType;
  typedef typename TensorialTraits<typename GT::ValType>::ValType ValType;
  typedef typename TensorialTraits<typename GT::ValType>::GradType GradType;
  typedef typename TensorialTraits<typename GT::ValType>::HessType HessType;
  typedef typename ContainerTraits<ValType,1>::Container ContainerValType;
  typedef typename ContainerTraits<GradType,1>::Container ContainerGradType;
  typedef typename ContainerTraits<HessType,1>::Container ContainerHessType;
  static const int Nsp=1;
protected :
  typename genFSpace<typename GT::ScalarType>::ConstHandle ScalarFS;
  std::vector<std::shared_ptr<GT> > multipliers;
  std::vector<int> comp;

public :
  ScalarToAnyFSpaceOriented(const ScalarToAnyFSpaceOriented<GT> &other) : ScalarFS(other.ScalarFS), multipliers(other.multipliers), comp(other.comp) {}

  ScalarToAnyFSpaceOriented(const typename genFSpace<ScalarType>::ConstHandle &SFS,
                    const std::shared_ptr<GT> &mult1, int comp1_) : ScalarFS(SFS)
  {
    multipliers.push_back(mult1); comp.push_back(comp1_);
  }
  ScalarToAnyFSpaceOriented(const typename genFSpace<ScalarType>::ConstHandle &SFS,
                    const std::shared_ptr<GT> &mult1, int comp1_, const std::shared_ptr<GT> &mult2, int comp2_) : ScalarFS(SFS)
  {
    multipliers.push_back(mult1); comp.push_back(comp1_);
    multipliers.push_back(mult2); comp.push_back(comp2_);
  }
  ScalarToAnyFSpaceOriented(const typename genFSpace<ScalarType>::ConstHandle &SFS,
                    const std::shared_ptr<GT> &mult1, int comp1_, const std::shared_ptr<GT> &mult2, int comp2_, const std::shared_ptr<GT> &mult3, int comp3_) : ScalarFS(SFS)
  {
    multipliers.push_back(mult1); comp.push_back(comp1_);
    multipliers.push_back(mult2); comp.push_back(comp2_);
    multipliers.push_back(mult3); comp.push_back(comp3_);

  }
//   virtual ~ScalarToAnyFSpaceOriented() {delete ScalarFS;}

  virtual int getNumKeys(MElement* ele, int k=0) const
  {
    return ScalarFS->getNumKeys(ele) * comp.size();
  }
  virtual void getKeys(MElement* ele, std::vector<Dof> &keys, int k=0) const
  {
    std::vector<Dof> bufk;
    ScalarFS->getKeys(ele, bufk);
    int nbdofs = bufk.size();
    int nbcomp = comp.size();
    int curpos = keys.size();
    keys.reserve(curpos + nbcomp * nbdofs);
    for (int j=0; j<nbcomp; ++j)
    {
      for (int k=0; k<nbdofs; ++k)
      {
        int vNum = bufk[k].getEntity();
        int iField = bufk[k].getType();
        int newType = DofIdManager::getInstance()->createTypeWithTwoInts(iField, comp[j]);
        keys.push_back(Dof(vNum,newType));
      }
    }
  }
  virtual int getIncidentSpaceTag(int k=0) const {return ScalarFS->getIncidentSpaceTag(k);}

  virtual void f(MElement* ele, double u, double v, double w, ContainerValType &vals) const
  {
    genFSpace<ValType>::f(ele, u, v, w, vals);
  }
  virtual void gradf(MElement* ele, double u, double v, double w, ContainerGradType &grads) const
  {
    genFSpace<ValType>::gradf(ele, u, v, w, grads);
  }
  virtual void gradfuvw(MElement* ele, double u, double v, double w, ContainerGradType &grads) const
  {
    genFSpace<ValType>::gradfuvw(ele, u, v, w, grads);
  }
  virtual void hessfuvw(MElement* ele, double u, double v, double w, ContainerHessType &hesss) const
  {
    genFSpace<ValType>::hessfuvw(ele, u, v, w, hesss);
  }

  // Fonction renvoyant un vecteur contenant, pour chaque point de Gauss, les vecteurs représentant la valeur de chaque FF
  virtual void get(MElement* ele, int npts, IntPt* GP, std::vector<ContainerValType> &vvals) const
  {
    std::vector< typename ContainerTraits<ScalarType,1>::Container > vvalsd(npts);
    ScalarFS->get(ele, npts, GP, vvalsd);

    int nbcomp = comp.size();
    std::vector< std::vector< genScalar<typename GT::ValType> > > v_vvalsm(nbcomp);
    for (int j=0; j<nbcomp; ++j)
    {
      v_vvalsm[j] = std::vector< genScalar<typename GT::ValType> >(npts);
      multipliers[j]->get(ele, npts, GP, v_vvalsm[j]);
    }

    int nbdofs = ScalarFS->getNumKeys(ele);
    int curpos;
    for (int i=0; i<npts; ++i)
    {
      curpos = vvals[i].size();
      vvals[i].resize(curpos + nbcomp * nbdofs, false);
      int l=0;
      for (int j=0; j<nbcomp; ++j)
      {
        for (int k=0; k<nbdofs; ++k)
        {
          vvals[i](curpos+l) = v_vvalsm[j][i]() * vvalsd[i](k);
          ++l;
        }
      }
    }
  }
  virtual void getgradf(MElement* ele, int npts, IntPt* GP, std::vector< ContainerGradType > &vgrads) const
  {
    std::vector< typename ContainerTraits<typename TensorialTraits<ScalarType>::GradType,1>::Container > vgradsd(npts);
    ScalarFS->getgradf(ele, npts, GP, vgradsd);

    int nbcomp = comp.size();
    std::vector< std::vector< genScalar<ValType> > > v_vvalsm(nbcomp);
    for (int j=0; j<nbcomp; ++j)
    {
      v_vvalsm[j] = std::vector< genScalar<ValType> >(npts);
      multipliers[j]->get(ele, npts, GP, v_vvalsm[j]);
    }

    int nbdofs = ScalarFS->getNumKeys(ele);
    int curpos;
    GradType grad;
    for (int i=0; i<npts; ++i)
    {
      curpos = vgrads[i].size();
      vgrads[i].resize(curpos + nbcomp * nbdofs, false);
      int l=0;
      for (int j=0; j<nbcomp; ++j)
      {
        for (int k=0; k<nbdofs; ++k)
        {
          tensprod(v_vvalsm[j][i](), vgradsd[i](k), grad);  // + v_vgradsm . vvalsd   if   GT=diffTerm
          vgrads[i](curpos+l) = grad;
          ++l;
        }
      }
    }
  }
  virtual void getgradfuvw(MElement* ele, int npts, IntPt* GP, std::vector< ContainerGradType > &vgrads) const
  {
    std::vector< typename ContainerTraits<typename TensorialTraits<ScalarType>::GradType,1>::Container > vgradsd(npts);
    ScalarFS->getgradfuvw(ele, npts, GP, vgradsd);

    int nbcomp = comp.size();
    std::vector< std::vector< genScalar<ValType> > > v_vvalsm(nbcomp);
    for (int j=0; j<nbcomp; ++j)
    {
      v_vvalsm[j] = std::vector< genScalar<ValType> >(npts);
      multipliers[j]->get(ele, npts, GP, v_vvalsm[j]);
    }

    int nbdofs = ScalarFS->getNumKeys(ele);
    int curpos;
    GradType grad;
    for (int i=0; i<npts; ++i)
    {
      curpos = vgrads[i].size();
      vgrads[i].resize(curpos + nbcomp * nbdofs, false);
      int l=0;
      for (int j=0; j<nbcomp; ++j)
      {
        for (int k=0; k<nbdofs; ++k)
        {
          tensprod(v_vvalsm[j][i](), vgradsd[i](k), grad);  // + v_vgradsm . vvalsd   if   GT=diffTerm
          vgrads[i](curpos+l) = grad;
          ++l;
        }
      }
    }
  }
  virtual void gethessfuvw(MElement* ele, int npts, IntPt* GP, std::vector< ContainerHessType > &vhesss) const {}

  virtual ScalarToAnyFSpaceOriented<GT>* clone () const {return new ScalarToAnyFSpaceOriented(*this);}
};

template<class scalar, int N=3> class VectorLagrangeFSpaceOriented : public ScalarToAnyFSpaceOriented<genTerm<genTensor1<scalar,N>,0> >
{
public:
  typedef typename TensorialTraits<genTensor1<scalar,N> >::ScalarType ScalarType;
  typedef typename TensorialTraits<genTensor1<scalar,N> >::ValType ValType;
  typedef typename TensorialTraits<genTensor1<scalar,N> >::GradType GradType;
  typedef typename TensorialTraits<genTensor1<scalar,N> >::HessType HessType;
  typedef typename ContainerTraits<ValType,1>::Container ContainerValType;
  typedef typename ContainerTraits<GradType,1>::Container ContainerGradType;
  typedef typename ContainerTraits<HessType,1>::Container ContainerHessType;
  static const int Nsp=1;
  enum Along { VECTOR_X = 0, VECTOR_Y = 1, VECTOR_Z = 2 };

public:
  VectorLagrangeFSpaceOriented(const VectorLagrangeFSpaceOriented<scalar,N> &other) : ScalarToAnyFSpaceOriented<genTerm<ValType,0> >::ScalarToAnyFSpaceOriented(other) {}
  VectorLagrangeFSpaceOriented() :
    ScalarToAnyFSpaceOriented<genTerm<ValType,0> >::ScalarToAnyFSpaceOriented(
      typename genFSpace<ScalarType>::ConstHandle (new ScalarLagrangeFSpace<ScalarType>()),
      typename genTerm<ValType,0>::Handle(new ConstantField<ValType>(ValType(1.,0.,0.))), VECTOR_X,
      typename genTerm<ValType,0>::Handle(new ConstantField<ValType>(ValType(0.,1.,0.))), VECTOR_Y,
      typename genTerm<ValType,0>::Handle(new ConstantField<ValType>(ValType(0.,0.,1.))), VECTOR_Z )
  {}
  VectorLagrangeFSpaceOriented(const std::shared_ptr<genTerm<ValType,0> > &mult1, Along comp1) :
          ScalarToAnyFSpaceOriented<genTerm<ValType,0> >::ScalarToAnyFSpaceOriented(
          typename genFSpace<ScalarType>::ConstHandle (new ScalarLagrangeFSpace<ScalarType>()), mult1, comp1)
  {}
  VectorLagrangeFSpaceOriented(const std::shared_ptr<genTerm<ValType,0> > &mult1, Along comp1,
                               const std::shared_ptr<genTerm<ValType,0> > &mult2, Along comp2) :
          ScalarToAnyFSpaceOriented<genTerm<ValType,0> >::ScalarToAnyFSpaceOriented(
          typename genFSpace<ScalarType>::ConstHandle (new ScalarLagrangeFSpace<ScalarType>()), mult1, comp1, mult2, comp2)
  {}
  VectorLagrangeFSpaceOriented(const std::shared_ptr<genTerm<ValType,0> > &mult1, Along comp1,
                               const std::shared_ptr<genTerm<ValType,0> > &mult2, Along comp2,
                               const std::shared_ptr<genTerm<ValType,0> > &mult3, Along comp3) :
          ScalarToAnyFSpaceOriented<genTerm<ValType,0> >::ScalarToAnyFSpaceOriented(
          typename genFSpace<ScalarType>::ConstHandle (new ScalarLagrangeFSpace<ScalarType>()), mult1, comp1, mult2, comp2, mult3, comp3)
  {}
  VectorLagrangeFSpaceOriented(int id) :
    ScalarToAnyFSpaceOriented<genTerm<ValType,0> >::ScalarToAnyFSpaceOriented(
      typename genFSpace<ScalarType>::ConstHandle (new ScalarLagrangeFSpace<ScalarType>(id)),
      typename genTerm<ValType,0>::Handle(new ConstantField<ValType>(ValType(1.,0.,0.))), VECTOR_X,
      typename genTerm<ValType,0>::Handle(new ConstantField<ValType>(ValType(0.,1.,0.))), VECTOR_Y,
      typename genTerm<ValType,0>::Handle(new ConstantField<ValType>(ValType(0.,0.,1.))), VECTOR_Z )
  {}
  VectorLagrangeFSpaceOriented(int id, const std::shared_ptr<genTerm<ValType,0> > &mult1, Along comp1) :
          ScalarToAnyFSpaceOriented<genTerm<ValType,0> >::ScalarToAnyFSpaceOriented(
          typename genFSpace<ScalarType>::ConstHandle (new ScalarLagrangeFSpace<ScalarType>(id)), mult1, comp1)
  {}
  VectorLagrangeFSpaceOriented(int id, const std::shared_ptr<genTerm<ValType,0> > &mult1, Along comp1,
                                       const std::shared_ptr<genTerm<ValType,0> > &mult2, Along comp2) :
          ScalarToAnyFSpaceOriented<genTerm<ValType,0> >::ScalarToAnyFSpaceOriented(
          typename genFSpace<ScalarType>::ConstHandle (new ScalarLagrangeFSpace<ScalarType>(id)), mult1, comp1, mult2, comp2)
  {}
  VectorLagrangeFSpaceOriented(int id, const std::shared_ptr<genTerm<ValType,0> > &mult1, Along comp1,
                                       const std::shared_ptr<genTerm<ValType,0> > &mult2, Along comp2,
                                       const std::shared_ptr<genTerm<ValType,0> > &mult3, Along comp3) :
          ScalarToAnyFSpaceOriented<genTerm<ValType,0> >::ScalarToAnyFSpaceOriented(
          typename genFSpace<ScalarType>::ConstHandle (new ScalarLagrangeFSpace<ScalarType>(id)), mult1, comp1, mult2, comp2, mult3, comp3)
  {}
  virtual ~VectorLagrangeFSpaceOriented() {}
  virtual VectorLagrangeFSpaceOriented<scalar,N>* clone () const {return new VectorLagrangeFSpaceOriented<scalar,N>(*this);}
};


#endif // _GENFUNCTION_SPACE_H_