On Thu, Nov 21, 2013 at 7:58 PM, Geoffrey Irving <[email protected]> wrote:
> On Thu, Nov 21, 2013 at 4:51 PM, Matthew Knepley <[email protected]> > wrote: > > 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. > > 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 > Yes, total number of nodes, accounting for dimension. > 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. A check would definitely go in PetscDTGaussJacobiQuadrature() if we want to disallow empty quadrature, and I just not sure we want to do that, since it makes things hard in parallel when you have empty processes. I think I would rather put the high level check in IntegrateResidual(). 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
