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

[Matthew Knepley]

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.

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