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

[Matthew Knepley]

On Tue, Jun 2, 2015 at 6:13 AM, Justin Chang <jychang48 at gmail.com> wrote:

> 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)
>

So here you would need 4 dual basis vectors, which I am assuming are:

  ev_(-1, -1), ev_(1, -1), ev_(-1, 1), ev_(-0.5, -0.5)

where ev_(x, y) is the point evaluation functional at (x, y). Then you need
some basis for the primal space, which naively is

  1, x, y, 1

As you can see, this basis in linearly dependent, so the Vandermonde matrix
that FIAT constructs will be singular. The construction of a nodal basis
will fail.

So Jed's question is, what are they actually doing internally?

  Thanks,

    Matt


> 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.
>>
>>
>


-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20150602/87dc5679/attachment.html>