[firedrake] Projecting CG1 function into DG0 function

[David Ham]

On Fri, 13 Nov 2015 at 16:29 Justin Chang <jychang48 at gmail.com> wrote:

> Yeah I found out about errornorm recently. Have been using that now.
>
> But as for the original issue... would D.assign(assemble(div(q) * e * dx)
> suffice?
>

I'm afraid not.

assemble(div(q)*e*dx)

is

\int_\Omega div(q) e dx  \forall e \in E

where E is the basis for DG0.

This is not the same thing as the projection of div(q) into DG0. In fact,
for the very specific case of DG0, this will differ from the correct answer
in each element by a factor of the volume of that element. This is why
dividing by

\int_\Omega e dx \forall e \in E

is necessary to achieve the correct answer.

To go into slightly more detail, note that the e_i \in E is:

             /   1 if x \in T_i (the i-th element)
e_i(x) = |
             \    0 otherwise.

Hence:

\int_\Omega \div(q) e_i dx   ==  \int_{T_i} div(q)) dx

and

\int_\Omega e_i  dx == \int_{T_i} dx
                                == vol(T_i)

but it must be the case that:

\int_T_i D_i e_i dx == \int_{T_i} div(q)) dx

and note:

\int_T_i D_i e_i dx == D_i \int_T_i dx

(since D_i is just a scalar, not a function of x) hence by rearranging:

D_i = \int_\Omega \div(q) e_i dx/ \int_\Omega e_i dx

and hence:

D = \int_\Omega \div(q) e_i dx/ \int_\Omega e_i dx \forall e_i \in E


>
> On Fri, Nov 13, 2015 at 8:27 AM, David Ham <David.Ham at imperial.ac.uk>
> wrote:
>
>>
>>
>> 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
>>>>
>>>>
>>>
>> _______________________________________________
>> firedrake mailing list
>> firedrake at imperial.ac.uk
>> https://mailman.ic.ac.uk/mailman/listinfo/firedrake
>>
>>
>
-------------- next part --------------
HTML attachment scrubbed and removed