<div dir="ltr"><div class="gmail_extra"><div class="gmail_quote">On Thu, Nov 21, 2013 at 7:42 PM, Geoffrey Irving <span dir="ltr"><<a href="mailto:irving@naml.us" target="_blank">irving@naml.us</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div class="HOEnZb"><div class="h5">On Thu, Nov 21, 2013 at 5:28 PM, Matthew Knepley <<a href="mailto:knepley@gmail.com">knepley@gmail.com</a>> wrote:<br>
> On Thu, Nov 21, 2013 at 5:01 PM, Geoffrey Irving <<a href="mailto:irving@naml.us">irving@naml.us</a>> wrote:<br>
>><br>
>> On Thu, Nov 21, 2013 at 1:31 PM, Geoffrey Irving <<a href="mailto:irving@naml.us">irving@naml.us</a>> wrote:<br>
>> > On Thu, Nov 21, 2013 at 1:27 PM, Matthew Knepley <<a href="mailto:knepley@gmail.com">knepley@gmail.com</a>><br>
>> > wrote:<br>
>> >> On Thu, Nov 21, 2013 at 3:23 PM, Geoffrey Irving <<a href="mailto:irving@naml.us">irving@naml.us</a>><br>
>> >> wrote:<br>
>> >>><br>
>> >>> On Wed, Nov 20, 2013 at 3:52 PM, Matthew Knepley <<a href="mailto:knepley@gmail.com">knepley@gmail.com</a>><br>
>> >>> wrote:<br>
>> >>> > On Tue, Nov 19, 2013 at 5:27 PM, Geoffrey Irving <<a href="mailto:irving@naml.us">irving@naml.us</a>><br>
>> >>> > wrote:<br>
>> >>> >><br>
>> >>> >> This question isn't about changing ex12, I'm just making sure I<br>
>> >>> >> understand what's going on:<br>
>> >>> >><br>
>> >>> >> SetupSection in snes ex12 calls DMPlexLabelComplete on the "marker"<br>
>> >>> >> label if Dirichlet conditions are chosen (not if Neumann is<br>
>> >>> >> chosen).<br>
>> >>> >> Specifically,<br>
>> >>> >><br>
>> >>> >> if (user->bcType == DIRICHLET) {ierr = DMPlexLabelComplete(dm,<br>
>> >>> >> label);CHKERRQ(ierr);}<br>
>> >>> >><br>
>> >>> >> It then uses the same label ("marker") to create an index set:<br>
>> >>> >><br>
>> >>> >> if (user->bcType == DIRICHLET) {ierr = DMPlexGetStratumIS(dm,<br>
>> >>> >> bdLabel, 1, &bcPoints[0]);CHKERRQ(ierr);}<br>
>> >>> >><br>
>> >>> >> I believe 1 is the codimension. If the index set is restricted to<br>
>> >>> >> codimension 1, and DMPlexLabelComplete labels vertices based on<br>
>> >>> >> edge<br>
>> >>> >> labels, isn't the DMPlexLabelComplete call unnecessary?<br>
>> >>> ><br>
>> >>> ><br>
>> >>> > 1) The value 1 is just a convention that I use for marking the<br>
>> >>> > boundary.<br>
>> >>> > In a real simulation,<br>
>> >>> > there are probably many makers and they are user specified.<br>
>> >>> > This is<br>
>> >>> > what happens in<br>
>> >>> > TS ex11.<br>
>> >>> ><br>
>> >>> > 2) I should probably complete the label in the Neumann case, but I<br>
>> >>> > just<br>
>> >>> > forgot. The completion<br>
>> >>> > is not necessary the way my defaults mesh generators work, but<br>
>> >>> > is<br>
>> >>> > sometimes necessary<br>
>> >>> > when reading stuff in.<br>
>> >>><br>
>> >>> Actually, it looks like calling DMPlexLabelComplete in the Neumann<br>
>> >>> case breaks with an error later. In plexfem.c, if feBd is nonzero,<br>
>> >>> DMPlexComputeCellGeometry is called on all "points", where points are<br>
>> >>> actually whatever the DMLabel "boundary" spits out.<br>
>> >>> DMPlexComputeGeometry bails with an exception if called on anything<br>
>> >>> with dimension 0:<br>
>> >>><br>
>> >>> [0]PETSC ERROR: DMPlexComputeCellGeometry() line 783 in<br>
>> >>> src/dm/impls/plex/plexgeometry.c Unsupported dimension 0 in cell 2 for<br>
>> >>> element geometry computation<br>
>> >>><br>
>> >>> Should the "points" variable be "faces"? Should I restrict the loops<br>
>> >>> in plexfem.c to only touch codimension 1 elements?<br>
>> >><br>
>> >> Definitely we should guard the feBd loop to only look at co-dimension 1<br>
>> >> things. I think<br>
>> >> we need to allow all points in boundary labels.<br>
>> ><br>
>> > Patch coming up.<br>
>><br>
>> Pushed as irving/fix-plexfem-boundary-loop. Let me know if there's a<br>
>> better way to do this. I can't just continue out of the loop since<br>
>> then PetscFEIntegrateBdResidual would see bad data.<br>
>><br>
>> If I complete unconditionally in ex12, plexfem crashes with an error<br>
>> without the patch. With the patch, most of the tests of ex12 pass,<br>
>> but some of them do not (thus BROKEN in the commit message). The<br>
>> output is attached. Any ideas where the mistake might be?<br>
><br>
><br>
> I could not easily find the bug, and I did not want to malloc, so I did it<br>
> again in<br>
><br>
> knepley/fix-fem-bd-integrate<br>
><br>
> and it passes the regression for ex12. I merged to next.<br>
<br>
</div></div>I verified that your change fails the regression in the same way as<br>
mine if you change ex12 to say<br>
<br>
if (!has) {<br>
ierr = DMPlexCreateLabel(dm, bdLabel);CHKERRQ(ierr);<br>
ierr = DMPlexGetLabel(dm, bdLabel, &label);CHKERRQ(ierr);<br>
ierr = DMPlexMarkBoundaryFaces(dm, label);CHKERRQ(ierr);<br>
}<br>
/* Complete even if has is true to test that fem boundary<br>
integration ignores codimension != 1 faces */<br>
ierr = DMPlexGetLabel(dm, bdLabel, &label);CHKERRQ(ierr);<br>
<div class="im"> ierr = DMPlexLabelComplete(dm, label);CHKERRQ(ierr);<br>
<br>
</div>Does this matter?</blockquote><div><br></div><div>Yes, it matters a lot. Thanks for finding this bug that I repeatedly missed. The tests that</div><div>were not passing were all 3D P2. I had marked the vertices and faces, but not boundary</div>
<div>edges (due to the way that CTetGen works). Thus, the midpoint dofs were unconstrained.</div><div>Still a legit solution I think, but not what the test was intended for, and not a match to the</div><div>provided exact solution. I have updates knepley/fix-fem-bd-integrate and merged to next.</div>
<div><br></div><div> Thanks,</div><div><br></div><div> Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><span class="HOEnZb"><font color="#888888"><br>
Geoffrey<br>
</font></span></blockquote></div><br><br clear="all"><div><br></div>-- <br>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>
-- Norbert Wiener
</div></div>