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

[Matthew Knepley]

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.


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

   Matt


>
> Geoffrey
>



-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20131121/01c5053d/attachment.html>