<div dir="ltr"><div class="gmail_extra"><div class="gmail_quote">On Thu, Nov 21, 2013 at 6:30 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="im">On Thu, Nov 21, 2013 at 4:11 PM, Matthew Knepley <<a href="mailto:knepley@gmail.com">knepley@gmail.com</a>> 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 == 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>
</div>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></blockquote><div><br></div><div>It appears that my use of the term 'order' is the problem. Jed does not like it either.</div><div>This is really the number of nodes.</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">

So arguably user error, but it would be nice if the defaults didn't<br>
cause integer division exceptions.<br></blockquote><div><br></div><div>Yes, we should have a check for numPoints > 0 in the integration function.</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">

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?</blockquote><div><br></div><div>I do not know. I tried to allow this corner case, but I am not sure its useful.</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>