Neumann / natural BCs in FEA

I'm trying to work out an as-general-as-possible 2D Laplace example for Finite Element Analysis. Starting from $\Delta u = 0$ for an unknown $u(x,y)$, I multiply both sides (well, in practice only the left-hand side) by an arbitrary test function $\varphi(x,y)$ and subsequently integrate over the domain $\Omega$. Then, using the divergence theorem (Green's first identity) I obtain

$$\int_\Omega \nabla \varphi \nabla u \; d\Omega = \int_\Gamma \varphi \nabla u \cdot \vec{n} \; d\Gamma,$$

where $\Gamma$ is the boundary of $\Omega$ and $\vec{n}$ the (unit) normal vector on that boundary. Next, I approximate $u$ and $\varphi$ using the same set of shape functions $N$, resulting in $u \approx u^h = \vec{u}^T N(x,y) = N^T(x,y) \vec{u}$, where $\vec{u}$ is the a column of unknown coefficients, and likewise $\varphi \approx \varphi^h = \vec{\varphi}^T N(x,y)$. Substituting this gives

$$\boxed{\vec{\varphi}^T} \int_\Omega \nabla N \nabla N^T \; d\Omega \;\boxed{\vec{u}} = \boxed{\vec{\varphi}} \int_\Gamma N \nabla u \cdot \vec{n} \; d\Gamma.$$

Since this has to be valid for arbitrary test functions $\varphi$ and therefore for arbitrary coefficients $\vec{\varphi}$, we can write this as

$$\int_\Omega \nabla N \nabla N^T \; d\Omega \;\boxed{\vec{u}} = \int_\Gamma N \nabla u \cdot \vec{n} \; d\Gamma.$$

Finally, we need to take the boundary conditions into account. In the case of Dirichlet / essential BCs, we already know some of the coefficients in $\vec{u}$, they are given values for some of the boundary nodes. We can then either partition the system, or use some kind of penalty method to solve the system $K\vec{u} = \vec{q}$, where $K$ is the left-hand side (matrix) and $\vec{q}$ the right-hand side (vector). However, what to do with Neumann / natural BCs?

Question:

Note that I'm referring to general Neumann BCs, so $\nabla u \cdot \vec{n} = g$. If $g = 0$, then it's rather trivial (unfortunately, this is the example almost all FEA books seem to use). However, what to do if $g \neq 0$? I could discretize $\nabla u$ as $\nabla N^T \vec{u}$, and I can evaluate $\vec{n}$ everywhere, but then

• What do I do with this new $\vec{u}$, combine it with the left-hand side $K$?
• How do I prescribe the Neumann BCs (in discretized form), is it just a value, a difference, something else?

Basically, how do I apply Neumann BCs (and how are they prescribed), and how do I subsequently solve the system $K\vec{u} = \vec{q}$?

On the portions of $\Gamma$ where $\nabla u \cdot \vec{n} = g$, you simply substitute that expression into the right-hand-side integral and evaluate it. I.e., you evaluate:
$$\int_\Gamma N \nabla u \cdot \vec{n} \; d\Gamma = \int_\Gamma N g \; d\Gamma.$$
Edited to add: You may need to represent $g$ in your finite-element basis in order to evaluate it sanely. It's possible to use exact integration if your $g$ functions are simple, but usually we interpolate $g$ as $g_h$ and then integrate using the same kind of numerical integration that we used to get $K$. There are some tricks here often because $\Gamma$ is one dimension lower than $\Omega$, so you have to be careful when you integrate on such imbedded surfaces.
• So to be sure, I don't need to discretize $g$ in the way described in my post (that is, $g \approx \nabla N^T \vec{u} \cdot \vec{n}$)? Ok, so a given/prescribed function $g$ can be evaluated in the Gauss points, that should be fine. But for Dirichlet, it's common to provide nodal values, not a function, correct? Following the same approach, how do I come up with nodal values for $g$ (if that's indeed the way to go)? – Ailurus May 5 '14 at 13:11
• Presumably you have the functional form of $g$, and can just evaluate it at the Gauss points (if your integration order is high enough to achieve the error you want), or evaluate it in your basis and integrate that. It won't help to discretize $g$ in terms of $u$, because you don't know $u$, but you can write $g \approx \vec{g} N$ and integrate that. – Bill Barth May 5 '14 at 13:17
• And this also means that I don't have to evaluate the normal $\vec{n}$ at the Gauss points, correct? – Ailurus May 5 '14 at 13:31
• If the basis is a nodal basis, then providing the nodal values of $g$ is sufficient to represent it in the basis. This is pretty typical. In my code, a user is also allowed to provide a function which is then sampled at the Gauss points for integration. In that case, it's up to the user to decide if the current integration rule is high enough order to integrate their function well. Either way you don't have to know the normal at the Gauss points or the nodes (the latter may be multiply-defined). I'll see if I can dig up a clear reference, but I think Oden, Carey, and Becker's book has it. – Bill Barth May 5 '14 at 13:34