On Thu, Nov 21, 2013 at 1:31 PM, Geoffrey Irving <[email protected]> wrote: > 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.
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
ex12-output
Description: Binary data
