#include "Quadrature.h"
#include "FiniteElemSpace.h"
#include "ElInfo.h"

namespace AMDiS {

  template<typename T>
  void SubAssembler::getVectorAtQPs(DOFVectorBase<T>* vec, 
				    const ElInfo* elInfo,
				    Quadrature *quad,
				    mtl::dense_vector<T>& vecAtQPs)
  {
    FUNCNAME("SubAssembler::getVectorAtQPs()");

    TEST_EXIT_DBG(vec)("No DOF vector!\n");
    TEST_EXIT_DBG(elInfo->getMesh() == vec->getFeSpace()->getMesh())
      ("Vector and FE space do not fit together!\n");

    Quadrature *localQuad = quad ? quad : quadrature;

    if (valuesAtQPs[vec] && valuesAtQPs[vec]->valid) {
      vecAtQPs = valuesAtQPs[vec]->values;
      return;
    }

    if (!valuesAtQPs[vec])
      valuesAtQPs[vec] = new ValuesAtQPs;
    valuesAtQPs[vec]->values.change_dim(localQuad->getNumPoints());
    vecAtQPs.change_dim(localQuad->getNumPoints());

    mtl::dense_vector<T>& values = valuesAtQPs[vec]->values;
    bool sameFeSpaces = 
      vec->getFeSpace() == rowFeSpace || vec->getFeSpace() == colFeSpace;

    if (opt && !quad && sameFeSpaces) {
      const BasisFunction *psi = rowFeSpace->getBasisFcts();
      const BasisFunction *phi = colFeSpace->getBasisFcts();
      if (vec->getFeSpace()->getBasisFcts() == psi)
	psiFast = updateFastQuadrature(psiFast, psi, INIT_PHI);
      else if(vec->getFeSpace()->getBasisFcts() == phi)
	phiFast = updateFastQuadrature(phiFast, phi, INIT_PHI);      
    }

    // calculate new values
    const BasisFunction *basFcts = vec->getFeSpace()->getBasisFcts();

    if (opt && !quad && sameFeSpaces) {
      if (psiFast->getBasisFunctions() == basFcts) {
	vec->getVecAtQPs(elInfo, NULL, psiFast, values);
      } else if (phiFast->getBasisFunctions() == basFcts) {
	vec->getVecAtQPs(elInfo, NULL, phiFast, values);
      } else {
	vec->getVecAtQPs(elInfo, localQuad, NULL, values);
      }
    } else {
      vec->getVecAtQPs(elInfo, localQuad, NULL, values);
    }

    valuesAtQPs[vec]->valid = true;
    vecAtQPs = values;
  }


  template<typename T>
  void SubAssembler::getVectorAtQPs(DOFVectorBase<T>* vec, 
				    const ElInfo* smallElInfo,
				    const ElInfo* largeElInfo,
				    Quadrature *quad,
				    mtl::dense_vector<T>& vecAtQPs)
  {
    FUNCNAME("SubAssembler::getVectorAtQPs()");

    TEST_EXIT_DBG(vec)("No DOF vector!\n");

    Quadrature *localQuad = quad ? quad : quadrature;

    if (!valuesAtQPs[vec])
      valuesAtQPs[vec] = new ValuesAtQPs;

    valuesAtQPs[vec]->values.change_dim(localQuad->getNumPoints());
    mtl::dense_vector<T>& values = valuesAtQPs[vec]->values;
    valuesAtQPs[vec]->valid = true;
    vec->getVecAtQPs(smallElInfo, largeElInfo, localQuad, NULL, values);
    vecAtQPs = values;
  }

}