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

[Geoffrey Irving]

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.

Geoffrey