# FEM current toy problem

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$.

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.

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, I get some nonsense. For example, consider the following circular domain and boundary values $u_0(x,y) = -y$

Then the flux along the boundary must coincide with $u_0(x,y)$

I am trying to compute the flux of with the equations used above. However I was wondering what functions to use for $h$? If I use piecewise constant on elements, I get that the $B$ matrix has rank 1 too small (because its a circle), but how can one know this beforehand? (I tried to adjoin this with that the total flux must be $0$ but it gives wrong results). If I try to use a continous at vertices and linear on elements function, I get a flux which misses a factor of $1.4$. If I try to use discontinous but linear elements the $B$ matrix has more columns then rows. For $w$ I used continous piecewise linear elements which are $1$ only at nodes on the boundary, is this right?

What is the precise elements to use for all the functions, and can someone write a pseudo code for this problem?

• You'd get more answers if you inlined the problem description, rather than link to some other site. – Wolfgang Bangerth Dec 6 '17 at 14:05
• This problem was also considered here scicomp.stackexchange.com/questions/27002/… – badmf Dec 6 '17 at 14:16
• I'm not sure your formulation is correct. Can you explain how you ended up with the weak formulation for $h$? Secondly, the test function $w$ lives everywhere, so it can't be in $H^{1/2}(\partial\Omega)$ but in fact needs to be in $H^1(\Omega)$. (Though this then implies that its restriction to the boundary is in $H^{1/2}(\partial\Omega)$.) – Wolfgang Bangerth Dec 6 '17 at 16:38
• The weak form is from 2nd answer here scicomp.stackexchange.com/questions/26774/… – TROLLHUNTER Dec 6 '17 at 17:51
• But at least the issue with the space of the test functions is wrong there already. – Wolfgang Bangerth Dec 6 '17 at 18:21

The formulation of this problem is tricky. Here is what you have in your original post:

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.$$ This already can't be quite right because you have $w$ on the right in a domain integral, so it cannot be a function $w \in H^{\frac 1 2}(\partial D)$ that only lives on the boundary. Rather, you need to consider all test functions $w\in H^1(D)$ for this.

Now, this means that on the left you have that the functions $h,w$ are from different spaces: you have "many more" test functions $w$ than solution functions $h$. If you discretized this, you'd get a rectangular matrix with more rows than columns.

Of course, in practice the "vast majority" of $H^1$ functions $w$ are in fact zero on the boundary. Indeed, the set of test functions $w\in H^1$ whose traces $w|_{\partial D}$ are linear independent (i.e.,for which the corresponding rows of the matrix on the left will be linearly independent) can all be written as the extension of $H^{1/2}(\partial D)$ functions into the interior. How exactly one extends into the interior does not actually matter; let us just say that $Ew \in H^1(D)$ is an extension of $w\in H^{1/2}(\partial D)$ so that $Ew|_{\partial D}=w$.

Then, the now correct problem formulation is as follows: Find $h \in H^{-1/2}(\partial D)$ such that for any $w \in H^{\frac 1 2}(\partial D)$, $$\int_{\partial D} hw \,dl = \int_D \nabla (Ew) \cdot \nabla u dS.$$

This is now the point where you can discretize. To this end, let $V_h \subset H^{1/2}(\partial D)$ be a finite dimensional subspace, for example consisting of piecewise constant functions on the boundary relative to the faces of finite element mesh of $D$. Then, if $\{\psi_j\}_{j=1}^N$ is a basis of $V_h$, you would get the following linear system for the finite element approximation $h_h = \sum_j H_j \psi_j$: $$\sum_j \left\{ \int_{\partial D} \psi_i \psi_j dl \right\} H_j = \int_D \nabla (E\psi_i) \nabla u_h dS.$$ But, importantly, the finite element approximation $u_h$ is defined with respect to a different set of shape functions, say $\varphi_k \in W_h \subset H^1(D)$ so that $u_h = \sum_k U_k \varphi_k$. It is important to realize that the $\varphi_k$ are not just extensions of the $\psi_k$ -- their number is entirely different, and they form an unrelated space. With this, you can then write the problem as you wanted to: $$\sum_j \left\{ \int_{\partial D} \psi_i \psi_j dl \right\} H_j = \sum_k \left\{\int_D \nabla (E\psi_k) \nabla \varphi_k dS\right\} U_k.$$ If you want to write this as a linear system with matrices, you'd get this: $$B_{ij} = \int_{\partial D} \psi_i \psi_j \, dl, \quad A_{ij} = \int_D \nabla (E\psi_i) \nabla \varphi_j dS,$$ and the linear system to find the coefficients $H_j$ of $h_h$ is given by $$B H = A U.$$ Here, $B$ is simply a mass matrix on the boundary, and consequently symmetric and positive definite. $A$, on the other hand, is a rectangular matrix.

Finally, we had not chosen an extension operator. For practical reasons, it would of course be nice if $Ew_i$ is a finite element function, for example one in $W_h$ (or at least a polynomial on each cell) because then you can just apply quadrature as always. This is easily achieved by making sure that $Ew_i$ is zero at all interior nodes. If you do this, then $A_{ij}$ is zero for all $j$ corresponding to nodes that are neither at the boundary or on cells that are adjacent to the boundary. (This would suggest more efficient ways to store the matrix $A$, but that's besides the point for the current post.)

• What is this $H^{-1/2}$ space? – badmf Dec 11 '17 at 12:00
• @WolfgangBangerth that is the great answer. Correct me if I am wrong, you restricting yourself to piecewise continuous in the domain and piecewise discontinuous functions on the boundary. You are taking advantage of that because $V_h$ has dual pace $V^*_h$, and space $V_h$ is pivot space, so those are the same and thus $V^*_h \subset H^{-1/2}$. As a result you have simple mass matrix. – likask Dec 11 '17 at 17:02
• @JustMe: $H^{-1/2}$ is the space of all functions that are Neumann values $\partial u/\partial n$ of $H^1$ functions $u$ restricted to the boundary. – Wolfgang Bangerth Dec 11 '17 at 18:30
• @likask: That's one way to see it. But also because $H^{1/2} \subset L_2 \subset H^{-1/2}$, you automatically have that $V_h \subset H^{-1/2}$ for any of the common finite element spaces you may choose. – Wolfgang Bangerth Dec 11 '17 at 18:32