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

Sun Jun 19 16:03:30 CDT 2016

And the files of course.

On 19 June 2016 at 23:02, Sander Arens <Sander.Arens at ugent.be> wrote:

> Thanks, Matt. That indeed made the problem for the test dissapear.
>
> However, there's still a problem when I try to work with a gmsh mesh.
>
> I modified snes/examples/tutorials/ex77 to read in this mesh (it's just a
> simple bar generated in gmsh) and I used DMPlexReverseCell because the
> jacobian determinants where negative. Snes doesn't converge and looking for
> example in the residual function f1_u_3d shows that x and u are not the
> same at the beginning (which should be though, as the coordinates where
> projected onto u). Please see the attached files ex77.c, makefile and
> bar.msh.
>
> I then tried working with a displacement field instead of deformation and
> this seemed to work (because the initial values where just zero), but the
> directions of the face normals are in the wrong direction because the bar
> was stretched instead of compressed. Please see the attached files
> ex77_2.c, makefile_2 and bar.msh.
>
> So did I forget something when using a read in mesh or is it something
> else?
>
> Also, /bin/petsc_gen_xdmf.py seemed to fail on the output ex77_2.h5 of
> make runex77_2_3 (when using ex77_2.c, makefile_2 and bar.msh).
>
> Thanks for all the help,
> Sander
>
> On 19 June 2016 at 02:17, Matthew Knepley <knepley at gmail.com> wrote:
>
>> 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/20160619/d2f608bc/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: ex77.c
Type: text/x-csrc
Size: 33264 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20160619/d2f608bc/attachment-0002.c>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: makefile
Type: application/octet-stream
Size: 1563 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20160619/d2f608bc/attachment-0002.obj>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: bar.msh
Type: model/mesh
Size: 42509 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20160619/d2f608bc/attachment-0001.msh>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: makefile_2
Type: application/octet-stream
Size: 2275 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20160619/d2f608bc/attachment-0003.obj>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: ex77_2.c
Type: text/x-csrc
Size: 33700 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20160619/d2f608bc/attachment-0003.c>