[petsc-users] Modified Taylor-Hood elements with piece-wise constant pressure for Stokes equation

[Justin Chang]

In FEniCS's Stokes example (example 19), one defines the Taylor-Hood
function spaces with these three lines:

V = VectorFunctionSpace(mesh, "CG", 2)
Q = FunctionSpace(mesh, "CG", 1)
W = V * Q

To implement P2/(P1+P0), all we gotta do is this:

V = VectorFunctionSpace(mesh, "CG", 2)
Q = FunctionSpace(mesh, "CG", 1)
P = FunctionSpace(mesh, "DG", 0)
W = V * (Q + P)

Everything else (boundary conditions, variational form, etc) remains the
same. The sum of element-wise divergence (via assemble(div(u)*dx)) for the
former is -0.000454806468878 whereas the latter its 3.91121126886e-14 hence
ensuring (better) local mass conservation.

I haven't looked into the exact framework on how this is done, but attached
are the UFL/FFC generated header files that implement the Taylor-Hood and
Modified Taylor-Hood discretization.

Not sure how it's done in Deal.II, but IIRC it should be possible.

Thanks,
Justin

On Tue, Jun 2, 2015 at 12:31 AM, Jed Brown <jed at jedbrown.org> wrote:

> Justin Chang <jychang48 at gmail.com> writes:
> > I am not quite sure what you're asking for. Are you asking for how people
> > actually implement this augmented TH? In other words, how the shape/basis
> > functions for this mixed function space would look?
>
> Ciarlet's finite element definition requires a local approximation space
> ("P1+P0" in your case -- though P0 is a subset of P1 so this is actually
> just P1) and a set of "nodes" (a dual basis, which will define the basis
> functions used for computation and the continuity between elements).
> The dual basis needs to be unisolvent, such that the Vandermonde matrix
> is invertible.
>
> Now what might that look like for "P1+P0"?
>
> And if you bypass Ciarlet's construction and just try to pick some
> "basis" functions, say {(1-x)/2, (1+x)/2, 1} for the 1D simplex [-1,1],
> then you can express the constant f(x)=1 as u=(1,1,0), v=(0,0,1), or t*u
> + (1-t)*v for any real value t.  Clearly it fails to be a basis because
> it's linearly dependent.
>
> Can you still compute with it?  The answer is probably yes if you pin
> some cell displacement, but that causes a bunch of somewhat unattractive
> side-effects that I asked about in my earlier message.
>
> Now you said it was "just a few additional lines of code", so maybe you
> can explain how people are implementing it.
>
> > I have only seen in some key note lectures and presentations at
> > conferences briefly mentioning this P2/(P1+P0) element, as if it's the
> > de facto discretization for Stokes flows.
>
> It's hardly "de facto".
>
> > That said, even I am not too sure how this would look.
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20150602/a99c97c1/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: TH.h
Type: text/x-chdr
Size: 409646 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20150602/a99c97c1/attachment-0002.h>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: MTH.h
Type: text/x-chdr
Size: 454613 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20150602/a99c97c1/attachment-0003.h>