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

[Sander Arens]

Ok, I found why it was failing. The gmsh reader didn't invert the
tetrahedra (or hexahedra) like it does in the exodus reader (both have the
same vertex ordering).
The hdf5 to xdmf issue was because the block size for the local coordinates
vector was never set.
I created a pull request.

Thanks,
Sander

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

> 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/20160622/d53a7f65/attachment-0001.html>