Assembler.cc

#include "Assembler.h"
#include "Operator.h"
#include "Element.h"
#include "QPsiPhi.h"
#include "ElementMatrix.h"
#include "ElementVector.h"
#include "DOFVector.h"
#include "OpenMP.h"
#include <vector>
#include <algorithm>

namespace AMDiS {

  std::vector<SubAssembler*> ZeroOrderAssembler::optimizedSubAssemblers;
  std::vector<SubAssembler*> FirstOrderAssembler::optimizedSubAssemblersGrdPhi;
  std::vector<SubAssembler*> FirstOrderAssembler::optimizedSubAssemblersGrdPsi;
  std::vector<SubAssembler*> SecondOrderAssembler::optimizedSubAssemblers;
  
  std::vector<SubAssembler*> ZeroOrderAssembler::standardSubAssemblers;
  std::vector<SubAssembler*> FirstOrderAssembler::standardSubAssemblersGrdPhi;
  std::vector<SubAssembler*> FirstOrderAssembler::standardSubAssemblersGrdPsi;
  std::vector<SubAssembler*> SecondOrderAssembler::standardSubAssemblers;

  SubAssembler::SubAssembler(Operator *op,
			     Assembler *assembler,
			     Quadrature *quadrat,
			     int order, 
			     bool optimized,
			     FirstOrderType type) 
    : nRow(0),
      nCol(0),
      coordsAtQPs(NULL),
      coordsNumAllocated(0),
      quadrature(quadrat),
      psiFast(NULL),
      phiFast(NULL),
      owner(assembler),
      symmetric(true),
      opt(optimized),
      firstCall(true)
  {
    const BasisFunction *psi = assembler->rowFESpace->getBasisFcts();
    const BasisFunction *phi = assembler->colFESpace->getBasisFcts();

    nRow = psi->getNumber();
    nCol = phi->getNumber();

    int maxThreads = omp_get_max_threads();
    terms.resize(maxThreads);

    switch (order) {
    case 0:
      terms = op->zeroOrder;
      break;
    case 1:
      if(type == GRD_PHI)
	terms = op->firstOrderGrdPhi;
      else 
	terms = op->firstOrderGrdPsi;
      break;
    case 2:
      terms = op->secondOrder;
      break;
    }

    // check if all terms are symmetric
    symmetric = true;
    for (int i = 0; i < static_cast<int>(terms[0].size()); i++) {
      if (!(terms[0][i])->isSymmetric()) {
	symmetric = false;
	break;
      }
    }  

    dim = assembler->rowFESpace->getMesh()->getDim();
  }

  FastQuadrature *SubAssembler::updateFastQuadrature(FastQuadrature *quadFast,
						     const BasisFunction *psi,
						     Flag updateFlag)
  {
    if (!quadFast) {
      quadFast = FastQuadrature::provideFastQuadrature(psi,
						       *quadrature,
						       updateFlag);
    } else {
      if (!quadFast->initialized(updateFlag))
	quadFast->init(updateFlag);
    }

    return quadFast;
  }

  void SubAssembler::initElement(const ElInfo* elInfo,
				 Quadrature *quad)
  {
    // set corrdsAtQPs invalid
    coordsValid = false;

    // set values at QPs invalid
    std::map<const DOFVectorBase<double>*, ValuesAtQPs*>::iterator it1;
    for (it1 = valuesAtQPs.begin(); it1 != valuesAtQPs.end(); ++it1) {
      ((*it1).second)->valid = false;
    }

    // set gradients at QPs invalid
    std::map<const DOFVectorBase<double>*, GradientsAtQPs*>::iterator it2;
    for (it2 = gradientsAtQPs.begin(); it2 != gradientsAtQPs.end(); ++it2) {
      ((*it2).second)->valid = false;
    }

    int myRank = omp_get_thread_num();
    // calls initElement of each term
    std::vector<OperatorTerm*>::iterator it;
    for (it = terms[myRank].begin(); it != terms[myRank].end(); ++it) {
      (*it)->initElement(elInfo, this, quad);
    }
  }

  WorldVector<double>* SubAssembler::getCoordsAtQPs(const ElInfo* elInfo,
						    Quadrature *quad)
  {
    Quadrature *localQuad = quad ? quad : quadrature;
  
    const int nPoints = localQuad->getNumPoints();

    // already calculated for this element ?
    if (coordsValid) {
      return coordsAtQPs;
    }
   
    if (coordsAtQPs)  {
      if (coordsNumAllocated != nPoints) {
	DELETE [] coordsAtQPs;
        coordsAtQPs = NEW WorldVector<double>[nPoints];
	coordsNumAllocated = nPoints;
      }
    } else {
      coordsAtQPs = NEW WorldVector<double>[nPoints];
      coordsNumAllocated = nPoints;
    }

    // set new values
    WorldVector<double>* k = &(coordsAtQPs[0]);
    for (int l = 0; k < &(coordsAtQPs[nPoints]); ++k, ++l) {
      elInfo->coordToWorld(localQuad->getLambda(l), k);
    }

    // mark values as valid
    coordsValid = true;

    return coordsAtQPs;
  }

  double* SubAssembler::getVectorAtQPs(DOFVectorBase<double>* dv, 
				       const ElInfo* elInfo,
				       Quadrature *quad)
  {
    FUNCNAME("SubAssembler::getVectorAtQPs()");

    const DOFVectorBase<double>* vec = dv ? dv : owner->operat->getUhOld();

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

    if (valuesAtQPs[vec] && valuesAtQPs[vec]->valid) 
      return valuesAtQPs[vec]->values.getValArray();

    Quadrature *localQuad = quad ? quad : quadrature;

    if (!valuesAtQPs[vec]) {
      valuesAtQPs[vec] = new ValuesAtQPs;
    }
    valuesAtQPs[vec]->values.resize(localQuad->getNumPoints());

    double *values = valuesAtQPs[vec]->values.getValArray();

    bool sameFESpaces = 
      (vec->getFESpace() == owner->rowFESpace) || 
      (vec->getFESpace() == owner->colFESpace);

    if (opt && !quad && sameFESpaces) {
      const BasisFunction *psi = owner->rowFESpace->getBasisFcts();
      const BasisFunction *phi = owner->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;

    return values;
  }

  WorldVector<double>* SubAssembler::getGradientsAtQPs(DOFVectorBase<double>* dv, 
						       const ElInfo* elInfo,
						       Quadrature *quad)
  {
    FUNCNAME("SubAssembler::getGradientsAtQPs()");

    const DOFVectorBase<double>* vec = dv ? dv : owner->operat->getUhOld();

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

    if (gradientsAtQPs[vec] && gradientsAtQPs[vec]->valid) 
      return gradientsAtQPs[vec]->values.getValArray();

    Quadrature *localQuad = quad ? quad : quadrature;

    if (!gradientsAtQPs[vec]) {
      gradientsAtQPs[vec] = new GradientsAtQPs;
    } 
    gradientsAtQPs[vec]->values.resize(localQuad->getNumPoints());

    WorldVector<double> *values = gradientsAtQPs[vec]->values.getValArray();

    const BasisFunction *psi = owner->rowFESpace->getBasisFcts();
    const BasisFunction *phi = owner->colFESpace->getBasisFcts();

    bool sameFESpaces = 
      (vec->getFESpace() == owner->rowFESpace) || 
      (vec->getFESpace() == owner->colFESpace);

    if (opt && !quad && sameFESpaces) {
      if (vec->getFESpace()->getBasisFcts() == psi) {
	psiFast = updateFastQuadrature(psiFast, psi, INIT_GRD_PHI);
      } else if(vec->getFESpace()->getBasisFcts() == phi) {
	phiFast = updateFastQuadrature(phiFast, phi, INIT_GRD_PHI);
      }
    }
  
    // calculate new values
    const BasisFunction *basFcts = vec->getFESpace()->getBasisFcts();

    if (opt && !quad && sameFESpaces) {
      if (psiFast->getBasisFunctions() == basFcts) {
	vec->getGrdAtQPs(elInfo, NULL, psiFast, values);
      } else {
	vec->getGrdAtQPs(elInfo, NULL, phiFast, values);
      }
    } else {
      vec->getGrdAtQPs(elInfo, localQuad, NULL, values);
    }
   
    gradientsAtQPs[vec]->valid = true;

    return values;
  }

  ZeroOrderAssembler::ZeroOrderAssembler(Operator *op,
					 Assembler *assembler,
					 Quadrature *quad,
					 bool optimized)
    : SubAssembler(op, assembler, quad, 0, optimized)
  {}

  FirstOrderAssembler::FirstOrderAssembler(Operator *op,
					   Assembler *assembler,
					   Quadrature *quad,
					   bool optimized,
					   FirstOrderType type)
    : SubAssembler(op, assembler, quad, 1, optimized, type)
  {}

  SecondOrderAssembler::SecondOrderAssembler(Operator *op,
					     Assembler *assembler,
					     Quadrature *quad,
					     bool optimized)
    : SubAssembler(op, assembler, quad, 2, optimized)
  {}

  ZeroOrderAssembler* 
  ZeroOrderAssembler::getSubAssembler(Operator* op,
				      Assembler *assembler,
				      Quadrature *quad,
				      bool optimized)
  {
    // check if a assembler is needed at all
    if (op->zeroOrder.size() == 0) {
      return NULL;
    }

    ZeroOrderAssembler *newAssembler;

    std::vector<SubAssembler*> *subAssemblers =
	optimized ?
	&optimizedSubAssemblers :
    &standardSubAssemblers;

    int myRank = omp_get_thread_num();
    std::vector<OperatorTerm*> opTerms  = op->zeroOrder[myRank];

    sort(opTerms.begin(), opTerms.end());

    // check if a new assembler is needed
    if (quad) {
      for (int i = 0; i < static_cast<int>( subAssemblers->size()); i++) {
	std::vector<OperatorTerm*> assTerms = *((*subAssemblers)[i]->getTerms());

	sort(assTerms.begin(), assTerms.end());

	if ((opTerms == assTerms) && 
	    ((*subAssemblers)[i]->getQuadrature() == quad)) {
	
	  return dynamic_cast<ZeroOrderAssembler*>((*subAssemblers)[i]);
	}
      }
    }
 
    // check if all terms are pw_const
    bool pwConst = true;
    for (int i = 0; i < static_cast<int>( op->zeroOrder[myRank].size()); i++) {
      if (!op->zeroOrder[myRank][i]->isPWConst()) {
	pwConst = false;
	break;
      }
    }  

    // create new assembler
    if (!optimized) {
      newAssembler = NEW Stand0(op, assembler, quad);
    } else {
      if (pwConst) {
	newAssembler = NEW Pre0(op, assembler, quad);
      } else {
	newAssembler = NEW Quad0(op, assembler, quad);
      }
    }

    subAssemblers->push_back(newAssembler);
    return newAssembler;
  }

  FirstOrderAssembler* 
  FirstOrderAssembler::getSubAssembler(Operator* op,
				       Assembler *assembler,
				       Quadrature *quad,
				       FirstOrderType type,
				       bool optimized)
  {
    std::vector<SubAssembler*> *subAssemblers =
	optimized ?
	(type == GRD_PSI ? 
	 &optimizedSubAssemblersGrdPsi : 
	 &optimizedSubAssemblersGrdPhi) :
    (type == GRD_PSI ? 
     &standardSubAssemblersGrdPsi :
     &standardSubAssemblersGrdPhi);

    int myRank = omp_get_thread_num();
    std::vector<OperatorTerm*> opTerms 
	= (type == GRD_PSI) ? 
	    op->firstOrderGrdPsi[myRank] : 
            op->firstOrderGrdPhi[myRank];

    // check if a assembler is needed at all
    if (opTerms.size() == 0) {
      return NULL;
    }

    sort(opTerms.begin(), opTerms.end());

    FirstOrderAssembler *newAssembler;

    // check if a new assembler is needed
    for (int i = 0; i < static_cast<int>( subAssemblers->size()); i++) {

      std::vector<OperatorTerm*> assTerms = *((*subAssemblers)[i]->getTerms());
    
      sort(assTerms.begin(), assTerms.end());   

      if ((opTerms == assTerms) && 
	  ((*subAssemblers)[i]->getQuadrature() == quad)) {

	return dynamic_cast<FirstOrderAssembler*>((*subAssemblers)[i]);
      }
    }

    // check if all terms are pw_const
    bool pwConst = true;
    for (int i = 0; i < static_cast<int>( opTerms.size()); i++) {
      if (!(opTerms[i])->isPWConst()) {
	pwConst = false;
	break;
      }
    }  

    // create new assembler
    if (!optimized) {
      newAssembler = 
	(type == GRD_PSI) ?
	dynamic_cast<FirstOrderAssembler*>(NEW Stand10(op, assembler, quad)) :
	dynamic_cast<FirstOrderAssembler*>(NEW Stand01(op, assembler, quad));    
    } else {
      if (pwConst) {
	newAssembler = 
	  (type == GRD_PSI) ?
	  dynamic_cast<FirstOrderAssembler*>(NEW Pre10(op, assembler, quad)) :
	  dynamic_cast<FirstOrderAssembler*>(NEW Pre01(op, assembler, quad));
      } else {
	newAssembler = 
	  (type == GRD_PSI) ?
	  dynamic_cast<FirstOrderAssembler*>( NEW Quad10(op, assembler, quad)) :
	  dynamic_cast<FirstOrderAssembler*>( NEW Quad01(op, assembler, quad));
      }
    }

    subAssemblers->push_back(newAssembler);
    return newAssembler;
  };

  SecondOrderAssembler* 
  SecondOrderAssembler::getSubAssembler(Operator* op,
					Assembler *assembler,
					Quadrature *quad,
					bool optimized)
  {
    int myRank = omp_get_thread_num();

    // check if a assembler is needed at all
    if (op->secondOrder[myRank].size() == 0) {
      return NULL;
    }

    SecondOrderAssembler *newAssembler;

    std::vector<SubAssembler*> *subAssemblers =
	optimized ?
	&optimizedSubAssemblers :
    &standardSubAssemblers;

    std::vector<OperatorTerm*> opTerms  = op->zeroOrder[myRank];

    sort(opTerms.begin(), opTerms.end());

    // check if a new assembler is needed
    for (int i = 0; i < static_cast<int>( subAssemblers->size()); i++) {
      std::vector<OperatorTerm*> assTerms = *((*subAssemblers)[i]->getTerms());
    
      sort(assTerms.begin(), assTerms.end());

      if ((opTerms == assTerms) && 
	  ((*subAssemblers)[i]->getQuadrature() == quad)) {
	
	return dynamic_cast<SecondOrderAssembler*>((*subAssemblers)[i]);
      }
    }

    // check if all terms are pw_const
    bool pwConst = true;
    for (int i = 0; i < static_cast<int>( op->secondOrder[myRank].size()); i++) {
      if (!op->secondOrder[myRank][i]->isPWConst()) {
	pwConst = false;
	break;
      }
    }  

    // create new assembler
    if (!optimized) {
      newAssembler = NEW Stand2(op, assembler, quad);
    } else {
      if (pwConst) {
	newAssembler = NEW Pre2(op, assembler, quad);
      } else {
	newAssembler = NEW Quad2(op, assembler, quad);
      }
    }

    subAssemblers->push_back(newAssembler);
    return newAssembler;
  }

  Stand0::Stand0(Operator *op, Assembler *assembler, Quadrature *quad)
    : ZeroOrderAssembler(op, assembler, quad, false)
  {
  }

  void Stand0::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
  {
    const BasisFunction *psi = owner->getRowFESpace()->getBasisFcts();
    const BasisFunction *phi = owner->getColFESpace()->getBasisFcts();

    double *phival = GET_MEMORY(double, nCol);
    int nPoints = quadrature->getNumPoints();
    double *c = GET_MEMORY(double, nPoints);

    for (int iq = 0; iq < nPoints; iq++) {
      c[iq] = 0.0;
    }

    int myRank = omp_get_thread_num();
    std::vector<OperatorTerm*>::iterator termIt;
    for (termIt = terms[myRank].begin(); termIt != terms[myRank].end(); ++termIt) {
      (static_cast<ZeroOrderTerm*>((*termIt)))->getC(elInfo, nPoints, c);
    }
      
    if (symmetric) {
      for (int iq = 0; iq < nPoints; iq++) {
	c[iq] *= elInfo->getDet();

	// calculate phi at QPs only once!
	for (int i = 0; i < nCol; i++) {
	  phival[i] = (*(phi->getPhi(i)))(quadrature->getLambda(iq));
	}

	for (int i = 0; i < nRow; i++) {
	  double psival = (*(psi->getPhi(i)))(quadrature->getLambda(iq));
	  (*mat)[i][i] += quadrature->getWeight(iq) * c[iq] * psival * phival[i];
	  for (int j = i + 1; j < nCol; j++) {
	    double val = quadrature->getWeight(iq) * c[iq] * psival * phival[j];
	    (*mat)[i][j] += val;
	    (*mat)[j][i] += val;
	  }
	}
      }
    } else {      //  non symmetric assembling 
      for (int iq = 0; iq < nPoints; iq++) {
	c[iq] *= elInfo->getDet();

	// calculate phi at QPs only once!
	for (int i = 0; i < nCol; i++) {
	  phival[i] = (*(phi->getPhi(i)))(quadrature->getLambda(iq));
	}

	for (int i = 0; i < nRow; i++) {
	  double psival = (*(psi->getPhi(i)))(quadrature->getLambda(iq));
	  for (int j = 0; j < nCol; j++) {
	    (*mat)[i][j] += quadrature->getWeight(iq) * c[iq] * psival * phival[j];
	  }
	}
      }
    }

    FREE_MEMORY(phival, double, nCol);
    FREE_MEMORY(c, double, nPoints);
  }

  void Stand0::calculateElementVector(const ElInfo *elInfo, ElementVector *vec)
  {
    int nPoints = quadrature->getNumPoints();

    double *c = GET_MEMORY(double, nPoints);

    for (int iq = 0; iq < nPoints; iq++) {
      c[iq] = 0.0;
    }

    int myRank = omp_get_thread_num();
    std::vector<OperatorTerm*>::iterator termIt;
    for (termIt = terms[myRank].begin(); termIt != terms[myRank].end(); ++termIt) {
      (static_cast<ZeroOrderTerm*>((*termIt)))->getC(elInfo, nPoints, c);
    }

    for (int iq = 0; iq < nPoints; iq++) {
      c[iq] *= elInfo->getDet();

      for (int i = 0; i < nRow; i++) {
	double psi = (*(owner->getRowFESpace()->getBasisFcts()->getPhi(i)))
	  (quadrature->getLambda(iq));
	(*vec)[i] += quadrature->getWeight(iq) * c[iq] * psi;
      }
    }
    
    FREE_MEMORY(c, double, nPoints);
  }

  Quad0::Quad0(Operator *op, Assembler *assembler, Quadrature *quad)
    : ZeroOrderAssembler(op, assembler, quad, true)
  {
    cPtrs.resize(omp_get_max_threads());
  }

  Quad0::~Quad0()
  {
    for (int i = 0; i < omp_get_max_threads(); i++) {
      FREE_MEMORY(cPtrs[i], double, quadrature->getNumPoints());
    }
  }

  void Quad0::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
  {
    int nPoints = quadrature->getNumPoints();
    int myRank = omp_get_thread_num();

    if (firstCall) {
      cPtrs[myRank] = GET_MEMORY(double, nPoints);
      const BasisFunction *basFcts = owner->getRowFESpace()->getBasisFcts();
      psiFast = updateFastQuadrature(psiFast, basFcts, INIT_PHI);
      basFcts = owner->getColFESpace()->getBasisFcts();
      phiFast = updateFastQuadrature(phiFast, basFcts, INIT_PHI);
      firstCall = false;
    }

    double *c = cPtrs[myRank];
    for (int iq = 0; iq < nPoints; iq++) {
      c[iq] = 0.0;
    }

    std::vector<OperatorTerm*>::iterator termIt;
    for (termIt = terms[myRank].begin(); termIt != terms[myRank].end(); ++termIt) {
      (static_cast<ZeroOrderTerm*>((*termIt)))->getC(elInfo, nPoints, c);
    }

    if (symmetric) {
      for (int iq = 0; iq < nPoints; iq++) {
	c[iq] *= elInfo->getDet();

	const double *psi = psiFast->getPhi(iq);
	const double *phi = phiFast->getPhi(iq);
	for (int i = 0; i < nRow; i++) {
	  (*mat)[i][i] += quadrature->getWeight(iq) * c[iq] * psi[i] * phi[i];
	  for (int j = i + 1; j < nCol; j++) {
	    double val = quadrature->getWeight(iq) * c[iq] * psi[i] * phi[j];
	    (*mat)[i][j] += val;
	    (*mat)[j][i] += val;
	  }
	}
      }
    } else {      /*  non symmetric assembling   */
      for (int iq = 0; iq < nPoints; iq++) {
	c[iq] *= elInfo->getDet();

	const double *psi = psiFast->getPhi(iq);
	const double *phi = phiFast->getPhi(iq);
	for (int i = 0; i < nRow; i++) {
	  for (int j = 0; j < nCol; j++) {
	    (*mat)[i][j] += quadrature->getWeight(iq) * c[iq] * psi[i] * phi[j];
	  }
	}
      }
    }
  }

  void Quad0::calculateElementVector(const ElInfo *elInfo, ElementVector *vec)
  {
    if (firstCall) {
      const BasisFunction *basFcts = owner->getRowFESpace()->getBasisFcts();
      psiFast = updateFastQuadrature(psiFast, basFcts, INIT_PHI);
      basFcts = owner->getColFESpace()->getBasisFcts();
      phiFast = updateFastQuadrature(phiFast, basFcts, INIT_PHI);
      firstCall = false;
    }

    int nPoints = quadrature->getNumPoints();
    double *c = GET_MEMORY(double, nPoints);

    for (int iq = 0; iq < nPoints; iq++) {
      c[iq] = 0.0;
    }

    int myRank = omp_get_thread_num();
    std::vector<OperatorTerm*>::iterator termIt;
    for (termIt = terms[myRank].begin(); termIt != terms[myRank].end(); ++termIt) {
      (static_cast<ZeroOrderTerm*>((*termIt)))->getC(elInfo, nPoints, c);
    }

    for (int iq = 0; iq < nPoints; iq++) {
      c[iq] *= elInfo->getDet();

      const double *psi = psiFast->getPhi(iq);
      for (int i = 0; i < nRow; i++) {
	(*vec)[i] += quadrature->getWeight(iq) * c[iq] * psi[i];
      }
    }
    FREE_MEMORY(c, double, nPoints);
  }

  Pre0::Pre0(Operator *op, Assembler *assembler, Quadrature *quad) 
    : ZeroOrderAssembler(op, assembler, quad, true)
  {
  }

  void Pre0::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
  {
    if (firstCall) {
      q00 = Q00PsiPhi::provideQ00PsiPhi(owner->getRowFESpace()->getBasisFcts(), 
					owner->getColFESpace()->getBasisFcts(), 
					quadrature);
      q0 = Q0Psi::provideQ0Psi(owner->getRowFESpace()->getBasisFcts(),
			       quadrature);
      firstCall = false;
    }

    //    c[0] = 0.0;
    double c = 0.0;
    int myRank = omp_get_thread_num();
    int size = static_cast<int>(terms[myRank].size());

    for (int i = 0; i < size; i++) {
      (static_cast<ZeroOrderTerm*>((terms[myRank][i])))->getC(elInfo, 1, &c);
    }

    c *= elInfo->getDet();

    if (symmetric) {
      for (int i = 0; i < nRow; i++) {
	(*mat)[i][i] += c * q00->getValue(i,i);
	for (int j = i + 1; j < nCol; j++) {
	  double val = c * q00->getValue(i, j);
	  (*mat)[i][j] += val;
	  (*mat)[j][i] += val;
	}
      }
    } else {
      for (int i = 0; i < nRow; i++)
	for (int j = 0; j < nCol; j++)
	  (*mat)[i][j] += c * q00->getValue(i, j);
    }
  }

  void Pre0::calculateElementVector(const ElInfo *elInfo, ElementVector *vec)
  {
    if (firstCall) {
      q00 = Q00PsiPhi::provideQ00PsiPhi(owner->getRowFESpace()->getBasisFcts(), 
					owner->getColFESpace()->getBasisFcts(), 
					quadrature);
      q0 = Q0Psi::provideQ0Psi(owner->getRowFESpace()->getBasisFcts(),
			       quadrature);
      firstCall = false;
    }

    std::vector<OperatorTerm*>::iterator termIt;

    int myRank = omp_get_thread_num();
    double c = 0.0;
    for (termIt = terms[myRank].begin(); termIt != terms[myRank].end(); ++termIt) {
      (static_cast<ZeroOrderTerm*>( *termIt))->getC(elInfo, 1, &c);
    }

    c *= elInfo->getDet();

    for (int i = 0; i < nRow; i++)
      (*vec)[i] += c * q0->getValue(i);
  }

  Stand10::Stand10(Operator *op, Assembler *assembler, Quadrature *quad) 
    : FirstOrderAssembler(op, assembler, quad, false, GRD_PSI)
  {}


  void Stand10::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
  {
    DimVec<double> grdPsi(dim, DEFAULT_VALUE, 0.0);
    double *phival = GET_MEMORY(double, nCol);

    const BasisFunction *psi = owner->getRowFESpace()->getBasisFcts();
    const BasisFunction *phi = owner->getColFESpace()->getBasisFcts();

    int nPoints = quadrature->getNumPoints();

    VectorOfFixVecs<DimVec<double> > Lb(dim, nPoints, NO_INIT);
    for (int iq = 0; iq < nPoints; iq++) {
      Lb[iq].set(0.0);
    }

    int myRank = omp_get_thread_num();
    for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) {
      (static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, nPoints, Lb);
    }
  
    for (int iq = 0; iq < nPoints; iq++) {
      Lb[iq] *= elInfo->getDet();

      for (int i = 0; i < nCol; i++) {
	phival[i] = (*(phi->getPhi(i)))(quadrature->getLambda(iq));
      }

      for (int i = 0; i < nRow; i++) {
	(*(psi->getGrdPhi(i)))(quadrature->getLambda(iq), grdPsi);
	for (int j = 0; j < nCol; j++) {
	  (*mat)[i][j] += quadrature->getWeight(iq) * (Lb[iq] * grdPsi) * phival[j];
	}
      }
    }
    FREE_MEMORY(phival, double, nCol);
  }


  Quad10::Quad10(Operator *op, Assembler *assembler, Quadrature *quad) 
    : FirstOrderAssembler(op, assembler, quad, true, GRD_PSI)
  {
  }


  void Quad10::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
  {
    VectorOfFixVecs<DimVec<double> > *grdPsi;
    const double *phi;

    if (firstCall) {
      const BasisFunction *basFcts = owner->getRowFESpace()->getBasisFcts();
      psiFast = updateFastQuadrature(psiFast, basFcts, INIT_GRD_PHI);
      basFcts = owner->getColFESpace()->getBasisFcts();
      phiFast = updateFastQuadrature(phiFast, basFcts, INIT_PHI);
      firstCall = false;
    }

    int nPoints = quadrature->getNumPoints();

    VectorOfFixVecs<DimVec<double> > Lb(dim,nPoints,NO_INIT);
    for (int iq = 0; iq < nPoints; iq++) {
      Lb[iq].set(0.0);
    }

    int myRank = omp_get_thread_num();
    for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) {
      (static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, nPoints, Lb);
    }
  
    for (int iq = 0; iq < nPoints; iq++) {
      Lb[iq] *= elInfo->getDet();

      phi    = phiFast->getPhi(iq);
      grdPsi = psiFast->getGradient(iq);

      for (int i = 0; i < nRow; i++) {
	for (int j = 0; j < nCol; j++)
	  (*mat)[i][j] += quadrature->getWeight(iq) * (Lb[iq] * (*grdPsi)[i]) * phi[j];
      }
    }
  }


  Pre10::Pre10(Operator *op, Assembler *assembler, Quadrature *quad) 
    : FirstOrderAssembler(op, assembler, quad, true, GRD_PSI)
  {
  }


  void Pre10::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
  {
    VectorOfFixVecs<DimVec<double> > Lb(dim, 1, NO_INIT);
    const int *k;
    const double *values;

    if (firstCall) {
      q10 = Q10PsiPhi::provideQ10PsiPhi(owner->getRowFESpace()->getBasisFcts(), 
					owner->getColFESpace()->getBasisFcts(), 
					quadrature);
      q1 = Q1Psi::provideQ1Psi(owner->getRowFESpace()->getBasisFcts(),
			       quadrature);
      firstCall = false;
    }

    const int **nEntries = q10->getNumberEntries();
    int myRank = omp_get_thread_num();
    Lb[0].set(0.0);

    for (int i = 0; i < static_cast<int>( terms[myRank].size()); i++) {
      (static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, 1, Lb);
    }

    Lb[0] *= elInfo->getDet();

    for (int i = 0; i < nRow; i++) {
      for (int j = 0; j < nCol; j++) {
	k = q10->getKVec(i, j);
	values = q10->getValVec(i, j);
	double val = 0.0;
	for (int m = 0; m < nEntries[i][j]; m++) {
	  val += values[m] * Lb[0][k[m]];
	}
	(*mat)[i][j] += val;
      }
    }
  }


  Stand01::Stand01(Operator *op, Assembler *assembler, Quadrature *quad) 
    : FirstOrderAssembler(op, assembler, quad, false, GRD_PHI)
  {}


  void Stand01::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
  {
    VectorOfFixVecs<DimVec<double> > grdPhi(dim, nCol, NO_INIT);

    const BasisFunction *psi = owner->getRowFESpace()->getBasisFcts();
    const BasisFunction *phi = owner->getColFESpace()->getBasisFcts();

    int nPoints = quadrature->getNumPoints();
    VectorOfFixVecs<DimVec<double> > Lb(dim, nPoints, NO_INIT);
    int myRank = omp_get_thread_num();

    for (int iq = 0; iq < nPoints; iq++) {
      Lb[iq].set(0.0);
    }
    for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) {
      (static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, nPoints, Lb);
    }
  
    for (int iq = 0; iq < nPoints; iq++) {
      Lb[iq] *= elInfo->getDet();

      for (int i = 0; i < nCol; i++) {
	(*(phi->getGrdPhi(i)))(quadrature->getLambda(iq), grdPhi[i]);
      }

      for (int i = 0; i < nRow; i++) {
	double psival = (*(psi->getPhi(i)))(quadrature->getLambda(iq));
	for (int j = 0; j < nCol; j++)
	  (*mat)[i][j] += quadrature->getWeight(iq) * ((Lb[iq] * psival) * grdPhi[j]);
      }
    } 
  }

  void Stand10::calculateElementVector(const ElInfo *elInfo, ElementVector *vec)
  {
    DimVec<double> grdPsi(dim, DEFAULT_VALUE, 0.0);
    const BasisFunction *psi = owner->getRowFESpace()->getBasisFcts();
    int nPoints = quadrature->getNumPoints();
    VectorOfFixVecs<DimVec<double> > Lb(dim,nPoints,NO_INIT);
    int myRank = omp_get_thread_num();

    for (int iq = 0; iq < nPoints; iq++) {
      Lb[iq].set(0.0);
    }
    for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) {
      (static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, nPoints, Lb);
    }
  
    for (int iq = 0; iq < nPoints; iq++) {
      Lb[iq] *= elInfo->getDet();

      for (int i = 0; i < nRow; i++) {
	(*(psi->getGrdPhi(i)))(quadrature->getLambda(iq), grdPsi);
	(*vec)[i] += quadrature->getWeight(iq) * (Lb[iq] * grdPsi);
      }
    }
  }

  Quad01::Quad01(Operator *op, Assembler *assembler, Quadrature *quad) 
    : FirstOrderAssembler(op, assembler, quad, true, GRD_PHI)
  {
  }

  void Quad01::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
  {
    VectorOfFixVecs<DimVec<double> > *grdPhi;

    if (firstCall) {
      const BasisFunction *basFcts = owner->getRowFESpace()->getBasisFcts();
      psiFast = updateFastQuadrature(psiFast, basFcts, INIT_PHI);
      basFcts = owner->getColFESpace()->getBasisFcts();
      phiFast = updateFastQuadrature(phiFast, basFcts, INIT_GRD_PHI);
      firstCall = false;
    }

    int nPoints = quadrature->getNumPoints();
    VectorOfFixVecs<DimVec<double> > Lb(dim,nPoints,NO_INIT);
    int myRank = omp_get_thread_num();

    for (int iq = 0; iq < nPoints; iq++) {
      Lb[iq].set(0.0);
    }
    for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) {
      (static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, nPoints, Lb);
    }
  
    for (int iq = 0; iq < nPoints; iq++) {
      Lb[iq] *= elInfo->getDet();

      const double *psi = psiFast->getPhi(iq);
      grdPhi = phiFast->getGradient(iq);

      for (int i = 0; i < nRow; i++) {
	for (int j = 0; j < nCol; j++)
	  (*mat)[i][j] += quadrature->getWeight(iq) * (Lb[iq] * (*grdPhi)[j]) * psi[i];
      }
    }
  }

  void Quad10::calculateElementVector(const ElInfo *elInfo, ElementVector *vec)
  {
    VectorOfFixVecs<DimVec<double> > *grdPsi;

    if (firstCall) {
      const BasisFunction *basFcts = owner->getRowFESpace()->getBasisFcts();
      psiFast = updateFastQuadrature(psiFast, basFcts, INIT_GRD_PHI);
      basFcts = owner->getColFESpace()->getBasisFcts();
      phiFast = updateFastQuadrature(phiFast, basFcts, INIT_PHI);
      firstCall = false;
    }