In this example, we solve the same problem as in Section \ref{s:vecellipt}, but here we treat the two equations as two coupled problems. The main difference is that the equations now aren't assembled into the same large system of equations, but into two separated systems of equations, that have to be solved separately. We define the two problems
\begin{eqnarray}
-\Delta u &=& f ~~~\mbox{in }\Omega\subset\mathbb{R}^{dim}\\
u &=& g ~~~ \mbox{on }\partial\Omega
\end{eqnarray}
and
\begin{equation}
v = u.
\end{equation}
We first solve the first problem and then use its solution to solve the second problem. This happens in every iteration of the adaptation loop. Both problems should use the same mesh. Mesh adaptation is done by the first problem. So one iteration now looks like illustrated in Figure \ref{f:coupled iteration}.
In the previous examples, the iteration was implemented always by the corresponding problem class. In this example, we want to couple two problems within one iteration, so the default implementation can't be used. For that reason, we define our own coupled iteration class\\\verb+MyCoupledIteration+ which implements the \verb+ProblemIterationInterface+:
\begin{lstlisting}{}
class MyCoupledIteration : public ProblemIterationInterface
{
public:
MyCoupledIteration(ProblemStatBase *prob1,
ProblemStatBase *prob2)
: problem1(prob1),
problem2(prob2)
{}
\end{lstlisting}
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.
MSG("end of iteration %d\n", adaptInfo->getSpaceIteration()+1);
MSG("=============================\n");
}
\end{lstlisting}
These two functions are called at the beginning and at the end of each iteration. Here, we only prompt some output.
The method \verb+oneIteration+ is the crucial part of our implementation:
\begin{lstlisting}{}
Flag oneIteration(AdaptInfo *adaptInfo, Flag toDo = FULL_ITERATION)
{
Flag flag, markFlag;
if (toDo.isSet(MARK)) markFlag = problem1->markElements(adaptInfo);
if (toDo.isSet(ADAPT) && markFlag.isSet(MESH_REFINED))
flag = problem1->refineMesh(adaptInfo);
if (toDo.isSet(BUILD)) problem1->buildAfterCoarsen(adaptInfo, markFlag);
if (toDo.isSet(SOLVE)) problem1->solve(adaptInfo);
if (toDo.isSet(BUILD)) problem2->buildAfterCoarsen(adaptInfo, markFlag);
if (toDo.isSet(SOLVE)) problem2->solve(adaptInfo);
if (toDo.isSet(ESTIMATE)) problem1->estimate(adaptInfo);
return flag;
}
\end{lstlisting}
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.
After \verb+problem1+ has solved its system of equations, \verb+problem2+ can assemble and solve its equations system using the solution of the first problem as right hand side. In the method \verb+oneIteration+, only the order of method calls is determined. The dependency to the solution of the first problem is created later when the operator for the right hand side of \verb+problem2+ is created.
After also the second problem computed its solution, \verb+problem1+ does the error estimation (remember: mesh adaptation is managed by \verb+problem1+).
Now, the access to the coupled problems is implemented and the member variables are defined:
\begin{lstlisting}{}
int getNumProblems()
{
return 2;
}
ProblemStatBase *getProblem(int number = 0)
{
FUNCNAME("CoupledIteration::getProblem()");
if (number == 0) return problem1;
if (number == 1) return problem2;
ERROR_EXIT("invalid problem number\n");
return NULL;
}
private:
ProblemStatBase *problem1;
ProblemStatBase *problem2;
};
\end{lstlisting}
The class \verb+MyCoupledIteration+ is finished now.
The next class, \verb+Identity+, implements the identity $I(x)=x$ for \verb+double+ variables. An arbitrary degree can be given to it. The class is used later to determine the quadrature degree used for the right hand side of \verb+problem2+.
\begin{lstlisting}{}
class Identity : public AbstractFunction<double, double>
// ===== add boundary conditions for problem1 =====
problem1.addDirichletBC(1, new G);
\end{lstlisting}
So far, we created and initialized \verb+problem1+ and its boundary conditions.
Now, we create \verb+problem2+. It should have its own finite element space, system, solver and file writer, but the mesh should be adopted from \verb+problem1+.
\begin{lstlisting}{}
// ===== create and init the second problem =====
Flag initFlag =
INIT_FE_SPACE |
INIT_SYSTEM |
INIT_SOLVER |
INIT_FILEWRITER;
Flag adoptFlag =
CREATE_MESH |
INIT_MESH;
ProblemScal problem2("problem2");
problem2.initialize(initFlag,
&problem1,
adoptFlag);
\end{lstlisting}
The operators for the first problem are defined like in Section
\ref{s:ellipt code}. Here, we only show the operators of
