[petsc-dev] Is DMPlexLabelComplete necessary in snes ex12?

[Geoffrey Irving]

On Thu, Nov 21, 2013 at 5:28 PM, Matthew Knepley <knepley at gmail.com> wrote:
> On Thu, Nov 21, 2013 at 5:01 PM, Geoffrey Irving <irving at naml.us> wrote:
>>
>> On Thu, Nov 21, 2013 at 1:31 PM, Geoffrey Irving <irving at naml.us> wrote:
>> > On Thu, Nov 21, 2013 at 1:27 PM, Matthew Knepley <knepley at gmail.com>
>> > wrote:
>> >> On Thu, Nov 21, 2013 at 3:23 PM, Geoffrey Irving <irving at naml.us>
>> >> wrote:
>> >>>
>> >>> On Wed, Nov 20, 2013 at 3:52 PM, Matthew Knepley <knepley at gmail.com>
>> >>> wrote:
>> >>> > On Tue, Nov 19, 2013 at 5:27 PM, Geoffrey Irving <irving at naml.us>
>> >>> > wrote:
>> >>> >>
>> >>> >> This question isn't about changing ex12, I'm just making sure I
>> >>> >> understand what's going on:
>> >>> >>
>> >>> >> SetupSection in snes ex12 calls DMPlexLabelComplete on the "marker"
>> >>> >> label if Dirichlet conditions are chosen (not if Neumann is
>> >>> >> chosen).
>> >>> >> Specifically,
>> >>> >>
>> >>> >>     if (user->bcType == DIRICHLET) {ierr = DMPlexLabelComplete(dm,
>> >>> >> label);CHKERRQ(ierr);}
>> >>> >>
>> >>> >> It then uses the same label ("marker") to create an index set:
>> >>> >>
>> >>> >>     if (user->bcType == DIRICHLET) {ierr  = DMPlexGetStratumIS(dm,
>> >>> >> bdLabel, 1, &bcPoints[0]);CHKERRQ(ierr);}
>> >>> >>
>> >>> >> I believe 1 is the codimension.  If the index set is restricted to
>> >>> >> codimension 1, and DMPlexLabelComplete labels vertices based on
>> >>> >> edge
>> >>> >> labels, isn't the DMPlexLabelComplete call unnecessary?
>> >>> >
>> >>> >
>> >>> > 1) The value 1 is just a convention that I use for marking the
>> >>> > boundary.
>> >>> > In a real simulation,
>> >>> >      there are probably many makers and they are user specified.
>> >>> > This is
>> >>> > what happens in
>> >>> >      TS ex11.
>> >>> >
>> >>> > 2) I should probably complete the label in the Neumann case, but I
>> >>> > just
>> >>> > forgot. The completion
>> >>> >      is not necessary the way my defaults mesh generators work, but
>> >>> > is
>> >>> > sometimes necessary
>> >>> >      when reading stuff in.
>> >>>
>> >>> Actually, it looks like calling DMPlexLabelComplete in the Neumann
>> >>> case breaks with an error later.  In plexfem.c, if feBd is nonzero,
>> >>> DMPlexComputeCellGeometry is called on all "points", where points are
>> >>> actually whatever the DMLabel "boundary" spits out.
>> >>> DMPlexComputeGeometry bails with an exception if called on anything
>> >>> with dimension 0:
>> >>>
>> >>> [0]PETSC ERROR: DMPlexComputeCellGeometry() line 783 in
>> >>> src/dm/impls/plex/plexgeometry.c Unsupported dimension 0 in cell 2 for
>> >>> element geometry computation
>> >>>
>> >>> Should the "points" variable be "faces"?  Should I restrict the loops
>> >>> in plexfem.c to only touch codimension 1 elements?
>> >>
>> >> Definitely we should guard the feBd loop to only look at co-dimension 1
>> >> things. I think
>> >> we need to allow all points in boundary labels.
>> >
>> > Patch coming up.
>>
>> Pushed as irving/fix-plexfem-boundary-loop.  Let me know if there's a
>> better way to do this.  I can't just continue out of the loop since
>> then PetscFEIntegrateBdResidual would see bad data.
>>
>> If I complete unconditionally in ex12, plexfem crashes with an error
>> without the patch.  With the patch, most of the tests of ex12 pass,
>> but some of them do not (thus BROKEN in the commit message).  The
>> output is attached.  Any ideas where the mistake might be?
>
>
> I could not easily find the bug, and I did not want to malloc, so I did it
> again in
>
>   knepley/fix-fem-bd-integrate
>
> and it passes the regression for ex12. I merged to next.

I verified that your change fails the regression in the same way as
mine if you change ex12 to say

  if (!has) {
    ierr = DMPlexCreateLabel(dm, bdLabel);CHKERRQ(ierr);
    ierr = DMPlexGetLabel(dm, bdLabel, &label);CHKERRQ(ierr);
    ierr = DMPlexMarkBoundaryFaces(dm, label);CHKERRQ(ierr);
  }
  /* Complete even if has is true to test that fem boundary
integration ignores codimension != 1 faces */
  ierr = DMPlexGetLabel(dm, bdLabel, &label);CHKERRQ(ierr);
  ierr = DMPlexLabelComplete(dm, label);CHKERRQ(ierr);

Does this matter?

Geoffrey