Neumann / natural BCs in FEA

Question

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}$?

Answer 1

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.