Commit 7b38ea10 authored by Praetorius, Simon's avatar Praetorius, Simon

some changes in PetscSolverGlobalMatrix

parent 923b0774
......@@ -29,7 +29,7 @@ namespace AMDiS {
VecAtQP_FOT::VecAtQP_FOT(DOFVectorBase<double> *dv,
AbstractFunction<double, double> *af,
WorldVector<double> *wv)
: FirstOrderTerm(af->getDegree()), vec(dv), f(af), b(wv)
: FirstOrderTerm(af ? af->getDegree() : dv->getFeSpace()->getBasisFcts()->getDegree()), vec(dv), f(af), b(wv)
{
TEST_EXIT(dv)("No vector!\n");
......@@ -39,7 +39,7 @@ namespace AMDiS {
VecAtQP_FOT::VecAtQP_FOT(DOFVectorBase<double> *dv,
AbstractFunction<double, double> *af,
int bIdx)
: FirstOrderTerm(af->getDegree()), vec(dv), f(af)
: FirstOrderTerm(af ? af->getDegree() : dv->getFeSpace()->getBasisFcts()->getDegree()), vec(dv), f(af)
{
TEST_EXIT(dv)("No vector!\n");
......@@ -69,14 +69,29 @@ namespace AMDiS {
const int nPoints = static_cast<int>(Lb.size());
if (bOne > -1) {
for (int iq = 0; iq < nPoints; iq++)
lb_one(grdLambda, Lb[iq], (*f)(vecAtQPs[iq]));
if (f) {
for (int iq = 0; iq < nPoints; iq++)
lb_one(grdLambda, Lb[iq], (*f)(vecAtQPs[iq]));
} else {
for (int iq = 0; iq < nPoints; iq++)
lb_one(grdLambda, Lb[iq], vecAtQPs[iq]);
}
} else if (b) {
for (int iq = 0; iq < nPoints; iq++)
lb(grdLambda, *b, Lb[iq], (*f)(vecAtQPs[iq]));
if (f) {
for (int iq = 0; iq < nPoints; iq++)
lb(grdLambda, *b, Lb[iq], (*f)(vecAtQPs[iq]));
} else {
for (int iq = 0; iq < nPoints; iq++)
lb(grdLambda, *b, Lb[iq], vecAtQPs[iq]);
}
} else {
for (int iq = 0; iq < nPoints; iq++)
l1(grdLambda, Lb[iq], (*f)(vecAtQPs[iq]));
if (f) {
for (int iq = 0; iq < nPoints; iq++)
l1(grdLambda, Lb[iq], (*f)(vecAtQPs[iq]));
} else {
for (int iq = 0; iq < nPoints; iq++)
l1(grdLambda, Lb[iq], vecAtQPs[iq]);
}
}
}
......@@ -91,7 +106,7 @@ namespace AMDiS {
if (num_rows(grdUhAtQP) > 0) {
for (int iq = 0; iq < nPoints; iq++) {
double factor = (*f)(vecAtQPs[iq]);
double factor = (f ? (*f)(vecAtQPs[iq]) : vecAtQPs[iq]);
double resultQP = 0.0;
for (int i = 0; i < dow; i++)
resultQP += grdUhAtQP[iq][i];
......@@ -522,7 +537,11 @@ namespace AMDiS {
Vec2AtQP_FOT::Vec2AtQP_FOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2,
BinaryAbstractFunction<double, double, double> *af,
WorldVector<double> *b_)
: FirstOrderTerm(af->getDegree()), vec1(dv1), vec2(dv2), f(af), b(b_)
: FirstOrderTerm(af ?
af->getDegree() :
dv1->getFeSpace()->getBasisFcts()->getDegree()
+ dv2->getFeSpace()->getBasisFcts()->getDegree()),
vec1(dv1), vec2(dv2), f(af), b(b_)
{
TEST_EXIT(dv1)("No first vector!\n");
TEST_EXIT(dv2)("No second vector!\n");
......@@ -534,7 +553,11 @@ namespace AMDiS {
Vec2AtQP_FOT::Vec2AtQP_FOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2,
BinaryAbstractFunction<double, double, double> *af,
int bIdx)
: FirstOrderTerm(af->getDegree()), vec1(dv1), vec2(dv2), f(af)
: FirstOrderTerm(af ?
af->getDegree() :
dv1->getFeSpace()->getBasisFcts()->getDegree()
+ dv2->getFeSpace()->getBasisFcts()->getDegree()),
vec1(dv1), vec2(dv2), f(af)
{
TEST_EXIT(dv1)("No first vector!\n");
TEST_EXIT(dv2)("No second vector!\n");
......@@ -559,14 +582,29 @@ namespace AMDiS {
const int nPoints = static_cast<int>(Lb.size());
if (bOne > -1) {
for (int iq = 0; iq < nPoints; iq++)
lb_one(grdLambda, Lb[iq], (*f)(vec1AtQPs[iq], vec2AtQPs[iq]));
if (f) {
for (int iq = 0; iq < nPoints; iq++)
lb_one(grdLambda, Lb[iq], (*f)(vec1AtQPs[iq], vec2AtQPs[iq]));
} else {
for (int iq = 0; iq < nPoints; iq++)
lb_one(grdLambda, Lb[iq], vec1AtQPs[iq] * vec2AtQPs[iq]);
}
} else if (b) {
for (int iq = 0; iq < nPoints; iq++)
lb(grdLambda, *b, Lb[iq], (*f)(vec1AtQPs[iq], vec2AtQPs[iq]));
if (f) {
for (int iq = 0; iq < nPoints; iq++)
lb(grdLambda, *b, Lb[iq], (*f)(vec1AtQPs[iq], vec2AtQPs[iq]));
} else {
for (int iq = 0; iq < nPoints; iq++)
lb(grdLambda, *b, Lb[iq], vec1AtQPs[iq] * vec2AtQPs[iq]);
}
} else {
for (int iq = 0; iq < nPoints; iq++)
l1(grdLambda, Lb[iq], (*f)(vec1AtQPs[iq], vec2AtQPs[iq]));
if (f) {
for (int iq = 0; iq < nPoints; iq++)
l1(grdLambda, Lb[iq], (*f)(vec1AtQPs[iq], vec2AtQPs[iq]));
} else {
for (int iq = 0; iq < nPoints; iq++)
l1(grdLambda, Lb[iq], vec1AtQPs[iq] * vec2AtQPs[iq]);
}
}
}
......@@ -577,9 +615,15 @@ namespace AMDiS {
mtl::dense_vector<double>& result,
double fac)
{
if (num_rows(grdUhAtQP) > 0)
for (int iq = 0; iq < nPoints; iq++)
result[iq] += fac * (*f)(vec1AtQPs[iq], vec2AtQPs[iq]) * ((*b) * grdUhAtQP[iq]);
if (num_rows(grdUhAtQP) > 0) {
if (f) {
for (int iq = 0; iq < nPoints; iq++)
result[iq] += fac * (*f)(vec1AtQPs[iq], vec2AtQPs[iq]) * ((*b) * grdUhAtQP[iq]);
} else {
for (int iq = 0; iq < nPoints; iq++)
result[iq] += fac * vec1AtQPs[iq] * vec2AtQPs[iq] * ((*b) * grdUhAtQP[iq]);
}
}
}
......@@ -590,7 +634,11 @@ namespace AMDiS {
DOFVectorBase<double> *dv3,
TertiaryAbstractFunction<double, double, double, double> *f_,
WorldVector<double> *bvec)
: FirstOrderTerm(f_->getDegree()),
: FirstOrderTerm(f_ ?
f_->getDegree() :
dv1->getFeSpace()->getBasisFcts()->getDegree()
+ dv2->getFeSpace()->getBasisFcts()->getDegree()
+ dv3->getFeSpace()->getBasisFcts()->getDegree()),
vec1(dv1),
vec2(dv2),
vec3(dv3),
......@@ -607,7 +655,11 @@ namespace AMDiS {
DOFVectorBase<double> *dv3,
TertiaryAbstractFunction<double, double, double, double> *f_,
int b)
: FirstOrderTerm(f_->getDegree()),
: FirstOrderTerm(f_ ?
f_->getDegree() :
dv1->getFeSpace()->getBasisFcts()->getDegree()
+ dv2->getFeSpace()->getBasisFcts()->getDegree()
+ dv3->getFeSpace()->getBasisFcts()->getDegree()),
vec1(dv1),
vec2(dv2),
vec3(dv3),
......@@ -635,14 +687,28 @@ namespace AMDiS {
const int nPoints = static_cast<int>(Lb.size());
if (bOne > -1) {
for (int iq = 0; iq < nPoints; iq++)
lb_one(grdLambda, Lb[iq], (*f)(vec1AtQPs[iq], vec2AtQPs[iq], vec3AtQPs[iq]));
} else {
for (int iq = 0; iq < nPoints; iq++) {
if (b)
lb(grdLambda, *b, Lb[iq], (*f)(vec1AtQPs[iq], vec2AtQPs[iq], vec3AtQPs[iq]));
else
l1(grdLambda, Lb[iq], (*f)(vec1AtQPs[iq], vec2AtQPs[iq], vec3AtQPs[iq]));
if (f) {
for (int iq = 0; iq < nPoints; iq++)
lb_one(grdLambda, Lb[iq], (*f)(vec1AtQPs[iq], vec2AtQPs[iq], vec3AtQPs[iq]));
} else {
for (int iq = 0; iq < nPoints; iq++)
lb_one(grdLambda, Lb[iq], vec1AtQPs[iq] * vec2AtQPs[iq] * vec3AtQPs[iq]);
}
} else {
if (f) {
for (int iq = 0; iq < nPoints; iq++) {
if (b)
lb(grdLambda, *b, Lb[iq], (*f)(vec1AtQPs[iq], vec2AtQPs[iq], vec3AtQPs[iq]));
else
l1(grdLambda, Lb[iq], (*f)(vec1AtQPs[iq], vec2AtQPs[iq], vec3AtQPs[iq]));
}
} else {
for (int iq = 0; iq < nPoints; iq++) {
if (b)
lb(grdLambda, *b, Lb[iq], vec1AtQPs[iq] * vec2AtQPs[iq] * vec3AtQPs[iq]);
else
l1(grdLambda, Lb[iq], vec1AtQPs[iq] * vec2AtQPs[iq] * vec3AtQPs[iq]);
}
}
}
}
......@@ -660,9 +726,15 @@ namespace AMDiS {
if (bOne > -1) {
ERROR_EXIT("Not yet implemented!\n");
} else {
for (int iq = 0; iq < nPoints; iq++)
result[iq] += fac * (*f)(vec1AtQPs[iq], vec2AtQPs[iq] ,vec3AtQPs[iq]) *
((*b) * grdUhAtQP[iq]);
if (f) {
for (int iq = 0; iq < nPoints; iq++)
result[iq] += fac * (*f)(vec1AtQPs[iq], vec2AtQPs[iq] ,vec3AtQPs[iq]) *
((*b) * grdUhAtQP[iq]);
} else {
for (int iq = 0; iq < nPoints; iq++)
result[iq] += fac * vec1AtQPs[iq] * vec2AtQPs[iq] * vec3AtQPs[iq] *
((*b) * grdUhAtQP[iq]);
}
}
}
......
......@@ -187,8 +187,11 @@ namespace AMDiS {
// ========== VecAtQP_SOT ==========
VecAtQP_SOT::VecAtQP_SOT(DOFVectorBase<double> *dv,
AbstractFunction<double, double> *af)
: SecondOrderTerm(af ? af->getDegree() : dv->getFeSpace()->getBasisFcts()->getDegree()), vec(dv), f(af)
AbstractFunction<double, double> *af,
double factor_
)
: SecondOrderTerm(af ? af->getDegree() : dv->getFeSpace()->getBasisFcts()->getDegree()),
vec(dv), f(af), factor(factor_)
{
setSymmetric(true);
......@@ -220,10 +223,10 @@ namespace AMDiS {
if (f) {
for (int iq = 0; iq < nPoints; iq++)
l1lt(grdLambda, LALt[iq], (*f)(vecAtQPs[iq]));
l1lt(grdLambda, LALt[iq], factor*(*f)(vecAtQPs[iq]));
} else {
for (int iq = 0; iq < nPoints; iq++)
l1lt(grdLambda, LALt[iq], vecAtQPs[iq]);
l1lt(grdLambda, LALt[iq], factor*vecAtQPs[iq]);
}
}
......@@ -237,19 +240,20 @@ namespace AMDiS {
int dow = Global::getGeo(WORLD);
if (num_rows(D2UhAtQP) > 0) {
double fac_ = factor * fac;
if (f) {
for (int iq = 0; iq < nPoints; iq++) {
double resultQP = 0.0;
for (int i = 0; i < dow; i++)
resultQP += D2UhAtQP[iq][i][i];
result[iq] += fac * (*f)(vecAtQPs[iq]) * resultQP;
result[iq] += fac_ * (*f)(vecAtQPs[iq]) * resultQP;
}
} else {
for (int iq = 0; iq < nPoints; iq++) {
double resultQP = 0.0;
for (int i = 0; i < dow; i++)
resultQP += D2UhAtQP[iq][i][i];
result[iq] += fac * vecAtQPs[iq] * resultQP;
result[iq] += fac_ * vecAtQPs[iq] * resultQP;
}
}
}
......@@ -261,13 +265,13 @@ namespace AMDiS {
int nPoints = grdUhAtQP.size();
if (f) {
for (int iq = 0; iq < nPoints; iq++) {
double factor = (*f)(vecAtQPs[iq]);
axpy(factor, grdUhAtQP[iq], result[iq]);
double fac_ = factor * (*f)(vecAtQPs[iq]);
axpy(fac_, grdUhAtQP[iq], result[iq]);
}
} else {
for (int iq = 0; iq < nPoints; iq++) {
double factor = vecAtQPs[iq];
axpy(factor, grdUhAtQP[iq], result[iq]);
double fac_ = factor * vecAtQPs[iq];
axpy(fac_, grdUhAtQP[iq], result[iq]);
}
}
}
......@@ -276,8 +280,11 @@ namespace AMDiS {
// ========== Vec2AtQP_SOT ==========
Vec2AtQP_SOT::Vec2AtQP_SOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2,
BinaryAbstractFunction<double, double, double> *af)
: SecondOrderTerm(af ? af->getDegree() : 2*dv1->getFeSpace()->getBasisFcts()->getDegree()), vec1(dv1), vec2(dv2), f(af)
BinaryAbstractFunction<double, double, double> *af,
double factor_
)
: SecondOrderTerm(af ? af->getDegree() : 2*dv1->getFeSpace()->getBasisFcts()->getDegree()),
vec1(dv1), vec2(dv2), f(af), factor(factor_)
{
setSymmetric(true);
......@@ -304,9 +311,9 @@ namespace AMDiS {
for (int iq = 0; iq < nPoints; iq++) {
if (f)
l1lt(grdLambda, LALt[iq], (*f)(vecAtQPs1[iq], vecAtQPs2[iq]));
l1lt(grdLambda, LALt[iq], factor * (*f)(vecAtQPs1[iq], vecAtQPs2[iq]));
else
l1lt(grdLambda, LALt[iq], vecAtQPs1[iq] * vecAtQPs2[iq]);
l1lt(grdLambda, LALt[iq], factor * vecAtQPs1[iq] * vecAtQPs2[iq]);
}
}
......@@ -321,15 +328,16 @@ namespace AMDiS {
if (num_rows(D2UhAtQP) > 0) {
for (int iq = 0; iq < nPoints; iq++) {
double factor = 1.0;
double fac_ = 1.0;
if (f)
factor = (*f)(vecAtQPs1[iq], vecAtQPs2[iq]);
fac_ = (*f)(vecAtQPs1[iq], vecAtQPs2[iq]);
else
factor = vecAtQPs1[iq] * vecAtQPs2[iq];
fac_ = vecAtQPs1[iq] * vecAtQPs2[iq];
fac_ *= fac * factor;
double resultQP = 0.0;
for (int i = 0; i < dow; i++)
resultQP += D2UhAtQP[iq][i][i];
result[iq] += fac * factor * resultQP;
result[iq] += fac_ * resultQP;
}
}
}
......@@ -339,12 +347,13 @@ namespace AMDiS {
{
int nPoints = grdUhAtQP.size();
for (int iq = 0; iq < nPoints; iq++) {
double factor = 1.0;
double fac_ = 1.0;
if (f)
factor = (*f)(vecAtQPs1[iq], vecAtQPs2[iq]);
fac_ = (*f)(vecAtQPs1[iq], vecAtQPs2[iq]);
else
factor = vecAtQPs1[iq] * vecAtQPs2[iq];
axpy(factor, grdUhAtQP[iq], result[iq]);
fac_ = vecAtQPs1[iq] * vecAtQPs2[iq];
fac_ *= factor;
axpy(fac_, grdUhAtQP[iq], result[iq]);
}
}
......
......@@ -362,7 +362,8 @@ namespace AMDiS {
{
public:
/// Constructor.
VecAtQP_SOT(DOFVectorBase<double> *dv, AbstractFunction<double, double> *af = NULL);
VecAtQP_SOT(DOFVectorBase<double> *dv, AbstractFunction<double, double> *af = NULL,
double factor_ = 1.0);
/// Implementation of \ref OperatorTerm::initElement().
void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
......@@ -398,6 +399,8 @@ namespace AMDiS {
/// Function evaluated at \ref vecAtQPs.
AbstractFunction<double, double> *f;
double factor;
};
/**
......@@ -414,7 +417,8 @@ namespace AMDiS {
/// Constructor.
Vec2AtQP_SOT(DOFVectorBase<double> *dv1,
DOFVectorBase<double> *dv2,
BinaryAbstractFunction<double, double, double> *af = NULL);
BinaryAbstractFunction<double, double, double> *af = NULL,
double factor_ = 1.0);
/// Implementation of \ref OperatorTerm::initElement().
void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
......@@ -444,6 +448,8 @@ namespace AMDiS {
/// Function evaluated at \ref vecAtQPs.
BinaryAbstractFunction<double, double, double> *f;
double factor;
};
......
......@@ -58,7 +58,7 @@ namespace AMDiS
}
};
PetscSolverGlobalMatrix(string name);
PetscSolverGlobalMatrix(string name, bool setOptions = true);
void fillPetscMatrix(Matrix<DOFMatrix*> *mat);
......
......@@ -186,7 +186,7 @@ namespace AMDiS { namespace Parallel {
KSPCreate(domainComm, &ksp);
KSPSetOperators(ksp, getMatInterior(), getMatInterior(), SAME_NONZERO_PATTERN);
KSPMonitorSet(ksp, KSPMonitorTrueResidualNorm, PETSC_NULL, PETSC_NULL);
petsc_helper::setSolver(ksp, "ch_", KSPFGMRES, PCSHELL, 1e-6, 1e-8, 10000);
petsc_helper::setSolver(ksp, "ch_", KSPFGMRES, PCSHELL, getRelative(), getTolerance(), getMaxIterations());
setConstantNullSpace(ksp, componentSpaces[0]->getMesh()->getDim() , true);
}
......
......@@ -42,36 +42,7 @@ namespace AMDiS { namespace Parallel {
{
public:
/// Creator class
class IdFct : public AbstractFunction<double, double>
{
public:
IdFct()
: AbstractFunction<double, double>(0)
{}
double operator()(const double& x) const
{
return x;
}
};
class MultConstFct : public AbstractFunction<double, double>
{
public:
MultConstFct(double c)
: AbstractFunction<double, double>(1),
mConst(c)
{}
double operator()(const double& x) const
{
return mConst * x;
}
private:
double mConst;
};
class LinearInterpolation : public AbstractFunction<double, double>
{
public:
......@@ -108,20 +79,6 @@ namespace AMDiS { namespace Parallel {
double a,b,cmin,cmax;
};
class Multiplier3 : public BinaryAbstractFunction<double, double, double>
{
public:
Multiplier3()
: BinaryAbstractFunction<double, double, double >(0)
{}
double operator()(const double& phi, const double& phase) const
{
return phase * phi;
}
};
class EinsMinus : public AbstractFunction<double, double>
{
......@@ -138,12 +95,6 @@ namespace AMDiS { namespace Parallel {
private:
double c;
};
struct Multiplication : public BinaryAbstractFunction<double, double, double>{
double operator()(const double &v1, const double &v2) const{
return v2 / v1;
}
};
......@@ -238,8 +189,6 @@ namespace AMDiS { namespace Parallel {
DOFVector<double>* phase;
IdFct idFct;
};
} }
......
......@@ -54,7 +54,7 @@ namespace AMDiS { namespace Parallel {
PetscSolverNavierStokes::PetscSolverNavierStokes(string name)
: PetscSolverGlobalMatrix(name),
: PetscSolverGlobalMatrix(name, false),
pressureComponent(-1),
pressureNullSpace(true),
useOldInitialGuess(false),
......@@ -120,7 +120,7 @@ namespace AMDiS { namespace Parallel {
KSPCreate(domainComm, &ksp);
KSPSetOperators(ksp, getMatInterior(), getMatInterior(), SAME_NONZERO_PATTERN);
KSPMonitorSet(ksp, KSPMonitorTrueResidualNorm, PETSC_NULL, PETSC_NULL);
petsc_helper::setSolver(ksp, "ns_", KSPFGMRES, PCNONE, 1e-6, 1e-8, 10000);
petsc_helper::setSolver(ksp, "ns_", KSPFGMRES, PCNONE, getRelative(), getTolerance(), getMaxIterations());
// Create null space information.
if (pressureNullSpace)
......@@ -132,8 +132,7 @@ namespace AMDiS { namespace Parallel {
{
FUNCNAME("PetscSolverNavierStokes::initPreconditioner()");
MPI::COMM_WORLD.Barrier();
double wtime = MPI::Wtime();
Timer t;
TEST_EXIT(nu)("nu pointer not set!\n");
TEST_EXIT(invTau)("invtau pointer not set!\n");
......@@ -142,10 +141,8 @@ namespace AMDiS { namespace Parallel {
int dim = componentSpaces[pressureComponent]->getMesh()->getDim();
vector<int> velocityComponents;
velocityComponents.push_back(0);
velocityComponents.push_back(1);
if (dim == 3)
velocityComponents.push_back(2);
for (size_t i = 0; i < dim; i++)
velocityComponents.push_back(i);
PCSetType(pc, PCFIELDSPLIT);
PCFieldSplitSetType(pc, PC_COMPOSITE_SCHUR);
......@@ -200,10 +197,10 @@ namespace AMDiS { namespace Parallel {
{
Operator massOp(pressureFeSpace, pressureFeSpace);
ZeroOrderTerm *massTerm = NULL;
if ((!phase)||(*nu==0.0))
if ((!phase) || (*nu == 0.0))
massTerm = new Simple_ZOT;
else
massTerm = new VecAtQP_ZOT(phase, &idFct);
massTerm = new VecAtQP_ZOT(phase);
massOp.addTerm(massTerm);
massMatrix.assembleOperator(massOp);
delete massTerm;
......@@ -218,10 +215,10 @@ namespace AMDiS { namespace Parallel {
{
Operator laplaceOp(pressureFeSpace, pressureFeSpace);
SecondOrderTerm *laplaceTerm = NULL;
if ((!phase)||(*nu==0.0))
if ((!phase) || (*nu == 0.0))
laplaceTerm = new Simple_SOT;
else
laplaceTerm = new VecAtQP_SOT(phase, &idFct);
laplaceTerm = new VecAtQP_SOT(phase);
laplaceOp.addTerm(laplaceTerm);
laplaceMatrix.assembleOperator(laplaceOp);
delete laplaceTerm;
......@@ -237,8 +234,9 @@ namespace AMDiS { namespace Parallel {
DOFVector<double> vz(pressureFeSpace, "vz");
DOFVector<double> vp(pressureFeSpace, "vp");
vx.interpol(solution->getDOFVector(0));
vy.interpol(solution->getDOFVector(1));
if (dim == 3)
if (dim >= 2)
vy.interpol(solution->getDOFVector(1));
if (dim >= 3)
vz.interpol(solution->getDOFVector(2));
DOFMatrix conDifMatrix(pressureFeSpace, pressureFeSpace);
......@@ -257,32 +255,34 @@ namespace AMDiS { namespace Parallel {
conDifOp.addTerm(conDif0);
conDif1 = new Simple_SOT(*nu);
conDifOp.addTerm(conDif1);
conDif2 = new VecAtQP_FOT(&vx, &idFct, 0);
conDif2 = new VecAtQP_FOT(&vx, NULL, 0);
conDifOp.addTerm(conDif2, GRD_PHI);
conDif3 = new VecAtQP_FOT(&vy, &idFct, 1);
conDifOp.addTerm(conDif3, GRD_PHI);
if (dim >= 2) {
conDif3 = new VecAtQP_FOT(&vy, NULL, 1);
conDifOp.addTerm(conDif3, GRD_PHI);
}
if (dim == 3) {
conDif4 = new VecAtQP_FOT(&vz, &idFct, 2);
conDif4 = new VecAtQP_FOT(&vz, NULL, 2);
conDifOp.addTerm(conDif4, GRD_PHI);
}
} else {
} else { // no phase given
vp.interpol(phase);
if (*nu>0.0) {
conDif0 = new VecAtQP_ZOT(&vp, new MultConstFct(*invTau));
if (*nu > 0.0) {
conDif0 = new VecAtQP_ZOT(&vp, NULL, *invTau);
conDifOp.addTerm(conDif0);
conDif1 = new VecAtQP_SOT(&vp, new MultConstFct(*nu));
conDif1 = new VecAtQP_SOT(&vp, NULL, *nu);
conDifOp.addTerm(conDif1);
conDif2 = new Vec2AtQP_FOT(&vx, &vp, new Multiplier3(), 0);
conDif2 = new Vec2AtQP_FOT(&vx, &vp, NULL, 0);
conDifOp.addTerm(conDif2, GRD_PHI);
conDif3 = new Vec2AtQP_FOT(&vy, &vp, new Multiplier3(), 1);
conDifOp.addTerm(conDif3, GRD_PHI);
if (dim >= 2) {
conDif3 = new Vec2AtQP_FOT(&vy, &vp, NULL, 1);
conDifOp.addTerm(conDif3, GRD_PHI);
}
if (dim == 3) {
conDif4 = new Vec2AtQP_FOT(&vz, &vp, new Multiplier3(), 2);
conDif4 = new Vec2AtQP_FOT(&vz, &vp, NULL, 2);
conDifOp.addTerm(conDif4, GRD_PHI);
}