<div dir="ltr"><div class="gmail_extra"><div class="gmail_quote">On Thu, Nov 21, 2013 at 7:58 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:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex"><div class=""><div class="h5">On Thu, Nov 21, 2013 at 4:51 PM, Matthew Knepley <<a href="mailto:knepley@gmail.com">knepley@gmail.com</a>> wrote:<br>

> On Thu, Nov 21, 2013 at 6:30 PM, Geoffrey Irving <<a href="mailto:irving@naml.us">irving@naml.us</a>> wrote:<br>
>><br>
>> On Thu, Nov 21, 2013 at 4:11 PM, Matthew Knepley <<a href="mailto:knepley@gmail.com">knepley@gmail.com</a>><br>
>> wrote:<br>
>> > On Thu, Nov 21, 2013 at 1:48 PM, Geoffrey Irving <<a href="mailto:irving@naml.us">irving@naml.us</a>> wrote:<br>
>> >><br>
>> >> I'm getting the following crash trying to use a two triangle mesh with<br>
>> >> DMPlex FEM.  The problem is that the PetscFEGetQuadrature call on line<br>
>> >> 546 of DMPlexComputeResidualFEM produces a quadrature object with<br>
>> >> q.numPoints = 0, which we then use for division.<br>
>> ><br>
>> > Yes, it causes batchSize = 0. I think we could error is q.numPoints ==<br>
>> > 0.<br>
>> ><br>
>> >> I'm using the default elements, which based on the "fe dofs" printout<br>
>> >> are piecewise constant.<br>
>> ><br>
>> > You are supposed to create the quadrature (like ex12.c:382). The default<br>
>> > thing is completely empty.<br>
>><br>
>> The code initializes quadrature just like ex12:<br>
>><br>
>>   /* Create quadrature (with specified order if given) */<br>
>>   PetscQuadrature q;<br>
>>   CHECK(PetscDTGaussJacobiQuadrature(dim,qorder>0?qorder:order,-1,1,&q));<br>
>>   CHECK(PetscFESetQuadrature(fe,q));<br>
>><br>
>> It works fine with "-petscspace_order 1" or "2", and computes the<br>
>> right answer to machine precision for 2 since the exact solution is<br>
>> quadratic.  The problem is that if you ask for zeroth order<br>
>> quadrature, you get no quadrature points.  Note that I'm using the<br>
>> same order quadrature as the polynomial space (I copied the code from<br>
>> ex12).<br>
><br>
> It appears that my use of the term 'order' is the problem. Jed does not like it either.<br>
> This is really the number of nodes.<br>
<br>
</div></div>If you're referring to the "order" parameter to<br>
PetscDTGaussJacobiQuadrature, there's no sense in which it could have<br>
been meant to mean number of nodes.  The first line is<br></blockquote><div><br></div><div>Yes, total number of nodes, accounting for dimension.</div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">

    PetscInt       npoints = dim > 1 ? dim > 2 ? order*PetscSqr(order)<br>
: PetscSqr(order) : order;<br>
<div class="im"><br>
>> So arguably user error, but it would be nice if the defaults didn't<br>
>> cause integer division exceptions.<br>
><br>
> Yes, we should have a check for numPoints > 0 in the integration function.<br>
><br>
>> Is silently producing no quadrature points the right thing for<br>
>> PetscFESetQuadrature to do?  If so, can I check somewhere else to<br>
>> avoid the zero division?<br>
><br>
> I do not know. I tried to allow this corner case, but I am not sure its<br>
> useful.<br>
<br>
</div>I can't think of a reason someone would want a quadrature rule with no<br>
points, other than to indicate that the quadrature rule simply isn't<br>
defined.  I suppose that's a valid reason, though.<br>
<br>
So, I started adding checks to the integration routines, and it seemed<br>
pretty ridiculous (PetscFEGetQuadrature is called a lot).  It seems<br>
like the check should go either in PetscDTGaussJacobiQuadrature and<br>
friends or PetscFESetQuadrature.  Or maybe PetscFEGetQuadrature?  I'm<br>
not sure what the standards are here, or whether we want to support<br>
setting and detecting explicitly uninitialized quadrature rules.</blockquote><div><br></div><div>A check would definitely go in PetscDTGaussJacobiQuadrature() if we</div><div>want to disallow empty quadrature, and I just not sure we want to do that,</div>
<div>since it makes things hard in parallel when you have empty processes. I</div><div>think I would rather put the high level check in IntegrateResidual().</div><div><br></div><div>   Matt</div><div><br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">
<span class=""><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>