## Re: [Libmesh-users] Bug in Tet10 4th and 5th order Gauss quadrature rules

 Re: [Libmesh-users] Bug in Tet10 4th and 5th order Gauss quadrature rules From: Benjamin S. Kirk - 2005-03-23 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, "Moderate-degree tetrahedral quadrature formulas", Comp Meth Appl Mech and Eng, vol 55:339-348 (1986). > I found this while implementing a Stokes solver and running it with Taylor-Hood Tet10 elements. > I was getting weird errors with a FOURTH and FIFTH Gauss quadrature order, which disappeared when I used the SIXTH order one (tensor-product 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 > > ```

 [Libmesh-users] Bug in Tet10 4th and 5th order Gauss quadrature rules From: - 2005-03-23 10:44:28 ```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, "Moderate-degree tetrahedral quadrature formulas", Comp Meth Appl Mech and Eng, vol 55:339-348 (1986). I found this while implementing a Stokes solver and running it with Taylor-Hood Tet10 elements. I was getting weird errors with a FOURTH and FIFTH Gauss quadrature order, which disappeared when I used the SIXTH order one (tensor-product 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 -- Luca Antiga, PhD ---------------------------------------------- Biomedical Technologies Laboratory Bioengineering Department, Mario Negri Institute Villa Camozzi, 24020, Ranica (BG), Italy ---------------------------------------------- phone: +39 035 4535-381 email: antiga@... web: http://villacamozzi.marionegri.it/~luca ---------------------------------------------- ```
 Re: [Libmesh-users] Bug in Tet10 4th and 5th order Gauss quadrature rules From: Benjamin S. Kirk - 2005-03-23 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, "Moderate-degree tetrahedral quadrature formulas", Comp Meth Appl Mech and Eng, vol 55:339-348 (1986). > I found this while implementing a Stokes solver and running it with Taylor-Hood Tet10 elements. > I was getting weird errors with a FOURTH and FIFTH Gauss quadrature order, which disappeared when I used the SIXTH order one (tensor-product 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 > > ```
 [Libmesh-users] Bug in Tet10 4th and 5th order Gauss quadrature rules From: John Peterson - 2005-03-23 16:39:29 ```Thanks! Open source strikes again. I wish we could start offering bounties for bugs like Gnome does... :) -John lantiga@... writes: > 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, "Moderate-degree tetrahedral quadrature > formulas", Comp Meth Appl Mech and Eng, vol 55:339-348 (1986). I > found this while implementing a Stokes solver and running it with > Taylor-Hood Tet10 elements. I was getting weird errors with a > FOURTH and FIFTH Gauss quadrature order, which disappeared when I > used the SIXTH order one (tensor-product rule). The errors finally > vanished when I introduced the corrections below. ```