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

[Sander Arens]

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/24d0fb14/attachment.html>