/******************************************************************************
 *
 * 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.
 *
 ******************************************************************************/


#include <list>
#include <algorithm>
#include <math.h>

#include <boost/numeric/mtl/mtl.hpp>
#include <boost/numeric/mtl/utility/tag.hpp>
#include <boost/numeric/mtl/utility/category.hpp>
#include <boost/numeric/linear_algebra/identity.hpp>

#include "FixVec.h"
#include "Boundary.h"
#include "DOFAdmin.h"
#include "ElInfo.h"
#include "Error.h"
#include "FiniteElemSpace.h"
#include "Global.h"
#include "Mesh.h"
#include "Quadrature.h"
#include "AbstractFunction.h"
#include "BoundaryManager.h"
#include "Assembler.h"
#include "Operator.h"
#include "Initfile.h"
#include "Traverse.h"
#include "DualTraverse.h"

#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
#include "parallel/MpiHelper.h"
#include "parallel/MeshDistributor.h"
#endif

// Defining the interface for MTL4
namespace mtl {

  // Let MTL4 know that DOFVector it is a column vector
  namespace traits {
    template <typename T>
    struct category< AMDiS::DOFVector<T> >
    {
      typedef tag::dense_col_vector type;
    };
  }

  namespace ashape {
    template <typename T>
    struct ashape< AMDiS::DOFVector<T> >
    {
      typedef cvec<typename ashape<T>::type> type;
    };
  }

  // Modelling Collection and MutableCollection
  template <typename T>
  struct Collection< AMDiS::DOFVector<T> >
  {
    typedef T               value_type;
    typedef const T&        const_reference;
    typedef std::size_t     size_type;
  };

  template <typename T>
  struct MutableCollection< AMDiS::DOFVector<T> >
    : public Collection< AMDiS::DOFVector<T> >
  {
    typedef T&              reference;
  };


} // namespace mtl



namespace AMDiS {

  template<typename T>
  DOFVectorBase<T>::DOFVectorBase(const FiniteElemSpace *f, std::string n)
    : feSpace(f),
      name(n),
      elementVector(f->getBasisFcts()->getNumber()),
      boundaryManager(NULL)
  {
    nBasFcts = feSpace->getBasisFcts()->getNumber();
    dim = feSpace->getMesh()->getDim();
  }


  template<typename T>
  DOFVectorBase<T>::~DOFVectorBase()
  {}


  template<typename T>
  DOFVector<T>::DOFVector(const FiniteElemSpace* f, std::string n, bool addToSynch)
    : DOFVectorBase<T>(f, n)
  {
    vec.resize(0);
    init(f, n, addToSynch);
  }


  template<typename T>
  void DOFVector<T>::init(const FiniteElemSpace* f, std::string n, bool addToSynch)
  {
    this->name = n;
    this->feSpace = f;
    
    if (this->feSpace && this->feSpace->getAdmin())
      (this->feSpace->getAdmin())->addDOFIndexed(this);
    
    this->boundaryManager = new BoundaryManager(f);
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    if (addToSynch && Parallel::MeshDistributor::globalMeshDistributor != NULL)
      Parallel::MeshDistributor::globalMeshDistributor->addInterchangeVector(this);
#endif
  }

  template<typename T>
  DOFVector<T>::~DOFVector()
  {
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    if ( Parallel::MeshDistributor::globalMeshDistributor != NULL)
      Parallel::MeshDistributor::globalMeshDistributor->removeInterchangeVector(this);
#endif
    
    #pragma omp critical
    if (this->feSpace && this->feSpace->getAdmin())
      (this->feSpace->getAdmin())->removeDOFIndexed(this);

    if (this->boundaryManager)
      delete this->boundaryManager;

    vec.clear();
  }

  template<typename T>
  void DOFVectorBase<T>::addElementVector(T factor,
					  const ElementVector &elVec,
					  const BoundaryType *bound,
					  ElInfo *elInfo,
					  bool add)
  {
    std::vector<DegreeOfFreedom> indices(nBasFcts);
    feSpace->getBasisFcts()->getLocalIndices(elInfo->getElement(),
					     feSpace->getAdmin(),
					     indices);

    for (int i = 0; i < nBasFcts; i++) {
      BoundaryCondition *condition =
	bound ? this->getBoundaryManager()->getBoundaryCondition(bound[i]) : NULL;

      if (!(condition && condition->isDirichlet())) {
	DegreeOfFreedom irow = indices[i];

	if (add)
	  (*this)[irow] += factor * elVec[i];
	else
	  (*this)[irow] = factor * elVec[i];
      }
    }
  }

  template<typename T>
  double DOFVector<T>::nrm2() const
  {
    checkFeSpace(this->feSpace, vec);

    double nrm = 0.0;
    Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(this)),
			 USED_DOFS);
    for (vecIterator.reset(); !vecIterator.end(); ++vecIterator)
      nrm += (*vecIterator) * (*vecIterator);

#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localNrm = nrm;
    MPI::COMM_WORLD.Allreduce(&localNrm, &nrm, 1, MPI_DOUBLE, MPI_SUM);
#endif

    return std::sqrt(nrm);
  }


  template<typename T>
  double DOFVector<T>::squareNrm2() const
  {
    checkFeSpace(this->feSpace, vec);

    double nrm = 0.0;
    Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(this)),
			 USED_DOFS);
    for (vecIterator.reset(); !vecIterator.end(); ++vecIterator)
      nrm += (*vecIterator) * (*vecIterator);

#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localNrm = nrm;
    MPI::COMM_WORLD.Allreduce(&localNrm, &nrm, 1, MPI_DOUBLE, MPI_SUM);
#endif

    return nrm;
  }


  template<typename T>
  T DOFVector<T>::asum() const
  {
    checkFeSpace(this->feSpace, vec);

    double nrm = 0.0;
    Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(this)),
			 USED_DOFS);
    for (vecIterator.reset(); !vecIterator.end(); ++vecIterator)
      nrm += abs(*vecIterator);

#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localNrm = nrm;
    MPI::COMM_WORLD.Allreduce(&localNrm, &nrm, 1, MPI_DOUBLE, MPI_SUM);
#endif

    return nrm;
  }


  template<typename T>
  T DOFVector<T>::sum() const
  {
    checkFeSpace(this->feSpace, vec);

    double nrm = 0.0;
    Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(this)),
			 USED_DOFS);
    for (vecIterator.reset(); !vecIterator.end(); ++vecIterator)
      nrm += *vecIterator;

#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localNrm = nrm;
    MPI::COMM_WORLD.Allreduce(&localNrm, &nrm, 1, MPI_DOUBLE, MPI_SUM);
#endif

    return nrm;
  }


  template<typename T>
  void DOFVector<T>::set(T alpha)
  {
    checkFeSpace(this->feSpace, vec);

    Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(this), USED_DOFS);
    for (vecIterator.reset(); !vecIterator.end(); ++vecIterator)
      *vecIterator = alpha ;
  }


  template<typename T>
  void DOFVector<T>::copy(const DOFVector<T>& x)
  {
    FUNCNAME_DBG("DOFVector<T>::copy()");

    checkFeSpace(this->feSpace, vec);

    TEST_EXIT_DBG(static_cast<int>(x.vec.size()) >=
		  this->feSpace->getAdmin()->getUsedSize())
      ("x.size = %d too small: admin->sizeUsed = %d\n", x.vec.size(),
       this->feSpace->getAdmin()->getUsedSize());

    Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(this), USED_DOFS);
    Iterator xIterator(dynamic_cast<DOFVector<T>*>(const_cast<DOFVector<T>*>(&x)),
		       USED_DOFS);
    for (vecIterator.reset(), xIterator.reset(); !vecIterator.end();
	 ++vecIterator, ++xIterator)
      *vecIterator = *xIterator;
  }


  template<typename T>
  T DOFVector<T>::min() const
  {
    checkFeSpace(this->feSpace, vec);

    T m;
    Iterator vecIterator(const_cast<DOFIndexed<T>*>(dynamic_cast<const DOFIndexed<T>*>(this)), USED_DOFS);
    for (vecIterator.reset(), m = *vecIterator; !vecIterator.end(); ++vecIterator)
      m = std::min(m, *vecIterator);

#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localMin = m;
    MPI::COMM_WORLD.Allreduce(&localMin, &m, 1, MPI_DOUBLE, MPI_MIN);
#endif

    return m;
  }


  template<typename T>
  T DOFVector<T>::max() const
  {
    checkFeSpace(this->feSpace, vec);

    T m;
    Iterator vecIterator(const_cast<DOFIndexed<T>*>(dynamic_cast<const DOFIndexed<T>*>(this)), USED_DOFS);
    for (vecIterator.reset(), m = *vecIterator; !vecIterator.end(); ++vecIterator)
      m = std::max(m, *vecIterator);

#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localMax = m;
    MPI::COMM_WORLD.Allreduce(&localMax, &m, 1, MPI_DOUBLE, MPI_MAX);
#endif

    return m;
  }


  template<typename T>
  T DOFVector<T>::absMax() const
  {
    return std::max(abs(max()), abs(min()));
  }


  template<typename T>
  T DOFVector<T>::average() const
  {
    checkFeSpace(this->feSpace, vec);

    int count = 0;
    T m = 0;
    Iterator vecIterator(const_cast<DOFIndexed<T>*>(dynamic_cast<const DOFIndexed<T>*>(this)), USED_DOFS);
    for (vecIterator.reset(); !vecIterator.end(); ++vecIterator) {
      m += *vecIterator;
      count++;
    }

#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localSum = m;
    int localCount = count;
    MPI::COMM_WORLD.Allreduce(&localSum, &m, 1, MPI_DOUBLE, MPI_SUM);
    MPI::COMM_WORLD.Allreduce(&localCount, &count, 1, MPI_INT, MPI_SUM);
#endif

    return m / count;
  }


  template<typename T>
  void DOFVector<T>::print() const
  {
    FUNCNAME("DOFVector<T>::print()");

    const DOFAdmin *admin = NULL;
    const char *format;

    if (this->feSpace)
      admin = this->feSpace->getAdmin();

    MSG("Vec `%s':\n", this->name.c_str());
    int j = 0;
    if (admin) {
      if (admin->getUsedSize() > 100)
	format = "%s(%3d,%10.5le)";
      else if (admin->getUsedSize() > 10)
	format = "%s(%2d,%10.5le)";
      else
	format = "%s(%1d,%10.5le)";

      Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(this)), USED_DOFS);
      for (vecIterator.reset(); !vecIterator.end(); ++vecIterator) {
	if ((j % 3) == 0) {
	  if (j)
	    Msg::print("\n");
	  MSG(format, "", vecIterator.getDOFIndex(), *vecIterator);
	} else {
	  Msg::print(format, " ", vecIterator.getDOFIndex(), *vecIterator);
	}
	j++;
      }
      Msg::print("\n");
    } else {
	MSG("no DOFAdmin, print whole vector.\n");

	for (int i = 0; i < static_cast<int>( vec.size()); i++) {
	  if ((j % 3) == 0) {
	    if (j)
	      Msg::print("\n");
	    MSG("(%d,%10.5e)",i,vec[i]);
	  } else {
	    Msg::print(" (%d,%10.5e)",i,vec[i]);
	  }
	  j++;
	}
	Msg::print("\n");
      }
    return;
  }


  template<typename T>
  size_t DOFVector<T>::calcMemoryUsage() const
  {
    size_t result = 0;
    result += sizeof(DOFVector<T>);
    result += vec.size() * sizeof(T);

    return result;
  }


  template<typename T>
  T DOFVectorBase<T>::evalUh(const DimVec<double>& lambda,
			     DegreeOfFreedom* dof_indices)
  {
    BasisFunction* phi = const_cast<BasisFunction*>(this->getFeSpace()->getBasisFcts());
    int nBasisFcts = phi->getNumber();
    T val = 0.0;

    for (int i = 0; i < nBasisFcts; i++)
      val += (*this)[dof_indices[i]]*(*phi->getPhi(i))(lambda);

#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localVal = val;
    MPI::COMM_WORLD.Allreduce(&localVal, &val, 1, MPI_DOUBLE, MPI_SUM);
#endif

    return val;
  }


  template<typename T>
  void DOFVector<T>::interpol(AbstractFunction<T, WorldVector<double> > *fct)
  {
    FUNCNAME("DOFVector::interpol()");

    TEST_EXIT_DBG(fct)("No function to interpolate!\n");

    if (!this->getFeSpace()) {
      MSG("no dof admin in vec %s, skipping interpolation\n",
	  this->getName().c_str());
      return;
    }

    if (!(this->getFeSpace()->getAdmin())) {
      MSG("no dof admin in feSpace %s, skipping interpolation\n",
	  this->getFeSpace()->getName().c_str());
      return;
    }

    if (!(this->getFeSpace()->getBasisFcts())) {
      MSG("no basis functions in admin of vec %s, skipping interpolation\n",
	  this->getName().c_str());
      return;
    }

    if (!(fct)) {
      MSG("function that should be interpolated only pointer to NULL, ");
      Msg::print("skipping interpolation\n");
      return;
    }

    const BasisFunction *basFct = this->getFeSpace()->getBasisFcts();
    const DOFAdmin* admin = this->getFeSpace()->getAdmin();
    int nBasFcts = basFct->getNumber();
    std::vector<DegreeOfFreedom> myLocalIndices(nBasFcts);
    mtl::dense_vector<T> fctInterpolValues(nBasFcts);

    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(this->getFeSpace()->getMesh(), -1,
					 Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
    while (elInfo) {
      basFct->interpol(elInfo, 0, NULL, fct, fctInterpolValues);
      basFct->getLocalIndices(const_cast<Element*>(elInfo->getElement()),
			      admin, myLocalIndices);
      for (int i = 0; i < nBasFcts; i++)
	vec[myLocalIndices[i]] = fctInterpolValues[i];

      elInfo = stack.traverseNext(elInfo);
    }
  }


  template<typename T>
  double DOFVector<T>::Int(int meshLevel, Quadrature* q) const
  {
    Mesh* mesh = this->feSpace->getMesh();

    if (!q) {
      int deg = 2 * this->feSpace->getBasisFcts()->getDegree();
      q = Quadrature::provideQuadrature(this->dim, deg);
    }

    FastQuadrature *quadFast =
      FastQuadrature::provideFastQuadrature(this->feSpace->getBasisFcts(), *q, INIT_PHI);

    Flag flag = Mesh::FILL_COORDS | Mesh::FILL_DET;
    if (meshLevel == -1)
      flag |= Mesh::CALL_LEAF_EL;
    else
      flag |= Mesh::CALL_EL_LEVEL;

    T result; nullify(result);
    int nPoints = quadFast->getNumPoints();
    mtl::dense_vector<T> uh_vec(nPoints);
    TraverseStack stack;
    ElInfo *elInfo =  stack.traverseFirst(mesh, meshLevel, flag);
    while (elInfo) {
      double det = elInfo->getDet();
      T normT; nullify(normT);
      this->getVecAtQPs(elInfo, NULL, quadFast, uh_vec);
      for (int iq = 0; iq < nPoints; iq++)
	normT += quadFast->getWeight(iq) * (uh_vec[iq]);
      result += det * normT;

      elInfo = stack.traverseNext(elInfo);
    }

#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    Parallel::mpi::globalAdd(result);
#endif

    return result;
  }


  template<typename TOut>
  TOut integrate_Coords(const FiniteElemSpace* feSpace,
		   AbstractFunction<TOut, WorldVector<double> > *fct)
  {
    FUNCNAME("integrate_Coords()");

    TEST_EXIT(fct)("No function defined!\n");

    int deg = std::max(fct->getDegree(), feSpace->getBasisFcts()->getDegree());
    Quadrature* quad =
      Quadrature::provideQuadrature(feSpace->getMesh()->getDim(), deg);
    FastQuadrature *fastQuad =
      FastQuadrature::provideFastQuadrature(feSpace->getBasisFcts(), *quad, INIT_PHI);

    WorldVector<double> coords;

    TOut value; nullify(value);
    Flag traverseFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET;
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(feSpace->getMesh(), -1, traverseFlag);
    while (elInfo) {
    TOut tmp; nullify(tmp);
      for (int iq = 0; iq < fastQuad->getNumPoints(); iq++) {
	elInfo->coordToWorld(fastQuad->getLambda(iq), coords);
 	tmp += (*fct)(coords) * fastQuad->getWeight(iq);
      }

      value += tmp * elInfo->getDet();

      elInfo = stack.traverseNext(elInfo);
    }

#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    Parallel::mpi::globalAdd(value);
#endif

    return value;
  }


  template<typename TOut, typename T>
  TOut integrate_Vec(const DOFVector<T> &vec,
		  AbstractFunction<TOut, T> *fct)
  {
    FUNCNAME("integrate_Vec()");

    TEST_EXIT(fct)("No function defined!\n");

    int deg = std::max(fct->getDegree(),
		       2 * vec.getFeSpace()->getBasisFcts()->getDegree());
    Quadrature* quad =
      Quadrature::provideQuadrature(vec.getFeSpace()->getMesh()->getDim(), deg);
    FastQuadrature *fastQuad =
      FastQuadrature::provideFastQuadrature(vec.getFeSpace()->getBasisFcts(), *quad, INIT_PHI);

    mtl::dense_vector<T> qp(fastQuad->getNumPoints());

    TOut value; nullify(value);

    Flag traverseFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET;
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(vec.getFeSpace()->getMesh(), -1, traverseFlag);
    while (elInfo) {
      vec.getVecAtQPs(elInfo, quad, fastQuad, qp);
      TOut tmp; nullify(tmp);
      for (int iq = 0; iq < fastQuad->getNumPoints(); iq++) {
 	tmp += fastQuad->getWeight(iq) * (*fct)(qp[iq]);
      }
      value += tmp * elInfo->getDet();

      elInfo = stack.traverseNext(elInfo);
    }

#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    Parallel::mpi::globalAdd(value);
#endif

    return value;
  }


  template<typename TOut, typename T1, typename T2>
  TOut integrate_Vec2(const DOFVector<T1> &vec1,
		   const DOFVector<T2> &vec2,
		   BinaryAbstractFunction<TOut, T1, T2> *fct)
  {
    if (vec1.getFeSpace()->getMesh() == vec2.getFeSpace()->getMesh())
      return int_Vec2_SingleMesh(vec1, vec2, fct);
    else
      return int_Vec2_DualMesh(vec1, vec2, fct);
  }


  template<typename TOut, typename T1, typename T2>
  TOut int_Vec2_SingleMesh(const DOFVector<T1> &vec1,
		   const DOFVector<T2> &vec2,
		   BinaryAbstractFunction<TOut, T1, T2> *fct)
  {
    FUNCNAME("intDualmesh()");

    TEST_EXIT(fct)("No function defined!\n");

    int deg = std::max(fct->getDegree(),
		       2 * vec1.getFeSpace()->getBasisFcts()->getDegree());
    Quadrature* quad =
      Quadrature::provideQuadrature(vec1.getFeSpace()->getMesh()->getDim(), deg);
    FastQuadrature *fastQuad =
      FastQuadrature::provideFastQuadrature(vec1.getFeSpace()->getBasisFcts(), *quad, INIT_PHI);

    mtl::dense_vector<T1> qp1(fastQuad->getNumPoints());
    mtl::dense_vector<T2> qp2(fastQuad->getNumPoints());

    TOut value; nullify(value);
    Flag traverseFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET;
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(vec1.getFeSpace()->getMesh(), -1, traverseFlag);
    while (elInfo) {
      vec1.getVecAtQPs(elInfo, quad, fastQuad, qp1);
      vec2.getVecAtQPs(elInfo, quad, fastQuad, qp2);

      TOut tmp; nullify(tmp);
      for (int iq = 0; iq < fastQuad->getNumPoints(); iq++)
 	tmp += fastQuad->getWeight(iq) * (*fct)(qp1[iq], qp2[iq]);
      value += tmp * elInfo->getDet();

      elInfo = stack.traverseNext(elInfo);
    }

#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    Parallel::mpi::globalAdd(value);
#endif

    return value;
  }


  template<typename TOut, typename T1, typename T2>
  TOut int_Vec2_DualMesh(const DOFVector<T1> &vec1,
		   const DOFVector<T2> &vec2,
		   BinaryAbstractFunction<TOut, T1, T2> *fct)
  {
    FUNCNAME("intDualmesh()");

    TEST_EXIT(fct)("No function defined!\n");

    int deg = std::max(fct->getDegree(),
		       2 * vec1.getFeSpace()->getBasisFcts()->getDegree());
    Quadrature* quad =
      Quadrature::provideQuadrature(vec1.getFeSpace()->getMesh()->getDim(), deg);
    FastQuadrature *fastQuad =
      FastQuadrature::provideFastQuadrature(vec1.getFeSpace()->getBasisFcts(), *quad, INIT_PHI);

    mtl::dense_vector<T1> qp1(fastQuad->getNumPoints());
    mtl::dense_vector<T2> qp2(fastQuad->getNumPoints());

    TOut value; nullify(value);
    Flag traverseFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET;
    DualTraverse dualTraverse;
    DualElInfo dualElInfo;
    bool cont = dualTraverse.traverseFirst(vec1.getFeSpace()->getMesh(),
					   vec2.getFeSpace()->getMesh(),
					   -1, -1, traverseFlag, traverseFlag,
					   dualElInfo);
    while (cont) {
      vec1.getVecAtQPs(dualElInfo.smallElInfo, dualElInfo.largeElInfo, quad, NULL, qp1);
      vec2.getVecAtQPs(dualElInfo.smallElInfo, dualElInfo.largeElInfo, quad, NULL, qp2);

      TOut tmp; nullify(tmp);
      for (int iq = 0; iq < quad->getNumPoints(); iq++)
 	tmp += quad->getWeight(iq) * (*fct)(qp1[iq], qp2[iq]);
      value += tmp * dualElInfo.smallElInfo->getDet();

      cont = dualTraverse.traverseNext(dualElInfo);
    }

#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    Parallel::mpi::globalAdd(value);
#endif

    return value;
  }


  template<typename T1, typename T2>
  typename traits::mult_type<T1,T2>::type
  integrate_VecTimesCoords(const DOFVector<T1> &vec,
			   AbstractFunction<T2, WorldVector<double> > *fct)
  {
    FUNCNAME("integrate_VecTimesCoords()");

    TEST_EXIT(fct)("No function defined!\n");

    typedef typename traits::mult_type<T1,T2>::type TOut;

    const FiniteElemSpace *feSpace = vec.getFeSpace();
    Mesh *mesh = feSpace->getMesh();

    int deg = std::max(fct->getDegree(), 2 * feSpace->getBasisFcts()->getDegree());
    Quadrature* quad = Quadrature::provideQuadrature(mesh->getDim(), deg);
    FastQuadrature *fastQuad =
      FastQuadrature::provideFastQuadrature(feSpace->getBasisFcts(), *quad, INIT_PHI);

    mtl::dense_vector<T1> qp(fastQuad->getNumPoints());
    WorldVector<double> coords;

    TOut value; nullify(value);
    Flag traverseFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET;
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, traverseFlag);
    while (elInfo) {
      vec.getVecAtQPs(elInfo, quad, fastQuad, qp);

      TOut tmp; nullify(tmp);
      for (int iq = 0; iq < fastQuad->getNumPoints(); iq++) {
	elInfo->coordToWorld(fastQuad->getLambda(iq), coords);
 	tmp += fastQuad->getWeight(iq) * (qp[iq] * (*fct)(coords));
      }

      value += tmp * elInfo->getDet();

      elInfo = stack.traverseNext(elInfo);
    }

#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    Parallel::mpi::globalAdd(value);
#endif

    return value;
  }


  template<typename TOut, typename T>
  TOut integrate_VecAndCoords(const DOFVector<T> &vec,
			   BinaryAbstractFunction<TOut, T, WorldVector<double> > *fct)
  {
    FUNCNAME("integrate_VecAndCoords()");

    TEST_EXIT(fct)("No function defined!\n");

    int deg = std::max(fct->getDegree(),
		       2 * vec.getFeSpace()->getBasisFcts()->getDegree());
    Quadrature* quad =
      Quadrature::provideQuadrature(vec.getFeSpace()->getMesh()->getDim(), deg);
    FastQuadrature *fastQuad =
      FastQuadrature::provideFastQuadrature(vec.getFeSpace()->getBasisFcts(), *quad, INIT_PHI);

    mtl::dense_vector<T> qp(fastQuad->getNumPoints());
    WorldVector<double> coordsAtQP;

    TOut value; nullify(value);

    Flag traverseFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET;
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(vec.getFeSpace()->getMesh(), -1, traverseFlag);
    while (elInfo) {
      vec.getVecAtQPs(elInfo, quad, fastQuad, qp);
      TOut tmp; nullify(tmp);
      for (int iq = 0; iq < fastQuad->getNumPoints(); iq++) {
	elInfo->coordToWorld(fastQuad->getLambda(iq), coordsAtQP);
 	tmp += fastQuad->getWeight(iq) * (*fct)(qp[iq], coordsAtQP);
      }
      value += tmp * elInfo->getDet();

      elInfo = stack.traverseNext(elInfo);
    }

#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    Parallel::mpi::globalAdd(value);
#endif

    return value;
  }


  template<typename T>
  double DOFVector<T>::L1Norm(Quadrature* q) const
  {
    Mesh* mesh = this->feSpace->getMesh();

    if (!q) {
      int deg = 2 * this->feSpace->getBasisFcts()->getDegree();
      q = Quadrature::provideQuadrature(this->dim, deg);
    }

    FastQuadrature *quadFast =
      FastQuadrature::provideFastQuadrature(this->feSpace->getBasisFcts(), *q, INIT_PHI);

    double result = 0.0;
    int nPoints = quadFast->getNumPoints();
    mtl::dense_vector<T> uh_vec(nPoints);
    TraverseStack stack;
    ElInfo *elInfo =
      stack.traverseFirst(mesh, -1,
			  Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET);
    while (elInfo) {
      double det = elInfo->getDet();
      double normT = 0.0;
      this->getVecAtQPs(elInfo, NULL, quadFast, uh_vec);
      for (int iq = 0; iq < nPoints; iq++)
	normT += quadFast->getWeight(iq) * abs(uh_vec[iq]);
      result += det * normT;

      elInfo = stack.traverseNext(elInfo);
    }

#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localResult = result;
    MPI::COMM_WORLD.Allreduce(&localResult, &result, 1, MPI_DOUBLE, MPI_SUM);
#endif

    return result;
  }


  template<typename T>
  double DOFVector<T>::L2NormSquare(Quadrature* q) const
  {
    Mesh* mesh = this->feSpace->getMesh();

    if (!q) {
      int deg = 2 * this->feSpace->getBasisFcts()->getDegree();
      q = Quadrature::provideQuadrature(this->dim, deg);
    }

    FastQuadrature *quadFast =
      FastQuadrature::provideFastQuadrature(this->feSpace->getBasisFcts(), *q, INIT_PHI);

    double result = 0.0;
    int nPoints = quadFast->getNumPoints();
    mtl::dense_vector<T> uh_vec(nPoints);
    TraverseStack stack;
    ElInfo *elInfo =
      stack.traverseFirst(mesh, -1,
			  Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET);
    while (elInfo) {
      double det = elInfo->getDet();
      double normT = 0.0;
      this->getVecAtQPs(elInfo, NULL, quadFast, uh_vec);
      for (int iq = 0; iq < nPoints; iq++)
	normT += quadFast->getWeight(iq) * sqr(uh_vec[iq]);
      result += det * normT;

      elInfo = stack.traverseNext(elInfo);
    }

#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localResult = result;
    MPI::COMM_WORLD.Allreduce(&localResult, &result, 1, MPI_DOUBLE, MPI_SUM);
#endif

    return result;
  }


  template<typename T>
  double DOFVector<T>::H1NormSquare(Quadrature *q) const
  {
    Mesh* mesh = this->feSpace->getMesh();

    if (!q) {
      int deg = 2 * this->feSpace->getBasisFcts()->getDegree() - 2;
      q = Quadrature::provideQuadrature(this->dim, deg);
    }

    FastQuadrature *quadFast =
      FastQuadrature::provideFastQuadrature(this->feSpace->getBasisFcts(), *q, INIT_GRD_PHI);

    double result = 0.0;
    int nPoints = quadFast->getNumPoints();
    int dimOfWorld = Global::getGeo(WORLD);
    mtl::dense_vector<WorldVector<T> > grduh_vec(nPoints);
    TraverseStack stack;
    ElInfo *elInfo =
      stack.traverseFirst(mesh, -1,
			  Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS |
			  Mesh::FILL_DET | Mesh::FILL_GRD_LAMBDA);
    while (elInfo) {
      double det = elInfo->getDet();
      double normT = 0.0;
      this->getGrdAtQPs(elInfo, NULL, quadFast, grduh_vec);

      for (int iq = 0; iq < nPoints; iq++) {
	double norm2 = 0.0;
	for (int j = 0; j < dimOfWorld; j++)
	  norm2 += sqr(grduh_vec[iq][j]);
	normT += quadFast->getWeight(iq) * norm2;
      }

      result += det * normT;

      elInfo = stack.traverseNext(elInfo);
    }

#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localResult = result;
    MPI::COMM_WORLD.Allreduce(&localResult, &result, 1, MPI_DOUBLE, MPI_SUM);
#endif

    return result;
  }



  template<typename T>
  bool DOFVector<T>::getDofIdxAtPoint(WorldVector<double> &p,
					    DegreeOfFreedom &idx,
					    ElInfo *oldElInfo,
					    bool useOldElInfo) const
  {
    Mesh *mesh = this->feSpace->getMesh();
    const BasisFunction *basFcts = this->feSpace->getBasisFcts();

    int dim = mesh->getDim();
    int numBasFcts = basFcts->getNumber();

    std::vector<DegreeOfFreedom> localIndices(numBasFcts);
    DimVec<double> lambda(dim, NO_INIT);

    ElInfo *elInfo = mesh->createNewElInfo();
    idx = 0;
    bool inside = false;

    if (oldElInfo && useOldElInfo && oldElInfo->getMacroElement()) {
      inside = mesh->findElInfoAtPoint(p, elInfo, lambda, oldElInfo->getMacroElement(), NULL, NULL);
      delete oldElInfo;
    } else {
      inside = mesh->findElInfoAtPoint(p, elInfo, lambda, NULL, NULL, NULL);
    }

    if (oldElInfo)
      oldElInfo = elInfo;

    if (!inside){
      delete elInfo;
      return false;
    }

    basFcts->getLocalIndices(elInfo->getElement(), this->feSpace->getAdmin(), localIndices);

    WorldVector<double> coord;
    int minIdx = -1;
    double minDist = 1.e15;

    for (int i = 0; i < numBasFcts; i++) {
      elInfo->coordToWorld(*(basFcts->getCoords(i)), coord);
      WorldVector<double> dist = coord - p;
      double newDist = norm(dist);
      if (newDist < minDist) {
	minIdx = i;
	minDist = newDist;
	if (minDist < 1.e-10)
	  break;
      }
    }

    if (minIdx >= 0)
      idx = localIndices[minIdx];

    if(!oldElInfo) delete elInfo;
    return inside;
  }


  template<typename T>
  void DOFVector<T>::compressDOFIndexed(int first, int last,
					std::vector<DegreeOfFreedom> &newDOF)
  {
    for (int i = first; i <= last; i++)
      if (newDOF[i] >= 0)
	vec[newDOF[i]] = vec[i];
  }


  template<typename T>
  Flag DOFVectorBase<T>::getAssembleFlag()
  {
    Flag fillFlag(0);

    for (std::vector<Operator*>::iterator op = this->operators.begin();
	 op != this->operators.end(); ++op)
      fillFlag |= (*op)->getFillFlag();

    return fillFlag;
  }


  template<typename T>
  void DOFVectorBase<T>::finishAssembling()
  {
    // call the operatos cleanup procedures
    for (std::vector<Operator*>::iterator it = this->operators.begin();
	 it != this->operators.end(); ++it)
      (*it)->finishAssembling();
  }




  template<typename T>
  void DOFVector<T>::refineInterpolImpl(RCNeighbourList& list, int n, id<double>)
  {
    switch (DOFIndexedBase::refineOperation) {
    case NO_OPERATION:
      return;
      break;
    case REFINE_INTERPOL:
    default:
      (const_cast<BasisFunction*>(this->feSpace->getBasisFcts()))->refineInter(this, &list, n);
      break;
    }
  }


  template<typename T>
  void DOFVector<T>::refineInterpolImpl(RCNeighbourList& list, int n, id<WorldVector<double> >)
  {
    if (DOFIndexedBase::refineOperation == NO_OPERATION)
      return;

    if (n < 1)
      return;

    Element *el = list.getElement(0);
    int n0 = this->feSpace->getAdmin()->getNumberOfPreDofs(VERTEX);
    DegreeOfFreedom dof0 = el->getDof(0, n0);
    DegreeOfFreedom dof1 = el->getDof(1, n0);
    DegreeOfFreedom dof_new = el->getChild(0)->getDof(this->feSpace->getMesh()->getDim(), n0);
    vec[dof_new] = vec[dof0];
    vec[dof_new] += vec[dof1];
    vec[dof_new] *= 0.5;
  }



  template<typename T>
  double DOFVector<T>::evalAtPointImpl(WorldVector<double> const& p,
					ElInfo *oldElInfo, id<double>) const
  {
    Mesh *mesh = this->feSpace->getMesh();
    const BasisFunction *basFcts = this->feSpace->getBasisFcts();

    int dim = mesh->getDim();
    int nBasFcts = basFcts->getNumber();

    std::vector<DegreeOfFreedom> localIndices(nBasFcts);
    DimVec<double> lambda(dim, NO_INIT);

    ElInfo *elInfo = mesh->createNewElInfo();
    double value = 0.0;
    bool inside = false;

    if (oldElInfo && oldElInfo->getMacroElement()) {
      inside = mesh->findElInfoAtPoint(p, elInfo, lambda, oldElInfo->getMacroElement(), NULL, NULL);
      delete oldElInfo;
    } else
      inside = mesh->findElInfoAtPoint(p, elInfo, lambda, NULL, NULL, NULL);

    if (oldElInfo)
      oldElInfo = elInfo;

    if (inside) {
      basFcts->getLocalIndices(elInfo->getElement(), this->feSpace->getAdmin(), localIndices);
      ElementVector uh(nBasFcts);
      for (int i = 0; i < nBasFcts; i++)
        uh[i] = operator[](localIndices[i]);
      value = basFcts->evalUh(lambda, uh);
    } else {
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
      value = 0.0;
#else
      ERROR_EXIT("Can not eval DOFVector at point p, because point is outside geometry.");
#endif
    }


    if (oldElInfo == NULL)
      delete elInfo;

#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    Parallel::mpi::globalAdd(value);
#endif
    return value;
  }


  template<typename T>
  WorldVector<double> DOFVector<T>::evalAtPointImpl(WorldVector<double> const& p,
								   ElInfo *oldElInfo, id<WorldVector<double> >) const
  {
    Mesh *mesh = this->feSpace->getMesh();
    const BasisFunction *basFcts = this->feSpace->getBasisFcts();

    int dim = mesh->getDim();
    int nBasFcts = basFcts->getNumber();

    std::vector<DegreeOfFreedom> localIndices(nBasFcts);
    DimVec<double> lambda(dim, NO_INIT);
    ElInfo *elInfo = mesh->createNewElInfo();
    WorldVector<double> value(DEFAULT_VALUE, 0.0);
    bool inside = false;

    if (oldElInfo && oldElInfo->getMacroElement()) {
      inside = mesh->findElInfoAtPoint(p, elInfo, lambda, oldElInfo->getMacroElement(), NULL, NULL);
      delete oldElInfo;
    } else
      inside = mesh->findElInfoAtPoint(p, elInfo, lambda, NULL, NULL, NULL);

    if (oldElInfo)
      oldElInfo = elInfo;

    if (inside) {
      basFcts->getLocalIndices(elInfo->getElement(), this->feSpace->getAdmin(), localIndices);
      mtl::dense_vector<WorldVector<double> > uh(nBasFcts);
      for (int i = 0; i < nBasFcts; i++)
        uh[i] = operator[](localIndices[i]);
      value = basFcts->evalUh(lambda, uh);
    } else {
      ERROR_EXIT("Can not eval DOFVector at point p, because point is outside geometry.");
    }

    if (oldElInfo == NULL)
      delete elInfo;

    return value;
  }

  template<typename T>
  DOFVector<T>& DOFVector<T>::operator=(const DOFVector<T>& rhs)
  {
    this->feSpace = rhs.feSpace;
    this->dim = rhs.dim;
    this->nBasFcts = rhs.nBasFcts;
    vec = rhs.vec;
    this->elementVector.change_dim(this->nBasFcts);
    this->operators = rhs.operators;
    this->operatorFactor = rhs.operatorFactor;

    if (rhs.boundaryManager) {
      if (this->boundaryManager)
	delete this->boundaryManager;

      this->boundaryManager = new BoundaryManager(*rhs.boundaryManager);
    } else {
      this->boundaryManager = NULL;
    }

    return *this;
  }


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

    TEST_EXIT_DBG(x.getFeSpace() && x.getFeSpace()->getAdmin())
      ("pointer is NULL: %8X, %8X\n", x.getFeSpace(), x.getFeSpace()->getAdmin());

    typename DOFVector<T>::Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(&x),
						USED_DOFS);
    for (vecIterator.reset(); !vecIterator.end(); ++vecIterator)
      (*vecIterator) *= scal;

    return x;
  }


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

    TEST_EXIT_DBG(x.getFeSpace() && y.getFeSpace())
      ("feSpace is NULL: %8X, %8X\n", x.getFeSpace(), y.getFeSpace());
    TEST_EXIT_DBG(x.getFeSpace()->getAdmin() &&
	      (x.getFeSpace()->getAdmin() == y.getFeSpace()->getAdmin()))
      ("no admin or different admins: %8X, %8X\n",
       x.getFeSpace()->getAdmin(), y.getFeSpace()->getAdmin());
    TEST_EXIT_DBG(x.getSize() == y.getSize())("different sizes\n");

    typename DOFVector<T>::Iterator xIterator(dynamic_cast<DOFIndexed<T>*>(&x), USED_DOFS);
    typename DOFVector<T>::Iterator yIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(&y)), USED_DOFS);
    for (xIterator.reset(), yIterator.reset(); !xIterator.end();
	 ++xIterator, ++yIterator)
      *xIterator += *yIterator;

    return x;
  }


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

    TEST_EXIT_DBG(x.getFeSpace() && y.getFeSpace())
      ("feSpace is NULL: %8X, %8X\n", x.getFeSpace(), y.getFeSpace());
    TEST_EXIT_DBG(x.getFeSpace()->getAdmin() &&
	      (x.getFeSpace()->getAdmin() == y.getFeSpace()->getAdmin()))
      ("no admin or different admins: %8X, %8X\n",
       x.getFeSpace()->getAdmin(), y.getFeSpace()->getAdmin());
    TEST_EXIT_DBG(x.getSize() == y.getSize())("different sizes\n");

    typename DOFVector<T>::Iterator xIterator(dynamic_cast<DOFIndexed<T>*>(&x), USED_DOFS);
    typename DOFVector<T>::Iterator yIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(&y)), USED_DOFS);
    for (xIterator.reset(), yIterator.reset(); !xIterator.end();
	 ++xIterator, ++yIterator)
      *xIterator -= *yIterator;

    return x;
  }


  template<typename T>
  const DOFVector<T>& operator*=(DOFVector<T>& x, const DOFVector<T>& y)
  {
    FUNCNAME_DBG("DOFVector<T>::operator*=(DOFVector<T>& x, const DOFVector<T>& y)");

    TEST_EXIT_DBG(x.getFeSpace() && y.getFeSpace())
      ("feSpace is NULL: %8X, %8X\n", x.getFeSpace(), y.getFeSpace());
    TEST_EXIT_DBG(x.getFeSpace()->getAdmin() &&
	      (x.getFeSpace()->getAdmin() == y.getFeSpace()->getAdmin()))
      ("no admin or different admins: %8X, %8X\n",
       x.getFeSpace()->getAdmin(), y.getFeSpace()->getAdmin());
    TEST_EXIT_DBG(x.getSize() == y.getSize())("different sizes\n");

    typename DOFVector<T>::Iterator xIterator(dynamic_cast<DOFIndexed<T>*>(&x), USED_DOFS);
    typename DOFVector<T>::Iterator yIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(&y)), USED_DOFS);
    for (xIterator.reset(), yIterator.reset(); !xIterator.end();
	 ++xIterator, ++yIterator)
      *xIterator *= *yIterator;

    return x;
  }


  template<typename T>
  T operator*(DOFVector<T>& x, DOFVector<T>& y)
  {
    FUNCNAME("DOFVector<T>::operator*(DOFVector<T>& x, DOFVector<T>& y)");
    const DOFAdmin *admin = NULL;

    TEST_EXIT(x.getFeSpace() && y.getFeSpace())
      ("feSpace is NULL: %8X, %8X\n", x.getFeSpace(), y.getFeSpace());
    TEST_EXIT((admin = x.getFeSpace()->getAdmin()) && (admin == y.getFeSpace()->getAdmin()))
      ("no admin or different admins: %8X, %8X\n",
       x.getFeSpace()->getAdmin(), y.getFeSpace()->getAdmin());
    TEST_EXIT(x.getSize() == y.getSize())("different sizes\n");

    T dot = 0;

    typename DOFVector<T>::Iterator xIterator(dynamic_cast<DOFIndexed<T>*>(&x), USED_DOFS);
    typename DOFVector<T>::Iterator yIterator(dynamic_cast<DOFIndexed<T>*>(&y), USED_DOFS);
    for (xIterator.reset(), yIterator.reset(); !xIterator.end();
	 ++xIterator, ++yIterator)
      dot += (*xIterator) * (*yIterator);

    return dot;
  }


  template<typename T>
  void mv(MatrixTranspose transpose, const DOFMatrix &a, const DOFVector<T>&x,
	  DOFVector<T> &result, bool add)
  {
    FUNCNAME("DOFVector<T>::mv()");

    TEST_EXIT(a.getRowFeSpace() && a.getColFeSpace() &&
	      x.getFeSpace() && result.getFeSpace())
      ("getFeSpace() is NULL: %8X, %8X, %8X, %8X\n",
       a.getRowFeSpace(), a.getColFeSpace(), x.getFeSpace(), result.getFeSpace());

    const DOFAdmin *rowAdmin = a.getRowFeSpace()->getAdmin();
    const DOFAdmin *colAdmin = a.getColFeSpace()->getAdmin();

    TEST_EXIT((rowAdmin && colAdmin) &&
	      (((transpose == NoTranspose) && (colAdmin == x.getFeSpace()->getAdmin()) &&
		(rowAdmin == result.getFeSpace()->getAdmin()))||
	       ((transpose == Transpose) && (rowAdmin == x.getFeSpace()->getAdmin()) &&
		(colAdmin == result.getFeSpace()->getAdmin()))))
      ("no admin or different admins: %8X, %8X, %8X, %8X\n",
       a.getRowFeSpace()->getAdmin(), a.getColFeSpace()->getAdmin(),
       x.getFeSpace()->getAdmin(), result.getFeSpace()->getAdmin());


    if (transpose == NoTranspose) {
      mtl::dense_vector<T> tmp_x(x.getUsedSize());
      mtl::dense_vector<T> tmp_result(result.getUsedSize());

      {
	int counter = 0;
	typename DOFVector<T>::Iterator it(const_cast<DOFVector<T>*>(&x), USED_DOFS);
	for (it.reset(); !it.end(); ++it)
	  tmp_x[counter++] = *it;
      }

      mult(a.getBaseMatrix(), tmp_x, tmp_result);

      {
	int counter = 0;
	typename DOFVector<T>::Iterator it(&result, USED_DOFS);
	for (it.reset(); !it.end(); ++it)
	  if (add)
	    *it += tmp_result[counter++];
	  else
	    *it = tmp_result[counter++];
      }


    } else if (transpose == Transpose) {
      ERROR_EXIT("Not yet implemented!\n");

//       mult(trans(const_cast<DOFMatrix::base_matrix_type&>(a.getBaseMatrix())),
// 	   x, result);
    }
    else
      ERROR_EXIT("transpose = %d\n", transpose);
  }


  template<typename T>
  void axpy(double alpha, const DOFVector<T>& x, DOFVector<T>& y)
  {
    FUNCNAME("DOFVector<T>::axpy()");

    TEST_EXIT(x.getFeSpace() && y.getFeSpace())
      ("feSpace is NULL: %8X, %8X\n", x.getFeSpace(), y.getFeSpace());

    const DOFAdmin *admin = x.getFeSpace()->getAdmin();

    TEST_EXIT((admin) && (admin == y.getFeSpace()->getAdmin()))
      ("no admin or different admins: %8X, %8X\n",
       x.getFeSpace()->getAdmin(), y.getFeSpace()->getAdmin());
    TEST_EXIT(static_cast<int>(x.getSize()) >= admin->getUsedSize())
      ("size = %d too small: admin->size = %d\n", x.getSize(),
       admin->getUsedSize());
    TEST_EXIT(static_cast<int>(y.getSize()) >= admin->getUsedSize())
      ("y.size = %d too small: admin->size = %d\n", y.getSize(),
       admin->getUsedSize());

    for (int i = 0; i < y.getSize(); i++)
      if (!admin->isDofFree(i))
	y[i] += alpha * x[i];
  }


  template<typename T>
  const DOFVector<T>& operator*(const DOFVector<T>& v, double d)
  {
    static DOFVector<T> result; // TODO: REMOVE STATIC
    return mult(d, v, result);
  }


  template<typename T>
  const DOFVector<T>& operator*(double d, const DOFVector<T>& v)
  {
    static DOFVector<T> result; // TODO: REMOVE STATIC
    return mult(d, v, result);
  }


  template<typename T>
  const DOFVector<T>& operator+(const DOFVector<T> &v1 , const DOFVector<T> &v2)
  {
    static DOFVector<T> result; // TODO: REMOVE STATIC
    return add(v1, v2, result);
  }


  template<typename T>
  void xpay(double alpha,
	    const DOFVector<T>& x,
	    DOFVector<T>& y)
  {
    FUNCNAME("DOFVector<T>::xpay()");

    TEST_EXIT(x.getFeSpace() && y.getFeSpace())
      ("feSpace is NULL: %8X, %8X\n", x.getFeSpace(), y.getFeSpace());

    const DOFAdmin *admin = x.getFeSpace()->getAdmin();

    TEST_EXIT(admin && (admin == y.getFeSpace()->getAdmin()))
      ("no admin or different admins: %8X, %8X\n",
       x.getFeSpace()->getAdmin(), y.getFeSpace()->getAdmin());
    TEST_EXIT(static_cast<int>(x.getSize()) >= admin->getUsedSize())
      ("size = %d too small: admin->size = %d\n", x.getSize(),
       admin->getUsedSize());
    TEST_EXIT(static_cast<int>(y.getSize()) >= admin->getUsedSize())
      ("y.size = %d too small: admin->size = %d\n", y.getSize(),
       admin->getUsedSize());

    for (int i = 0; i < y.getSize(); i++)
      if (!admin->isDofFree(i))
	y[i] = alpha * y[i] + x[i];
  }


  template<typename T>
  inline const DOFVector<T>& mult(double scal,
				  const DOFVector<T>& v,
				  DOFVector<T>& result)
  {
    typename DOFVector<T>::Iterator vIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(&v)), USED_DOFS);
    typename DOFVector<T>::Iterator rIterator(dynamic_cast<DOFIndexed<T>*>(&result), USED_DOFS);
    for (vIterator.reset(), rIterator.reset();
	 !vIterator.end(); ++vIterator, ++rIterator)
      *rIterator = scal * (*vIterator);

    return result;
  }


  template<typename T>
  inline const DOFVector<T>& add(const DOFVector<T>& v,
				 double scal,
				 DOFVector<T>& result)
  {
    typename DOFVector<T>::Iterator vIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(&v)), USED_DOFS);
    typename DOFVector<T>::Iterator rIterator(dynamic_cast<DOFIndexed<T>*>(&result), USED_DOFS);
    for (vIterator.reset(), rIterator.reset();
	!vIterator.end(); ++vIterator, ++rIterator)
      *rIterator = (*vIterator) + scal;

    return result;
  }


  template<typename T>
  inline const DOFVector<T>& add(const DOFVector<T>& v1,
				 const DOFVector<T>& v2,
				 DOFVector<T>& result)
  {
    typename DOFVector<T>::Iterator v1Iterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(&v1)), USED_DOFS);
    typename DOFVector<T>::Iterator v2Iterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(&v2)), USED_DOFS);
    typename DOFVector<T>::Iterator rIterator(dynamic_cast<DOFIndexed<T>*>(&result), USED_DOFS);
    for (v1Iterator.reset(), v2Iterator.reset(), rIterator.reset();
	 !v1Iterator.end(); ++v1Iterator, ++v2Iterator, ++rIterator)
      *rIterator = (*v1Iterator) + (*v2Iterator);

    return result;
  }


  template<typename T>
  void DOFVectorBase<T>::getLocalVector(const Element *el,
					mtl::dense_vector<T>& d) const
  {
    FUNCNAME_DBG("DOFVectorBase<T>::getLocalVector()");

    TEST_EXIT_DBG(feSpace->getMesh() == el->getMesh())
      ("Element is defined on a different mesh than the DOF vector!\n");

    std::vector<DegreeOfFreedom> localIndices(nBasFcts);
    const DOFAdmin* admin = feSpace->getAdmin();
    feSpace->getBasisFcts()->getLocalIndices(el, admin, localIndices);

    for (int i = 0; i < nBasFcts; i++)
      d[i] = (*this)[localIndices[i]];
  }


  template<typename T>
  void DOFVectorBase<T>::getVecAtQPs(const ElInfo *elInfo,
				     const Quadrature *quad,
				     const FastQuadrature *quadFast,
				     mtl::dense_vector<T>& vecAtQPs) const
  {
    FUNCNAME_DBG("DOFVector<T>::getVecAtQPs()");

    TEST_EXIT_DBG(quad || quadFast)
      ("Neither quad nor quadFast defined!\n");
    TEST_EXIT_DBG(!(quad && quadFast) || quad == quadFast->getQuadrature())
      ("quad != quadFast->quadrature\n");
    TEST_EXIT_DBG(!quadFast || quadFast->getBasisFunctions() == feSpace->getBasisFcts())
      ("Invalid basis functions!");

    const BasisFunction *basFcts = feSpace->getBasisFcts();
    int nBasFcts  = basFcts->getNumber();
    mtl::dense_vector<T> localVec(nBasFcts);
    getLocalVector(elInfo->getElement(), localVec);

    if (quadFast) {
      // using precalculated values at QPs
      const mtl::dense2D<double>& phi = quadFast->getPhi();
      vecAtQPs.change_dim(num_rows(phi)); 	// = quadrature->getNumPoints()
      vecAtQPs = phi * localVec; 		// Matrix<double> * Vector<T>
    } else {
      // evaluate basisFunctions at QPs
      int nPoints = quad->getNumPoints();
      vecAtQPs.change_dim(nPoints);

      for (int iq = 0; iq < nPoints; iq++) {
       	nullify(vecAtQPs[iq]);
       	for (int j = 0; j < nBasFcts; j++)
       	  vecAtQPs[iq] +=
       	    localVec[j] * (*(basFcts->getPhi(j)))(quad->getLambda(iq));
      }
    }
  }


  template<typename T>
  void DOFVectorBase<T>::getVecAtQPs(const ElInfo *smallElInfo,
					  const ElInfo *largeElInfo,
					  const Quadrature *quad,
					  const FastQuadrature *quadFast,
					  mtl::dense_vector<T>& vecAtQPs) const
  {
    FUNCNAME_DBG("DOFVector<T>::getVecAtQPs()");

    TEST_EXIT_DBG(quad || quadFast)("neither quad nor quadFast defined\n");
    TEST_EXIT_DBG(!(quad && quadFast) || quad == quadFast->getQuadrature())
      ("quad != quadFast->quadrature\n");
    TEST_EXIT_DBG(!quadFast || quadFast->getBasisFunctions() == feSpace->getBasisFcts())
      ("invalid basis functions");

    if (smallElInfo->getMesh() == feSpace->getMesh())
      return getVecAtQPs(smallElInfo, quad, quadFast, vecAtQPs);

    const BasisFunction *basFcts = feSpace->getBasisFcts();
    int nBasFcts  = basFcts->getNumber();
    int nPoints = quad->getNumPoints();

    vecAtQPs.change_dim(nPoints);

    mtl::dense_vector<T> localVec(nBasFcts);
    this->getLocalVector(largeElInfo->getElement(), localVec);
    mtl::dense2D<double> &m = smallElInfo->getSubElemCoordsMat(basFcts->getDegree());

    for (int iq = 0; iq < nPoints; iq++) {
      nullify(vecAtQPs[iq]);
      for (int j = 0; j < nBasFcts; j++) {
	double val = 0.0;
	for (int k = 0; k < nBasFcts; k++)
	  val += m[j][k] * (*(basFcts->getPhi(k)))(quad->getLambda(iq));

	vecAtQPs[iq] += localVec[j] * val;
      }
    }
  }


  template<typename T>
  void DOFVectorBase<T>::getGrdAtQPs(const ElInfo *elInfo,
				     const Quadrature *quad,
				     const FastQuadrature *quadFast,
				     mtl::dense_vector<typename GradientType<T>::type> &grdAtQPs) const
  {
    FUNCNAME_DBG("DOFVector<T>::getGrdAtQPs()");

    TEST_EXIT_DBG(quad || quadFast)("neither quad nor quadFast defined\n");
    TEST_EXIT_DBG(!(quad && quadFast) || quad == quadFast->getQuadrature())
      ("quad != quadFast->quadrature\n");
    TEST_EXIT_DBG(!quadFast || quadFast->getBasisFunctions() == feSpace->getBasisFcts())
      ("invalid basis functions");

    const BasisFunction *basFcts = feSpace->getBasisFcts();
    int nBasFcts  = basFcts->getNumber();
    int dow = Global::getGeo(WORLD);
    int nPoints = quadFast ? quadFast->getQuadrature()->getNumPoints() : quad->getNumPoints();

    mtl::dense_vector<T> localVec(nBasFcts);
    this->getLocalVector(elInfo->getElement(), localVec);

    mtl::dense_vector<T> grd1(dim + 1);
    int parts = Global::getGeo(PARTS, dim);
    const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();

    grdAtQPs.change_dim(nPoints);
    if (quadFast) {
      for (int iq = 0; iq < nPoints; iq++) {
	nullify(grd1);

	for (int j = 0; j < nBasFcts; j++) { // #BasisFunctions
	  for (int k = 0; k < parts; k++)  // #edges (2d) or #faces (3d)
	    grd1[k] += quadFast->getGradient(iq, j, k) * localVec[j];
	}

	for (int l = 0; l < dow; l++) {
	  nullify(grdAtQPs[iq][l]);
	  for (int k = 0; k < parts; k++)
	    grdAtQPs[iq][l] += grdLambda[k][l] * grd1[k];
	}
      }

    } else {
      mtl::dense_vector<double> grdPhi(dim + 1);

      for (int iq = 0; iq < nPoints; iq++) {
	nullify(grd1);

	for (int j = 0; j < nBasFcts; j++) {
	  (*(basFcts->getGrdPhi(j)))(quad->getLambda(iq), grdPhi);
	  for (int k = 0; k < parts; k++)
	    grd1[k] += grdPhi[k] * localVec[j];
	}

	for (int l = 0; l < dow; l++) {
	  nullify(grdAtQPs[iq][l]);
	  for (int k = 0; k < parts; k++)
	    grdAtQPs[iq][l] += grdLambda[k][l] * grd1[k];
	}
      }
    }
  }


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

    TEST_EXIT_DBG(quad || quadFast)("neither quad nor quadFast defined\n");
    TEST_EXIT_DBG(!(quad && quadFast) || quad == quadFast->getQuadrature())
      ("quad != quadFast->quadrature\n");
    TEST_EXIT_DBG(!quadFast || quadFast->getBasisFunctions() == feSpace->getBasisFcts())
      ("invalid basis functions");

    if (smallElInfo->getMesh() == feSpace->getMesh())
      return getGrdAtQPs(smallElInfo, quad, quadFast, grdAtQPs);

    const BasisFunction *basFcts = feSpace->getBasisFcts();
    int nBasFcts  = basFcts->getNumber();
    int dow = Global::getGeo(WORLD);
    int nPoints = quadFast ? quadFast->getQuadrature()->getNumPoints() : quad->getNumPoints();

    mtl::dense_vector<T> localVec(nBasFcts);
    this->getLocalVector(largeElInfo->getElement(), localVec);

    mtl::dense2D<double> &m = smallElInfo->getSubElemCoordsMat(basFcts->getDegree());

    int parts = Global::getGeo(PARTS, dim);
    const DimVec<WorldVector<double> > &grdLambda = largeElInfo->getGrdLambda();

    mtl::dense_vector<T> grd1(dim + 1);
    mtl::dense_vector<double> grdPhi(dim + 1);
    mtl::dense_vector<double> tmp(dim + 1);

    grdAtQPs.change_dim(nPoints);
    for (int iq = 0; iq < nPoints; iq++) {
      nullify(grd1);

      for (int j = 0; j < nBasFcts; j++) {
	grdPhi = 0.0;

	for (int k = 0; k < nBasFcts; k++) {
	  (*(basFcts->getGrdPhi(k)))(quad->getLambda(iq), tmp);
	  tmp *= m[j][k];
	  grdPhi += tmp;
	}

	for (int k = 0; k < parts; k++)
	  grd1[k] += grdPhi[k] * localVec[j];
      }

      for (int l = 0; l < dow; l++) {
	nullify(grdAtQPs[iq][l]);
	for (int k = 0; k < parts; k++)
	  grdAtQPs[iq][l] += grdLambda[k][l] * grd1[k];
      }
    }
  }


  template<typename T>
  void DOFVectorBase<T>::getDerivativeAtQPs(const ElInfo *elInfo,
					    const Quadrature *quad,
					    const FastQuadrature *quadFast,
					    int comp,
					    mtl::dense_vector<T> &derivativeAtQPs) const
  {
    FUNCNAME_DBG("DOFVector<T>::getGrdAtQPs()");

    TEST_EXIT_DBG(quad || quadFast)("neither quad nor quadFast defined\n");
    TEST_EXIT_DBG(!(quad && quadFast) || quad == quadFast->getQuadrature())
      ("quad != quadFast->quadrature\n");
    TEST_EXIT_DBG(!quadFast || quadFast->getBasisFunctions() == feSpace->getBasisFcts())
      ("invalid basis functions");

    const BasisFunction *basFcts = feSpace->getBasisFcts();
    int nBasFcts  = basFcts->getNumber();
    int nPoints = quadFast ? quadFast->getQuadrature()->getNumPoints() : quad->getNumPoints();

    mtl::dense_vector<T> localVec(nBasFcts);
    this->getLocalVector(elInfo->getElement(), localVec);

    mtl::dense_vector<T> grd1(dim + 1);
    int parts = Global::getGeo(PARTS, dim);
    const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();

    derivativeAtQPs.change_dim(nPoints);
    if (quadFast) {
      for (int iq = 0; iq < nPoints; iq++) {
	nullify(grd1);

	for (int j = 0; j < nBasFcts; j++) // #BasisFunctions
	  for (int k = 0; k < parts; k++)  // #edges (2d) or #faces (3d)
	    grd1[k] += quadFast->getGradient(iq, j, k) * localVec[j];

	  nullify(derivativeAtQPs[iq]);
	  for (int k = 0; k < parts; k++)
	    derivativeAtQPs[iq] += grdLambda[k][comp] * grd1[k];
      }

    } else {
      mtl::dense_vector<double> grdPhi(dim + 1);

      for (int iq = 0; iq < nPoints; iq++) {
	nullify(grd1);

	for (int j = 0; j < nBasFcts; j++) {
	  (*(basFcts->getGrdPhi(j)))(quad->getLambda(iq), grdPhi);
	  for (int k = 0; k < parts; k++)
	    grd1[k] += grdPhi[k] * localVec[j];
	}

	  nullify(derivativeAtQPs[iq]);
	  for (int k = 0; k < parts; k++)
	    derivativeAtQPs[iq] += grdLambda[k][comp] * grd1[k];
      }
    }
  }


  template<typename T>
  void DOFVectorBase<T>::getDerivativeAtQPs(const ElInfo *smallElInfo,
					    const ElInfo *largeElInfo,
					    const Quadrature *quad,
					    const FastQuadrature *quadFast,
					    int comp,
					    mtl::dense_vector<T> &derivativeAtQPs) const
  {
    FUNCNAME_DBG("DOFVector<T>::getGrdAtQPs()");

    TEST_EXIT_DBG(quad || quadFast)("neither quad nor quadFast defined\n");
    TEST_EXIT_DBG(!(quad && quadFast) || quad == quadFast->getQuadrature())
      ("quad != quadFast->quadrature\n");
    TEST_EXIT_DBG(!quadFast || quadFast->getBasisFunctions() == feSpace->getBasisFcts())
      ("invalid basis functions");

    if (smallElInfo->getMesh() == feSpace->getMesh())
      return getDerivativeAtQPs(smallElInfo, quad, quadFast, comp, derivativeAtQPs);

    const BasisFunction *basFcts = feSpace->getBasisFcts();
    int nBasFcts  = basFcts->getNumber();
    int nPoints = quadFast ? quadFast->getQuadrature()->getNumPoints() : quad->getNumPoints();

    mtl::dense_vector<T> localVec(nBasFcts);
    this->getLocalVector(largeElInfo->getElement(), localVec);

    mtl::dense2D<double> &m = smallElInfo->getSubElemCoordsMat(basFcts->getDegree());

    int parts = Global::getGeo(PARTS, dim);
    const DimVec<WorldVector<double> > &grdLambda = largeElInfo->getGrdLambda();

    mtl::dense_vector<T> grd1(dim + 1);
    mtl::dense_vector<double> grdPhi(dim + 1);
    mtl::dense_vector<double> tmp(dim + 1);

    derivativeAtQPs.change_dim(nPoints);
    for (int iq = 0; iq < nPoints; iq++) {
      nullify(grd1);

      for (int j = 0; j < nBasFcts; j++) {
	grdPhi = 0.0;

	for (int k = 0; k < nBasFcts; k++) {
	  (*(basFcts->getGrdPhi(k)))(quad->getLambda(iq), tmp);
	  tmp *= m[j][k];
	  grdPhi += tmp;
	}

	for (int k = 0; k < parts; k++)
	  grd1[k] += grdPhi[k] * localVec[j];
      }

      nullify(derivativeAtQPs[iq]);
      for (int k = 0; k < parts; k++)
	derivativeAtQPs[iq] += grdLambda[k][comp] * grd1[k];
    }
  }

  template <typename T>
  inline void set_to_zero(AMDiS::DOFVector<T>& v)
  {
    using math::zero;
    T  ref, my_zero(zero(ref));

    std::fill(v.begin(), v.end(), my_zero);
  }


  template<typename T>
  double DOFVector<T>::DoubleWell(Quadrature* q) const
  {
    Mesh* mesh = this->feSpace->getMesh();

    if (!q) {
      int deg = 2 * this->feSpace->getBasisFcts()->getDegree();
      q = Quadrature::provideQuadrature(this->dim, deg);
    }

    FastQuadrature *quadFast =
      FastQuadrature::provideFastQuadrature(this->feSpace->getBasisFcts(), *q, INIT_PHI);

    double result = 0.0;
    int nPoints = quadFast->getNumPoints();
    mtl::dense_vector<T> uh_vec(nPoints);
    TraverseStack stack;
    ElInfo *elInfo =
      stack.traverseFirst(mesh, -1,
			  Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET);
    while (elInfo) {
      double det = elInfo->getDet();
      double normT = 0.0;
      this->getVecAtQPs(elInfo, NULL, quadFast, uh_vec);
      for (int iq = 0; iq < nPoints; iq++)
	normT += quadFast->getWeight(iq) * sqr(uh_vec[iq]) * sqr(1.0 - uh_vec[iq]);
      result += det * normT;

      elInfo = stack.traverseNext(elInfo);
    }

#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localResult = result;
    MPI::COMM_WORLD.Allreduce(&localResult, &result, 1, MPI_DOUBLE, MPI_SUM);
#endif

    return result;
  }


  template<typename T>
  DOFVector<typename GradientType<T>::type>*
  DOFVector<T>::getGradient(DOFVector<typename GradientType<T>::type> *grad) const
  {
    FUNCNAME_DBG("DOFVector<T>::getGradient()");
    const FiniteElemSpace *feSpace = DOFVector<T>::feSpace;

    // define result vector
    static DOFVector<typename GradientType<T>::type> *result = NULL; // TODO: REMOVE STATIC

    if (grad) {
      result = grad;
    } else {
      if(result && result->getFeSpace() != feSpace) {
	delete result;
	result = new DOFVector<typename GradientType<T>::type>(feSpace, "gradient");
      }
    }

    Mesh *mesh = feSpace->getMesh();
    int dim = mesh->getDim();
    const BasisFunction *basFcts = feSpace->getBasisFcts();
    DOFAdmin *admin = feSpace->getAdmin();

    // count number of nodes and dofs per node
    std::vector<int> nNodeDOFs;
    std::vector<int> nNodePreDofs;
    std::vector<DimVec<double>*> bary;

    int nNodes = 0;
    int nDofs = 0;

    for (int i = 0; i < dim + 1; i++) {
      GeoIndex geoIndex = INDEX_OF_DIM(i, dim);
      int nPositions = mesh->getGeo(geoIndex);
      int numPreDofs = admin->getNumberOfPreDofs(i);
      for (int j = 0; j < nPositions; j++) {
	int dofs = basFcts->getNumberOfDofs(geoIndex);
	nNodeDOFs.push_back(dofs);
	nDofs += dofs;
	nNodePreDofs.push_back(numPreDofs);
      }
      nNodes += nPositions;
    }

    TEST_EXIT_DBG(nDofs == basFcts->getNumber())
      ("number of dofs != number of basis functions\n");

    for (int i = 0; i < nDofs; i++)
      bary.push_back(basFcts->getCoords(i));

    mtl::dense_vector<T> localUh(basFcts->getNumber());

    // traverse mesh
    std::vector<bool> visited(getUsedSize(), false);
    TraverseStack stack;
    Flag fillFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_GRD_LAMBDA | Mesh::FILL_COORDS;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, fillFlag);

    while (elInfo) {
      const DegreeOfFreedom **dof = elInfo->getElement()->getDof();
      const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
      this->getLocalVector(elInfo->getElement(), localUh);

      int localDOFNr = 0;
      for (int i = 0; i < nNodes; i++) { // for all nodes
	for (int j = 0; j < nNodeDOFs[i]; j++) { // for all dofs at this node
	  DegreeOfFreedom dofIndex = dof[i][nNodePreDofs[i] + j];
	  if (!visited[dofIndex]) {
	    basFcts->evalGrdUh(*(bary[localDOFNr]), grdLambda,
			       localUh, ((*result)[dofIndex]));

	    visited[dofIndex] = true;
	  }
	  localDOFNr++;
	}
      }

      elInfo = stack.traverseNext(elInfo);
    }

    return result;
  }


  template<typename T>
  DOFVector<typename GradientType<T>::type>*
  DOFVector<T>::getRecoveryGradient(DOFVector<typename GradientType<T>::type> *grad) const
  {
    const FiniteElemSpace *feSpace = DOFVector<T>::feSpace;
    int dim = DOFVector<T>::dim;

    // define result vector
    static DOFVector<typename GradientType<T>::type> *vec = NULL; // TODO: REMOVE STATIC

    DOFVector<typename GradientType<T>::type> *result = grad;

    if (!result) {
      if (vec && vec->getFeSpace() != feSpace) {
	delete vec;
	vec = NULL;
      }
      if (!vec)
	vec = new DOFVector<typename GradientType<T>::type>(feSpace, "gradient");

      result = vec;
    }

    typename GradientType<T>::type grd;
    nullify(grd);

    result->set(grd);

    DOFVector<double> volume(feSpace, "volume");
    volume.set(0.0);

    const BasisFunction *basFcts = feSpace->getBasisFcts();
    int nBasisFcts = basFcts->getNumber();
    DimVec<double> bary(dim, DEFAULT_VALUE, (1.0 / (dim + 1.0)));

    // traverse mesh
    Mesh *mesh = feSpace->getMesh();
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1,
					 Mesh::CALL_LEAF_EL | Mesh::FILL_DET |
					 Mesh::FILL_GRD_LAMBDA | Mesh::FILL_COORDS);

    mtl::dense_vector<T> localUh(nBasisFcts);
    std::vector<DegreeOfFreedom> localIndices(nBasisFcts);

    while (elInfo) {
      double det = elInfo->getDet();
      const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
      this->getLocalVector(elInfo->getElement(), localUh);
      basFcts->evalGrdUh(bary, grdLambda, localUh, grd);
      basFcts->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(), localIndices);

      for (int i = 0; i < nBasisFcts; i++) {
	(*result)[localIndices[i]] += grd * det;
	volume[localIndices[i]] += det;
      }

      elInfo = stack.traverseNext(elInfo);

    }

    // NOTE: We have to synchronize the vectors in PARALLEL_DOMAIN_AMDIS mode
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    Parallel::MeshDistributor::globalMeshDistributor->checkMeshChange(false);
    Parallel::MeshDistributor::globalMeshDistributor->synchAddVector(*result);
    Parallel::MeshDistributor::globalMeshDistributor->synchAddVector(volume);
#endif

    DOFVector<double>::Iterator volIt(&volume, USED_DOFS);
    typename DOFVector<typename GradientType<T>::type>::Iterator grdIt(result, USED_DOFS);

    for (volIt.reset(), grdIt.reset(); !volIt.end(); ++volIt, ++grdIt)
      if (*volIt != 0.0)
	*grdIt *= 1.0 / (*volIt);

    return result;
  }


  template<typename T>
  std::vector<DOFVector<double>*> *transform(DOFVector<typename GradientType<T>::type> *vec,
					     std::vector<DOFVector<double>*> *res)
  {
    FUNCNAME_DBG("DOFVector<T>::transform()");

    TEST_EXIT_DBG(vec)("no vector\n");

    static std::vector<DOFVector<double>*> *result = NULL; // TODO: REMOVE STATIC

    int len = num_rows(GradientType<T>::getValues((*vec)[0]));

    if (!res && !result) {
      result = new std::vector<DOFVector<double>*>(len);
      for (int i = 0; i < len; i++)
	(*result)[i] = new DOFVector<double>(vec->getFeSpace(), "transform");
    } else if (res->size() == 0 || (*res)[0] == NULL) {
      res->resize(len, NULL);
      for (int i = 0; i < len; i++)
	(*res)[i] = new DOFVector<double>(vec->getFeSpace(), "transform");
    }

    std::vector<DOFVector<double>*> *r = res ? res : result;

    int vecSize = vec->getSize();
    for (int i = 0; i < vecSize; i++)
      for (int j = 0; j < len; j++)
	(*((*r)[j]))[i] = (GradientType<T>::getValues((*vec)[i]))[j];

    return r;
  }
}