Commit 802ccc5d authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files

* Operator VecGrad_FOT added

* AdaptInfo::reset() does not set timestep to 0.0, instead it reads the parameter from init file again
* ParaReal-stuff changed
parent 4ce5b65b
......@@ -213,6 +213,8 @@ namespace AMDiS {
timestepNumber = 0;
solverIterations = 0;
solverResidual = 0.0;
GET_PARAMETER(0, name + "->timestep", "%f", &timestep);
};
/** \brief
......
......@@ -407,8 +407,6 @@ namespace AMDiS {
coordsAtQPs = subAssembler->getCoordsAtQPs(elInfo, quad);
}
// bis hierhin erstmal
void VecAndGradAtQP_ZOT::initElement(const ElInfo* elInfo,
SubAssembler* subAssembler,
Quadrature *quad)
......@@ -2157,4 +2155,34 @@ namespace AMDiS {
}
}
void VecGrad_FOT::initElement(const ElInfo* elInfo,
SubAssembler* subAssembler,
Quadrature *quad)
{
vecAtQPs = subAssembler->getVectorAtQPs(vec1, elInfo, quad);
gradAtQPs = subAssembler->getGradientsAtQPs(vec2, elInfo, quad);
}
void VecGrad_FOT::getLb(const ElInfo *elInfo, int numPoints, VectorOfFixVecs<DimVec<double> >& Lb) const {
const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
for (int iq = 0; iq < numPoints; iq++) {
lb(Lambda, (*vecFct)(vecAtQPs[iq], gradAtQPs[iq]), Lb[iq], 1.0);
}
}
void VecGrad_FOT::eval(int numPoints,
const double *uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
double *result,
double factor) const
{
if (grdUhAtQP) {
for (int iq = 0; iq < numPoints; iq++) {
WorldVector<double> b = (*vecFct)(vecAtQPs[iq], gradAtQPs[iq]);
result[iq] += b * grdUhAtQP[iq] * factor;
}
}
}
}
......@@ -1900,6 +1900,67 @@ namespace AMDiS {
AbstractFunction<WorldVector<double>, WorldVector<double> > *g;
};
/**
* \ingroup Assembler
*
* \brief
* First order term: \f$ b(v(\vec{x}), \nabla w(\vec{x})) \cdot \nabla u(\vec{x})\f$.
*/
class VecGrad_FOT : public FirstOrderTerm
{
public:
/** \brief
* Constructor.
*/
VecGrad_FOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2,
BinaryAbstractFunction<WorldVector<double>, double,
WorldVector<double> > *vecFct_)
: FirstOrderTerm(vecFct_->getDegree()), vec1(dv1), vec2(dv2), vecFct(vecFct_)
{};
/** \brief
* Implementation of \ref OperatorTerm::initElement().
*/
void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
Quadrature *quad = NULL);
/** \brief
* Implements FirstOrderTerm::getLb().
*/
void getLb(const ElInfo *elInfo, int qPoint,
VectorOfFixVecs<DimVec<double> >& Lb) const;
/** \brief
* Implements FirstOrderTerm::eval().
*/
void eval(int numPoints,
const double *uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
double *result,
double factor) const;
protected:
/** \brief
* DOFVector to be evaluated at quadrature points.
*/
DOFVectorBase<double>* vec1;
DOFVectorBase<double>* vec2;
/** \brief
* Vector v at quadrature points.
*/
double *vecAtQPs;
WorldVector<double> *gradAtQPs;
/** \brief
* Function for b.
*/
BinaryAbstractFunction<WorldVector<double>, double, WorldVector<double> > *vecFct;
};
// ============================================================================
class General_FOT : public FirstOrderTerm
......
......@@ -12,6 +12,7 @@ namespace AMDiS {
adaptInfo_->setTimestep(coarseTimestep);
pararealProb->setStoreSolutions(false, true);
pararealProb->setStoreInitSolution(true);
pararealProb->setInitSolutionVec(NULL);
AdaptInstationary::adapt();
// ParaReal iterations
......
......@@ -52,7 +52,7 @@ namespace AMDiS {
TEST_EXIT(fineTimestep > 0.0)("fineTimestep must be greater than zero!\n");
TEST_EXIT(coarseTimestep > 0.0)("coarseTimestep must be greater than zero!\n");
TEST_EXIT(fineTimestep < coarseTimestep)("fineTimestep must be smaller than coarseTimestep!\n");
TEST_EXIT(pararealIter > 0)("ParaReal Iterations must be greater than zero!\n");
TEST_EXIT(pararealIter >= 0)("ParaReal Iterations must be >= 0!\n");
}
/** \brief
......
......@@ -165,7 +165,11 @@ namespace AMDiS {
}
virtual void solveInitialProblem(AdaptInfo *adaptInfo) {
ProblemInstatScal::solveInitialProblem(adaptInfo);
if (initSolutionVec == NULL) {
ProblemInstatScal::solveInitialProblem(adaptInfo);
} else {
*(problemStat->getSolution()) = *initSolutionVec;
}
if (storeInitSolution) {
storeSolution(problemStat->getSolution());
......@@ -194,7 +198,11 @@ namespace AMDiS {
}
virtual void solveInitialProblem(AdaptInfo *adaptInfo) {
ProblemInstatVec::solveInitialProblem(adaptInfo);
if (initSolutionVec == NULL) {
ProblemInstatVec::solveInitialProblem(adaptInfo);
} else {
*(problemStat->getSolution()) = *initSolutionVec;
}
if (storeInitSolution) {
storeSolution(problemStat->getSolution());
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment