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):
Then the flux along the boundary must coincide with $u_0(x,y)$ but I get the following noisy picture:
What am I doing wrong?
