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

[Matthew Knepley]

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.

  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/20131121/c7b58f95/attachment.html>