POperators.cc

/******************************************************************************
 *
 * Extension of 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 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.
 *
 *
 * See also license.opensource.txt in the distribution.
 * 
 ******************************************************************************/


#include "POperators.h"
#include "Helpers.h"

using namespace std;
using namespace AMDiS;

Phase_SOT::Phase_SOT(DOFVectorBase<double>* phaseDV_, double fac_)
: SecondOrderTerm(phaseDV_->getFeSpace()->getBasisFcts()->getDegree()), 
  phaseDV(phaseDV_), 
  fac(fac_)
{
  TEST_EXIT(phaseDV_)("no vector phase!\n");
  auxFeSpaces.insert(phaseDV_->getFeSpace());
}

void Phase_SOT::initElement(const ElInfo* elInfo, 
			SubAssembler* subAssembler,
			Quadrature *quad)
{
  getVectorAtQPs(phaseDV, elInfo, subAssembler, quad, phase);
}

void Phase_SOT::initElement(const ElInfo* largeElInfo, const ElInfo* smallElInfo,
			SubAssembler* subAssembler,
			Quadrature *quad)
{
  getVectorAtQPs(phaseDV, smallElInfo, largeElInfo, subAssembler, quad, phase);
}

void Phase_SOT::getLALt(const ElInfo *elInfo,
			vector<mtl::dense2D<double> > &LALt) const
{
  const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
  const int nPoints = static_cast<int>(LALt.size());
  
  for (int iq = 0; iq < nPoints; iq++)
    l1lt(Lambda, LALt[iq], f(iq) * fac);
}

void Phase_SOT::eval(int nPoints,
		    const mtl::dense_vector<double>& uhAtQP,
		    const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
		    const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
		    mtl::dense_vector<double>& result,
		    double opFactor)
{	
  if (num_rows(D2UhAtQP) > 0) {
    for (int iq = 0; iq < nPoints; iq++) {
      double feval = f(iq) * opFactor * fac;
      double resultQP = 0.0;
      for (int i = 0; i < dimOfWorld; i++)
	resultQP += D2UhAtQP[iq][i][i];	
      result[iq] += feval * resultQP;
    }
  }
}

void Phase_SOT::weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
			std::vector<WorldVector<double> > &result)
{
  int nPoints = grdUhAtQP.size();
  for (int iq = 0; iq < nPoints; iq++) {
    double factor = f(iq) * fac;
    axpy(factor, grdUhAtQP[iq], result[iq]);
  }
}

double Phase_SOT::f(const int iq) const
{
  return std::max(0.0, std::min( 1.0, phase[iq] ) );
}

/* ----------------------------------------------------------- */

Phase_FOT::Phase_FOT(DOFVectorBase<double>* phaseDV_, double fac_)
: FirstOrderTerm(phaseDV_->getFeSpace()->getBasisFcts()->getDegree()), 
  phaseDV(phaseDV_), 
  fac(fac_)
{
  TEST_EXIT(phaseDV_)("no vector phase!\n");
  auxFeSpaces.insert(phaseDV_->getFeSpace());
}

void Phase_FOT::initElement(const ElInfo* elInfo, 
      SubAssembler* subAssembler,
      Quadrature *quad)
{
  getVectorAtQPs(phaseDV, elInfo, subAssembler, quad, phase);
}

void Phase_FOT::initElement(const ElInfo* largeElInfo, const ElInfo* smallElInfo,
      SubAssembler* subAssembler,
      Quadrature *quad)
{
  getVectorAtQPs(phaseDV, smallElInfo, largeElInfo, subAssembler, quad, phase);
}

void Phase_FOT::getLb(const ElInfo *elInfo, 
            vector<mtl::dense_vector<double> >& result) const
{
  const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
  const int nPoints = static_cast<int>(result.size());
  
  for (int iq = 0; iq < nPoints; iq++)
    lb(Lambda, f(iq), result[iq], phase[iq]*fac);
}

void Phase_FOT::eval(int nPoints,
		    const mtl::dense_vector<double>& uhAtQP,
		    const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
		    const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
		    mtl::dense_vector<double>& result,
		    double opFactor)
{   
  double factor = fac * opFactor;
  for (int iq = 0; iq < nPoints; iq++) {
    double resultQP = grdUhAtQP[iq] * f(iq);
    result[iq] += phase[iq] * factor * resultQP;
  }
}

WorldVector<double> Phase_FOT::f(const int iq) const
{
  WorldVector<double> result;
  result.set(1.0);
  return result;
}

/* ----------------------------------------------------------- */

Phase_ZOT::Phase_ZOT(DOFVectorBase<double>* phaseDV_, double fac_)
: ZeroOrderTerm(phaseDV_->getFeSpace()->getBasisFcts()->getDegree()),
  phaseDV(phaseDV_), 
  fac(fac_)
{	
  TEST_EXIT(phaseDV_)("No vector Phase!\n");
	
  auxFeSpaces.insert(phaseDV_->getFeSpace());
}

void Phase_ZOT::initElement(const ElInfo* elInfo, 
			SubAssembler* subAssembler,
			Quadrature *quad) 
{
  getVectorAtQPs(phaseDV, elInfo, subAssembler, quad, phase);
}

void Phase_ZOT::initElement(const ElInfo* largeElInfo,
			const ElInfo* smallElInfo,
			SubAssembler* subAssembler,
			Quadrature *quad)
{
  getVectorAtQPs(phaseDV, smallElInfo, largeElInfo, subAssembler, quad, phase);
}

void Phase_ZOT::getC(const ElInfo *elInfo, int nPoints, ElementVector &C)
{ 
  for (int iq = 0; iq < nPoints; iq++)
    C[iq] += f(iq) * phase[iq] * fac;
}

void Phase_ZOT::eval(int nPoints,
		    const mtl::dense_vector<double>& uhAtQP,
		    const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
		    const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
		    mtl::dense_vector<double>& result,
		    double opFactor)
{
  double factor = opFactor * fac;
  for (int iq = 0; iq < nPoints; iq++)
    result[iq] += factor * f(iq) * phase[iq] * uhAtQP[iq];
}

double Phase_ZOT::f(const int iq) const
{
  return 1.0;
}

/* ----------------------------------------------------------- */

PhaseInverse_ZOT::PhaseInverse_ZOT(DOFVectorBase<double>* phaseDV_, double fac_)
: ZeroOrderTerm(phaseDV_->getFeSpace()->getBasisFcts()->getDegree()),
  phaseDV(phaseDV_), 
  fac(fac_), 
  component(0)
{	
  TEST_EXIT(phaseDV_)("No vector Phase!\n");
  facVec = new WorldVector<double>();
  facVec->set(1.0);

  auxFeSpaces.insert(phaseDV_->getFeSpace());
}

PhaseInverse_ZOT::PhaseInverse_ZOT(DOFVectorBase<double>* phaseDV_, double fac_, WorldVector<double> *facVec_, int component_)
: ZeroOrderTerm(phaseDV_->getFeSpace()->getBasisFcts()->getDegree()),
  phaseDV(phaseDV_), 
  fac(fac_), 
  component(component_), 
  facVec(facVec_)
{	
  TEST_EXIT(phaseDV_)("No vector Phase!\n");
  
  auxFeSpaces.insert(phaseDV_->getFeSpace());
}

void PhaseInverse_ZOT::initElement(const ElInfo* elInfo, 
			SubAssembler* subAssembler,
			Quadrature *quad) 
{
  getVectorAtQPs(phaseDV, elInfo, subAssembler, quad, phase);
}

void PhaseInverse_ZOT::initElement(const ElInfo* largeElInfo,
			const ElInfo* smallElInfo,
			SubAssembler* subAssembler,
			Quadrature *quad)
{
  getVectorAtQPs(phaseDV, smallElInfo, largeElInfo, subAssembler, quad, phase);
}

void PhaseInverse_ZOT::getC(const ElInfo *, int nPoints, ElementVector &C)
{ 
  double factor = fac * (*facVec)[component];
  for (int iq = 0; iq < nPoints; iq++)
    C[iq] += f(iq) * factor;
}

void PhaseInverse_ZOT::eval(int nPoints,
			    const mtl::dense_vector<double>& uhAtQP,
			    const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
			    const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
			    mtl::dense_vector<double>& result,
			    double opFactor) 
{
  double factor = opFactor * fac * (*facVec)[component];
  for (int iq = 0; iq < nPoints; iq++)
    result[iq] += factor * f(iq) * uhAtQP[iq];
}

double PhaseInverse_ZOT::f(const int iq) const
{
  return std::max(0.0, std::min( 1.0, (1.0-phase[iq]) ) );
}

/* ----------------------------------------------------------- */

WorldVecAndVecFct_FOT::WorldVecAndVecFct_FOT(WorldVector<DOFVectorBase<double>*> vecs_, DOFVectorBase<double> *phaseDV_, double fac_)
: FirstOrderTerm(vecs_[0]->getFeSpace()->getBasisFcts()->getDegree() + phaseDV_->getFeSpace()->getBasisFcts()->getDegree()), 
  phaseDV(phaseDV_), 
  fac(fac_)
{
  numVecs=vecs_.size();
  TEST_EXIT(numVecs==2 || numVecs==3)("Only Dim=2 or Dim=3 possible\n");
  
  TEST_EXIT(phaseDV_)("phaseDV is NULL!\n");
  for (int i = 0; i < numVecs; i++) {
    TEST_EXIT(vecs_[i])("One vector is NULL!\n");

    auxFeSpaces.insert(vecs_[i]->getFeSpace());
  }
  
  vec0DV=vecs_[0];
  vec1DV=vecs_[1];
  if(numVecs>=3) vec2DV=vecs_[2];
}

void WorldVecAndVecFct_FOT::initElement(const ElInfo* elInfo, 
					SubAssembler* subAssembler,
					Quadrature *quad) 
{
  getVectorAtQPs(vec0DV, elInfo, subAssembler, quad, vec0);
  getVectorAtQPs(vec1DV, elInfo, subAssembler, quad, vec1);
  if(numVecs>=3) getVectorAtQPs(vec2DV, elInfo, subAssembler, quad, vec2);
  getVectorAtQPs(phaseDV, elInfo, subAssembler, quad, phase);
}

void WorldVecAndVecFct_FOT::initElement(const ElInfo* largeElInfo, const ElInfo* smallElInfo,
					SubAssembler* subAssembler,
					Quadrature *quad) 
{
  getVectorAtQPs(vec0DV, smallElInfo, largeElInfo, subAssembler, quad, vec0);
  getVectorAtQPs(vec1DV, smallElInfo, largeElInfo, subAssembler, quad, vec1);
  if(numVecs>=3) getVectorAtQPs(vec2DV, smallElInfo, largeElInfo, subAssembler, quad, vec2);
  getVectorAtQPs(phaseDV, smallElInfo, largeElInfo, subAssembler, quad, phase);
}

void WorldVecAndVecFct_FOT::getLb(const ElInfo *elInfo, 
			vector<mtl::dense_vector<double> >& result) const
{
  const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
  const int nPoints = static_cast<int>(result.size());
	
  for (int iq = 0; iq < nPoints; iq++) {		
    WorldVector<double> vec;
    vec[0]=vec0[iq];
    vec[1]=vec1[iq];
    if(numVecs>=3) vec[2]=vec2[iq];
    lb(Lambda, vec, result[iq], phase[iq]*fac);
  }
}

void WorldVecAndVecFct_FOT::eval(int nPoints,
			const mtl::dense_vector<double>& uhAtQP,
			const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
			const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
			mtl::dense_vector<double>& result,
			double opFactor)
{	
  double factor = fac * opFactor;
  for (int iq = 0; iq < nPoints; iq++) {
    double resultQP = 0.0;
		    
    resultQP += grdUhAtQP[iq][0] * vec0[iq];
    resultQP += grdUhAtQP[iq][1] * vec1[iq];
    if(numVecs>=3) resultQP += grdUhAtQP[iq][2] * vec2[iq];
    
    result[iq] += phase[iq] * factor * resultQP;
  }
}

/* ----------------------------------------------------------- */

WorldVec_FOT::WorldVec_FOT(WorldVector<DOFVector<double>*> vecs_, double fac_)
: FirstOrderTerm(vecs_[0]->getFeSpace()->getBasisFcts()->getDegree()), 
  fac(fac_)
{
  numVecs=vecs_.size();
  
  for (int i = 0; i < numVecs; i++) {
    TEST_EXIT(vecs_[i])("One vector is NULL!\n");	
    auxFeSpaces.insert(vecs_[i]->getFeSpace());
  }
  
  vec0DV=vecs_[0];
  if(numVecs>=2) vec1DV=vecs_[1];
  if(numVecs>=3) vec2DV=vecs_[2];
}

void WorldVec_FOT::initElement(const ElInfo* elInfo, 
					SubAssembler* subAssembler,
					Quadrature *quad) 
{
  getVectorAtQPs(vec0DV, elInfo, subAssembler, quad, vec0);
  if(numVecs>=2) getVectorAtQPs(vec1DV, elInfo, subAssembler, quad, vec1);
  if(numVecs>=3) getVectorAtQPs(vec2DV, elInfo, subAssembler, quad, vec2);
}

void WorldVec_FOT::initElement(const ElInfo* largeElInfo, const ElInfo* smallElInfo,
					SubAssembler* subAssembler,
					Quadrature *quad) 
{
  getVectorAtQPs(vec0DV, smallElInfo, largeElInfo, subAssembler, quad, vec0);
  if(numVecs>=2) getVectorAtQPs(vec1DV, smallElInfo, largeElInfo, subAssembler, quad, vec1);
  if(numVecs>=3) getVectorAtQPs(vec2DV, smallElInfo, largeElInfo, subAssembler, quad, vec2);
}

void WorldVec_FOT::getLb(const ElInfo *elInfo, 
			vector<mtl::dense_vector<double> >& result) const
{
  const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
  const int nPoints = static_cast<int>(result.size());
	
  for (int iq = 0; iq < nPoints; iq++) {
    WorldVector<double> vec;
    vec[0]=vec0[iq];
    if(numVecs>=2) vec[1]=vec1[iq];
    if(numVecs>=3) vec[2]=vec2[iq];
    lb(Lambda, vec, result[iq], fac);
  }
}

void WorldVec_FOT::eval(int nPoints,
			const mtl::dense_vector<double>& uhAtQP,
			const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
			const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
			mtl::dense_vector<double>& result,
			double opFactor)
{	
  double factor = fac * opFactor;
  for (int iq = 0; iq < nPoints; iq++) {
    double resultQP = 0.0;
		    
    resultQP += grdUhAtQP[iq][0] * vec0[iq];
    if(numVecs>=2) resultQP += grdUhAtQP[iq][1] * vec1[iq];
    if(numVecs>=3) resultQP += grdUhAtQP[iq][2] * vec2[iq];
    
    result[iq] += factor * resultQP;
  }
}

/* ----------------------------------------------------------- */

WorldVector_FOT::WorldVector_FOT(DOFVector<WorldVector<double> >* dv, double fac_)
: FirstOrderTerm(dv->getFeSpace()->getBasisFcts()->getDegree()),
  vec(dv),
  fac(fac_)
{
  TEST_EXIT(dv)("No vector!\n");

  auxFeSpaces.insert(dv->getFeSpace());
}

void WorldVector_FOT::initElement(const ElInfo* elInfo,
				  SubAssembler* subAssembler,
				  Quadrature *quad)
{
  getVectorAtQPs(vec, elInfo, subAssembler, quad, vecAtQPs);
}

void WorldVector_FOT::initElement(const ElInfo* largeElInfo,
				  const ElInfo* smallElInfo,
				  SubAssembler* subAssembler,
				  Quadrature *quad)
{
  getVectorAtQPs(vec, smallElInfo, largeElInfo, subAssembler, quad, vecAtQPs);
}

void WorldVector_FOT::getLb(const ElInfo *elInfo,
			    vector<mtl::dense_vector<double> >& Lb) const
{
  const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
  const int nPoints = static_cast<int>(Lb.size());

  for (int iq = 0; iq < nPoints; iq++)
    lb(grdLambda, vecAtQPs[iq], Lb[iq], fac);
}

void WorldVector_FOT::eval(int nPoints,
			   const mtl::dense_vector<double>&,
			   const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
			   const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
			   mtl::dense_vector<double>& result,
		 	   double opFactor)
{
  if (num_rows(grdUhAtQP) > 0) {
    double fac_ = fac * opFactor;
    for (int iq = 0; iq < nPoints; iq++)
      result[iq] += (vecAtQPs[iq] * grdUhAtQP[iq]) * fac_;
  }
}

/* ----------------------------------------------------------- */

WorldVecPhase_FOT::WorldVecPhase_FOT(DOFVectorBase<double> *phaseDV_, WorldVector<DOFVector<double>*> vecs_,double fac_)
: FirstOrderTerm(phaseDV_->getFeSpace()->getBasisFcts()->getDegree() + vecs_[0]->getFeSpace()->getBasisFcts()->getDegree()), 
  phaseDV(phaseDV_), 
  fac(fac_)
{
  numVecs=vecs_.size();
  TEST_EXIT(numVecs==2 || numVecs==3)("Only Dim=2 or Dim=3 possible\n");
  
  TEST_EXIT(phaseDV_)("phaseDV is NULL!\n");
  for (int i = 0; i < numVecs; i++) {
    TEST_EXIT(vecs_[i])("One vector is NULL!\n");

    auxFeSpaces.insert(vecs_[i]->getFeSpace());
  }
  auxFeSpaces.insert(phaseDV_->getFeSpace());
  
  vec0DV=vecs_[0];
  vec1DV=vecs_[1];
  if(numVecs>=3) vec2DV=vecs_[2];
}

void WorldVecPhase_FOT::initElement(const ElInfo* elInfo, 
                    SubAssembler* subAssembler,
                    Quadrature *quad) 
{
  getVectorAtQPs(vec0DV, elInfo, subAssembler, quad, vec0);
  getVectorAtQPs(vec1DV, elInfo, subAssembler, quad, vec1);
  if(numVecs>=3) getVectorAtQPs(vec2DV, elInfo, subAssembler, quad, vec2);
  getVectorAtQPs(phaseDV, elInfo, subAssembler, quad, phase);
}

void WorldVecPhase_FOT::initElement(const ElInfo* largeElInfo, const ElInfo* smallElInfo,
                    SubAssembler* subAssembler,
                    Quadrature *quad) 
{
  getVectorAtQPs(vec0DV, smallElInfo, largeElInfo, subAssembler, quad, vec0);
  getVectorAtQPs(vec1DV, smallElInfo, largeElInfo, subAssembler, quad, vec1);
  if(numVecs>=3) getVectorAtQPs(vec2DV, smallElInfo, largeElInfo, subAssembler, quad, vec2);
  getVectorAtQPs(phaseDV, smallElInfo, largeElInfo, subAssembler, quad, phase);
}

void WorldVecPhase_FOT::getLb(const ElInfo *elInfo, 
            vector<mtl::dense_vector<double> >& result) const
{
  const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
  const int nPoints = static_cast<int>(result.size());
  
  for (int iq = 0; iq < nPoints; iq++) {      
    WorldVector<double> vec;
    vec[0]=vec0[iq];
    vec[1]=vec1[iq];
    if(numVecs>=3) vec[2]=vec2[iq];
    lb(Lambda, vec, result[iq], phase[iq]*fac);
  }
}

void WorldVecPhase_FOT::eval(int nPoints,
            const mtl::dense_vector<double>& uhAtQP,
            const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
            const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
            mtl::dense_vector<double>& result,
            double opFactor)
{   
  double factor = fac * opFactor;
  for (int iq = 0; iq < nPoints; iq++) {
    double resultQP = 0.0;
	    
    resultQP += grdUhAtQP[iq][0] * vec0[iq];
    resultQP += grdUhAtQP[iq][1] * vec1[iq];
    if(numVecs>=3) resultQP += grdUhAtQP[iq][2] * vec2[iq];
    
    result[iq] += phase[iq] * factor * resultQP;
  }
}

/* ----------------------------------------------------------- */

VecAndWorldVec_FOT::VecAndWorldVec_FOT(DOFVector<double>* vec0_, WorldVector<DOFVector<double>*> vecs_, AbstractFunction<double, double>* fct_, double fac_)
: FirstOrderTerm(fct_->getDegree() + vecs_[0]->getFeSpace()->getBasisFcts()->getDegree()), 
  vDV(vec0_), 
  fct(fct_), 
  fac(fac_)
{
  numVecs = vecs_.size();
  TEST_EXIT(numVecs == 2 || numVecs == 3)
	("Only Dim=2 or Dim=3 possible\n");
  
  for (int i = 0; i < numVecs; i++) {
    TEST_EXIT(vecs_[i])("One vector is NULL!\n");	
    auxFeSpaces.insert(vecs_[i]->getFeSpace());
  }
  
  vec0DV = vecs_[0];
  vec1DV = vecs_[1];
  if (numVecs >= 3)
    vec2DV = vecs_[2];
}

void VecAndWorldVec_FOT::initElement(const ElInfo* elInfo, 
					SubAssembler* subAssembler,
					Quadrature *quad) 
{
  getVectorAtQPs(vDV, elInfo, subAssembler, quad, v);
  getVectorAtQPs(vec0DV, elInfo, subAssembler, quad, vec0);
  getVectorAtQPs(vec1DV, elInfo, subAssembler, quad, vec1);
  if (numVecs >= 3)
    getVectorAtQPs(vec2DV, elInfo, subAssembler, quad, vec2);
}

void VecAndWorldVec_FOT::initElement(const ElInfo* largeElInfo, const ElInfo* smallElInfo,
					SubAssembler* subAssembler,
					Quadrature *quad) 
{
  getVectorAtQPs(vDV, smallElInfo, largeElInfo, subAssembler, quad, v);
  getVectorAtQPs(vec0DV, smallElInfo, largeElInfo, subAssembler, quad, vec0);
  getVectorAtQPs(vec1DV, smallElInfo, largeElInfo, subAssembler, quad, vec1);
  if (numVecs >= 3)
    getVectorAtQPs(vec2DV, smallElInfo, largeElInfo, subAssembler, quad, vec2);
}

void VecAndWorldVec_FOT::getLb(const ElInfo *elInfo, 
			vector<mtl::dense_vector<double> >& result) const
{
  const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
  const int nPoints = static_cast<int>(result.size());
  
  for (int iq = 0; iq < nPoints; iq++) {
    double f = (*fct)(v[iq]);
    WorldVector<double> vec;
    vec[0] = vec0[iq];
    vec[1] = vec1[iq];
    if (numVecs >= 3)
      vec[2] = vec2[iq];
    lb(Lambda, f*vec, result[iq], fac);
  }
}

void VecAndWorldVec_FOT::eval(int nPoints,
			const mtl::dense_vector<double>& uhAtQP,
			const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
			const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
			mtl::dense_vector<double>& result,
			double opFactor)
{	
  double factor = fac * opFactor;
  for (int iq = 0; iq < nPoints; iq++) {
    double resultQP = 0.0;
		    
    resultQP += grdUhAtQP[iq][0] * vec0[iq];
    resultQP += grdUhAtQP[iq][1] * vec1[iq];
    if (numVecs >= 3)
      resultQP += grdUhAtQP[iq][2] * vec2[iq];
    
    result[iq] += factor * resultQP * (*fct)(v[iq]);
  }
}

/* ----------------------------------------------------------- */

// vec1*vec2*d/dxi(...)
VecAndPartialDerivative_FOT::VecAndPartialDerivative_FOT(DOFVectorBase<double> *dv1,
      int component_, double fac_)
: FirstOrderTerm(dv1->getFeSpace()->getBasisFcts()->getDegree()), 
  vec1DV(dv1), 
  fac(fac_), 
  component(component_)
{
  auxFeSpaces.insert(vec1DV->getFeSpace());

  setB(component);
}

void VecAndPartialDerivative_FOT::initElement(const ElInfo* elInfo, 
        SubAssembler* subAssembler,
        Quadrature *quad) 
{ 
  getVectorAtQPs(vec1DV, elInfo, subAssembler, quad, vec1);
}

void VecAndPartialDerivative_FOT::initElement(const ElInfo* largeElInfo, const ElInfo* smallElInfo,
  SubAssembler* subAssembler,
  Quadrature *quad)
{ 
  getVectorAtQPs(vec1DV, smallElInfo, largeElInfo, subAssembler, quad, vec1);
}

void VecAndPartialDerivative_FOT::getLb(const ElInfo *elInfo,
      vector<mtl::dense_vector<double> >& result) const
{
  const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
    const int nPoints = static_cast<int>(result.size());
    
  for(int iq = 0; iq < nPoints; iq++) {
    lb_one(Lambda, result[iq], vec1[iq] * fac);
  }
}

void VecAndPartialDerivative_FOT::eval(int nPoints,
      const mtl::dense_vector<double>& uhAtQP,
      const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
      const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
      mtl::dense_vector<double>& result,
      double opFactor)
{
  for(int iq = 0; iq < nPoints; iq++) 
    result[iq] += opFactor * vec1[iq] * grdUhAtQP[iq][component] * fac;
}

/* ----------------------------------------------------------- */

// vec1*vec2*d/dxi(...)
VecAndPartialDerivativeIJ_FOT::VecAndPartialDerivativeIJ_FOT(DOFVectorBase<double> *dv1, 
							    DOFVectorBase<double> *dv2,
							    int i_,
							    int j_,
							    double fac_)
: FirstOrderTerm(dv1->getFeSpace()->getBasisFcts()->getDegree() + dv2->getFeSpace()->getBasisFcts()->getDegree() - 1), 
  vec1DV(dv1), 
  vec2DV(dv2), 
  fct(NULL), 
  fac(fac_), 
  xi(i_), 
  xj(j_)
{
  auxFeSpaces.insert(vec1DV->getFeSpace());
  auxFeSpaces.insert(vec2DV->getFeSpace());

  setB(xi);
}

VecAndPartialDerivativeIJ_FOT::VecAndPartialDerivativeIJ_FOT(DOFVectorBase<double> *dv1, 
							  DOFVectorBase<double> *dv2,
							  int i_,
							  int j_,
							  BinaryAbstractFunction<double, double, double>* fct_, // = f(v1, d_j(v2))
							  double fac_)
: FirstOrderTerm(fct_->getDegree()), 
  vec1DV(dv1), 
  vec2DV(dv2), 
  fct(fct_), 
  fac(fac_), 
  xi(i_), 
  xj(j_)
{
  auxFeSpaces.insert(vec1DV->getFeSpace());
  auxFeSpaces.insert(vec2DV->getFeSpace());

  setB(xi);
}

void VecAndPartialDerivativeIJ_FOT::initElement(const ElInfo* elInfo, 
        SubAssembler* subAssembler,
        Quadrature *quad) 
{ 
  getVectorAtQPs(vec1DV, elInfo, subAssembler, quad, vec1);
  getGradientsAtQPs(vec2DV, elInfo, subAssembler, quad, grad);
}

void VecAndPartialDerivativeIJ_FOT::initElement(const ElInfo* largeElInfo, const ElInfo* smallElInfo,
  SubAssembler* subAssembler,
  Quadrature *quad)
{ 
  getVectorAtQPs(vec1DV, smallElInfo, largeElInfo, subAssembler, quad, vec1);
  getGradientsAtQPs(vec1DV, smallElInfo, largeElInfo, subAssembler, quad, grad);
}

void VecAndPartialDerivativeIJ_FOT::getLb(const ElInfo *elInfo,
      vector<mtl::dense_vector<double> >& result) const
{
  const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
  const int nPoints = static_cast<int>(result.size());
  
  if (fct == NULL) {
    for(int iq = 0; iq < nPoints; iq++) {
      lb_one(Lambda, result[iq], vec1[iq] * grad[iq][xj] * fac);
    }
  } else {
    for(int iq = 0; iq < nPoints; iq++) {
      lb_one(Lambda, result[iq], (*fct)(vec1[iq], grad[iq][xj]) * fac);
    }
  }
}

void VecAndPartialDerivativeIJ_FOT::eval(int nPoints,
      const mtl::dense_vector<double>& uhAtQP,
      const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
      const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
      mtl::dense_vector<double>& result,
      double opFactor)
{
  if (fct == NULL) {
    for(int iq = 0; iq < nPoints; iq++) 
      result[iq] += opFactor * vec1[iq] * grad[iq][xj] * grdUhAtQP[iq][xi] * fac;
  } else {
    for(int iq = 0; iq < nPoints; iq++) 
      result[iq] += opFactor * (*fct)(vec1[iq], grad[iq][xj]) * grdUhAtQP[iq][xi] * fac;
  }
}

/* ----------------------------------------------------------- */

// vec1*vec2*d/dxi(...)
Vec2AndPartialDerivative_FOT::Vec2AndPartialDerivative_FOT(
			       DOFVectorBase<double> *dv1, 
			       DOFVectorBase<double> *dv2,
			       int component_, 
			       double fac_)
: FirstOrderTerm(dv1->getFeSpace()->getBasisFcts()->getDegree() + dv2->getFeSpace()->getBasisFcts()->getDegree()), 
  vec1DV(dv1), 
  vec2DV(dv2), 
  fct(NULL), 
  fac(fac_), 
  component(component_)
{
  auxFeSpaces.insert(vec1DV->getFeSpace());
  auxFeSpaces.insert(vec2DV->getFeSpace());

  setB(component);
}


Vec2AndPartialDerivative_FOT::Vec2AndPartialDerivative_FOT(
			       DOFVectorBase<double> *dv1, 
			       DOFVectorBase<double> *dv2,
			       int component_,
			       BinaryAbstractFunction<double, double, double>* fct_, // = f(v1, v2)
			       double fac_)
: FirstOrderTerm(fct_->getDegree()), 
  vec1DV(dv1), 
  vec2DV(dv2), 
  fct(fct_), 
  fac(fac_), 
  component(component_)
{
  auxFeSpaces.insert(vec1DV->getFeSpace());
  auxFeSpaces.insert(vec2DV->getFeSpace());

  setB(component);
}

void Vec2AndPartialDerivative_FOT::initElement(const ElInfo* elInfo, 
				SubAssembler* subAssembler,
				Quadrature *quad) 
{ 
  getVectorAtQPs(vec1DV, elInfo, subAssembler, quad, vec1);
  getVectorAtQPs(vec2DV, elInfo, subAssembler, quad, vec2);
}

void Vec2AndPartialDerivative_FOT::initElement(const ElInfo* largeElInfo, const ElInfo* smallElInfo,
	SubAssembler* subAssembler,
	Quadrature *quad)
{ 
  getVectorAtQPs(vec1DV, smallElInfo, largeElInfo, subAssembler, quad, vec1);
  getVectorAtQPs(vec2DV, smallElInfo, largeElInfo, subAssembler, quad, vec2);
}

void Vec2AndPartialDerivative_FOT::getLb(const ElInfo *elInfo,
			vector<mtl::dense_vector<double> >& result) const
{
  const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
  const int nPoints = static_cast<int>(result.size());
  
  if (fct == NULL) {
    for(int iq = 0; iq < nPoints; iq++)
      lb_one(Lambda, result[iq], vec1[iq] * vec2[iq] * fac);
  } else {
    for(int iq = 0; iq < nPoints; iq++)
      lb_one(Lambda, result[iq], (*fct)(vec1[iq], vec2[iq]) * fac);
  }
}

void Vec2AndPartialDerivative_FOT::eval(int nPoints,
			const mtl::dense_vector<double>& uhAtQP,
			const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
			const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
			mtl::dense_vector<double>& result,
			double opFactor)
{
  if (fct == NULL) {
    for(int iq = 0; iq < nPoints; iq++) 
      result[iq] += opFactor * vec1[iq] * vec2[iq] * grdUhAtQP[iq][component] * fac;
  } else {
    for(int iq = 0; iq < nPoints; iq++) 
      result[iq] += opFactor * (*fct)(vec1[iq], vec2[iq]) * grdUhAtQP[iq][component] * fac;
  }
}

/* ----------------------------------------------------------- */

PartialDerivative_FOT::PartialDerivative_FOT(int component_, double fac_)
: FirstOrderTerm(0), 
  component(component_), 
  fac(fac_)
{
  setB(component);
}

void PartialDerivative_FOT::initElement(const ElInfo* elInfo, 
					SubAssembler* subAssembler,
					Quadrature *quad) {}
					
void PartialDerivative_FOT::initElement(const ElInfo* largeElInfo, const ElInfo* smallElInfo,
					SubAssembler* subAssembler,
					Quadrature *quad) {}
					
void PartialDerivative_FOT::getLb(const ElInfo *elInfo, 
			vector<mtl::dense_vector<double> >& result) const
{
  const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
  const int nPoints = static_cast<int>(result.size());
	
  for (int iq = 0; iq < nPoints; iq++) {	
    lb_one(Lambda, result[iq], fac);
  }
}

void PartialDerivative_FOT::eval(int nPoints,
			const mtl::dense_vector<double>& uhAtQP,
			const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
			const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
			mtl::dense_vector<double>& result,
			double opFactor)
{	
  for (int iq = 0; iq < nPoints; iq++)	
    result[iq] += opFactor * grdUhAtQP[iq][component] * fac;
}

/* ----------------------------------------------------------- */

PartialDerivative_ZOT::PartialDerivative_ZOT(DOFVectorBase<double> *vecDV_, int component_, double fac_)
: ZeroOrderTerm(vecDV_->getFeSpace()->getBasisFcts()->getDegree() - 1), 
  vecDV(vecDV_), 
  fac(fac_), 
  component(component_)
{
  TEST_EXIT(vecDV_)("No first vector!\n");
  auxFeSpaces.insert(vecDV_->getFeSpace());

  facVec = new WorldVector<double>(); facVec->set(1.0);
  component2 = 0;
}

PartialDerivative_ZOT::PartialDerivative_ZOT(DOFVectorBase<double> *vecDV_, int component_, double fac_, WorldVector<double> *facVec_, int component2_)
: ZeroOrderTerm(vecDV_->getFeSpace()->getBasisFcts()->getDegree() - 1), 
  vecDV(vecDV_), 
  fac(fac_), 
  facVec(facVec_), 
  component(component_), 
  component2(component2_)
{
  TEST_EXIT(vecDV_)("No first vector!\n");
  auxFeSpaces.insert(vecDV_->getFeSpace());
}

void PartialDerivative_ZOT::initElement(const ElInfo* elInfo, 
					SubAssembler* subAssembler,
					Quadrature *quad)
{
  getGradientsAtQPs(vecDV, elInfo, subAssembler, quad, grad);
}

void PartialDerivative_ZOT::initElement(const ElInfo* largeElInfo, const ElInfo* smallElInfo,
					SubAssembler* subAssembler,
					Quadrature *quad)
{
  getGradientsAtQPs(vecDV, smallElInfo, largeElInfo, subAssembler, quad, grad);
}

void PartialDerivative_ZOT::getC(const ElInfo *, int nPoints, ElementVector &C)
{ 	
	double factor = fac * (*facVec)[component2];
	for (int iq = 0; iq < nPoints; iq++)
		C[iq] += grad[iq][component]*factor;
}

void PartialDerivative_ZOT::eval(int nPoints,
				const mtl::dense_vector<double>& uhAtQP,
				const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
				const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
				mtl::dense_vector<double>& result,
				double opFactor)
{	
  double factor = opFactor * fac * (*facVec)[component2];
  for (int iq = 0; iq < nPoints; iq++)		
    result[iq] += grad[iq][component] * uhAtQP[iq] * factor;
}

/* ----------------------------------------------------------- */

VecAndPartialDerivative_ZOT::VecAndPartialDerivative_ZOT(DOFVectorBase<double> *vecDV_, DOFVectorBase<double> *gradDV_, int component_,  double fac_)
: ZeroOrderTerm(vecDV_->getFeSpace()->getBasisFcts()->getDegree() + gradDV_->getFeSpace()->getBasisFcts()->getDegree() - 1),
  vecDV(vecDV_),
  gradDV(gradDV_),
  component(component_),
  fct(NULL),
  fac(fac_)
{
  TEST_EXIT(vecDV_)("No value vector!\n");
  TEST_EXIT(gradDV_)("No gradient vector!\n");

  auxFeSpaces.insert(vecDV_->getFeSpace());
  auxFeSpaces.insert(gradDV_->getFeSpace());
}

VecAndPartialDerivative_ZOT::VecAndPartialDerivative_ZOT(DOFVectorBase<double> *vecDV_, 
							DOFVectorBase<double> *gradDV_, 
							int component_,
							AbstractFunction<double, double>* fct_,
							double fac_)
: ZeroOrderTerm(fct_->getDegree() + gradDV_->getFeSpace()->getBasisFcts()->getDegree() - 1),
  vecDV(vecDV_),
  gradDV(gradDV_),
  component(component_),
  fct(fct_),
  fac(fac_)
{
  TEST_EXIT(vecDV_)("No value vector!\n");
  TEST_EXIT(gradDV_)("No gradient vector!\n");

  auxFeSpaces.insert(vecDV_->getFeSpace());
  auxFeSpaces.insert(gradDV_->getFeSpace());
}

void VecAndPartialDerivative_ZOT::initElement(const ElInfo* elInfo,
					      SubAssembler* subAssembler,
					      Quadrature *quad)
{
  getVectorAtQPs(vecDV, elInfo, subAssembler, quad, vec);
  getGradientsAtQPs(gradDV, elInfo, subAssembler, quad, grad);
}

void VecAndPartialDerivative_ZOT::initElement(const ElInfo* largeElInfo,
					      const ElInfo* smallElInfo,
					      SubAssembler* subAssembler,
					      Quadrature *quad)
{
  getVectorAtQPs(vecDV, smallElInfo, largeElInfo, subAssembler, quad, vec);
  getGradientsAtQPs(gradDV, smallElInfo, largeElInfo, subAssembler, quad, grad);
}

void VecAndPartialDerivative_ZOT::getC(const ElInfo *, int nPoints, ElementVector &C)
{
  if (fct == NULL) {
    for (int iq = 0; iq < nPoints; iq++)
      C[iq] += vec[iq] * grad[iq][component] * fac;
  } else {
    for (int iq = 0; iq < nPoints; iq++)
      C[iq] += (*fct)(vec[iq]) * grad[iq][component] * fac;
  }	
}

void VecAndPartialDerivative_ZOT::eval(int nPoints,
				       const mtl::dense_vector<double>& uhAtQP,
				       const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
				       const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
				       mtl::dense_vector<double>& result,
				       double opFactor)
{
  double fac_ = fac * opFactor;
  if (fct == NULL) {
    for (int iq = 0; iq < nPoints; iq++)
      result[iq] += vec[iq] * grad[iq][component] * uhAtQP[iq] * fac_;
  } else {
    for (int iq = 0; iq < nPoints; iq++)
      result[iq] += (*fct)(vec[iq]) * grad[iq][component] * uhAtQP[iq] * fac_;
  }
}

/* ----------------------------------------------------------- */

Vec2AndPartialDerivative_ZOT::Vec2AndPartialDerivative_ZOT(DOFVectorBase<double> *vec1DV_, 
							    DOFVectorBase<double> *vec2DV_, 
							    DOFVectorBase<double> *gradDV_,
							    int component_,  
							    double fac_)
: ZeroOrderTerm(vec1DV_->getFeSpace()->getBasisFcts()->getDegree() + vec2DV_->getFeSpace()->getBasisFcts()->getDegree() + gradDV_->getFeSpace()->getBasisFcts()->getDegree() - 1), 
  vec1DV(vec1DV_), 
  vec2DV(vec2DV_), 
  gradDV(gradDV_), 
  fct(NULL), 
  fac(fac_), 
  component(component_)
{
  TEST_EXIT(vec1DV_)("No first vector!\n");
  TEST_EXIT(vec2DV_)("No second vector!\n");
  TEST_EXIT(gradDV_)("No gradient vector!\n");
  auxFeSpaces.insert(vec1DV_->getFeSpace());
  auxFeSpaces.insert(vec2DV_->getFeSpace());
  auxFeSpaces.insert(gradDV_->getFeSpace());
}

Vec2AndPartialDerivative_ZOT::Vec2AndPartialDerivative_ZOT(DOFVectorBase<double> *vec1DV_, 
							  DOFVectorBase<double> *vec2DV_, 
							  DOFVectorBase<double> *gradDV_,
							  int component_,  
							  BinaryAbstractFunction<double, double, double>* fct_,
							  double fac_)
: ZeroOrderTerm(fct_->getDegree() + gradDV_->getFeSpace()->getBasisFcts()->getDegree() - 1), 
  vec1DV(vec1DV_), 
  vec2DV(vec2DV_), 
  gradDV(gradDV_), 
  fct(fct_), 
  fac(fac_), 
  component(component_)
{
  TEST_EXIT(vec1DV_)("No first vector!\n");
  TEST_EXIT(vec2DV_)("No second vector!\n");
  TEST_EXIT(gradDV_)("No gradient vector!\n");
  auxFeSpaces.insert(vec1DV_->getFeSpace());
  auxFeSpaces.insert(vec2DV_->getFeSpace());
  auxFeSpaces.insert(gradDV_->getFeSpace());
}

void Vec2AndPartialDerivative_ZOT::initElement(const ElInfo* elInfo,
					SubAssembler* subAssembler,
					Quadrature *quad)
{
  getVectorAtQPs(vec1DV, elInfo, subAssembler, quad, vec1);
  getVectorAtQPs(vec2DV, elInfo, subAssembler, quad, vec2);
  getGradientsAtQPs(gradDV, elInfo, subAssembler, quad, grad);
}

void Vec2AndPartialDerivative_ZOT::initElement(const ElInfo* largeElInfo, const ElInfo* smallElInfo,
					SubAssembler* subAssembler,
					Quadrature *quad)
{
  getVectorAtQPs(vec1DV, smallElInfo, largeElInfo, subAssembler, quad, vec1);
  getVectorAtQPs(vec2DV, smallElInfo, largeElInfo, subAssembler, quad, vec2);
  getGradientsAtQPs(gradDV, smallElInfo, largeElInfo, subAssembler, quad, grad);
}

void Vec2AndPartialDerivative_ZOT::getC(const ElInfo *, int nPoints, ElementVector &C)
{
  if (fct == NULL) {
    for (int iq = 0; iq < nPoints; iq++)
      C[iq] += vec1[iq]*vec2[iq]*grad[iq][component]*fac;
  } else {
    for (int iq = 0; iq < nPoints; iq++)
      C[iq] += (*fct)(vec1[iq],vec2[iq])*grad[iq][component]*fac;
  }
}

void Vec2AndPartialDerivative_ZOT::eval(int nPoints,
				const mtl::dense_vector<double>& uhAtQP,
				const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
				const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
				mtl::dense_vector<double>& result,
				double opFactor)
{
  if (fct == NULL) {
    for (int iq = 0; iq < nPoints; iq++) {
      result[iq] += (vec1[iq] * vec2[iq]) * grad[iq][component] * uhAtQP[iq] * opFactor * fac;
    }
  } else {
    for (int iq = 0; iq < nPoints; iq++) {
      result[iq] += (*fct)(vec1[iq], vec2[iq]) * grad[iq][component] * uhAtQP[iq] * opFactor * fac;
    }
  }
}

/* ----------------------------------------------------------- */

MatrixPhase_SOT::MatrixPhase_SOT(DOFVectorBase<double>* phaseDV_, int comp_, double fac_) :