#ifndef FILE_FESPACE
#define FILE_FESPACE

/*********************************************************************/
/* File:   fespace.hh                                                */
/* Author: Joachim Schoeberl                                         */
/* Date:   25. Mar. 2000                                             */
/*********************************************************************/

/*
   Finite Element Space
*/


/// transformation from local to global orientation
enum TRANSFORM_TYPE { TRANSFORM_MAT_LEFT = 1,
		      TRANSFORM_MAT_RIGHT = 2,
		      TRANSFORM_MAT_LEFT_RIGHT = 3,
		      TRANSFORM_RHS = 4,
		      TRANSFORM_SOL = 8 };
		     

/**
   Base class for finite element space.
   Provides finite elements, global degrees of freedom, 
   and transformations of element-matrices and element-vectors
 */
class FESpace : public NGS_Object
{
protected:
  /// order of finite elements
  int order;
  /// how many components
  int dimension;
  /// complex space
  bool iscomplex;

  ///
  bool eliminate_internal;

  /// prolongation from coarser level to finer level
  ngmg::Prolongation *prol;

  /// stores matrix graph. attention: either sym or non-sym ... obsolete !!!
  ARRAY<MatrixGraph*> graphs;

  /// on which subdomains is the space defined ?
  ARRAY<int> definedon;
  /// on which boundaries is the space defined ?
  ARRAY<int> definedonbound;
  ///
  ARRAY<int> BEMboundary;

  BitArray dirichlet_boundaries;
  BitArray dirichlet_dofs;

  ///
  FiniteElement * tet;
  ///
  FiniteElement * prism;
  ///
  FiniteElement * pyramid;
  ///
  FiniteElement * hex;
  ///
  FiniteElement * trig;
  ///
  FiniteElement * quad;
  ///
  FiniteElement * segm;

  ///
  BilinearFormIntegrator * evaluator;
  ///
  BilinearFormIntegrator * boundary_evaluator;


  /// if non-zero, pointer to low order space
  FESpace * low_order_space;

  /// if directsolverclustered[i] is true, then the unknowns of domain i are clustered
  ARRAY<bool> directsolverclustered;
  
public:
  ///
  FESpace (const MeshAccess & ama, int aorder,
	   int adim, bool acomplex);
  FESpace (const MeshAccess & ama, const Flags & flags);
  ///
  virtual ~FESpace ();
  
  ///
  virtual void Update();

  ///
  virtual void Update(LocalHeap & lh);

  ///
  virtual void PrintReport (ostream & ost);

  ///
  int GetOrder () const { return order; }

  ///
  int GetDimension () const { return dimension; }

  ///
  bool IsComplex () const { return iscomplex; }

  ///
  void SetBEM (bool abem);

  /// will be replaced by getclassname !
  virtual const char * GetType() 
  { return GetClassName().c_str(); }


  ///
  virtual int GetNDof () const = 0;
  ///
  virtual int GetNDofLevel (int level) const;
  
  ///
  virtual const FiniteElement & GetFE (int elnr, LocalHeap & lh) const;
  ///
  virtual void GetDofNrs (int elnr, ARRAY<int> & dnums) const = 0;
  /// 
  virtual void GetExternalDofNrs (int elnr, ARRAY<int> & dnums) const;

  ///
  virtual const FiniteElement & GetSFE (int selnr, LocalHeap & lh) const;
  ///
  virtual void GetSDofNrs (int selnr, ARRAY<int> & dnums) const = 0;
  ///
  virtual void GetBEMDofNrs (ARRAY<int> & dnums) const;

  ///
  bool DefinedOn (int elnr) const
  { return !definedon.Size() || definedon[elnr]; }
  ///
  bool DefinedOnBoundary (int belnr) const
  {return !definedonbound.Size() || definedonbound[belnr]; }
  ///
  void SetDefinedOn (const BitArray & defon);
  ///
  void SetDefinedOnBoundary (const BitArray & defon);
  ///
  void SetDirichletBoundaries (const BitArray & dirbnds);
  ///
  const FiniteElement & GetFE (ELEMENT_TYPE type) const;

  ///
  FESpace & LowOrderFESpace () { return *low_order_space; }
  ///
  const FESpace & LowOrderFESpace () const { return *low_order_space; }



  //  class GridFunction * CreateGridFunction (const string & name) const;

  ///
  virtual void LockSomeDofs (BaseMatrix & mat) const { };
  ///
  virtual Table<int> * CreateSmoothingBlocks (int type = 0) const;
  /// for anisotropic plane smoothing
  virtual BitArray * CreateIntermediatePlanes (int type = 0) const
  { return 0; }
  virtual ARRAY<int> * CreateDirectSolverClusters (int type = 0) const
  { return 0; }



  void TransformMat (int elnr, bool boundary,
		     FlatMatrix<double> & mat, TRANSFORM_TYPE type) const
  {
    VTransformMR (elnr, boundary, mat, type);
  }
  
  void TransformMat (int elnr, bool boundary,
		     FlatMatrix<Complex> & mat, TRANSFORM_TYPE type) const
  {
    VTransformMC (elnr, boundary, mat, type);
  }
  
  void TransformVec (int elnr, bool boundary,
		     FlatVector<double> & vec, TRANSFORM_TYPE type) const
  {
    VTransformVR (elnr, boundary, vec, type);
  }
  
  void TransformVec (int elnr, bool boundary,
		     FlatVector<Complex> & vec, TRANSFORM_TYPE type) const
  {
    VTransformVC (elnr, boundary, vec, type);
  }
  

  virtual void VTransformMR (int elnr, bool boundary,
			     FlatMatrix<double> & mat, TRANSFORM_TYPE type) const
    { ; }
  virtual void VTransformMC (int elnr, bool boundary,
			     FlatMatrix<Complex> & mat, TRANSFORM_TYPE type) const
    { ; }
  virtual void VTransformVR (int elnr, bool boundary,
			     FlatVector<double> & vec, TRANSFORM_TYPE type) const
    { ; }
  virtual void VTransformVC (int elnr, bool boundary,
			     FlatVector<Complex> & vec, TRANSFORM_TYPE type) const
    { ; }
  
  
  ///
  virtual const ngmg::Prolongation * GetProlongation () const
  { return prol; }
  ///
  void SetProlongation (ngmg::Prolongation * aprol)
  { prol = aprol; }

  
  const BilinearFormIntegrator * GetEvaluator () const
  { return evaluator; }
  const BilinearFormIntegrator * GetBoundaryEvaluator () const
  { return boundary_evaluator; }


  ///
  MatrixGraph * GetGraph (int level, bool symmetric);

  ARRAY<SpecialElement*> specialelements;
};







/*


template <class T>
class T_FESpace : public FESpace
{
public:
  T_FESpace (const MeshAccess & ama, int aorder,
	     int adim, bool acomplex) 
    : FESpace(ama, aorder, adim, acomplex) { ; }
  
  const T & Spec() const { return static_cast<const T&> (*this);  }
  T & Spec() { return static_cast<T&> (*this); }

  template <class MAT>
  void TransformMat (int elnr, bool boundary,
		     MAT & mat, TRANSFORM_TYPE type) const
  {
    ;
  }
  template <class VEC>
  void TransformVec (int elnr, bool boundary,
		     VEC & vec, TRANSFORM_TYPE type) const
  {
    ;
  }
    

  virtual void VTransformMR (int elnr, bool boundary,
			     FlatMatrix<double> & mat, TRANSFORM_TYPE type) const 
  {
    Spec().TransformMat (elnr, boundary, mat, type);
  }

  virtual void VTransformMC (int elnr, bool boundary,
			     FlatMatrix<Complex> & mat, TRANSFORM_TYPE type) const
  {
    Spec().TransformMat (elnr, boundary, mat, type);
  }

  virtual void VTransformVR (int elnr, bool boundary,
			     FlatVector<double> & vec, TRANSFORM_TYPE type) const 
  {
    Spec().TransformVec (elnr, boundary, vec, type);
  }

  virtual void VTransformVC (int elnr, bool boundary,
			     FlatVector<Complex> & vec, TRANSFORM_TYPE type) const 
  {
    Spec().TransformVec (elnr, boundary, vec, type);
  }
};
*/





///
class NodalFESpace : public FESpace
{
  ///
  ARRAY<int> ndlevel;

public:

  NodalFESpace (const MeshAccess & ama, const Flags & flags);
  //		int aorder, int adim, bool acomplex,
  //		bool hierarchical = 1);
  ///
  virtual ~NodalFESpace ();

  virtual string GetClassName () const
  {
    return "NodalFESpace";
  }

  ///
  virtual void Update();
  ///
  virtual int GetNDof () const;
  ///
  virtual int GetNDofLevel (int level) const;
  ///
  virtual void GetDofNrs (int elnr, ARRAY<int> & dnums) const;
  ///
  virtual void GetSDofNrs (int selnr, ARRAY<int> & dnums) const;
};






///
class NonconformingFESpace : public FESpace
{
  ///
  ARRAY<int> ndlevel;

public:
  NonconformingFESpace (const MeshAccess & ama, const Flags & flags);
  virtual ~NonconformingFESpace ();

  virtual string GetClassName () const
  { return "Nonconforming FESpace"; }

  static FESpace * Create (const MeshAccess & ma, const Flags & flags)
  { return new NonconformingFESpace (ma, flags); }

  ///
  virtual void Update();
  ///
  virtual int GetNDof () const;
  ///
  virtual void GetDofNrs (int elnr, ARRAY<int> & dnums) const;
  ///
  virtual void GetSDofNrs (int selnr, ARRAY<int> & dnums) const;
};




#ifdef OLD
/// use edge and face tables
class NodalFESpaceAlt : public FESpace
{
  ///
  ARRAY<int> ndlevel;
  ///
  int nv, ned, nfa;
public:

  ///
  NodalFESpaceAlt (const MeshAccess & ama,
		   int aorder, int adim, bool acomplex);

  ///
  ~NodalFESpaceAlt ();


  virtual string GetClassName () const
  {
    return "NodalFESpaceAlt";
  }

  ///
  virtual void Update();
  ///
  virtual int GetNDof () const;
  ///
  virtual int GetNDofLevel (int level) const;

  /*
  ///
  virtual const NodalFiniteElement & GetFE (int elnr, LocalHeap & lh) const
  {
    return static_cast<const NodalFiniteElement&> (FESpace::GetFE(elnr));
  }

  ///
  virtual const FiniteElement & GetSFE (int selnr) const
  {
    return static_cast<const NodalFiniteElement&> (FESpace::GetSFE(selnr));
  }
  */
  ///
  virtual void GetDofNrs (int elnr, ARRAY<int> & dnums) const;
  ///
  virtual void GetSDofNrs (int selnr, ARRAY<int> & dnums) const;


  template <class MAT>
  void TransformMat (int elnr, bool boundary,
		     MAT & mat, TRANSFORM_TYPE tt) const;

  template <class VEC>
  void TransformVec (int elnr, bool boundary,
		     VEC & vec, TRANSFORM_TYPE tt) const;


  virtual void VTransformMR (int elnr, bool boundary,
			     FlatMatrix<double> & mat, TRANSFORM_TYPE tt) const 
  {
    TransformMat (elnr, boundary, mat, tt);
  }

  virtual void VTransformMC (int elnr, bool boundary,
			     FlatMatrix<Complex> & mat, TRANSFORM_TYPE tt) const
  {
    TransformMat (elnr, boundary, mat, tt);
  }

  virtual void VTransformVR (int elnr, bool boundary,
			     FlatVector<double> & vec, TRANSFORM_TYPE tt) const 
  {
    TransformVec (elnr, boundary, vec, tt);
  }

  virtual void VTransformVC (int elnr, bool boundary,
			     FlatVector<Complex> & vec, TRANSFORM_TYPE tt) const 
  {
    TransformVec (elnr, boundary, vec, tt);
  }

  /*
  ///
  virtual void TransformMatrix (int elnr, DenseMatrix & mat) const;
  ///
  virtual void TransformSurfMatrix (int elnr, DenseMatrix & mat) const;
  */
  ///
  virtual void LockSomeDofs (BaseMatrix & mat) const;
  ///
  virtual Table<int> * CreateSmoothingBlocks (int type = 0) const;
};
#endif



///
class ElementFESpace : public FESpace
{
  ///  ARRAY<int> startelement;
  ARRAY<int> ndlevel;
  int n_el_dofs;
public:

  ///
  ElementFESpace (const MeshAccess & ama,
		  int aorder, int adim, bool acomplex);

  ///
  ~ElementFESpace ();

  virtual string GetClassName () const
  {
    return "ElementFESpace";
  }

  ///
  virtual void Update();

  ///
  virtual int GetNDof () const
    { return ndlevel.Last(); }
  
  ///
  //  virtual const FiniteElement & GetFE (int elnr) const;
  ///
  virtual void GetDofNrs (int elnr, ARRAY<int> & dnums) const;

  ///
  virtual int GetNDofLevel (int level) const;

  ///
  // virtual const FiniteElement & GetSFE (int selnr) const;
  ///
  virtual void GetSDofNrs (int selnr, ARRAY<int> & dnums) const;
};





/// Non-continous fe space on boundary
class SurfaceElementFESpace : public FESpace
{
  ///
  ARRAY<int> ndlevel;
  int n_el_dofs;
public:

  ///
  SurfaceElementFESpace (const MeshAccess & ama,
			 int aorder, int adim, bool acomplex);

  ///
  ~SurfaceElementFESpace ();

  ///
  virtual string GetClassName() const
    { return "SurfaceElement"; }

  ///
  virtual void Update();

  ///
  virtual int GetNDof () const
    { return ndlevel.Last(); }

  ///
  virtual const FiniteElement & GetFE (int elnr, LocalHeap & lh) const;
  ///
  virtual void GetDofNrs (int elnr, ARRAY<int> & dnums) const;

  ///
  virtual int GetNDofLevel (int level) const;

  ///
  virtual void GetSDofNrs (int selnr, ARRAY<int> & dnums) const;
};








///
class NonConformingFESpace : public FESpace
{
  ///
  FE_NcSegm1 segm1;
  ///
  FE_NcTrig1 trig1;
  ///
  FE_NcTet1 tet1;

  ///
  HashTable<ngstd::INT<2>,int> *node2face2d;
  ///
  HashTable<ngstd::INT<3>,int> *node2face3d;
  ///
  ARRAY<ngstd::INT<2> > faces;
  ///
  ARRAY<int[4]> elementfaces;
  ///
  ARRAY<int> surfelementfaces;

  ///
  ARRAY<int[5]> parentfaces;

  ///
  ARRAY<short int> finelevelofedge;
  ///
  ARRAY<int> nflevel;
  
public:
  ///
  NonConformingFESpace (const MeshAccess & ama,
			int aorder, int adim, bool acomplex);

  ///
  ~NonConformingFESpace ();

  ///
  virtual string GetClassName() const
    { return "Non-conforming"; }

  ///
  virtual void Update();

  ///
  virtual int GetNDof () const;
  ///
  virtual int GetNDofLevel (int level) const;

  ///
  virtual const FiniteElement & GetFE (int elnr, LocalHeap & lh) const;
  ///
  virtual void GetDofNrs (int elnr, ARRAY<int> & dnums) const;

  ///
  virtual const FiniteElement & GetSFE (int selnr, LocalHeap & lh) const;
  ///
  virtual void GetSDofNrs (int selnr, ARRAY<int> & dnums) const;

  ///
  int GetFacePoint1 (int fnr) const { return faces[fnr][0]; }
  ///
  int GetFacePoint2 (int fnr) const { return faces[fnr][1]; }

  ///
  int GetParentFace1 (int fnr) const { return parentfaces[fnr][0]; }
  ///
  int GetParentFace2 (int fnr) const { return parentfaces[fnr][1]; }
  ///
  int GetParentFace3 (int fnr) const { return parentfaces[fnr][2]; }
  ///
  int GetParentFace4 (int fnr) const { return parentfaces[fnr][3]; }
  ///
  int GetParentFace5 (int fnr) const { return parentfaces[fnr][4]; }
  ///
  int GetFineLevelOfFace (int ednr) const { return finelevelofedge[ednr]; }
};












/// A combination of fe-spaces
class CompoundFESpace : public FESpace
{
protected:
  /// pointer to components
  ARRAY<const FESpace*> spaces;
  /// cummlated #dofs of components
  ARRAY<int> cummulative_nd;
  /// 
  ARRAY<int> ndlevel;
public:
  CompoundFESpace (const MeshAccess & ama,
		   const ARRAY<const FESpace*> & aspaces);
  ///
  virtual ~CompoundFESpace ();

  virtual string GetClassName () const
  {
    return "CompoundFESpace";
  }

  ///
  virtual void Update();

  ///
  virtual int GetNDof () const
  { return ndlevel.Last(); }
  ///
  virtual int GetNDofLevel (int level) const
  { return ndlevel[level]; }

  // returns start and end points of dofs corresponding to space "spacenr"
  // first space: spacenr = 0
  int GetStorageStart(int spacenr) const
  { return cummulative_nd[spacenr]; }

  int GetStorageEnd(int spacenr) const
  { return cummulative_nd[spacenr+1]; }



  const FESpace * operator[] (int i) const
  { return spaces[i]; }

  virtual const FiniteElement & GetFE (int elnr, LocalHeap & lh) const;
  ///
  virtual void GetDofNrs (int elnr, ARRAY<int> & dnums) const;
  ///
  virtual const FiniteElement & GetSFE (int selnr, LocalHeap & lh) const;
  ///
  virtual void GetSDofNrs (int selnr, ARRAY<int> & dnums) const;


  template <class MAT>
  void TransformMat (int elnr, bool boundary,
		     MAT & mat, TRANSFORM_TYPE tt) const;

  template <class VEC>
  void TransformVec (int elnr, bool boundary,
		     VEC & vec, TRANSFORM_TYPE tt) const;


  virtual void VTransformMR (int elnr, bool boundary,
			     FlatMatrix<double> & mat, TRANSFORM_TYPE tt) const 
  {
    TransformMat (elnr, boundary, mat, tt);
  }

  virtual void VTransformMC (int elnr, bool boundary,
			     FlatMatrix<Complex> & mat, TRANSFORM_TYPE tt) const
  {
    TransformMat (elnr, boundary, mat, tt);
  }

  virtual void VTransformVR (int elnr, bool boundary,
			     FlatVector<double> & vec, TRANSFORM_TYPE tt) const 
  {
    TransformVec (elnr, boundary, vec, tt);
  }

  virtual void VTransformVC (int elnr, bool boundary,
			     FlatVector<Complex> & vec, TRANSFORM_TYPE tt) const 
  {
    TransformVec (elnr, boundary, vec, tt);
  }

};









/// Registered FESpace classes
class FESpaceClasses
{
public:
  struct FESpaceInfo
  {
    string name;
    FESpace* (*creator)(const MeshAccess & ma, const Flags & flags);
    FESpaceInfo (const string & aname,
		 FESpace* (*acreator)(const MeshAccess & ma, const Flags & flags));
  };

  ARRAY<FESpaceInfo*> fesa;
public:
  FESpaceClasses();
  ~FESpaceClasses();  
  void AddFESpace (const string & aname, 
		   FESpace* (*acreator)(const MeshAccess & ma, const Flags & flags));
  
  const ARRAY<FESpaceInfo*> & GetFESpaces() { return fesa; }
  const FESpaceInfo * GetFESpace(const string & name);

  void Print (ostream & ost) const;
};
 
extern FESpaceClasses & GetFESpaceClasses ();






#endif