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

  • $\begingroup$ 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 ? $\endgroup$
    – cfdlab
    Commented May 3, 2017 at 12:46
  • $\begingroup$ 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$). $\endgroup$ Commented May 3, 2017 at 13:02
  • $\begingroup$ 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. $\endgroup$
    – cfdlab
    Commented May 3, 2017 at 13:09
  • $\begingroup$ 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. $\endgroup$
    – nicoguaro
    Commented May 3, 2017 at 13:28
  • $\begingroup$ 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) $\endgroup$ Commented May 3, 2017 at 13:42

3 Answers 3


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.

  • $\begingroup$ In the problem for $T^*$, the test function $w$ belongs to a trace space, not $H^1(\Omega)$. $\endgroup$
    – cfdlab
    Commented May 4, 2017 at 15:54
  • $\begingroup$ @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 $\endgroup$
    – DanielRch
    Commented May 7, 2017 at 8:16
  • $\begingroup$ Yes, would be nice to fix your answer. I have a Fenics implementation of this approach here github.com/cpraveen/fenics/blob/master/2d/… $\endgroup$
    – cfdlab
    Commented May 7, 2017 at 15:48
  • $\begingroup$ My problem is $\nabla \cdot \sigma (x) \nabla u =0$ how do you do that? $\endgroup$
    – badmf
    Commented Oct 16, 2017 at 10:17
  • $\begingroup$ @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. $\endgroup$
    – DanielRch
    Commented 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.

  • $\begingroup$ With more degrees of freedom will this take much longer to solve? $\endgroup$
    – badmf
    Commented Jun 2, 2017 at 8:42
  • $\begingroup$ 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. $\endgroup$
    – 56th
    Commented Jun 2, 2017 at 10:16
  • $\begingroup$ 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 $\endgroup$
    – 56th
    Commented 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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.