Projection of div-free vector field on crossed mesh is not div-free

mesh
1

Hi all,

I feel like this should have come up already, but I was not able to find anything through the search function. Here is my MWE (2019.1.0):

from dolfin import *

n = 16
mesh_L = UnitSquareMesh( n, n, 'left')
mesh_C = UnitSquareMesh( n, n, 'crossed')

V_L = VectorFunctionSpace(mesh_L, 'CG', 1)
V_C = VectorFunctionSpace(mesh_C, 'CG', 1)

ic = Expression( ('sin(x[1])', 'sin(x[0])'), degree=1)

v_L = project(ic, V_L)
v_C = project(ic, V_C)

print("div v_L_L2 =", assemble(div(v_L)*div(v_L)*dx) )
# = 4.77822736512426e-29
print("div v_C_L2 =", assemble(div(v_C)*div(v_C)*dx) )
# = 6.6524224789526e-05

Is there a way to achieve “exact” projection on crossed meshes? I would like to do some calculations on locally refined meshes and the same problem arises there. Many thanks in advance.

Edit: The same behavior effect is there for interpolation.

2

The zero divergence with mesh_L is a red herring, as L^2 projection (as implemented by the project function) does not in general preserve the divergence-free property of a vector field. The result you’re seeing is a coincidence, due to the fact that ic is of the form f(x_2)\mathbf{e}_1 + g(x_1)\mathbf{e}_2, and is being interpolated into the CG1 space prior to projection (so project is not actually doing anything, because ic is being projected into the space it is already in). If you project the same vector field but using SpatialCoordinates directly in UFL instead of interpolating first as an Expression, you will not have zero divergence on mesh_L.

To robustly obtain H^1-conforming div-free vector fields on simplicial meshes, you might look into Scott–Vogelius elements, although the resulting vector field will need to have higher polynomial order and projection would need to be done in a mixed space, including a pressure-like Lagrange multiplier, e.g., find \Pi\mathbf{u}\in\mathcal{V}^h and p^h\in\mathcal{Q}^h s.t. \forall\mathbf{v}^h\in\mathcal{V}^h and q^h\in\mathcal{Q}^h

\left(\Pi\mathbf{u} - \mathbf{u},\mathbf{v}^h\right)_{L^2} + \left(q^h,\nabla\cdot\Pi\mathbf{u}\right)_{L^2} - \left(p^h,\nabla\cdot\mathbf{v}^h\right)_{L^2} = 0\text{ ,}

so you’d need to set up the variational problem yourself, since project only includes the first term on the left-hand side. However, you can actually get away with never creating a mixed space or setting up a saddle-point problem, if you instead use the “iterated penalty” method described (with FEniCS code) in these slides, but for the Stokes problem instead of projection:

Alternatively, if maintaining a low-order approximation is important and you only really need H(\text{div})-conforming fields, you could use Raviart–Thomas elements.

3

Thank you very much for your reply. I see now that i stripped down my MWE a little bit too much, I am in fact using Scott-Vogelius elements. I noticed the problem because my solution “jumped” at the first timestep (from u0 to u1), which is due to the fact that there the proper saddle point problem is solved and u1 is exactly divergence-free while u0, the interpolation, was not. I now do the proper “manual” projection as you suggested as a “0th” timestep and the behavior looks fine.

I am still not sure if I understand the mesh-dependence though. Especially in a simple case where the x-component is a function of y only and vice-versa, shouldn’t the interpolation preserve this fact? What is so special about the “standard” simplical meshes (“left” and “right”) that it holds in this case?

4

It’s a little hard to explain clearly without drawing a picture, but imagine what the linear interpolant of f(x_2) will look like on the “left” or “right” mesh: It is a linear interpolant in the x_2 direction that is extruded in a constant way along the x_1 direction, so \partial_{x_1}I_h(f) = 0, where I_h is interpolation on linears. Interpolating also at the center node of each “crossed” cell breaks this pattern.