On Thu, Nov 21, 2013 at 5:01 PM, Geoffrey Irving <[email protected]> wrote:
> 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? I could not easily find the bug, and I did not want to malloc, so I did it again in knepley/fix-fem-bd-integrate and it passes the regression for ex12. I merged to next. Thanks, 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
