\section{Coupled problems} \label{s:couple} 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 $$v = u.$$ 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}. \begin{figure} \center \includegraphics[width=0.9\textwidth]{fig/coupled_iteration.pdf} \caption{State diagram of the coupled iteration.} \label{f:coupled iteration} \end{figure} \begin{table} \center \begin{tabular}{|cc|c|c|c|} \hline & & {\bf 1d} & {\bf 2d} & {\bf 3d} \\ \hline {\bf source code} & \tt src/ & \multicolumn{3}{|c|}{\tt couple.cc} \\ \hline {\bf parameter file} & \tt init/ & \tt couple.dat.1d & \tt couple.dat.2d & \tt couple.dat.3d \\ \hline {\bf macro file} & \tt macro/ & \tt macro.stand.1d & \tt macro.stand.2d & \tt macro.stand.3d \\ \hline {\bf output files} & \tt output/ & \multicolumn{3}{|c|}{\tt couple.mesh, couple.dat} \\ \hline \end{tabular} \caption{Files of the {\tt couple} example.} \end{table} \subsection{Source code} 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. Now, we implement the needed interface methods: \begin{lstlisting}{} void beginIteration(AdaptInfo *adaptInfo) { FUNCNAME("StandardProblemIteration::beginIteration()"); MSG("\n"); MSG("begin of iteration %d\n", adaptInfo->getSpaceIteration()+1); MSG("=============================\n"); } void endIteration(AdaptInfo *adaptInfo) { FUNCNAME("StandardProblemIteration::endIteration()"); MSG("\n"); 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+). \begin{figure} \center \includegraphics[width=0.9\textwidth]{fig/stationary_loop.pdf} \caption{Stationary adaptation loop.} \label{f:stationary loop} \end{figure} 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 { public: Identity(int degree) : AbstractFunction(degree) {} double operator()(const double& x) const { return x; } }; \end{lstlisting} Now, we start with the main program: \begin{lstlisting}{} int main(int argc, char* argv[]) { FUNCNAME("main"); TEST_EXIT(argc == 2)("usage: couple initfile\n"); Parameters::init(true, argv[1]); // ===== create and init the first problem ===== ProblemScal problem1("problem1"); problem1.initialize(INIT_ALL); // ===== 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 \verb+problem2+. \begin{lstlisting}{} // ===== create operators for problem2 ===== Operator matrixOperator2(Operator::MATRIX_OPERATOR, problem2.getFESpace()); matrixOperator2.addZeroOrderTerm(new Simple_ZOT); problem2.addMatrixOperator(&matrixOperator2); Operator rhsOperator2(Operator::VECTOR_OPERATOR, problem2.getFESpace()); rhsOperator2.addZeroOrderTerm(new VecAtQP_ZOT(problem1.getSolution(), new Identity(degree))); problem2.addVectorOperator(&rhsOperator2); \end{lstlisting} At the left hand side, we have just an ordinary \verb+Simple_ZOT+. At the right hand side, we have a zero order term of the form $f(u)$ with $f=I$ the identity. $u$ is given by the solution DOF vector of \verb+problem1+. $I$ maps the values of the DOF vector evaluated at quadrature points to itself. The function degree is used to determine the needed quadrature degree in the assembler. Now, the adaptation loop is created: \begin{lstlisting}{} // ===== create adaptation loop and iteration interface ===== AdaptInfo *adaptInfo = new AdaptInfo("couple->adapt", 1); MyCoupledIteration coupledIteration(&problem1, &problem2); AdaptStationary *adapt = new AdaptStationary("couple->adapt", &coupledIteration, adaptInfo); \end{lstlisting} Note that not a pointer to one of the problems is passed to the adaptation loop, but a pointer to the \verb+coupledIteration+ object, which in turn knows both problems. The adaptation loop is now started. After it is finished, the solutions of both problems are written. \begin{lstlisting}{} // ===== start adaptation loop ===== adapt->adapt(); // ===== write solution ===== problem1.writeFiles(adaptInfo, true); problem2.writeFiles(adaptInfo, true); } \end{lstlisting} \subsection{Parameter file} We have one adaptation loop called \verb+couple->adapt+: \begin{lstlisting}{} couple->adapt->tolerance: 1e-8 couple->adapt->max iteration: 10 couple->adapt->refine bisections: 2 \end{lstlisting} The coupled problem consits of two sub problems. The first problem creates the mesh, solves its linear system of equations, estimates the error, adapts the mesh, and finally writes its output: \begin{lstlisting}{} coupleMesh->macro file name: ./macro/macro.stand.2d coupleMesh->global refinements: 0 problem1->mesh: coupleMesh problem1->dim: 2 problem1->polynomial degree: 1 problem1->solver: cg % no, bicgstab, cg, gmres, odir, ores problem1->solver->max iteration: 1000 problem1->solver->tolerance: 1.e-8 problem1->solver->left precon: diag problem1->estimator: residual problem1->estimator->C0: 0.1 % constant of element residual problem1->estimator->C1: 0.1 % constant of jump residual problem1->marker->strategy: 2 % 0: no 1: GR 2: MS 3: ES 4:GERS problem1->marker->MSGamma: 0.5 problem1->output->filename: output/problem1 problem1->output->AMDiS format: 1 problem1->output->AMDiS mesh ext: .mesh problem1->output->AMDiS data ext: .dat \end{lstlisting} The second problem uses the mesh of \verb+problem1+. So it creates no mesh, no estimator, and no marker. But a solver is needed to solve \verb+problem2+s linear system of equations, and a file writer to write the solution: \begin{lstlisting}{} problem2->dim: 2 problem2->polynomial degree: 1 problem2->solver: cg % no, bicgstab, cg, gmres, odir, ores problem2->solver->max iteration: 1000 problem2->solver->tolerance: 1.e-8 problem2->solver->left precon: diag problem2->estimator: no problem2->marker->strategy: 0 problem2->output->filename: output/problem2 problem2->output->AMDiS format: 1 problem2->output->AMDiS mesh ext: .mesh problem2->output->AMDiS data ext: .dat \end{lstlisting} \subsection{Macro file} We again use the macro file \verb+macro/macro.stand.2d+, which was described in Section \ref{s:ellipt macro}. \subsection{Output} The solution of the first problem is written to the files \verb+output/problem1.mesh+ and\\ \verb+output/problem1.dat+. The solution of the second problem is written to the files\\ \verb+output/problem2.mesh+ and \verb+output/problem2.dat+. We don't visualize the results here, because they conform with the results showed in Section \ref{s:vecellipt output}.