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

[Geoffrey Irving]

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).

So arguably user error, but it would be nice if the defaults didn't
cause integer division exceptions.

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?

Geoffrey