[petsc-users] Problem with interpolating projected field of quadratic order
Matthew Knepley
knepley at gmail.com
Wed Jun 22 10:14:32 CDT 2016
On Wed, Jun 22, 2016 at 7:52 AM, Sander Arens <Sander.Arens at ugent.be> wrote:
> 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 for finding this. I come back from the conference tomorrow, so I can
get it integrated before the end of the week.
Matt
> 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
>>>>
>>>
>>>
>>
>
--
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/543618d7/attachment-0001.html>
More information about the petsc-users
mailing list