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

[Geoffrey Irving]

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?

Geoffrey
-------------- next part --------------
A non-text attachment was scrubbed...
Name: ex12-output
Type: application/octet-stream
Size: 37106 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20131121/42f20436/attachment.obj>