On Thu, Nov 21, 2013 at 6:30 PM, Geoffrey Irving <[email protected]> wrote:
> On Thu, Nov 21, 2013 at 4:11 PM, Matthew Knepley <[email protected]> > wrote: > > On Thu, Nov 21, 2013 at 1:48 PM, Geoffrey Irving <[email protected]> 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
