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? Geoffrey
