8
$\begingroup$

I am solving the Dirichlet problem $$ \begin{cases} \Delta u = 0, \\ u|_{\partial D} = f, \end{cases} $$ in a $2d$ domain $D$ using the finite element method. What I want to get is the normal derivative of the solution $u$ on the boundary: $\tfrac{\partial u}{\partial \nu}$. In other words, I'm interested in FEM implementation of the Dirichlet-to-Neumann (Poincare-Steklov) map $\Lambda$. In Wikipedia I've found the following:

When the partial differential equation is discretized, for example by finite elements or finite differences, the discretization of the Poincaré–Steklov operator is the Schur complement obtained by eliminating all degrees of freedom inside the domain.

Hovewer, I could not find an explicit answer to this problem in Google. I'm programming in Matlab and I would extremely appreciate an optimized solution.

Update (30 May 2017). As it was pointed out by Praveen Chandrashekar in comments, I tried to solve this problem by formulating it in a weak form: find $h \in L^2(\partial D)$ such that for any $w \in H^{\frac 1 2}(\partial D)$, $$ \int_{\partial D} hw \,dl = \int_D \nabla w \nabla u dS. $$ As usually, I represent $h = \sum' h_j \psi_j$, $w = \sum w_i \psi_i$, $u = \sum u_l \psi_l$, where $\psi_i$ are the basis elements and $\sum'$ denotes summation over the boundary elements only, and I get $$ \sum_i' \sum_j' h_j w_i \int_{\partial D} \psi_i \psi_j dl = \sum_k \sum_l w_k u_l \int_D \nabla \psi_k \nabla \psi_l dS. $$ Since $\{w_j\}$ are arbitrary, we get the system $$ \sum_j' \int_{\partial D} \psi_i \psi_j dl \cdot h_j = \sum_l\int_D \nabla \psi_k \nabla \psi_l dS \, \cdot u_l. $$ Introducing the matrices $$ B_{ij} = \int_{\partial D} \psi_i \psi_j \, dl, \quad A_{ij} = \int_D \nabla \psi_i \nabla \psi_j dS, $$ where in the first matrix $i$, $j$ correspond to boundary elements and in the second matrix $i$ corresponds to the boundary elements and $j$ corresponds to all elements, I get the matrix equation $$ B h = Au, $$ for finding $h$.

However, if I implement this approach, the results are far from being precise. For example, consider the following circular domain and boundary values $u_0(x,y) = -y$ (triangulation created using DistMesh):

enter image description here

Then the flux along the boundary must coincide with $u_0(x,y)$ but I get the following noisy picture: enter image description here

What am I doing wrong?

Improve this question
$\endgroup$
11
  • 1
    $\begingroup$ I think this answers your question scicomp.stackexchange.com/a/26779/3970 $\endgroup$ – cfdlab May 28 '17 at 14:08
  • $\begingroup$ It is not clear if you are asking about DtN elements which are applied for example to take into account boundary conditions in infinity or about the calculation of fluxes on the boundary. Or you like to statically condense DOFs inside the body and express fluxes explicitly by boundary conditions data. Can you clarify? $\endgroup$ – likask May 29 '17 at 21:50
  • $\begingroup$ @likask I've updated the question $\endgroup$ – Appliqué May 30 '17 at 15:42
  • $\begingroup$ @Appliqué Note that your fluxes should be in L2, but you have piecewise continuous on the boundary. I could be wrong here, but those oscillations do not look good. One way to use P2 elements (6-node triangles) in the volume and use lower order to approximate fluxes or use constant base on each boundary element. Another way is to use mix formulation, like here mofem.eng.gla.ac.uk/mofem/html/tutorials.html $\endgroup$ – likask May 30 '17 at 17:35
  • 1
    $\begingroup$ Could you plot for comparison the normal derivatives without this post-processing stage? I suspect this is the best you can obtain with your mesh. Try to refine it and see if there's any improvement. $\endgroup$ – DanielRch May 30 '17 at 17:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.