# Computing accurate fluxes with FEM

I have solved Poisson equation on a 3d domain with neumann and dirichlet boundary condition. I get the potential, take the gradient for each element and integrate on a surface of an element, I do this for all elements on a surface enclosing an electrode I am interested in knowing the current through. This does not give very accurate results, upto 50% error. And commercial FEM software gets very accurate fluxes.
Does anyone know how they do it or know a better way to compute flux? I am unable to find any references on this

• Can you give more details on the problem preferably with equations ? Are you trying to compute the flux on the Dirichlet part of the boundary ? May 3, 2017 at 12:46
• Its kind of a general problem with FEM, but in my example Im interested in Poissons equation $\nabla \cdot \sigma \nabla u =0$ and yes on dirichlet part of boundary ($u=p_i$ on the boundary_i, rest is insulators $du/dn=0$). May 3, 2017 at 13:02
• You can use Lagrange multipliers in a post-processing step to compute the fluxes. See Gunzburger and Hou, "TREATING INHOMOGENEOUS ESSENTIAL BOUNDARY CONDITIONS IN FINITE ELEMENT METHODS AND THE CALCULATION OF BOUNDARY STRESSES", SIAM J. Num. Anal., vol. 20, no. 2, pp. 390-424, April 1992. May 3, 2017 at 13:09
• What is the discretization level in each case? Is the mesh the same? Are both calculations using the same type of elements? How are you computing the integral? Those are all good questions to answer (add to your question) in order to get better answers. May 3, 2017 at 13:28
• What is discretization level? What are you talking about? There is only one computation. The integral is just the sum of the gradient of the potential times normal vector multiplied by surface area. (Using linear elements) May 3, 2017 at 13:42

If averaging the element fluxes, as suggested by @Bill Greene, does not produce accurate enough fluxes, there are ways to improve the accuracy of the fluxes by post-processing. A standard and quite straightforward method, which is mentioned as an exercise in the book by Hughes T.J.R., "The Finite Element Method", is as follows:

If your problem is: \begin{align*} & \nabla^2 u = f \quad \text{in} \quad \Omega \\ & \frac{\partial u}{\partial n} = T \quad \text{on} \quad \Gamma_T && \text {Neumann data} \\ & \frac{\partial u}{\partial n} = T^* \quad \text{on} \quad \Gamma_g && \text {sought flux on the Dirichlet boundary} \\ & u = g \quad \text{on} \quad \Gamma_g && \text {Dirichlet data} \end{align*}

Then you find the unknown fluxes by solving the corresponding weak problem: Find $T^*\in L_2(\Gamma_g)$ such that for all $w\in H^{\frac{1}{2}}(\Gamma_g)$ there holds:

\begin{equation*} \int_{\Gamma_g} wT^* d\Gamma = \int_\Omega w_{,i}u_{,i} d\Omega - \int_{\Omega} wf d\Omega - \int_{\Gamma_T} wT d\Gamma, \end{equation*}

where the comma followed by index $i$ denotes differentiation with respect to the coordinates (the summation convention is implied). Note that the Dirichlet data are not used here.

You proceed to discretize this additional weak problem as you did for the main problem.

• In the problem for $T^*$, the test function $w$ belongs to a trace space, not $H^1(\Omega)$. May 4, 2017 at 15:54
• @PraveenChandrashekar, True. That would be $H^{\frac{1}{2}}(\Gamma_g)$, right? If I'm not mistaken, this does not change the actual computation here May 7, 2017 at 8:16
• Yes, would be nice to fix your answer. I have a Fenics implementation of this approach here github.com/cpraveen/fenics/blob/master/2d/… May 7, 2017 at 15:48
• My problem is $\nabla \cdot \sigma (x) \nabla u =0$ how do you do that? Oct 16, 2017 at 10:17
• @Just Me It should be quite similar. I suggest you write down the standard weak for for the Laplace case and see how the form above is different from that. Then apply the same process to your case. That is no proof but it should suffice. Oct 16, 2017 at 17:31

If $\nabla u$ is of more interest than $u$ itself, it is reasonable to set $\mathbb v := -\sigma \, \nabla u$ and convert your Poisson problem $$-\nabla \cdot (\sigma \, \nabla u) = f \quad \text{in }\Omega, \\ u = g_E \quad \text{on }\partial\Omega_E, \\ \sigma \, \nabla u \cdot \hat{\mathbb n} = g_N \quad \text{on }\partial\Omega_N$$ to a mixed one: $$\frac{1}{\sigma} \, \mathbb v + \nabla u = 0 \quad \text{in }\Omega, \\ \nabla \cdot \mathbb v = f \quad \text{in }\Omega, \\ u = g_E \quad \text{on }\partial\Omega_E, \\ \mathbb v \cdot \hat{\mathbb n} = -g_N \quad \text{on }\partial\Omega_N.$$ (Note that for the mixed problem essential BCs are natural and natural BCs are essential.)

Here you build an approximation for $\mathbb v$ directly since it is one of the unknowns. Hence you do not lose accuracy due to numerical differentiation.

This mixed Poisson formulation is pretty well discussed in literature.

• With more degrees of freedom will this take much longer to solve? Jun 2, 2017 at 8:42
• Saddle point problems are generally harder to solve than elliptic ones, but it is still possible to solve them efficiently with a suitable preconditioner. Please check the book from this answer: scicomp.stackexchange.com/a/27026/21916 Chapter 5 discusses block preconditioners.
– 56th
Jun 2, 2017 at 10:16
• It may be useful to look through this deal.II step: dealii.org/8.4.1/doxygen/deal.II/step_20.html Please also check “Better linear solvers” section
– 56th
Jun 2, 2017 at 10:30

One simple way to improve the accuracy of the fluxes is to average the discontinuous values of the element fluxes at each node. This gives you a continuous, nodal flux function over the surface that can be integrated.

If all the elements connected to a node are (roughly) the same size, simply averaging the element fluxes at the node is sufficiently accurate. If the element sizes are significantly different, the contribution of each element flux to the nodal value should be weighted by the size of the element.

There are more sophisticated methods for improving the values of fluxes at the nodes but the fundamental idea is to "smooth" the discontinuous values across inter-element boundaries. In structural analysis, this is commonly referred to as "stress smoothing."

Of course, the fluxes are going to be less accurate than the fundamental nodal unknowns because they are derivatives of these quantities. If you are using the simple 4-node tetrahedron, the flux is constant in the element. Particularly, when flux calculations are important, it is common to use the 10-node tetrahedron element instead.