After doing some mathematics related to the stability of elements in 3D Stokes problem I was slightly shocked to realize that $P_2-P_1$ is not stable for an arbitrary tetrahedral mesh. More precisely, in case you have an element where all nodes and three out of four facets lie on the boundary of the domain with a Dirichlet condition, you end up getting a singular matrix. This is in fact fairly trivial to conclude from the weak form of the Stokes system.

I tested the only commercial Stokes code I have an access to (COMSOL) and it allowed me to create such a mesh. Upon clicking solve I get 'Error: Singular matrix' as expected. (I am under the impression that COMSOL uses $P_2-P_1$ for its creeping flow module.)

The bad mesh.

To further test that the problem was not related to other configurations I tried the following mesh and everything works as expected.

enter image description here

Questions: Is this kind of constraint taken into account in (adaptive or non-adaptive) mesh generators? I see from various research papers that this element seems to be quite popular. Are these kind of boundary instabilities generally disregarded as insignificant when choosing a method to use? More importantly, what does it really mean to have a stable finite element, i.e., what kind of mesh-dependent instabilities are too much so that we conclude that the method is bad?

  • 5
    $\begingroup$ Interesting question! As far as I see, these elements typically result from structured tetrahedral mesh generation on cubes and such and play only a minor role in real applications where you have unstructured nodalization algorithms. I have tried a bit around a while ago and was not able to produce such a mesh with a mesh generator producing fully unstructured meshes. I suspect they employ a mechanism to avoid such over-constrained elements. I have no access to COMSOL though but I guess that for most solvers these element do not pose a significant problem. $\endgroup$ Mar 1, 2016 at 19:30
  • $\begingroup$ I wonder if this is also a problem with the MINI element? $\endgroup$ Mar 4, 2016 at 0:45
  • 1
    $\begingroup$ I think MINI does not have this problem. The single interior DOF will save the situation. Uniqueness condition for Stokes is $(\nabla\cdot v, p)=0~\forall v \in V_h \Rightarrow p =\text{global const.}$. Let $p(x,y)=a+bx+cy$. Choose $v=(b\phi,c\phi)$ where $\phi$ is the bubble in some tetra. This gives that $p$ is local constant and continuity takes care of the rest. $\endgroup$
    – knl
    Mar 4, 2016 at 8:41

1 Answer 1


This problem also arises in 2D, when all the vertices of a triangle lie on the boundary. And Stokes problem is not the only problem which may fail with such meshes, a $p$-Laplacian problem for some values of $p$ will find a zero gradient in the element and fail.

Mesh generators generally have an option to handle this, e.g. the 2D mesh generator bamg of freefem++ has a -splitpbedge option that adds a node in the middle of any edge having both ends on the boundary. According to bamg documentation, unstructured mesh generation can return such triangles.

  • $\begingroup$ Are you sure that this is the case with e.g. Taylor-Hood in 2D Stokes? My intuition tells me that the DOF related to the edge saves the situation there. In 3D Taylor-Hood, there is no DOF related to the facet and hence the instability occurs. $\endgroup$
    – knl
    Nov 16, 2016 at 18:21
  • $\begingroup$ You're right, it may be the case. I think the Verfuhrt proof of the inf-sup condition for Taylor-Hood is constructive enough to check this, but no time just now. $\endgroup$
    – Joce
    Nov 17, 2016 at 8:34

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.