// ============================================================================
// ==                                                                        ==
// == AMDiS - Adaptive multidimensional simulations                          ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
// ==  crystal growth group                                                  ==
// ==                                                                        ==
// ==  Stiftung caesar                                                       ==
// ==  Ludwig-Erhard-Allee 2                                                 ==
// ==  53175 Bonn                                                            ==
// ==  germany                                                               ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
// ==  http://www.caesar.de/cg/AMDiS                                         ==
// ==                                                                        ==
// ============================================================================

/** \file DOFVector.h */

#ifndef AMDIS_DOFVECTOR_H 
#define AMDIS_DOFVECTOR_H 
 
// ===========================================================================
// ===== includes ============================================================
// ===========================================================================

#include "FixVec.h"
#include "Global.h" 
#include "Flag.h" 
#include "RCNeighbourList.h" 
#include "DOFIterator.h"
#include "DOFIndexed.h"
#include "DOFContainer.h"
#include "MemoryManager.h"
#include "Boundary.h"
#include "CreatorInterface.h"
#include "Serializable.h"
#include "DOFMatrix.h" 
#include <vector> 
#include <memory> 
#include <list> 
 
namespace AMDiS {

  // ===========================================================================
  // ===== forward declarations ================================================
  // ===========================================================================

  class Mesh; 
  class FiniteElemSpace; 
  class ElInfo; 
  class DOFAdmin; 
  class BasisFunction; 
  class FillInfo; 
  class Quadrature; 
  class FastQuadrature;
  class DOFMatrix; 
  class MultiGridSortSolver; 
  class Operator;
  class ElementVector;
  class BoundaryManager;

  template<typename ReturnType, typename ArgumentType> class AbstractFunction;


  template<typename T> 
  class DOFVectorBase : public DOFIndexed<T>
  {
  public:
    DOFVectorBase() 
      : feSpace(NULL),
	elementVector(NULL),
	boundaryManager(NULL)
    {};

    DOFVectorBase(const FiniteElemSpace *f, ::std::string n)
      : feSpace(f),
	name(n),
	elementVector(NULL),
	boundaryManager(NULL)
    {};

    virtual const T *getLocalVector(const Element *el, T* localVec) const;

    const T *getVecAtQPs(const ElInfo         *elInfo, 
			 const Quadrature     *quad,
			 const FastQuadrature *quadFast,
			 T                    *vecAtQPs) const;

    const WorldVector<T> *getGrdAtQPs(const ElInfo         *elInfo, 
				      const Quadrature     *quad,
				      const FastQuadrature *quadFast,
				      WorldVector<T>       *grdAtQPs) const;

    const WorldMatrix<T> *getD2AtQPs(const ElInfo         *elInfo, 
				     const Quadrature     *quad,
				     const FastQuadrature *quadFast,
				     WorldMatrix<T>       *d2AtQPs) const;

    virtual const FiniteElemSpace* getFESpace() const {
      return feSpace;
    };

    ElementVector *assemble(T factor, ElInfo *elInfo,
			    const BoundaryType *bound, 
			    Operator *op = NULL);

    void addElementVector(T sign,
			  const ElementVector &elVec, 
			  const BoundaryType *bound,
			  bool add = true); 
 
    inline void addOperator(Operator* op, 
			    double *factor = NULL,
			    double *estFactor = NULL) 
    {
      operators.push_back(op);
      operatorFactor.push_back(factor);
      operatorEstFactor.push_back(estFactor);
    };

    inline ::std::vector<double*>::iterator getOperatorFactorBegin() {
      return operatorFactor.begin();
    };

    inline ::std::vector<double*>::iterator getOperatorFactorEnd() {
      return operatorFactor.end();
    };

    inline ::std::vector<double*>::iterator getOperatorEstFactorBegin() {
      return operatorEstFactor.begin();
    };

    inline ::std::vector<double*>::iterator getOperatorEstFactorEnd() {
      return operatorEstFactor.end();
    };


    inline ::std::vector<Operator*>::iterator getOperatorsBegin() {
      return operators.begin();
    };

    inline ::std::vector<Operator*>::iterator getOperatorsEnd() {
      return operators.end();
    };

    Flag getAssembleFlag();

    /** \brief
     * Evaluates \f[ u_h(x(\lambda)) = \sum_{i=0}^{m-1} vec[ind[i]] * 
     * \varphi^i(\lambda) \f] where \f$ \varphi^i \f$ is the i-th basis function,
     * \f$ x(\lambda) \f$ are the world coordinates of lambda and
     * \f$ m \f$ is the number of basis functions
     */
    T evalUh(const DimVec<double>& lambda, DegreeOfFreedom* ind);

    inline ::std::vector<Operator*>& getOperators() { return operators; };

    inline ::std::vector<double*>& getOperatorFactor() { return operatorFactor; };
    inline ::std::vector<double*>& getOperatorEstFactor() { 
      return operatorEstFactor; 
    };

    /** \brief
     * Returns \ref name
     */
    inline const ::std::string& getName() const { return name; }; 

    inline BoundaryManager* getBoundaryManager() const { return boundaryManager; };

    inline void setBoundaryManager(BoundaryManager *bm) {
      boundaryManager = bm;
    };

  protected:
    const FiniteElemSpace *feSpace;

    ::std::string name;

    ElementVector *elementVector;

    ::std::vector<Operator*> operators;

    ::std::vector<double*> operatorFactor;

    ::std::vector<double*> operatorEstFactor;

    BoundaryManager *boundaryManager;
  };


  // template<typename T>
  // class QPEvaluatable : public DOFVectorBase<T>
  // {
  // public:
  //   QPEvaluatable() 
  //     : DOFVectorBase<T>()
  //   {};

  //   QPEvaluatable(const FiniteElemSpace *f, ::std::string n)
  //     : DOFVectorBase<T>(f, n)
  //   {};

  // protected:
  // };

  // ===========================================================================
  // ===== defs ================================================================
  // ===========================================================================

  /** \brief
   * Specifies which operation should be done after coarsening
   */
  typedef enum{
    NO_OPERATION = 0,   
    COARSE_RESTRICT = 1,
    COARSE_INTERPOL = 2 
  } CoarsenOperation;
 
  // ===========================================================================
  // ===== class DOFVector =====================================================
  // ===========================================================================

  /** \ingroup DOFAdministration 
   * \brief
   * The DOFs described above are just integers that can be used as indices into 
   * vectors and matrices. During refinement and coarsening of the mesh, the 
   * number of used DOFs, the meaning of one integer index, and even the total 
   * range of DOFs change. To be able to handle these changes automatically for 
   * all vectors, which are indexed by the DOFs, special data structures are 
   * used which contain such vector data. Lists of these structures are kept in 
   * DOFAdmin, so that all vectors in the lists can be resized together with the
   * range of DOFs. During refinement and coarsening of elements, values can be
   * interpolated automatically to new DOFs, and restricted from old DOFs.
   */
  template<typename T> 
  class DOFVector : public DOFVectorBase<T>, public Serializable
  {  
  public:
    MEMORY_MANAGED(DOFVector<T>);

    /** \ingroup DOFAdministration
     * \brief
     * Enables the access of DOFVector<T>::Iterator. Alias for DOFIterator<T>
     */
    class Iterator : public DOFIterator<T> {
    public:
      Iterator(DOFIndexed<T> *c, DOFIteratorType type)
	: DOFIterator<T>(c, type)
      {};   

      Iterator(DOFAdmin *admin, DOFIndexed<T> *c, DOFIteratorType type)
	: DOFIterator<T>(admin, c, type)
      {};   
    };

    class Creator : public CreatorInterface<DOFVector<T> > {
    public:
      MEMORY_MANAGED(Creator);

      Creator(FiniteElemSpace *feSpace_) : feSpace(feSpace_) {};

      DOFVector<T> *create() { 
	return NEW DOFVector<T>(feSpace, ""); 
      };

      void free(DOFVector<T> *vec) { 
	DELETE vec; 
      };

    private:
      FiniteElemSpace *feSpace;
    };

  public:
    /** \brief
     * Empty constructor. No initialization!
     */
    DOFVector() 
      : DOFVectorBase<T>(), 
	//elementVector(NULL), 
	refineInter(false), 
	feSpace(NULL),
	coarsenOperation(NO_OPERATION)
    {};

    /** \brief
     * Constructs a DOFVector with name n belonging to FiniteElemSpace f
     */
    DOFVector(const FiniteElemSpace* f,::std::string n); 

    /** \brief
     * Initialization.
     */
    void init(const FiniteElemSpace* f, ::std::string n);

    /** \brief
     * Copy Constructor
     */
    DOFVector(const DOFVector& rhs) {
      *this=rhs;   
      name=rhs.name+"copy";
      if(feSpace && feSpace->getAdmin()) {
	(dynamic_cast<DOFAdmin*>(feSpace->getAdmin()))->addDOFIndexed(this);
      }
    }; 

    /** \brief
     * Destructor
     */
    virtual ~DOFVector();

    /** \brief
     * Returns iterator to the begin of \ref vec
     */
    typename ::std::vector<T>::iterator begin() { 
      return vec.begin(); 
    };

    /** \brief
     * Returns iterator to the end of \ref vec
     */
    typename ::std::vector<T>::iterator end() { 
      return vec.end(); 
    };
  
    /** \brief
     * Used by DOFAdmin to compress this DOFVector. Implementation of
     * DOFIndexedBase::compress()
     */
    virtual void compressDOFIndexed(int first, int last,
				    ::std::vector<DegreeOfFreedom> &newDof);

    /** \brief
     * Sets \ref refineInter to b
     */
    inline void refineInterpol(bool b) { 
      refineInter = b; 
    };

    /** \brief
     * Sets \ref coarsenOperation to op
     */
    inline void setCoarsenOperation(CoarsenOperation op) { 
      coarsenOperation = op; 
    };

    /** \brief
     * Returns \ref coarsenOperation
     */
    inline CoarsenOperation getCoarsenOperation() { return coarsenOperation; };

    /** \brief
     * Restriction after coarsening. Implemented for DOFVector<double>
     */
    inline void coarseRestrict(RCNeighbourList&, int) {};

    /** \brief
     * Returns \ref refineInter
     */
    inline bool refineInterpol() { return refineInter; }; 

    /** \brief
     * Interpolation after refinement.
     */
    inline void refineInterpol(RCNeighbourList&, int) {}; 

    /** \brief
     * Returns \ref vec
     */
    //  ::std::vector<T>& getVector() const {return (::std::vector<T>&) vec;} ; 
    ::std::vector<T>& getVector() { return  vec;} ; 

    //   /** \brief
    //    * Returns \ref feSpace
    //    */
    //   inline const FiniteElemSpace* getFESpace() const { return feSpace; }; 

    /** \brief
     * Returns size of \ref vec
     */
    inline int getSize() const { 
      return vec.size();
    }; 

    /** \brief
     * Returns used size of the vector.
     */
    inline int getUsedSize() const { 
      return feSpace->getAdmin()->getUsedSize(); 
    };

    /** \brief
     * Resizes \ref vec to n
     */
    inline void resize(int n) { 
      FUNCNAME("DOFVector<T>::resize");
      TEST_EXIT((n>=0))("Can't resize DOFVector to negative size\n"); 
      vec.resize(n);
    }; 

    /** \brief
     * Resizes \ref vec to n and inits new values with init
     */
    inline void resize(int n, T init) { 
      FUNCNAME("DOFVector<T>::resize");
      TEST_EXIT((n>=0))("Can't resize DOFVector to negative size\n"); 
      vec.resize(n, init);
    }; 

    /** \brief
     * Returns \ref vec[i]
     */
    inline const T& operator[](DegreeOfFreedom i) const {
      FUNCNAME("DOFVector<T>::operator[]");
      TEST_EXIT((i>=0)&&(i<static_cast<int>(vec.size())))("Illegal vector index %d.\n",i);
      return vec[i];
    }; 

    inline int getSize() { 
      return vec.size(); 
    };

    /** \brief
     * Returns \ref vec[i]
     */
    inline T& operator[](DegreeOfFreedom i) {
      FUNCNAME("DOFVector<T>::operator[]");
      TEST_EXIT((i>=0)&&(i<static_cast<int>(vec.size())))("Illegal vector index %d.\n",i);
      return vec[i];
    };

    // Andreas Ergaenzungen:
  
    /** \brief
     * Calculates Integral of this DOFVector
     */
    double Int(Quadrature* q = NULL) const;

    /** \brief
     * Calculates L1 norm of this DOFVector
     */
    double L1Norm(Quadrature* q = NULL) const;

    // Ende Andreas Ergaenzungen (hier) unten geht es weiter

  
    /** \brief
     * Calculates L2 norm of this DOFVector
     */
    inline double L2Norm(Quadrature* q = NULL) const {
      return sqrt(L2NormSquare());
    };

    /** \brief
     * Calculates square of L2 norm of this DOFVector
     */
    double L2NormSquare(Quadrature* q = NULL) const;

    /** \brief
     * Calculates H1 norm of this DOFVector
     */
    inline double H1Norm(Quadrature* q = NULL) const {
      return sqrt(H1NormSquare());
    };  

    /** \brief
     * Calculates square of H1 norm of this DOFVector
     */
    double H1NormSquare(Quadrature* q = NULL) const;  

    /** \brief
     * Calculates euclidian norm of this DOFVector
     */
    double nrm2() const; 

    /** \brief
     * Returns square of the euclidian norm.
     */
    double squareNrm2() const;

    /** \brief
     * Calculates l2 norm of this DOFVector
     */
    inline double l2norm() const { 
      return nrm2();
    }; 

    /** \brief
     * Calculates the absolute sum of this DOFVector
     */
    T asum() const; 

    /** \brief
     * Calculates the l1 norm of this DOFVector
     */
    inline double l1norm() const { 
      return asum();
    }; 

    // hier Andreas Ergaenzung
  
    /** \brief
     * Calculates the sum of this DOFVector
     */
    T sum() const; 

    // bis hier Andreas E.
  
    /** \brief
     * Sets \ref vec[i] = val, i=0 , ... , size
     */
    void set(T val); 

    /** \brief
     * Assignment operator for setting complete vector to a certain value d
     */
    inline DOFVector<T>& operator=(T d) {
      set(d); 
      return *this;
    }; 

    /** \brief
     * Assignment operator between two vectors
     */
    DOFVector<T>& operator=(const DOFVector<T>& ); 


    //virtual DOFVector<T>& operator-=(const DOFVector<T>& rhs) { 
    //  axpy(-1.0, rhs); return *this;
    //};

    //virtual DOFVector<T>& operator+=(const DOFVector<T>& rhs) {
    //  axpy(1.0, rhs); return *this;
    //};

    /** \brief
     * Multiplies all entries of \ref vec with s
     */
    //virtual void scal(T s); 

    /** \brief
     * Multiplication with a scalar
     */
    //virtual const DOFVector<T>& operator*(T d) const {
    //  static DOFVector<T> vec(this->feSpace, "operator*-vector");
    //  vec = *this;
    //  vec.scal(d);
    //  return vec;
    //};
 
    /** \brief
     * \ref vec[i] = d * \ref vec[i]
     */
    //virtual const DOFVector<T>& operator*=(T d) { 
    //  scal(d); 
    //  return (*this);
    //}; 

    /** \brief
     * vec[i] = v.vec[i]
     */
    void copy(const DOFVector<T>& v); 

    /** \brief
     * vec[i] = a * x.vec[i] + vec[i]
     */
    //void axpy(T a,const DOFVector<T>& x); 

    /** \brief
     * vec[i] = a * x.vec[i] + y.vec[i]
     */
    //void axpy(T a,const DOFVector<T>& x, const DOFVector<T>& y);

    /** \brief
     * vec[i] = x.vec[i] + a * vec[i]
     */
    //void xpay(T a,const DOFVector<T>& x); 

    /** \brief
     * Returns minimum of DOFVector
     */
    T min() const; 

    /** \brief
     * Returns maximum of DOFVector
     */
    T max() const; 

    /** \brief
     * Used by interpol while mesh traversal
     */
    static int interpolFct(ElInfo* elinfo);

    /** \brief
     * Generalized matrix-vector multiplication
     */
    //virtual void gemv(MatrixTranspose transpose, T alpha, 
    //		    const DOFMatrix &a, const DOFVector<T>& x, 
    //		    T beta); 

    /** \brief
     * Matrix-vector multiplication
     */
    //virtual void mv(MatrixTranspose transpose, 
    //		  const DOFMatrix&a, const DOFVector<T>&x); 
   
    /** \brief
     * Prints \ref vec to stdout
     */
    void print() const; 

    //   ElementVector *assemble(T factor, ElInfo *elInfo,
    // 			  const BoundaryType *bound, 
    // 			  Operator *op = NULL);

    //   /** \brief
    //    * Adds element contribution to the vector
    //    */
    //   void addElementVector(T sign,
    // 			const ElementVector &elVec, 
    // 			const BoundaryType *bound,
    // 			bool add = true); 
 


    //   inline void addOperator(Operator* op, 
    // 			  double *factor = NULL,
    // 			  double *estFactor = NULL) 
    //   {
    //     operators.push_back(op);
    //     operatorFactor.push_back(factor);
    //     operatorEstFactor.push_back(estFactor);
    //   };

    //   inline ::std::vector<double*>::iterator getOperatorFactorBegin() {
    //     return operatorFactor.begin();
    //   };

    //   inline ::std::vector<double*>::iterator getOperatorFactorEnd() {
    //     return operatorFactor.end();
    //   };

    //   inline ::std::vector<double*>::iterator getOperatorEstFactorBegin() {
    //     return operatorEstFactor.begin();
    //   };

    //   inline ::std::vector<double*>::iterator getOperatorEstFactorEnd() {
    //     return operatorEstFactor.end();
    //   };


    //   inline ::std::vector<Operator*>::iterator getOperatorsBegin() {
    //     return operators.begin();
    //   };

    //   inline ::std::vector<Operator*>::iterator getOperatorsEnd() {
    //     return operators.end();
    //   };

    //   Flag getAssembleFlag();

    //   /** \brief
    //    * Evaluates \f[ u_h(x(\lambda)) = \sum_{i=0}^{m-1} vec[ind[i]] * 
    //    * \varphi^i(\lambda) \f] where \f$ \varphi^i \f$ is the i-th basis function,
    //    * \f$ x(\lambda) \f$ are the world coordinates of lambda and
    //    * \f$ m \f$ is the number of basis functions
    //    */
    //   T evalUh(const DimVec<double>& lambda, DegreeOfFreedom* ind);

    /** \brief
     * Computes the coefficients of the interpolant of the function fct and
     * stores these in the DOFVector
     */
    void interpol(AbstractFunction<T, WorldVector<double> > *fct);

    void interpol(DOFVector<T> *v, double factor);

    //   inline ::std::vector<Operator*>& getOperators() { return operators; };

    //   inline ::std::vector<double*>& getOperatorFactor() { return operatorFactor; };
    //   inline ::std::vector<double*>& getOperatorEstFactor() { 
    //     return operatorEstFactor; 
    //   };


    // ===== Serializable implementation =====

    void serialize(::std::ostream &out) 
    {
      unsigned int size = vec.size();
      out.write(reinterpret_cast<const char*>(&size), sizeof(unsigned int));
      out.write(reinterpret_cast<const char*>(&(vec[0])), size * sizeof(T));
    };

    void deserialize(::std::istream &in) {
      unsigned int size;
      in.read(reinterpret_cast<char*>(&size), sizeof(unsigned int));
      vec.resize(size);
      in.read(reinterpret_cast<char*>(&(vec[0])), size * sizeof(T));
    };

    DOFVector<WorldVector<T> > *getGradient(DOFVector<WorldVector<T> >*) const;
    //   {
    //     FUNCNAME("DOFVector<T>::getGradient()");
    //     ERROR_EXIT("only specialized for DOFVector<double>\n");
    //     return NULL;
    //   };

    WorldVector<DOFVector<T>*> *getGradient(WorldVector<DOFVector<T>*> *grad) const;

    DOFVector<WorldVector<T> >*
    getRecoveryGradient(DOFVector<WorldVector<T> >*) const;
    //   {
    //     FUNCNAME("DOFVector<T>::getRecoveryGradient()");
    //     ERROR_EXIT("only specialized for DOFVector<double>\n");
    //     return NULL;
    //   };

    //const T *getLocalVector(const Element *el, T* localVec);

    //   const T *getVecAtQPs(const ElInfo         *elInfo, 
    // 		       const Quadrature     *quad,
    // 		       const FastQuadrature *quadFast,
    // 		       T                    *vecAtQPs) const;

    //   const WorldVector<T> *getGrdAtQPs(const ElInfo         *elInfo, 
    // 				    const Quadrature     *quad,
    // 				    const FastQuadrature *quadFast,
    // 				    WorldVector<T>       *grdAtQPs) const;
    //   {
    //     FUNCNAME("DOFVector<T>::getGrdAtQPs()");
    //     ERROR_EXIT("only specialized for DOFVector<double>\n");
    //     return NULL;
    //   };

    //   const WorldMatrix<T> *getD2AtQPs(const ElInfo         *elInfo, 
    // 				   const Quadrature     *quad,
    // 				   const FastQuadrature *quadFast,
    // 				   WorldMatrix<T>       *d2AtQPs) const;
    //   {
    //     FUNCNAME("DOFVector<T>::getD2AtQPs()");
    //     ERROR_EXIT("only specialized for DOFVector<double>\n");
    //     return NULL;
    //   };


  protected:

    //Hier auch Ergaenzungen (Andreas)
  
    /** \brief
     * Used by Int while mesh traversal
     */
    static int Int_fct(ElInfo* elinfo);

    /** \brief
     * Used by L1Norm while mesh traversal
     */
    static int L1Norm_fct(ElInfo* elinfo);

    // Letzte Erganzung in diesem File (Andreas)


  
    /** \brief
     * Used by L2Norm while mesh traversal
     */
    static int L2Norm_fct(ElInfo* elinfo);

    /** \brief
     * Used by H1Norm while mesh traversal
     */
    static int H1Norm_fct(ElInfo* elinfo);

  protected: 
    /** \brief
     * Name of this DOFVector
     */
    ::std::string          name; 

    /** \brief
     * FiniteElemSpace of the vector
     */
    const FiniteElemSpace *feSpace; 

    /** \brief
     * Data container
     */
    ::std::vector<T> vec; 
 
    //   ElementVector *elementVector;

    /** \brief
     * Specifies whether interpolation should be performed after refinement
     */
    bool refineInter;

    /** \brief
     * Specifies what operation should be performed after coarsening
     */
    CoarsenOperation coarsenOperation;

    //   ::std::vector<Operator*> operators;

    //   ::std::vector<double*> operatorFactor;

    //   ::std::vector<double*> operatorEstFactor;

    /** \brief
     * Used by \ref interpol
     */
    AbstractFunction<T, WorldVector<double> > *interFct;

    /** \brief
     * Used for mesh traversal
     */
    static DOFVector<T> *traverseVector;

  protected:
    /** \brief
     * Used for mesh traversal
     */
    //   static DOFVector<T> *traverseVector;

    /** \brief
     * Used while calculating vector norms
     */
    static FastQuadrature *quad_fast;

    /** \brief
     * Stores the last calculated vector norm
     */
    static double norm;

    /** \brief
     * Dimension of the mesh this DOFVector belongs to
     */
    static int dim;
  }; 

  template<>
  void DOFVector<double>::refineInterpol(RCNeighbourList&, int); 

  template<>
  void DOFVector<double>::coarseRestrict(RCNeighbourList&, int); 

  inline double min(const DOFVector<double>& v) {return v.min();}; 
  inline double max(const DOFVector<double>& v) {return v.max();}; 

  // ===========================================================================
  // ===== class DOFVectorDOF ==================================================
  // ===========================================================================

  /** \ingroup DOFAdministration
   * \brief
   * A DOFVector that stores DOF indices.
   */
  class DOFVectorDOF : public DOFVector<DegreeOfFreedom>,
		       public DOFContainer
  {
  public:  
    MEMORY_MANAGED(DOFVectorDOF);

    /** \brief
     * Calls constructor of DOFVector<DegreeOfFreedom> and registers itself
     * as DOFContainer at DOFAdmin
     */
    DOFVectorDOF(const FiniteElemSpace* feSpace_, ::std::string name_)
      : DOFVector<DegreeOfFreedom>(feSpace_, name_)
    {
      feSpace->getAdmin()->addDOFContainer(this);
    };
  
    /** \brief
     * Deregisters itself at DOFAdmin.
     */
    ~DOFVectorDOF() {
      feSpace->getAdmin()->removeDOFContainer(this);
    };

    /** \brief
     * Implements DOFContainer::operator[]() by calling 
     * DOFVEctor<DegreeOfFreedom>::operator[]()
     */
    DegreeOfFreedom& operator[](DegreeOfFreedom i) {
      return DOFVector<DegreeOfFreedom>::operator[](i);
    };

    /** \brief
     * Implements DOFIndexedBase::getSize()
     */
    int getSize() const {
      return DOFVector<DegreeOfFreedom>::getSize();
    };

    /** \brief
     * Implements DOFIndexedBase::resize()
     */
    void resize(int size) {
      DOFVector<DegreeOfFreedom>::resize(size);
    }

    void freeDOFContent(DegreeOfFreedom dof);

  protected:
    DOFVectorDOF();
  };

  template<typename T>
  double norm(DOFVector<T> *vec) {
    return vec->nrm2();
  };

  template<typename T>
  double L2Norm(DOFVector<T> *vec) {
    return vec->L2Norm();
  };

  template<typename T>
  double H1Norm(DOFVector<T> *vec) {
    return vec->H1Norm();
  };

  template<typename T>
  void print(DOFVector<T> *vec) {
    vec->print();
  };

  // point wise multiplication
  template<typename T>
  const DOFVector<T>& operator*=(DOFVector<T>& x, const DOFVector<T>& y);

  // multiplication with scalar

  template<typename T>
  const DOFVector<T>& operator*=(DOFVector<T>& x, T scal);

  // scalar product

  template<typename T>
  T operator*(DOFVector<T>& x, DOFVector<T>& y);

  // addition

  template<typename T>
  const DOFVector<T>& operator+=(DOFVector<T>& x, const DOFVector<T>& y);

  // subtraction

  template<typename T>
  const DOFVector<T>& operator-=(DOFVector<T>& x, const DOFVector<T>& y);

  template<typename T>
  const DOFVector<T>& operator*(const DOFVector<T>& v, double d);

  template<typename T>
  const DOFVector<T>& operator*(double d, const DOFVector<T>& v);

  template<typename T>
  const DOFVector<T>& operator+(const DOFVector<T>&v1 , const DOFVector<T>& v2);



  // y = a*x + y

  template<typename T>
  void axpy(double a,const DOFVector<T>& x, DOFVector<T>& y);

  // matrix vector product
  template<typename T>
  void mv(MatrixTranspose transpose, 
	  const DOFMatrix &a, 
	  const DOFVector<T> &x,
	  DOFVector<T> &result,
	  bool add = false); 

  template<typename T>
  void xpay(double a,const DOFVector<T>& x,DOFVector<T>& y);

  template<typename T>
  inline void scal(T a, DOFVector<T>& y) {y*=a;};

  template<typename T>
  inline const DOFVector<T>& mult(double scal,
				  const DOFVector<T>& v,
				  DOFVector<T>& result);

  template<typename T>
  inline const DOFVector<T>& add(const DOFVector<T>& v,
				 double scal,
				 DOFVector<T>& result);

  template<typename T>
  inline const DOFVector<T>& add(const DOFVector<T>& v1,
				 const DOFVector<T>& v2,
				 DOFVector<T>& result);

  template<typename T>
  inline void set(DOFVector<T>& vec, T d) 
  {
    vec.set(d);
  };

  template<typename T>
  inline void setValue(DOFVector<T>& vec, T d) 
  {
    vec.set(d);
  };

  template<typename T>
  inline int size(DOFVector<T> *vec) {
    return vec->getUsedSize();
  };

  //template<typename T>
  //void interpolGrd(DOFVector<T> *v, 
  //		 WorldVector<DOFVector<T>*> *grad,
  //		 double factor, bool add = false);

  WorldVector<DOFVector<double>*> *transform(DOFVector<WorldVector<double> > *vec,
					     WorldVector<DOFVector<double>*> *result);
#if 0
  void transform2(DOFVector<WorldVector<double> > *vec,
		  ::std::vector<DOFVector<double>*> *result);
#endif

}

#include "DOFVector.hh"

#endif  // !_DOFVECTOR_H_