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

[Matthew Knepley]

On Thu, Nov 21, 2013 at 7:42 PM, Geoffrey Irving <irving at naml.us> wrote:

> 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?


Yes, it matters a lot. Thanks for finding this bug that I repeatedly
missed. The tests that
were not passing were all 3D P2. I had marked the vertices and faces, but
not boundary
edges (due to the way that CTetGen works). Thus, the midpoint dofs were
unconstrained.
Still a legit solution I think, but not what the test was intended for, and
not a match to the
provided exact solution. I have updates knepley/fix-fem-bd-integrate and
merged to next.

  Thanks,

      Matt


>
> Geoffrey
>



-- 
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-dev/attachments/20131122/3527ffee/attachment.html>