+/**
+ @page mse Minimal surface equation
+
<a name="Intro"></a>
<h1>Introduction</h1>
<h3>Foreword</h3>
This programm deals with an example of a non-linear elliptic pde, the minimal
-surface equation. You can imagine the solution as a soap bubble inside a
-closed wire, where the wire isn't smooth, but curved. The soap bubble will
+surface equation. You can imagine the solution as a soap bubble inside a
+closed wire, where the wire isn't smooth, but curved. The soap bubble will
take a shape with minimal surface. The solution of the minimal surface equation
-describes this shape with the wire as a boundary condition.
+describes this shape with the wire as a boundary condition.
-Because the equation is non-linear, we can't solve it directly, but have to use
+Because the equation is non-linear, we can't solve it directly, but have to use
the newton-method to compute the solution iterativly.
@f[
u=g \qquad\qquad on ~ \partial \Omega
@f]
-
+
In this example, we choose the unitball as our domain $\Omega$.
-As described above, we have to formulate the Newton-method for this problem
+As described above, we have to formulate the Newton-method for this problem
with a damping parameter $\lambda$ to have a better global convergence behaviour:
@par
@f[
u^{n+1}=u^{n}+\delta u^{n}
@f]
-
+
with:
@f[
F(u):= -\nabla \cdot \left( \frac{1}{\sqrt{1+|\nabla u|^{2}}}\nabla u \right)
@f]
-
+
and $F'(u,\delta u)$ the derivative of F in direction of $\delta u$:
@f[
the solution of:
@f[
- - \nabla \cdot \left( \frac{1}{(1+|\nabla u^{n}|^{2})^{\frac{1}{2}}}\nabla
- \delta u^{n} \right) +
- \nabla \cdot \left( \frac{\nabla u^{n} \cdot
- \nabla \delta u^{n}}{(1+|\nabla u^{n}|^{2})^{\frac{3}{2}}} \nabla u^{n}
- \right) =
+ - \nabla \cdot \left( \frac{1}{(1+|\nabla u^{n}|^{2})^{\frac{1}{2}}}\nabla
+ \delta u^{n} \right) +
+ \nabla \cdot \left( \frac{\nabla u^{n} \cdot
+ \nabla \delta u^{n}}{(1+|\nabla u^{n}|^{2})^{\frac{3}{2}}} \nabla u^{n}
+ \right) =
-\left( - \nabla \cdot \left( \frac{1}{(1+|\nabla u^{n}|^{2})^{\frac{1}{2}}}
\nabla u^{n} \right) \right)
@f]
In order to solve the minimal surface equation, we have to solve this equation in every
-Newton step. To solve this, we have to take a look at the boundary condition of this
-problem. Assuming that $u^{n}$ already has the right boundary values, the Newton update
+Newton step. To solve this, we have to take a look at the boundary condition of this
+problem. Assuming that $u^{n}$ already has the right boundary values, the Newton update
$\delta u^{n}$ should have zero boundary conditions, in order to have the right boundary
condition after adding both.
-In the first Newton step, we are starting with the solution $u^{0}\equiv 0$, the Newton
-update still has to deliever the right boundary condition to the solution $u^{1}$.
+In the first Newton step, we are starting with the solution $u^{0}\equiv 0$, the Newton
+update still has to deliever the right boundary condition to the solution $u^{1}$.
@par
-Summing up, we have to solve the pde above with the boundary condition $\delta u^{0}=g$
+Summing up, we have to solve the pde above with the boundary condition $\delta u^{0}=g$
in the first step and with $\delta u^{n}=0$ in all the other steps.
both sides of the pde with a testfunction $\varphi$ and integrating by parts on both sides:
@f[
- \left( \nabla \varphi , \frac{1}{(1+|\nabla u^{n}|^{2})^{\frac{1}{2}}}\nabla
- \delta u^{n} \right)-\left(\nabla \varphi ,\frac{\nabla u^{n} \cdot \nabla
- \delta u^{n}}{(1+|\nabla u^{n}|^{2})^{\frac{3}{2}}}\nabla u^{n} \right)
+ \left( \nabla \varphi , \frac{1}{(1+|\nabla u^{n}|^{2})^{\frac{1}{2}}}\nabla
+ \delta u^{n} \right)-\left(\nabla \varphi ,\frac{\nabla u^{n} \cdot \nabla
+ \delta u^{n}}{(1+|\nabla u^{n}|^{2})^{\frac{3}{2}}}\nabla u^{n} \right)
= -\left(\nabla \varphi , \frac{1}{(1+|\nabla u^{n}|^{2})^{\frac{1}{2}}} \nabla u^{n}
\right)
@f]
-
-Where the solution $\delta u^{n}$ is a function in the infinte space $H^{1}(\Omega)$.
+
+Where the solution $\delta u^{n}$ is a function in the infinte space $H^{1}(\Omega)$.
Reducing this space to a finite space with basis $\left\{ \varphi_{0},\dots ,
\varphi_{N-1}\right\}$, we can write the solution:
{\sqrt{1+|\nabla u^{n}|^{2}}}$, we can rewrite the weak formualtion:
@f[
- \sum_{j=0}^{N-1}\left[ \left( \nabla \varphi_{i} , a_{n} \nabla \varphi_{j} \right) -
- \left(\nabla u^{n}\cdot \nabla \varphi_{i} , a_{n}^{3} \nabla u^{n} \cdot \nabla
- \varphi_{j} \right) \right] \cdot U_{j}=\left( \nabla \varphi_{i} , a_{n}
- \nabla u^{n}\right) \qquad \forall i=0,\dots ,N-1
+ \sum_{j=0}^{N-1}\left[ \left( \nabla \varphi_{i} , a_{n} \nabla \varphi_{j} \right) -
+ \left(\nabla u^{n}\cdot \nabla \varphi_{i} , a_{n}^{3} \nabla u^{n} \cdot \nabla
+ \varphi_{j} \right) \right] \cdot U_{j}=\left( \nabla \varphi_{i} , a_{n}
+ \nabla u^{n}\right) \qquad \forall i=0,\dots ,N-1
@f]
where the solution $\delta u^{n}$ is given by the coefficents $\delta U^{n}_{j}$.
where the entries of the matrix $A^{n}$ are given by:
@f[
- A^{n}_{ij}:= \left( \nabla \varphi_{i} , a_{n} \nabla \varphi_{j} \right) -
- \left(\nabla u^{n}\cdot \nabla \varphi_{i} , a_{n}^{3} \nabla u^{n} \cdot \nabla
+ A^{n}_{ij}:= \left( \nabla \varphi_{i} , a_{n} \nabla \varphi_{j} \right) -
+ \left(\nabla u^{n}\cdot \nabla \varphi_{i} , a_{n}^{3} \nabla u^{n} \cdot \nabla
\varphi_{j} \right)
@f]
b^{n}_{i}:=\left( \nabla \varphi_{i} , a_{n} \nabla u^{n}\right)
@f]
-The matrix A is symmetric, but it is indefinite. So we have to take a better look
-at the solver we choose for this linear system. The CG-method needs
+The matrix A is symmetric, but it is indefinite. So we have to take a better look
+at the solver we choose for this linear system. The CG-method needs
positive-definiteness of the matrix A, which is not given, so it can't be used.
Using the symmetry of the matrix we can choose the minimal residual method as a
solver, which needs symmetry but no definiteness.
<h3>Summary</h3>
-Starting with the function $u^{0}\equiv 0$, the first Newton update is computed by
+Starting with the function $u^{0}\equiv 0$, the first Newton update is computed by
solving the system $A^{0}U^{0}=b^{0}$ with boundary condition $\delta u^{0}=g$ on
- $\partial \Omega$. The new approximation of the solution is given by
- $u^{1}=u^{0}+\delta u^{0}$. The next updates are given as solution of
- the linear system $A^{n}U^{n}=b^{n}$ with boundary condition $\delta u^{n}=0$ on
+ $\partial \Omega$. The new approximation of the solution is given by
+ $u^{1}=u^{0}+\delta u^{0}$. The next updates are given as solution of
+ the linear system $A^{n}U^{n}=b^{n}$ with boundary condition $\delta u^{n}=0$ on
$\partial \Omega$ and the new approximation given by $u^{n+1}=u^{n}+\delta u^{n}$.
+
+*/
\ No newline at end of file
-The final solution, as written by the program at the end of the
+The last solution, as written by the program at the end of the
<code>%run()</code> function, looks as follows:
-In each cycle, the program furthermore writes the grid in EPS
-format. These are shown in the following:
+After each refinement, the solution after the fifth Newton-
+iteration is shown in the following pictures. On the left
+the solution is visualized, on the right side the solution
+with the mesh is shown:
</tr>
</table>
+It is clearly visible, that the solution minimizes the surface
+after each refinement. The solution converges to a picture one
+would imagine a soapbubble inside a wire, which is bended like
+the boundary. Also it is visible, how the boundary
+is smoothed out after each refinement. On the coarse mesh,
+the boundary doesn't look like a sine, whereas it does the
+finer the mesh gets.
+The mesh is mostly refined near the boundary, where the solution
+increases or decreases strongly, whereas it is coarsened on
+the inside of the domain, where nothing interesting happens,
+because there isn't much change in the solution.
-It is clearly visible that the region where the solution has a kink,
-i.e. the circle at radial distance 0.5 from the center, is
-refined most. Furthermore, the central region where the solution is
-very smooth and almost flat, is almost not refined at all, but this
-results from the fact that we did not take into account that the
-coefficient is large there. The region outside is refined rather
-randomly, since the second derivative is constant there and refinement
-is therefore mostly based on the size of the cells and their deviation
-from the optimal square.
while(first_step || (res>1e-3))
{
+ // In the first step, we compute the solution on the two times globally
+ // refined mesh. After that the mesh will be refined
+ // adaptively, in order to not get too many cells. The refinement
+ // is the first thing done every time we restart the process in the while-loop.
+ if(!first_step)
+ {
+ refine_grid();
+
+ std::cout<<"********mesh-refinement:"<<refinement+1<<" ********"<<std::endl;
+ refinement++;
+ }
+
+
// First thing to do after refining the mesh, is to setup the vectors,
// matrices, etc., which is done in the <code>setup_system</code>
// function.
std::ofstream output (filename.c_str());
data_out.write_vtk (output);
- // last thing to do is to refine the mesh:
-
- refine_grid();
-
- std::cout<<"********mesh-refinement:"<<refinement+1<<" ********"<<std::endl;
- refinement++;
-
}
-
- // After the residual is less than $10^{-3}$, the final-solution
- // is written in a vtk-file:
-
- DataOutBase::EpsFlags vtk_flags;
-
- DataOut<dim> data_out;
- data_out.set_flags (vtk_flags);
-
- data_out.attach_dof_handler (dof_handler);
- data_out.add_data_vector (present_solution, "solution");
- data_out.build_patches (6);
-
- std::ofstream output ("final-solution.vtk");
- data_out.write_vtk (output);
}
// @ sect4{The main function}
}
return 0;
}
+