DOFVector.h

/******************************************************************************
 *
 * AMDiS - Adaptive multidimensional simulations
 *
 * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
 * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
 *
 * Authors: 
 * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
 *
 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 *
 *
 * This file is part of AMDiS
 *
 * See also license.opensource.txt in the distribution.
 * 
 ******************************************************************************/



/** \file DOFVector.h */

#ifndef AMDIS_DOFVECTOR_H 
#define AMDIS_DOFVECTOR_H 
 
#include <vector> 
#include <memory> 
#include <map> 

#include "AMDiS_fwd.h"
#include "FixVec.h"
#include "Global.h" 
#include "Flag.h" 
#include "RCNeighbourList.h" 
#include "DOFIterator.h"
#include "DOFIndexed.h"
#include "DOFContainer.h"
#include "Boundary.h"
#include "CreatorInterface.h"
#include "Serializable.h"
#include "DOFMatrix.h" 
#include "BasisFunction.h"
#include "FiniteElemSpace.h"
#include "SurfaceQuadrature.h"
#include "Traits.h"

namespace AMDiS {

  template<typename T> 
  class DOFVectorBase : public DOFIndexed<T>
  {
  public:

    DOFVectorBase() 
      : feSpace(nullptr),
	elementVector(3),
        boundaryManager(nullptr),
        nBasFcts(0)
    {}
    
    DOFVectorBase(const FiniteElemSpace *f, std::string n);

    virtual ~DOFVectorBase();

    /// For the given element, this function returns an array of all DOFs of 
    /// this DOFVector that are defined on this element.
    virtual void getLocalVector(const Element *el, 
				mtl::dense_vector<T>& localVec) const;

    /// Evaluates the DOF vector at a set of quadrature points defined on the 
    /// given element.
    void getVecAtQPs(const ElInfo *elInfo, 
		     const Quadrature *quad,
		     const FastQuadrature *quadFast,
		     mtl::dense_vector<T>& vecAtQPs) const;

    void getVecAtQPs(const ElInfo *smallElInfo, 
		     const ElInfo *largeElInfo,
		     const Quadrature *quad,
		     const FastQuadrature *quadFast,
		     mtl::dense_vector<T>& vecAtQPs) const;

    /// Evaluates the gradient of a DOF vector at a set of quadrature points
    /// defined on the given element.
    void getGrdAtQPs( const ElInfo *elInfo,
		      const Quadrature *quad,
		      const FastQuadrature *quadFast,
		      mtl::dense_vector<typename GradientType<T>::type> &grdAtQPs) const;

    void getGrdAtQPs( const ElInfo *smallElInfo,
		      const ElInfo *largeElInfo,
		      const Quadrature *quad,
		      const FastQuadrature *quadFast,
		      mtl::dense_vector<typename GradientType<T>::type> &grdAtQPs) const;

    /// Evaluates the comp'th component of the derivative of a DOF vector at a
    /// set of quadrature points defined on the given element.
    void getDerivativeAtQPs( const ElInfo *elInfo,
		      const Quadrature *quad,
		      const FastQuadrature *quadFast,
		      int comp,
		      mtl::dense_vector<T> &derivativeAtQPs) const;

    void getDerivativeAtQPs( const ElInfo *smallElInfo,
		      const ElInfo *largeElInfo,
		      const Quadrature *quad,
		      const FastQuadrature *quadFast,
		      int comp,
		      mtl::dense_vector<T> &derivativeAtQPs) const;

    /// Evaluates the jacobian of a DOF vector at a set of quadrature points
    ///  defined on the given element.
    void getD2AtQPs(const ElInfo *elInfo,
		    const Quadrature *quad,
		    const FastQuadrature *quadFast,
		    mtl::dense_vector<typename D2Type<T>::type> &d2AtQPs) const;

    /// Returns the FE space the DOF vector is defined on.
    inline const FiniteElemSpace* getFeSpace() const 
    {
      return feSpace;
    }

    /// Assembles the element vector for the given ellement and adds the
    /// element matrix to the current DOF vector.
    void assemble(T factor, ElInfo *elInfo,			    
		  const BoundaryType *bound, 
		  Operator *op = nullptr);

    void assemble2(T factor, 
		   ElInfo *mainElInfo, ElInfo *auxElInfo,
		   ElInfo *smallElInfo, ElInfo *largeElInfo,
		   const BoundaryType *bound, 
		   Operator *op = nullptr);

    void addElementVector(T sign,
			  const ElementVector& elVec, 
			  const BoundaryType *bound,
			  ElInfo *elInfo,
			  bool add = true); 

    ///
    void assembleOperator(Operator &op);
    
    /// That function must be called after the matrix assembling has been
    /// finished. This makes it possible to start some cleanup or matrix
    /// data compressing procedures.
    void finishAssembling();
 
    inline void addOperator(Operator* op, 
			    double *factor = nullptr,
			    double *estFactor = nullptr) 
    {
      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();

    /// 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; 
    }

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

    inline void setName(std::string n)
    {
      name = n;
    }

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

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

    inline void setDirichletDofValue(DegreeOfFreedom dof,
				     T value)
    {
      dirichletDofValues[dof] = value;
    }

    std::map<DegreeOfFreedom, T>& getDirichletValues()
    {
      return dirichletDofValues;
    }

  protected:
    ///
    const FiniteElemSpace *feSpace;

    ///
    std::string name;

    ///
    ElementVector elementVector;

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

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

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

    ///
    BoundaryManager *boundaryManager;

    /// Number of basis functions of the used finite element space.
    int nBasFcts;

    /// Dimension of the mesh this DOFVectorBase belongs to
    int dim;

    std::map<DegreeOfFreedom, T> dirichletDofValues;
  };


  /// Specifies which operation should be done after coarsening
  typedef enum{
    NO_OPERATION = 0,   
    COARSE_RESTRICT = 1,
    COARSE_INTERPOL = 2, 
    REFINE_INTERPOL = 4
  } RefineCoarsenOperation;
 

  /** \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:
    typedef typename DOFVectorBase<T>::value_type	value_type;
    typedef typename DOFVectorBase<T>::size_type		size_type;
    typedef typename DOFVectorBase<T>::reference		reference;
    typedef typename DOFVectorBase<T>::const_reference	const_reference;
  
  public:
    /** \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:
      Creator(FiniteElemSpace *feSpace_) 
        : feSpace(feSpace_) 
      {}

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

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

    private:
      FiniteElemSpace *feSpace;
    };

  public:
    /// Empty constructor. No initialization!
    DOFVector() 
      : DOFVectorBase<T>(), 
	coarsenOperation(NO_OPERATION)
    {}

    /// Constructs a DOFVector with name n belonging to FiniteElemSpace f
    DOFVector(const FiniteElemSpace* f, std::string n, bool addToSynch = true); 

    /// Initialization.
    void init(const FiniteElemSpace* f, std::string n, bool addToSynch = true);

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

    /// Destructor
    virtual ~DOFVector();

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

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

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

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

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

    /// Returns \ref coarsenOperation
    inline RefineCoarsenOperation getCoarsenOperation() 
    { 
      return coarsenOperation; 
    }

    /// Sets \ref refineOperation to op
    inline void setRefineOperation(RefineCoarsenOperation op) 
    { 
      refineOperation = op; 
    }

    /// Returns \ref refineOperation
    inline RefineCoarsenOperation getRefineOperation() 
    { 
      return refineOperation; 
    }
    
    /// Restriction after coarsening. Implemented for DOFVector<double>
    inline void coarseRestrict(RCNeighbourList&, int) {}

    /// Interpolation after refinement. Implemented for DOFVector<double>
    inline void refineInterpol(RCNeighbourList&, int) {}

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

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

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

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

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

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

    /// Returns \ref vec[i]
    inline reference operator[](DegreeOfFreedom i) 
    {
      FUNCNAME_DBG("DOFVector<T>::operator[]");

      TEST_EXIT_DBG(i >= 0 && i < static_cast<int>(vec.size())) 
 	("Illegal vector index %d.\n", i); 

      return vec[i];
    }
 
    /// Calculates Integral of this DOFVector
    double Int(Quadrature* q = nullptr) const
    {
      return Int(-1, q);
    }

    /** \brief
     * Calculates Integral of this DOFVector. 
     *
     * \param[in]  meshlevel  Then mesh level on which the integral should be
     *                        calculated. Usually, only -1 for the leaf level
     *                        makes sence, but other values can be used, e.g.,
     *                        for debugging.
     * \param[in]  q          Quadrature object. If not specified, the function
     *                        creates a new quadrature object.
     */
    double Int(int meshLevel, Quadrature* q = nullptr) const;


    /// Calculates Integral of this DOFVector over parts of the domain
    /// boundary, indicated by boundaryType. Implemented for DOFVector<double>
    T IntOnBoundary(BoundaryType boundary, Quadrature* q = nullptr) const
    {
      FUNCNAME("DOFVector::IntOnBoundary())");
      TEST_EXIT(false)("Please implement your integration\n");
      return 0.0;
    }

    /// Calculates Integral of this DOFVector times normal vector over parts 
    /// of the domain boundary, indicated by boundaryType. Implemented for 
    /// DOFVector<WorldVector<double> >
    double IntOnBoundaryNormal(BoundaryType boundary, Quadrature* q = nullptr) const
    {
      FUNCNAME("DOFVector::IntOnBoundaryNormal())");
      TEST_EXIT(false)("Please implement your integration\n");
      return 0.0;
    }

    /// Calculates L1 norm of this DOFVector
    double L1Norm(Quadrature* q = nullptr) const;
 
    /// Calculates L2 norm of this DOFVector
    inline double L2Norm(Quadrature* q = nullptr) const 
    {
      return sqrt(L2NormSquare());
    }

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

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

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

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

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

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

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

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

    /// Calculates doublewell of this DOFVector
    double DoubleWell(Quadrature* q = nullptr) const;
 
    /// Calculates the sum of this DOFVector
    T sum() const; 
 
    /// Sets \ref vec[i] = val, i=0 , ... , size
    void set(T val); 

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

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

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

    /// Returns minimum of DOFVector.
    T min() const; 

    /// Returns maximum of DOFVector.
    T max() const;

    /// Returns absolute maximum of DOFVector.
    T absMax() const;

    /// Returns the average value of the DOFVector.
    T average() const;

    /// Prints \ref vec to stdout.
    void print() const; 

    ///
    size_t calcMemoryUsage() const;

    /// 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 = 1.0);


    /// Eval DOFVector at given point p. If oldElInfo != nullptr the search for 
    /// the element, where p is inside, starts from oldElInfo. implemented for:
    /// double, WorldVector< double >
    T evalAtPoint(WorldVector<double> &p, 
		  ElInfo *oldElInfo = nullptr) const 
    {
      FUNCNAME("DOFVector::evalAtPoint())");
      TEST_EXIT(false)("Please implement your evaluation\n");
    }

    /// Determine the DegreeOfFreedom that has coords with minimal euclidean 
    /// distance to WorldVector p. return true if DOF is found, and false 
    /// otherwise.
    bool getDofIdxAtPoint(WorldVector<double> &p, 
				DegreeOfFreedom &idx, 
				ElInfo *oldElInfo = nullptr, 
				bool useOldElInfo = false) const;

    template<typename S>
    void serialize(std::ostream &out, boost::mpl::identity<S>) 
    {
      unsigned int size = static_cast<unsigned int>(vec.size());
      out.write(reinterpret_cast<const char*>(&size), sizeof(unsigned int));
      out.write(reinterpret_cast<const char*>(&(vec[0])), size * sizeof(T));
    }
    
    // no serialization for bool vectors
    void serialize(std::ostream &out, boost::mpl::identity<bool>) {}
    
    /// Writes the data of the DOFVector to an output stream.
    void serialize(std::ostream &out) 
    {
      serialize(out, boost::mpl::identity<T>());
    }
    
    template<typename S>
    void deserialize(std::istream &in, boost::mpl::identity<S>) 
    {
      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));
    }
    
    // no de-serialization for bool vectors
    void deserialize(std::istream &in, boost::mpl::identity<bool>) {}

    /// Reads data of an DOFVector from an input stream.
    void deserialize(std::istream &in) 
    {
      deserialize(in, boost::mpl::identity<T>());
    }

    DOFVector<typename GradientType<T>::type>*
      getGradient(DOFVector<typename GradientType<T>::type> *grad) const;

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

    DOFVector<typename GradientType<T>::type>*
      getRecoveryGradient(DOFVector<typename GradientType<T>::type> *grad) const;

  protected: 

    /// Data container
    std::vector<T> vec; 
 
    /// Specifies what operation should be performed after coarsening
    RefineCoarsenOperation coarsenOperation;
    RefineCoarsenOperation refineOperation;
  }; 


  template<>
  double DOFVector<double>::IntOnBoundary(
    BoundaryType boundaryType, Quadrature* q) const;

  template<>
  double DOFVector<WorldVector<double> >::IntOnBoundaryNormal(
    BoundaryType boundaryType, Quadrature* q) const;

  template<>
  double DOFVector<double>::evalAtPoint(WorldVector<double> &p, 
					 ElInfo *oldElInfo) const;

  template<>
  WorldVector<double> DOFVector<WorldVector<double> >::evalAtPoint(WorldVector<double> &p, 
								   ElInfo *oldElInfo) const;

  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();
  }


  /** \ingroup DOFAdministration
   * \brief
   * A DOFVector that stores DOF indices.
   */
  class DOFVectorDOF : public DOFVector<DegreeOfFreedom>,
		       public DOFContainer
  {
  public:  
    /// 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);
    }
  
    /// Deregisters itself at DOFAdmin.
    ~DOFVectorDOF() 
    {
      feSpace->getAdmin()->removeDOFContainer(this);
    }

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

    const DegreeOfFreedom& operator[](DegreeOfFreedom i) const 
    {
      return DOFVector<DegreeOfFreedom>::operator[](i);
    }

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

    /// 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 const DOFVector<T>& mod(const DOFVector<T>& v,
				 DOFVector<T>& result);

  template<typename T>
  inline const DOFVector<T>& Tanh(const DOFVector<T>& v,
				  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>
//   inline int size(const DOFVector<T>& vec) 
//   {
//     return vec.getUsedSize();
//   }