[petsc-dev] integer divide by zero when using a two triangle DMPlex

[Geoffrey Irving]

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

If you're referring to the "order" parameter to
PetscDTGaussJacobiQuadrature, there's no sense in which it could have
been meant to mean number of nodes.  The first line is

    PetscInt       npoints = dim > 1 ? dim > 2 ? order*PetscSqr(order)
: PetscSqr(order) : order;

>> So arguably user error, but it would be nice if the defaults didn't
>> cause integer division exceptions.
>
> Yes, we should have a check for numPoints > 0 in the integration function.
>
>> Is silently producing no quadrature points the right thing for
>> PetscFESetQuadrature to do?  If so, can I check somewhere else to
>> avoid the zero division?
>
> I do not know. I tried to allow this corner case, but I am not sure its
> useful.

I can't think of a reason someone would want a quadrature rule with no
points, other than to indicate that the quadrature rule simply isn't
defined.  I suppose that's a valid reason, though.

So, I started adding checks to the integration routines, and it seemed
pretty ridiculous (PetscFEGetQuadrature is called a lot).  It seems
like the check should go either in PetscDTGaussJacobiQuadrature and
friends or PetscFESetQuadrature.  Or maybe PetscFEGetQuadrature?  I'm
not sure what the standards are here, or whether we want to support
setting and detecting explicitly uninitialized quadrature rules.

Geoffrey