Skip to content

GitLab

A
amdis
Commit 6ff14ecb authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files
Options

Navier-Stokes corrected

parent 515af748
Show whitespace changes
Inline Side-by-side
Showing with 255 additions and 290 deletions
extensions/ExtendedProblemStat.h
View file @ 6ff14ecb
......@@ -129,7 +129,7 @@ public:
if (asmMatrix && (singularDirichletBC.size() > 0 || manualPeriodicBC.size() > 0)) {
solverMatrix.setMatrix(*getSystemMatrix());
}
// writeMatrix("matrix.mat");
#endif
}
......
extensions/POperators_ZOT.h
View file @ 6ff14ecb
......@@ -22,9 +22,9 @@ void initElement(const ElInfo* largeElInfo, const ElInfo* smallElInfo,
SubAssembler* subAssembler,
Quadrature *quad=NULL);
inline void getC(const ElInfo *, int nPoints, ElementVector &C);
void getC(const ElInfo *, int nPoints, ElementVector &C);
inline void eval(int nPoints,
void eval(int nPoints,
const mtl::dense_vector<double>&,
const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
......@@ -76,16 +76,16 @@ void initElement(const ElInfo* largeElInfo, const ElInfo* smallElInfo,
SubAssembler* subAssembler,
Quadrature *quad=NULL);
inline void getC(const ElInfo *, int nPoints, ElementVector &C);
void getC(const ElInfo *, int nPoints, ElementVector &C);
inline void eval(int nPoints,
void 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);
inline double f(const int iq) const;
double f(const int iq) const;
void setFactor(const double fac_) { fac=fac_; }
protected:
......@@ -109,15 +109,15 @@ public:
void initElement(const ElInfo* largeElInfo, const ElInfo* smallElInfo,
SubAssembler*,
Quadrature *quad = NULL);
inline void getC(const ElInfo *, int nPoints, ElementVector &C);
void getC(const ElInfo *, int nPoints, ElementVector &C);
inline void eval(int nPoints,
void 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);
inline double f(const int iq) const;
double f(const int iq) const;
protected:
DOFVectorBase<double> *rhoDV;
mtl::dense_vector<double> rho;
......@@ -137,7 +137,7 @@ public:
void initElement(const ElInfo* largeElInfo, const ElInfo* smallElInfo,
SubAssembler*,
Quadrature *quad = NULL);
inline double f(const int iq) const;
double f(const int iq) const;
protected:
DOFVectorBase<double> *rhoDV;
mtl::dense_vector<double> rho;
......@@ -155,15 +155,15 @@ public:
void initElement(const ElInfo* largeElInfo, const ElInfo* smallElInfo,
SubAssembler*,
Quadrature *quad = NULL);
inline void getC(const ElInfo *, int nPoints, ElementVector &C);
void getC(const ElInfo *, int nPoints, ElementVector &C);
inline void eval(int nPoints,
void 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);
inline double f(const int iq) const;
double f(const int iq) const;
protected:
DOFVectorBase<double> *rhoDV;
mtl::dense_vector<double> rho;
......@@ -183,7 +183,7 @@ public:
void initElement(const ElInfo* largeElInfo, const ElInfo* smallElInfo,
SubAssembler*,
Quadrature *quad = NULL);
inline double f(const int iq) const;
double f(const int iq) const;
protected:
DOFVectorBase<double> *rhoDV;
mtl::dense_vector<double> rho;
......
extensions/SignedDistFunctors.h
View file @ 6ff14ecb
......@@ -457,10 +457,15 @@ public:
double operator()(const WorldVector<double>& x) const
{
double result=0.0;
PerturbedCircleRadius perturbedCircle(radius, 1.0, 3, rotation);
PerturbedCircleRadius perturbedCircle(radius, 0.7, 3, rotation);
result = signedDist2D(x, midPoint, &perturbedCircle, 1.e-6);
return -result;
};
}
void setRotation(double rotation_)
{
rotation = rotation_;
}
private:
......
extensions/SingularDirichletBC.h
View file @ 6ff14ecb
......@@ -185,88 +185,30 @@ namespace details {
* The 'nr' and 'meshIndicator' describe one boundary of the mesh and to each DegreeOfFreedom
* the 'periodicMap' associates a second DegreeOfFreedom.
**/
inline void getPeriodicAssociation(const FiniteElemSpace* feSpace,
BoundaryType nr,
AbstractFunction<bool, WorldVector<double> >* meshIndicator,
AbstractFunction<WorldVector<double>, WorldVector<double> >* periodicMap,
std::vector<std::pair<DegreeOfFreedom, DegreeOfFreedom> > &association)
{
std::set<DegreeOfFreedom> indices;
details::getBoundaryIndices(feSpace, nr, meshIndicator, indices);
association.clear();
DOFVector<WorldVector<double> > coords(feSpace, "coords");
feSpace->getMesh()->getDofIndexCoords(coords);
std::set<DegreeOfFreedom>::iterator it;
for (it = indices.begin(); it != indices.end(); it++)
{
DegreeOfFreedom idx2;
WorldVector<double> p;
p = (*periodicMap)(coords[*it]);
TEST_EXIT(coords.getDofIdxAtPoint(p,idx2))("periodic association not found!\n");
association.push_back(std::make_pair(*it, idx2));
}
}
inline bool getDofIdxAtPoint(const FiniteElemSpace* feSpace,
WorldVector<double> &p,
DegreeOfFreedom &idx,
ElInfo *oldElInfo = NULL,
bool useOldElInfo = false)
{
FUNCNAME("DOFVector::getDofIdxAtPoint()");
Mesh *mesh = feSpace->getMesh();
const BasisFunction *basFcts = feSpace->getBasisFcts();
int dim = mesh->getDim();
int numBasFcts = basFcts->getNumber();
std::vector<DegreeOfFreedom> localIndices(numBasFcts);
DimVec<double> lambda(dim, NO_INIT);
ElInfo *elInfo = mesh->createNewElInfo();
idx = 0;
bool inside = false;
if (oldElInfo && useOldElInfo && oldElInfo->getMacroElement()) {
inside = mesh->findElInfoAtPoint(p, elInfo, lambda, oldElInfo->getMacroElement(), NULL, NULL);
delete oldElInfo;
} else {
inside = mesh->findElInfoAtPoint(p, elInfo, lambda, NULL, NULL, NULL);
}
if (oldElInfo)
oldElInfo = elInfo;
if (!inside)
return false;
basFcts->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(), localIndices);
WorldVector<double> coord;
int minIdx = -1;
double minDist = 1.e15;
for (int i = 0; i < numBasFcts; i++) {
elInfo->coordToWorld(*(basFcts->getCoords(i)), coord);
WorldVector<double> dist = coord - p;
double newDist = norm(dist);
if (newDist < minDist) {
minIdx = i;
minDist = newDist;
}
}
if (minIdx >= 0)
idx = localIndices[minIdx];
if(!oldElInfo) delete elInfo;
return inside;
}
// inline void getPeriodicAssociation(const FiniteElemSpace* feSpace,
// BoundaryType nr,
// AbstractFunction<bool, WorldVector<double> >* meshIndicator,
// AbstractFunction<WorldVector<double>, WorldVector<double> >* periodicMap,
// std::vector<std::pair<DegreeOfFreedom, DegreeOfFreedom> > &association)
// {
// std::set<DegreeOfFreedom> indices;
// details::getBoundaryIndices(feSpace, nr, meshIndicator, indices);
//
// association.clear();
//
// DOFVector<WorldVector<double> > coords(feSpace, "coords");
// feSpace->getMesh()->getDofIndexCoords(coords);
// std::set<DegreeOfFreedom>::iterator it;
// for (it = indices.begin(); it != indices.end(); it++)
// {
// DegreeOfFreedom idx2;
// WorldVector<double> p;
// p = (*periodicMap)(coords[*it]);
// TEST_EXIT(coords.getDofIdxAtPoint(p,idx2))("periodic association not found!\n");
//
// association.push_back(std::make_pair(min(*it, idx2), max(*it, idx2)));
// }
// }
/**
* get al list of pairs that describes wich DegreeOfFreedom are assiciated in the given mesh
......@@ -274,28 +216,51 @@ namespace details {
* The 'meshIndicator' describes one boundary of the mesh and to each DegreeOfFreedom
* the 'periodicMap' associates a second DegreeOfFreedom.
**/
// inline void getPeriodicAssociation(const FiniteElemSpace* feSpace,
// AbstractFunction<bool, WorldVector<double> >* meshIndicator,
// AbstractFunction<WorldVector<double>, WorldVector<double> >* periodicMap,
// std::vector<std::pair<DegreeOfFreedom, DegreeOfFreedom> > &association)
// {
// std::set<DegreeOfFreedom> indices;
// details::getBoundaryIndices(feSpace, meshIndicator, indices);
// MSG_DBG("Nr of indices: %d\n", indices.size());
//
// association.clear();
//
// DOFVector<WorldVector<double> > coords(feSpace, "coords");
// feSpace->getMesh()->getDofIndexCoords(coords);
// std::set<DegreeOfFreedom>::iterator it;
// for (it = indices.begin(); it != indices.end(); it++)
// {
// DegreeOfFreedom idx2;
// WorldVector<double> p;
// p = (*periodicMap)(coords[*it]);
// TEST_EXIT(coords.getDofIdxAtPoint(p,idx2))("periodic association not found!\n");
//
// association.push_back(std::make_pair(min(*it, idx2), max(*it, idx2));
// }
//
// MSG_DBG("Nr of associations: %d\n", association.size());
// }
inline void getPeriodicAssociation(const FiniteElemSpace* feSpace,
AbstractFunction<bool, WorldVector<double> >* meshIndicator,
std::set<DegreeOfFreedom> indices,
AbstractFunction<WorldVector<double>, WorldVector<double> >* periodicMap,
std::vector<std::pair<DegreeOfFreedom, DegreeOfFreedom> > &association)
{
std::set<DegreeOfFreedom> indices;
details::getBoundaryIndices(feSpace, meshIndicator, indices);
MSG_DBG("Nr of indices: %d\n", indices.size());
association.clear();
DOFVector<WorldVector<double> > coords(feSpace, "coords");
feSpace->getMesh()->getDofIndexCoords(coords);
std::set<DegreeOfFreedom>::iterator it;
DegreeOfFreedom idx_;
WorldVector<double> p;
for (it = indices.begin(); it != indices.end(); it++)
{
DegreeOfFreedom idx2;
WorldVector<double> p;
p = (*periodicMap)(coords[*it]);
TEST_EXIT(coords.getDofIdxAtPoint(p,idx2))("periodic association not found!\n");
association.push_back(std::make_pair(*it, idx2));
TEST_EXIT(coords.getDofIdxAtPoint(p,idx_))("periodic association not found!\n");
association.push_back(std::make_pair(std::min(*it, idx_), std::max(*it, idx_)));
}
MSG_DBG("Nr of associations: %d\n", association.size());
......@@ -547,11 +512,14 @@ struct PeriodicBcData {
void addToList(const FiniteElemSpace *feSpace, std::vector<ManualPeriodicBC> &list)
{
std::vector<std::pair<DegreeOfFreedom, DegreeOfFreedom> > associations;
std::set<DegreeOfFreedom> indices;
if (nr == 0)
details::getPeriodicAssociation(feSpace, meshIndicator, periodicMap, associations);
details::getBoundaryIndices(feSpace, meshIndicator, indices);
else
details::getPeriodicAssociation(feSpace, nr, meshIndicator, periodicMap, associations);
details::getBoundaryIndices(feSpace, nr, meshIndicator, indices);
std::vector<std::pair<DegreeOfFreedom, DegreeOfFreedom> > associations;
details::getPeriodicAssociation(feSpace, indices, periodicMap, associations);
list.push_back(ManualPeriodicBC(row, associations));
}
int row;
......
extensions/base_problems/NavierStokesPhase_TaylorHood.cc
View file @ 6ff14ecb
......@@ -10,7 +10,8 @@ NavierStokesPhase_TaylorHood::NavierStokesPhase_TaylorHood(const std::string &na
beta(1.0),
epsilon(1.e-1),
alpha(3.0),
fileWriterPhase(NULL)
fileWriterPhase(NULL),
phaseOld(NULL)
{
for (int i = 0; i < bcDOF.getSize(); ++i)
bcDOF[i] = NULL;
......@@ -24,8 +25,14 @@ NavierStokesPhase_TaylorHood::NavierStokesPhase_TaylorHood(const std::string &na
NavierStokesPhase_TaylorHood::~NavierStokesPhase_TaylorHood()
{
if (fileWriterPhase != NULL)
if (fileWriterPhase != NULL) {
delete fileWriterPhase;
fileWriterPhase = NULL;
}
if (phaseOld != NULL) {
delete phaseOld;
phaseOld = NULL;
}
}
......@@ -44,29 +51,13 @@ void NavierStokesPhase_TaylorHood::transferInitialSolution(AdaptInfo *adaptInfo)
fileWriterPhase = new FileWriter(name + "->phase->output", getFeSpace()->getMesh(), phase);
for (int i = 0; i < dow; i++)
transformDOF(prob->getSolution()->getDOFVector(0), phase, prob->getSolution()->getDOFVector(0), new AMDiS::Mult<double>);
transformDOF(prob->getSolution()->getDOFVector(i), phase, prob->getSolution()->getDOFVector(i), new AMDiS::Mult<double>);
super::transferInitialSolution(adaptInfo);
phaseOld->interpol(phase);
}
void NavierStokesPhase_TaylorHood::closeTimestep(AdaptInfo *adaptInfo)
{
super::closeTimestep(adaptInfo);
phaseOld->interpol(phase);
writeFiles(adaptInfo, false);
}
void NavierStokesPhase_TaylorHood::writeFiles(AdaptInfo *adaptInfo, bool force)
{ FUNCNAME("NavierStokesPhase_TaylorHood::closeTimestep()");
super::writeFiles(adaptInfo, force);
fileWriterPhase->writeFiles(adaptInfo, false);
}
void NavierStokesPhase_TaylorHood::fillOperators()
{ FUNCNAME("NavierStokesPhase_TaylorHood::fillOperators()");
......@@ -82,40 +73,34 @@ void NavierStokesPhase_TaylorHood::fillOperators()
Operator *opTime = new Operator(prob->getFeSpace(i), prob->getFeSpace(i));
opTime->addTerm(new Phase_ZOT(phase, density));
prob->addMatrixOperator(*opTime, i, i, getInvTau());
/// < (1/tau)*u_i^old , psi >
Operator *opTimeOld = new Operator(prob->getFeSpace(i), prob->getFeSpace(i));
opTimeOld->addTerm(new Phase_ZOT(phaseOld, density));
opTimeOld->setUhOld(getOldSolution(i));
opTimeOld->addTerm(new FactorPhase_ZOT(phaseOld, getOldSolution(i), density));
prob->addVectorOperator(*opTimeOld, i, getInvTau());
/// < u^old*grad(u_i^old) , psi >
Operator *opUGradU0 = new Operator(prob->getFeSpace(i),prob->getFeSpace(i));
opUGradU0->addTerm(new WorldVecPhase_FOT(phase,vel, -density), GRD_PHI);
Operator *opUGradU0 = new Operator(getFeSpace(i), getFeSpace(i));
opUGradU0->addTerm(new WorldVecPhase_FOT(phase, vel, density), GRD_PHI);
opUGradU0->setUhOld(prob->getSolution()->getDOFVector(i));
if (nonLinTerm == 0) {
prob->addVectorOperator(*opUGradU0, i);
} else {
prob->addVectorOperator(*opUGradU0, i, &theta1);
}
if (nonLinTerm == 1) {
/// < u'*grad(u_i^old) , psi >
for (unsigned j = 0; j < dow; ++j) {
Operator *opUGradU1 = new Operator(prob->getFeSpace(i),prob->getFeSpace(i));
/// < u*grad(u_i^old) , psi >
for (unsigned k = 0; k < dow; ++k) {
Operator *opUGradU1 = new Operator(getFeSpace(i), getFeSpace(k));
opUGradU1->addTerm(new VecAndPartialDerivative_ZOT(
phase,
prob->getSolution()->getDOFVector(i), j, density));
prob->addMatrixOperator(*opUGradU1, i, j, &theta);
prob->getSolution()->getDOFVector(i), k, density));
prob->addMatrixOperator(*opUGradU1, i, k);
}
} else if (nonLinTerm == 2) {
/// < u^old*grad(u'_i) , psi >
for(unsigned j = 0; j < dow; ++j) {
Operator *opUGradU2 = new Operator(prob->getFeSpace(i),prob->getFeSpace(i));
/// < u^old*grad(u_i) , psi >
for(unsigned k = 0; k < dow; ++k) {
Operator *opUGradU2 = new Operator(getFeSpace(i),getFeSpace(i));
opUGradU2->addTerm(new Vec2ProductPartial_FOT(
phase,
prob->getSolution()->getDOFVector(j), j, density), GRD_PHI);
prob->addMatrixOperator(*opUGradU2, i, i, &theta);
}
vel[k], k, density), GRD_PHI);
prob->addMatrixOperator(*opUGradU2, i, i);
}
/// Diffusion-Operator (including Stress-Tensor for space-dependent viscosity
......@@ -147,6 +132,10 @@ void NavierStokesPhase_TaylorHood::fillOperators()
opDivU2->addTerm(new PartialDerivative_ZOT(phase,i));
prob->addMatrixOperator(*opDivU2, dow, i);
}
Operator *opNull = new Operator(getFeSpace(dow), getFeSpace(dow));
opNull->addTerm(new Simple_ZOT(0.0));
prob->addMatrixOperator(*opNull, dow, dow);
}
......@@ -226,3 +215,18 @@ void NavierStokesPhase_TaylorHood::fillBoundaryConditions()
TEST_EXIT(!bcCond)("implicit boundary condition specified twice, by bcFct and bcDOF. This is not allowed!\n");
}
}
void NavierStokesPhase_TaylorHood::closeTimestep(AdaptInfo *adaptInfo)
{ FUNCNAME("NavierStokesPhase_TaylorHood::closeTimestep()");
phaseOld->interpol(phase);
super::closeTimestep(adaptInfo);
}
void NavierStokesPhase_TaylorHood::writeFiles(AdaptInfo *adaptInfo, bool force)
{ FUNCNAME("NavierStokesPhase_TaylorHood::writeFiles()");
super::writeFiles(adaptInfo, force);
fileWriterPhase->writeFiles(adaptInfo, false);
}
extensions/base_problems/NavierStokes_TaylorHood.cc
View file @ 6ff14ecb
......@@ -5,33 +5,17 @@ using namespace AMDiS;
NavierStokes_TaylorHood::NavierStokes_TaylorHood(const std::string &name_) :
super(name_),
forceDBC(false),
initialVelocityIsSet(false),
laplaceType(0),
nonLinTerm(2),
oldTimestep(0.0),
viscosity(1.0),
density(1.0),
c(0.0),
theta(0.5),
theta1(0.5),
minusTheta1(-0.5),
theta(0.5+1.e-2),
theta1(1-theta),
minusTheta1(-theta1),
velocity(NULL),
fileWriter(NULL),
initialX(NULL),
initialY(NULL)
fileWriter(NULL)
{
// force the homogeniouse dirichlet BC if combination of dirichlet and neumann BC
// are set and AMDiS can not handle this combination automatically
Initfile::get(name + "->force dirichlet bc", forceDBC);
if (forceDBC) {
int numDirichletPoints = 0;
Initfile::get("number of dirichlet points", numDirichletPoints);
dirichletPoints.resize(numDirichletPoints);
for (int i = 0; i < numDirichletPoints; i++)
Initfile::get("dirichlet point[" + Helpers::toString(i) + "]",dirichletPoints[i]);
}
// parameters for navier-stokes
Initfile::get(name + "->viscosity", viscosity);
Initfile::get(name + "->density", density);
......@@ -54,7 +38,7 @@ NavierStokes_TaylorHood::NavierStokes_TaylorHood(const std::string &name_) :
} else
dimension.set(1.0);
oldSolution.resize(dow);
for (size_t i = 0; i < dow; i++)
oldSolution[i] = NULL;
}
......@@ -63,11 +47,6 @@ NavierStokes_TaylorHood::NavierStokes_TaylorHood(const std::string &name_) :
NavierStokes_TaylorHood::~NavierStokes_TaylorHood()
{ FUNCNAME("NavierStokes_TaylorHood::~NavierStokes_TaylorHood()");
if (initialVelocityIsSet) {
delete initialX;
delete initialY;
}
if (velocity != NULL) {
delete velocity;
velocity = NULL;
......@@ -101,6 +80,8 @@ void NavierStokes_TaylorHood::initData()
void NavierStokes_TaylorHood::solveInitialProblem(AdaptInfo *adaptInfo)
{ FUNCNAME("NavierStokes_TaylorHood::solveInitialProblem()");
double c = 1.0;
int initialVelocity = 0, flowDirection = 1;
Initfile::get(name + "->initial velocity", initialVelocity);
Initfile::get(name + "->initial velocity->value", c);
......@@ -123,10 +104,6 @@ void NavierStokes_TaylorHood::solveInitialProblem(AdaptInfo *adaptInfo)
} else
throw(std::runtime_error("Unknown initial velocity."));
initialX = initialFcts[0];
initialY = initialFcts[1];
initialVelocityIsSet = true;
MSG("solve initial problem...\n");
for (unsigned i = 0; i < dow; ++i) {
prob->getSolution()->getDOFVector(i)->interpol(initialFcts[i]);
......@@ -138,74 +115,64 @@ void NavierStokes_TaylorHood::solveInitialProblem(AdaptInfo *adaptInfo)
void NavierStokes_TaylorHood::transferInitialSolution(AdaptInfo *adaptInfo)
{ FUNCNAME("NavierStokes_TaylorHood::transferInitialSolution()");
calcVelocity();
for (int i = 0; i < dow; i++)
prob->setExactSolution(prob->getSolution()->getDOFVector(i), i);
calcVelocity();
for (size_t i = 0; i < dow; i++)
oldSolution[i]->copy(*prob->getSolution()->getDOFVector(i));
fileWriter->writeFiles(adaptInfo, false);
writeFiles(adaptInfo, false);
// initial parameters for detecting mesh changes
oldMeshChangeIdx= getMesh()->getChangeIndex();
super::transferInitialSolution(adaptInfo);
}
void NavierStokes_TaylorHood::fillOperators()
{ FUNCNAME("NavierStokes_TaylorHood::fillOperators()");
WorldVector<DOFVector<double>* > vel;
for (unsigned k = 0; k < dow; k++)
vel[k] = prob->getSolution()->getDOFVector(k);
// fill operators for prob
for (unsigned i = 0; i < dow; ++i) {
/// < (1/tau)*u'_i , psi >
Operator *opTime = new Operator(getFeSpace(i), getFeSpace(i));
opTime->addTerm(new Simple_ZOT(density));
prob->addMatrixOperator(*opTime, i, i, getInvTau(), getInvTau());
/// < (1/tau)*u_i^old , psi >
opTime->setUhOld(prob->getSolution()->getDOFVector(i));
prob->addVectorOperator(*opTime, i, getInvTau(), getInvTau());
Operator *opTimeOld = new Operator(getFeSpace(i), getFeSpace(i));
opTimeOld->addTerm(new Phase_ZOT(getOldSolution(i), density));
prob->addVectorOperator(*opTimeOld, i, getInvTau(), getInvTau());
/// < u^old*grad(u_i^old) , psi >
Operator *opUGradU0 = new Operator(getFeSpace(i), getFeSpace(i));
opUGradU0->addTerm(new WorldVector_FOT(velocity, -density), GRD_PHI);
opUGradU0->addTerm(new WorldVec_FOT(vel, density), GRD_PHI);
opUGradU0->setUhOld(prob->getSolution()->getDOFVector(i));
if (nonLinTerm == 0) {
prob->addVectorOperator(*opUGradU0, i);
} else if (abs(theta1) > DBL_TOL) {
prob->addVectorOperator(*opUGradU0, i, &theta1, &theta1);
}
if (nonLinTerm == 1) {
/// < u'*grad(u_i^old) , psi >
for (unsigned j = 0; j < dow; ++j) {
Operator *opUGradU1 = new Operator(getFeSpace(i),getFeSpace(i));
/// < u*grad(u_i^old) , psi >
for (unsigned k = 0; k < dow; ++k) {
Operator *opUGradU1 = new Operator(getFeSpace(i), getFeSpace(k));
opUGradU1->addTerm(new PartialDerivative_ZOT(
prob->getSolution()->getDOFVector(i), j, density));
prob->addMatrixOperator(*opUGradU1, i,