[firedrake] Projecting CG1 function into DG0 function

[David Ham]

On Thu, 12 Nov 2015 at 21:48 Justin Chang <jychang48 at gmail.com> wrote:

> Hi guys,
>
> So I was attempting to do something similar:
>
> L2_error_norm = norm(assemble(v*(u-u_exact)*dx)/assemble(v*dx))
>
> where v is test function, u is FE solution, and u_exact is the analytical
> solution. All of which are CG1 space. I get an error saying
> "fl.log.UFLException: Division by non-scalar is undefined"
>
> Know what's up?
>

Yes, that's either a missing feature or a bug depending on how you look at
it. However I think you are just confused in what you are trying to do.

What you probably meant was:

errornorm(u, u_exact)

(errornorm is a Firedrake-provided function).

What errornorm is basically doing is:

assemble((u-u_exact)**2*dx))**0.5

Note that this is a zero-form (no test or trial function; the answer is a
scalar) so there are no solves, or divisions by lumped mass.

errornorm actually does somewhat more sophisticated things such as
increasing the polynomial order of the functions in order to avoid certain
sorts of error,  so it's typically better to use errornorm than to do it
yourself.

Regards,

David



>
> Thanks,
> Justin
>
> On Thu, Nov 12, 2015 at 4:33 AM, Lawrence Mitchell <
> lawrence.mitchell at imperial.ac.uk> wrote:
>
>> On 12/11/15 11:27, Justin Chang wrote:
>> > David,
>> >
>> > So if D.assign(assemble(div(q) * e * dx)/assemble(e * dx)) returns
>> > cell-wise div(q), what's the denominator "/assemble(e * dx)" for?
>>
>>
>> It's just the normal FE L2 projection:
>>
>> you want:
>>
>> u = div(q)
>>
>> So you hit both sides with a test function and integrate:
>>
>> u*e*dx = div(q)*e*dx
>>
>> Where u is a trial function in DG0, e is a test function in DG and q
>> is in whereever.
>>
>> But, as David points out, the DG0 mass-matrix is completely diagonal.
>>  The values on the diagonal are just obtained by assembling the lone
>> test function:
>>
>> diag = assemble(e*dx)
>>
>> But now, because the matrix is diagonal, you can solve the linear
>> system pointwise by doing a pointwise division:
>>
>> rhs = assemble(div(q)*e*dx)
>>
>> solution = assemble(rhs/diag)
>>
>>
>> Cheers,
>>
>> Lawrence
>>
>>
>> _______________________________________________
>> firedrake mailing list
>> firedrake at imperial.ac.uk
>> https://mailman.ic.ac.uk/mailman/listinfo/firedrake
>>
>>
>
-------------- next part --------------
HTML attachment scrubbed and removed