In the constructor pointers to the two problems are assigned to the private members \verb+problem1+ and \verb+problem2+. Note that the pointers point to the interface \verb+ProblemStatBase+ and not to\\\verb+ProblemScal+. This leads to a more general implementation. If e.g. two vector valued problems should be coupled in the future, we could use our iteration class without modifications.
...
...
@@ -60,14 +60,14 @@ Now, we implement the needed interface methods:
MSG("\n");
MSG("begin of iteration %d\n", adaptInfo->getSpaceIteration()+1);
The \verb+toDo+ flag is used by the adaptation loop to determine which parts of the iteration should be performed. The first iteration is always an iteration without mesh adaptation (see Figure \ref{f:stationary loop}). So we start our iteration by marking and adapting the mesh. The mesh and its adaptation is managed by the first problem. So we call \verb+markElements+ and \verb+refineMesh+ of \verb+problem1+. Note that no mesh coarsenings have to be performed in our example. Afterwards, \verb+problem1+ assembles its system of equations by \verb+buildAfterCoarsen+. Assemblage and mesh adaptation are nested operations in AMDiS (\verb+buildBeforeRefine+, \verb+refineMesh+, \verb+buildBeforeCoarsen+,\\\verb+coarsenMesh+,\verb+buildAfterCoarsen+). Here, we implement a simplified version.
...
...
@@ -111,16 +111,16 @@ Now, the access to the coupled problems is implemented and the member variables
int getNumProblems()
{
return 2;
};
}
ProblemStatBase *getProblem(int number = 0)
{
FUNCNAME("CoupledIteration::getProblem()");
if(number == 0) return problem1;
if(number == 1) return problem2;
if(number == 0) return problem1;
if(number == 1) return problem2;
ERROR_EXIT("invalid problem number\n");
return NULL;
};
}
private:
ProblemStatBase *problem1;
...
...
@@ -136,15 +136,12 @@ The next class, \verb+Identity+, implements the identity $I(x)=x$ for \verb+doub
class Identity : public AbstractFunction<double, double>
@@ -136,10 +126,10 @@ The next steps are the creation of the adaptation loop and the corresponding \ve
\begin{lstlisting}{}
// === create adapt info ===
AdaptInfo *adaptInfo = NEW AdaptInfo("ellipt->adapt", 1);
AdaptInfo *adaptInfo = new AdaptInfo("ellipt->adapt", 1);
// === create adapt ===
AdaptStationary *adapt = NEW AdaptStationary("ellipt->adapt", &ellipt,
AdaptStationary *adapt = new AdaptStationary("ellipt->adapt", &ellipt,
adaptInfo);
\end{lstlisting}
...
...
@@ -150,7 +140,7 @@ Now, we define boundary conditions:
\begin{lstlisting}{}
// ===== add boundary conditions =====
ellipt.addDirichletBC(1, NEW G);
ellipt.addDirichletBC(1, new G);
\end{lstlisting}
We have one Dirichlet boundary condition associated with identifier $1$. All nodes belonging to this boundary are set to the value of function \verb+G+ at the corresponding coordinates. In the macro file (see Section \ref{s:ellipt macro}) the Dirichlet boundary is marked with identifier $1$, too. So the nodes can be uniquely determined.
...
...
@@ -160,13 +150,13 @@ The operators now are defined as follows:
@@ -87,8 +85,6 @@ The first lines of class \verb+Heat+ are:
class Heat : public ProblemInstatScal
{
public:
MEMORY_MANAGED(Heat);
Heat(ProblemScal *heatSpace)
: ProblemInstatScal("heat", heatSpace)
{
...
...
@@ -96,7 +92,7 @@ public:
GET_PARAMETER(0, name + "->theta", "%f", &theta);
TEST_EXIT(theta >= 0)("theta not set!\n");
theta1 = theta - 1;
};
}
\end{lstlisting}
The argument \verb+heatSpace+ is a pointer to the stationary problem which is solved each timestep. It is directly handed to the base class constructor of \verb+ProblemInstatScal+. In the body of the constructor, $\theta$ is read from the parameter file and stored in a member variable. The member variable \verb+theta1+ stores the value of $\theta-1$. A pointer to this value is used later as factor in the $\theta$-scheme.
...
...
@@ -104,18 +100,20 @@ The argument \verb+heatSpace+ is a pointer to the stationary problem which is so
The next lines show the implementation of the time interface.
\begin{lstlisting}{}
void setTime(AdaptInfo *adaptInfo) {
void setTime(AdaptInfo *adaptInfo)
{
rhsTime =
adaptInfo->getTime()-
(1-theta)*adaptInfo->getTimestep();
boundaryTime = adaptInfo->getTime();
tau1 = 1.0 / adaptInfo->getTimestep();
};
}
void closeTimestep(AdaptInfo *adaptInfo) {
void closeTimestep(AdaptInfo *adaptInfo)
{
ProblemInstatScal::closeTimestep(adaptInfo);
WAIT;
};
}
\end{lstlisting}
The method \verb+setTime+ is called by the adaptation loop to inform the problem about the current time. The right hand side function $f$ will be evaluated at $t^{old}+\theta\tau= t^{new}-(1-\theta)\tau$, the Dirichlet boundary function $g$ at $t^{new}$. $t^{new}$ is the current time, $\tau$ is the current timestep, both set by the adaptation loop and stored in \verb+adaptInfo+. \verb+tau1+ stores the value of $\frac{1}{\tau}$, which is used later as factor for the zero order time discretization terms.
...
...
@@ -128,7 +126,7 @@ Now, the implementation of the \verb+ProblemStatBase+ interface begins. As menti
@@ -138,7 +136,7 @@ Now, the implementation of the \verb+ProblemStatBase+ interface begins. As menti
0, &errMax, false);
adaptInfo->setEstSum(errSum, 0);
adaptInfo->setEstMax(errMax, 0);
};
}
\end{lstlisting}
Here, only the solve and the estimate step are overloaded. For the other steps, there are empty default implementations in \verb+ProblemInstatScal+. Since the mesh is not adapted in the initial problem, the initial adaptation loop will stop after one iteration. In the solve step, the exact solution is interpolated on the macro mesh and stored in the solution vector of the stationary problem. In the estimate step, the L2 error is computed. The maximal element error and the sum over all element errors are stored in \verb+adaptInfo+. To make the exact solution known to the problem, we need a setting function:
...
...
@@ -181,7 +179,7 @@ int main(int argc, char** argv)
Parameters::init(false, argv[1]);
// ===== create and init stationary problem =====
ProblemScal *heatSpace = NEW ProblemScal("heat->space");
ProblemScal *heatSpace = new ProblemScal("heat->space");
heatSpace->initialize(INIT_ALL);
// ===== create instationary problem =====
...
...
@@ -195,13 +193,13 @@ The next step is the creation of the needed \verb+AdaptInfo+ objects and of the
\begin{lstlisting}{}
// create adapt info for heat
AdaptInfo *adaptInfo = NEW AdaptInfo("heat->adapt");
AdaptInfo *adaptInfo = new AdaptInfo("heat->adapt");
// create initial adapt info
AdaptInfo *adaptInfoInitial = NEW AdaptInfo("heat->initial->adapt");
AdaptInfo *adaptInfoInitial = new AdaptInfo("heat->initial->adapt");
// create instationary adapt
AdaptInstationary *adaptInstat = NEW AdaptInstationary("heat->adapt",
AdaptInstationary *adaptInstat = new AdaptInstationary("heat->adapt",
heatSpace,
adaptInfo,
heat,
...
...
@@ -214,7 +212,7 @@ The definitions of functions $f$ and $g$ are:
@@ -109,14 +103,11 @@ The class \verb+X3+ implements the function $u^3(x)$ used within the Newton step
class X3 : public AbstractFunction<double, double>
{
public:
MEMORY_MANAGED(X3);
X3() : AbstractFunction<double, double>(3) {}
X3() : AbstractFunction<double, double>(3) {};
/** \brief
* Implementation of AbstractFunction::operator().
*/
double operator()(const double& x) const {
/// Implementation of AbstractFunction::operator().
double operator()(const double& x) const
{
return x * x * x;
}
};
...
...
@@ -158,7 +149,7 @@ public:
&newtonMaxIter);
solution = problemNonlin->getSolution();
correction = newtonStep->getCorrection();
};
}
\end{lstlisting}
In the constructor, pointers to the nonlinear problem and to the Newton-step object are stored to the class members \verb+problemNonlin+ and \verb+newtonStep+. Furthermore, the parameters \verb+newtonTolerance+ and \verb+newtonMaxIter+ are initialized, and pointers to the nonlinear solution and to the correction vector are stored.
...
...
@@ -178,13 +169,13 @@ The following methods define one iteration in the adaptation loop.
@@ -282,54 +273,60 @@ In the destructor, the allocated memory is freed.
\begin{lstlisting}{}
~Nonlin()
{
DELETE correction;
DELETE dirichletZero;
DELETE dirichletG;
};
delete correction;
delete dirichletZero;
delete dirichletG;
}
\end{lstlisting}
Now, we implement the Newton step functionality. First, in \verb+initNewtonStep+, we fill the solution vector with boundary values. This will not be done automatically because we let the \verb+solution_+ pointer point to \verb+correction+. The address of \verb+solution_+ is stored in \verb+tmp+. After the Newton method is finished, the \verb+solution_+ pointer is reset to its original value in \verb+exitNewtonStep+.
The implementation of \verb+assembleNewtonStep+ and \verb+solveNewtonStep+ just delegates the calls to the base class implementations in \verb+ProblemScal+.
\begin{lstlisting}{}
void assembleNewtonStep(AdaptInfo *adaptInfo, Flag flag) {
void assembleNewtonStep(AdaptInfo *adaptInfo, Flag flag)
{
ProblemScal::buildAfterCoarsen(adaptInfo, flag);
};
}
void solveNewtonStep(AdaptInfo *adaptInfo) {
void solveNewtonStep(AdaptInfo *adaptInfo)
{
ProblemScal::solve(adaptInfo);
};
}
\end{lstlisting}
Finally, the \verb+getCorrection+ method is implemented and private class members are defined.
