On Thu, Nov 21, 2013 at 1:27 PM, Matthew Knepley <[email protected]> wrote: > On Thu, Nov 21, 2013 at 3:23 PM, Geoffrey Irving <[email protected]> wrote: >> >> On Wed, Nov 20, 2013 at 3:52 PM, Matthew Knepley <[email protected]> >> wrote: >> > On Tue, Nov 19, 2013 at 5:27 PM, Geoffrey Irving <[email protected]> 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
