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(6), 
  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) * phase[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) * phase[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) * phase[iq] * fac;
    axpy(factor, grdUhAtQP[iq], result[iq]);
  }
}

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

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

Phase_FOT::Phase_FOT(DOFVectorBase<double>* phaseDV_, double fac_)
: FirstOrderTerm(6), 
  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(6),
  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(6),
  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(6),
  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]) ) );
}

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

Pow3_ZOT::Pow3_ZOT(DOFVectorBase<double> *rhoDV_, double fac_)
: ZeroOrderTerm(3), 
  rhoDV(rhoDV_), 
  fac(fac_)
{
  TEST_EXIT(rhoDV_)("No vector rho!\n");
  auxFeSpaces.insert(rhoDV_->getFeSpace());
}

void Pow3_ZOT::initElement(const ElInfo* elInfo, 
			SubAssembler* subAssembler,
			Quadrature *quad) 
{
  getVectorAtQPs(rhoDV, elInfo, subAssembler, quad, rho);
}

void Pow3_ZOT::initElement(const ElInfo* largeElInfo,
			const ElInfo* smallElInfo,
			SubAssembler* subAssembler,
			Quadrature *quad)
{
  getVectorAtQPs(rhoDV, smallElInfo, largeElInfo, subAssembler, quad, rho);
}

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

void Pow3_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) * uhAtQP[iq];
}

double Pow3_ZOT::f(const int iq) const
{
  return rho[iq]*sqr(rho[iq]);
}

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

Pow3Phase_ZOT::Pow3Phase_ZOT(DOFVectorBase<double>* phaseDV, DOFVectorBase<double> *rhoDV_, double fac_)
: Phase_ZOT(phaseDV, fac_),
  rhoDV(rhoDV_)
{
  TEST_EXIT(rhoDV_)("No vector rho!\n");
  auxFeSpaces.insert(rhoDV_->getFeSpace());
}

void Pow3Phase_ZOT::initElement(const ElInfo* elInfo, 
			SubAssembler* subAssembler,
			Quadrature *quad) 
{
  Phase_ZOT::initElement(elInfo,subAssembler,quad);
  getVectorAtQPs(rhoDV, elInfo, subAssembler, quad, rho);
}

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

double Pow3Phase_ZOT::f(const int iq) const
{
  return rho[iq]*sqr(rho[iq]);
}

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

Pow2_ZOT::Pow2_ZOT(DOFVectorBase<double> *rhoDV_, double fac_)
: ZeroOrderTerm(2), 
  rhoDV(rhoDV_), 
  fac(fac_)
{
  TEST_EXIT(rhoDV_)("No vector rho!\n");
  auxFeSpaces.insert(rhoDV_->getFeSpace());
}

void Pow2_ZOT::initElement(const ElInfo* elInfo, 
			SubAssembler* subAssembler,
			Quadrature *quad) 
{
  getVectorAtQPs(rhoDV, elInfo, subAssembler, quad, rho);
}

void Pow2_ZOT::initElement(const ElInfo* largeElInfo,
			const ElInfo* smallElInfo,
			SubAssembler* subAssembler,
			Quadrature *quad)
{
  getVectorAtQPs(rhoDV, smallElInfo, largeElInfo, subAssembler, quad, rho);
}

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

void Pow2_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) * uhAtQP[iq];
}

double Pow2_ZOT::f(const int iq) const
{
  return sqr(rho[iq]);
}

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

Pow2Phase_ZOT::Pow2Phase_ZOT(DOFVectorBase<double>* phaseDV, DOFVectorBase<double> *rhoDV_, double fac_)
: Phase_ZOT(phaseDV, fac_), 
  rhoDV(rhoDV_)
{
  TEST_EXIT(rhoDV_)("No vector rho!\n");
  auxFeSpaces.insert(rhoDV_->getFeSpace());
}

void Pow2Phase_ZOT::initElement(const ElInfo* elInfo, 
			SubAssembler* subAssembler,
			Quadrature *quad) 
{
  Phase_ZOT::initElement(elInfo,subAssembler,quad);
  getVectorAtQPs(rhoDV, elInfo, subAssembler, quad, rho);
}

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

double Pow2Phase_ZOT::f(const int iq) const
{
  return sqr(rho[iq]);
}

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

ConstrainedFrac_ZOT::ConstrainedFrac_ZOT(DOFVectorBase<double> *phaseDV_, DOFVectorBase<double> *rhoDV_, double fac_)
: Phase_ZOT(phaseDV_,fac_), 
  rhoDV(rhoDV_)
{
  TEST_EXIT(rhoDV_)("No vector rho!\n");
  auxFeSpaces.insert(rhoDV_->getFeSpace());
}

void ConstrainedFrac_ZOT::initElement(const ElInfo* elInfo, 
			SubAssembler* subAssembler,
			Quadrature *quad) 
{
  Phase_ZOT::initElement(elInfo, subAssembler,quad);
  getVectorAtQPs(rhoDV, elInfo, subAssembler, quad, rho);
}

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

double ConstrainedFrac_ZOT::f(const int iq) const
{
  return 1.0/std::max(1.e-6, rho[iq]+0.9);
}

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

ConstrainedFracSqr_ZOT::ConstrainedFracSqr_ZOT(DOFVectorBase<double> *phaseDV_, DOFVectorBase<double> *rhoDV_, double factor_)
: Phase_ZOT(phaseDV_,factor_), 
  rhoDV(rhoDV_)
{
  TEST_EXIT(rhoDV_)("No vector rho!\n");
  auxFeSpaces.insert(rhoDV_->getFeSpace());
}

void ConstrainedFracSqr_ZOT::initElement(const ElInfo* elInfo, 
			SubAssembler* subAssembler,
			Quadrature *quad) 
{
  Phase_ZOT::initElement(elInfo, subAssembler,quad);
  getVectorAtQPs(rhoDV, elInfo, subAssembler, quad, rho);
}

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

double ConstrainedFracSqr_ZOT::f(const int iq) const
{
  return -1.0/std::max(1.e-6, sqr(rho[iq]+0.9));
}

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

ConstrainedPowImpl_ZOT::ConstrainedPowImpl_ZOT(DOFVectorBase<double> *rhoDV_, double p_, double rhoMin_, double fac_)
: ZeroOrderTerm(2), 
  rhoDV(rhoDV_), 
  p(p_), 
  rhoMin(rhoMin_), 
  fac(fac_)
{
  TEST_EXIT(rhoDV_)("No vector rho!\n");
  auxFeSpaces.insert(rhoDV_->getFeSpace());
}

void ConstrainedPowImpl_ZOT::initElement(const ElInfo* elInfo, 
            SubAssembler* subAssembler,
            Quadrature *quad) 
{
  getVectorAtQPs(rhoDV, elInfo, subAssembler, quad, rho);
}

void ConstrainedPowImpl_ZOT::initElement(const ElInfo* largeElInfo,
            const ElInfo* smallElInfo,
            SubAssembler* subAssembler,
            Quadrature *quad)
{
  getVectorAtQPs(rhoDV, smallElInfo, largeElInfo, subAssembler, quad, rho);
}

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

void ConstrainedPowImpl_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) * uhAtQP[iq];
}

double ConstrainedPowImpl_ZOT::f(const int iq) const
{
  return 2.0*p*Helpers::toInt(0.0<=rhoMin-rho[iq]);
}

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

ConstrainedPowExpl_ZOT::ConstrainedPowExpl_ZOT(DOFVectorBase<double> *rhoDV_, double p_, double rhoMin_, double fac_) 
: ZeroOrderTerm(2), 
  rhoDV(rhoDV_), 
  p(p_), 
  rhoMin(rhoMin_), 
  fac(fac_)
{
  TEST_EXIT(rhoDV_)("No vector rho!\n");
  auxFeSpaces.insert(rhoDV_->getFeSpace());
}

void ConstrainedPowExpl_ZOT::initElement(const ElInfo* elInfo, 
            SubAssembler* subAssembler,
            Quadrature *quad) 
{
  getVectorAtQPs(rhoDV, elInfo, subAssembler, quad, rho);
}

void ConstrainedPowExpl_ZOT::initElement(const ElInfo* largeElInfo,
            const ElInfo* smallElInfo,
            SubAssembler* subAssembler,
            Quadrature *quad)
{
  getVectorAtQPs(rhoDV, smallElInfo, largeElInfo, subAssembler, quad, rho);
}

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

void ConstrainedPowExpl_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) * uhAtQP[iq];
}

double ConstrainedPowExpl_ZOT::f(const int iq) const
{
  return -2.0*p*std::max(0.0, rhoMin-rho[iq]);
}

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

ConstrainedPow_ZOT::ConstrainedPow_ZOT(DOFVectorBase<double> *rhoDV_, double p_, double rhoMin_, double fac_) 
: ZeroOrderTerm(2), 
  rhoDV(rhoDV_), 
  p(p_), 
  rhoMin(rhoMin_), 
  fac(fac_)
{
  TEST_EXIT(rhoDV_)("No vector rho!\n");
  auxFeSpaces.insert(rhoDV_->getFeSpace());
}

void ConstrainedPow_ZOT::initElement(const ElInfo* elInfo, 
            SubAssembler* subAssembler,
            Quadrature *quad) 
{
  getVectorAtQPs(rhoDV, elInfo, subAssembler, quad, rho);
}

void ConstrainedPow_ZOT::initElement(const ElInfo* largeElInfo,
            const ElInfo* smallElInfo,
            SubAssembler* subAssembler,
            Quadrature *quad)
{
  getVectorAtQPs(rhoDV, smallElInfo, largeElInfo, subAssembler, quad, rho);
}

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

void ConstrainedPow_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) * uhAtQP[iq];
}

double ConstrainedPow_ZOT::f(const int iq) const
{
  return -2.0*p*std::max(0.0, rhoMin-rho[iq]);
}

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

FactorPhase_ZOT::FactorPhase_ZOT(DOFVectorBase<double> *phaseDV_, DOFVectorBase<double> *rhoDV_, double fac_)
: Phase_ZOT(phaseDV_,fac_), 
  rhoDV(rhoDV_)
{
  TEST_EXIT(rhoDV_)("No vector rho!\n");
  auxFeSpaces.insert(rhoDV_->getFeSpace());
}

void FactorPhase_ZOT::initElement(const ElInfo* elInfo, 
			SubAssembler* subAssembler,
			Quadrature *quad) 
{
  Phase_ZOT::initElement(elInfo, subAssembler,quad);
  getVectorAtQPs(rhoDV, elInfo, subAssembler, quad, rho);
}

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

double FactorPhase_ZOT::f(const int iq) const
{
  return rho[iq];
}

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

WorldVecAndVecFct_FOT::WorldVecAndVecFct_FOT(WorldVector<DOFVectorBase<double>*> vecs_, DOFVectorBase<double> *phaseDV_, double fac_)
: FirstOrderTerm(2), 
  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(2), 
  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 WorldVec_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);
}

void WorldVec_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);
}

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];
    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];
    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(1),
  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(2), 
  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(2), 
  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++) {