On Thu, Nov 21, 2013 at 7:42 PM, Geoffrey Irving <[email protected]> wrote:
> On Thu, Nov 21, 2013 at 5:28 PM, Matthew Knepley <[email protected]> > wrote: > > 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. > > I verified that your change fails the regression in the same way as > mine if you change ex12 to say > > if (!has) { > ierr = DMPlexCreateLabel(dm, bdLabel);CHKERRQ(ierr); > ierr = DMPlexGetLabel(dm, bdLabel, &label);CHKERRQ(ierr); > ierr = DMPlexMarkBoundaryFaces(dm, label);CHKERRQ(ierr); > } > /* Complete even if has is true to test that fem boundary > integration ignores codimension != 1 faces */ > ierr = DMPlexGetLabel(dm, bdLabel, &label);CHKERRQ(ierr); > ierr = DMPlexLabelComplete(dm, label);CHKERRQ(ierr); > > Does this matter? Yes, it matters a lot. Thanks for finding this bug that I repeatedly missed. The tests that were not passing were all 3D P2. I had marked the vertices and faces, but not boundary edges (due to the way that CTetGen works). Thus, the midpoint dofs were unconstrained. Still a legit solution I think, but not what the test was intended for, and not a match to the provided exact solution. I have updates knepley/fix-fem-bd-integrate and 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
