[petsc-users] Problem with interpolating projected field of quadratic order

[Matthew Knepley]

On Sat, Jun 18, 2016 at 12:49 PM, Sander Arens <Sander.Arens at ugent.be>
wrote:

> 1) Does the test run for you?
>
>
> I ran the tests and they all pass.
>
> 2) If so, then could it be that you assume something about the order of
>> functions?
>
>
> I'm not really sure what you mean with this.
>
> So the values in 'u' are the FEM coefficients for the coordinate function
>> {x, y}. Why are you calling 'x' projected and
>> 'y' interpolated?
>>
>
> No, u is a vector with three fields which are obtained in the following
> way:
>
> 1) Project the coordinates onto field nr 0 so that I have the coordinates
> for all the dofs
> 2) Project the coordinates onto field nr 1, so that I get the coordinates
> at the centroid of the cell
> 3) Project field nr 0 onto field nr 2, so that I get the coordinates at
> the centroid of the cell
>
> I then output the values of the different fields, so that I can compare
> the values of field 1 and 2.
>
> What I would expect is that the values of field 1 and 2 would be the same.
> This is the case for field 0 being P1, but not for P2.
> I also projected the coordinates onto a P2 field in
> snes/examples/tutorials/ex77.c, but this was done using DMPlexCreateBoxMesh
> and gave no problems.
> Maybe it's some assumption on the vertex numbering in
> DMPlexCreateFromCellList or DMPlexCreateFromFile that I'm not aware of?
>
> I hope I made the problem clearer now?
>

Okay, I see what happened now. It is indeed a shortcoming. The way I have
things setup now, all fields share the quadrature
since that way I can evaluate all fields in the same loop at a quadrature
point. This could of course be changed at some cost.
I have changed your example to run.

  Thanks,

    Matt


> Thanks,
> Sander
>
> On 18 June 2016 at 16:33, Matthew Knepley <knepley at gmail.com> wrote:
>
>> On Sat, Jun 18, 2016 at 9:26 AM, Matthew Knepley <knepley at gmail.com>
>> wrote:
>>
>>> On Sat, Jun 18, 2016 at 8:56 AM, Sander Arens <Sander.Arens at ugent.be>
>>> wrote:
>>>
>>>> Hello all,
>>>>
>>>> Recently I've been trying to run a fem problem using DMPlex, PetscFE,
>>>> etc. using a mesh created with gmsh. I used DMPlexReverseCell to get rid of
>>>> the negative jacobian determinants. However, I noticed another problem:
>>>> when I projected the coordinates on a quadratic fem-field and interpolated
>>>> them they didn't match the interpolated values of the original coordinate
>>>> field anymore.
>>>>
>>>> I attached a simple test to reproduce this. With linear interpolation I
>>>> see no problem, it's only when I start using quadratic interpolation.
>>>>
>>>> I have no idea why this is giving problems. Am I missing something
>>>> simple here or is this a bug?
>>>>
>>>
>>> Tests for this already exist
>>>
>>>
>>> https://bitbucket.org/petsc/petsc/src/996a00f1a1693e2d85044ed633fe1848cf357f2a/src/dm/impls/plex/examples/tests/ex3.c?at=master&fileviewer=file-view-default
>>>
>>> with run parameters here:
>>>
>>>
>>> https://bitbucket.org/petsc/petsc/src/996a00f1a1693e2d85044ed633fe1848cf357f2a/config/builder.py?at=master&fileviewer=file-view-default#builder.py-58
>>>
>>> and specifically for P2 triangles
>>>
>>>
>>> https://bitbucket.org/petsc/petsc/src/996a00f1a1693e2d85044ed633fe1848cf357f2a/config/builder.py?at=master&fileviewer=file-view-default#builder.py-72
>>>
>>> 1) Does the test run for you?
>>>
>>>   ./config/builder2.py check src/dm/impls/plex/examples/tests/ex3.c
>>>
>>> 2) If so, then could it be that you assume something about the order of
>>> functions?
>>>
>>
>> I am having a hard time understanding what you are doing here:
>>
>>   /* View the projected coordinates field */
>>   ierr = VecGetSubVector(u, fields[0], &subvec);CHKERRQ(ierr);
>>   ierr = PetscObjectSetName((PetscObject) subvec, fieldNames[0]);CHKERRQ(ierr);
>>   ierr = VecViewFromOptions(subvec, NULL, "-projectedcoordinates_vec_view");CHKERRQ(ierr);
>>   ierr = VecRestoreSubVector(u, fields[0], &subvec);CHKERRQ(ierr);
>>
>>   /* View the interpolated coordinates field */
>>   ierr = VecGetSubVector(u, fields[1], &subvec);CHKERRQ(ierr);
>>   ierr = PetscObjectSetName((PetscObject) subvec, fieldNames[1]);CHKERRQ(ierr);
>>   ierr = VecViewFromOptions(subvec, NULL, "-interpolatedcoordinates_vec_view");CHKERRQ(ierr);
>>   ierr = VecRestoreSubVector(u, fields[1], &subvec);CHKERRQ(ierr);
>>
>>
>> So the values in 'u' are the FEM coefficients for the coordinate function
>> {x, y}. Why are you calling 'x' projected and
>> 'y' interpolated?
>>
>>   Thanks,
>>
>>      Matt
>>
>>
>>>   Thanks,
>>>
>>>      Matt
>>>
>>>
>>>> Thanks,
>>>> Sander
>>>>
>>>
>>>
>>>
>>> --
>>> 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
>>>
>>
>>
>>
>> --
>> 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
>>
>
>


-- 
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/20160618/13422427/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: makefile
Type: application/octet-stream
Size: 3079 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20160618/13422427/attachment-0001.obj>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: test.c
Type: text/x-csrc
Size: 9261 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20160618/13422427/attachment-0001.c>