From: Benjamin S. Kirk <benkirk@cf...>  20050323 13:21:58

Thank you very much! I either caused the error in the original implementation or the book has an error. Either way, good catch. Let me apologize for wasting your time, this is a bug you shouldn't have had to find. If I was rigorous about writing unit tests I'd put together some code assuring that all our quadrature rules can integrate polynomials of the required degrees accurately, but that would leave no time to implement new features! Does the paper you reference happen to have any rules that aren't implemented in the library? I'm not on campus right now, otherwise I'd head over to the library and check. Ben lantiga@... wrote: > Dear libmesh people, > I found what seems to be a bug in 4th and 5th order Gauss quadrature rules for Tet10 elements. > I don't have the Zienkiewicz book, so I can't check on the reference given in the code, however I'm referring to P. Keast, "Moderatedegree tetrahedral quadrature formulas", Comp Meth Appl Mech and Eng, vol 55:339348 (1986). > I found this while implementing a Stokes solver and running it with TaylorHood Tet10 elements. > I was getting weird errors with a FOURTH and FIFTH Gauss quadrature order, which disappeared when I used the SIXTH order one (tensorproduct rule). The errors finally vanished when I introduced the corrections below. > > In particular, in file src/quadrature/quadrature_3D.C, > > case FOURTH: > { > ... > points[5](0) = a; > points[5](1) = a; > points[5](2) = a; > ... > points[10](0) = b; > points[10](1) = b; > points[10](2) = b; > ... > } > > should be instead > > case FOURTH: > { > ... > points[5](0) = b; > points[5](1) = a; > points[5](2) = a; > ... > points[10](0) = a; > points[10](1) = b; > points[10](2) = a; > ... > } > > Similarly, > > case FIFTH: > { > ... > points[9](0) = a; > points[9](1) = a; > points[9](2) = a; > ... > points[14](0) = b; > points[14](1) = b; > points[14](2) = b; > ... > } > > should be instead > > case FIFTH: > { > ... > points[9](0) = b; > points[9](1) = a; > points[9](2) = a; > ... > points[14](0) = a; > points[14](1) = b; > points[14](2) = a; > ... > } > > > Luca > > 