function fem2d_pack_test24 ( ) %*****************************************************************************80 % %% TEST24 tests TRIANGLE_UNIT_SET. % % Licensing: % % This code is distributed under the GNU LGPL license. % % Modified: % % 14 February 2003 % % Author: % % John Burkardt % order_max = 64; rule_max = 20; fprintf ( 1, '\n' ); fprintf ( 1, 'TEST24\n' ); fprintf ( 1, ' TRIANGLE_UNIT_SET sets up a quadrature\n' ); fprintf ( 1, ' in the unit triangle,\n' ); fprintf ( 1, '\n' ); for a = 0 : 10 for b = 0 : 10 - a coef = ( a + b + 2 ) * ( a + b + 1 ); for i = 1 : b coef = coef * ( a + i ) / i; end fprintf ( 1, '\n' ); fprintf ( 1, ' A = %f B = %f\n', a, b ); fprintf ( 1, '\n' ); fprintf ( 1, ' Rule QUAD ERROR\n' ); fprintf ( 1, '\n' ); for rule = 1 : rule_max order = triangle_unit_size ( rule ); [ xtab, ytab, weight ] = triangle_unit_set ( rule ); quad = 0.0; for i = 1 : order x = xtab(i); y = ytab(i); if ( a == 0 && b == 0 ) value = coef; elseif ( a == 0 && b ~= 0 ) value = coef * ytab(i)^b; elseif ( a ~= 0 && b == 0 ) value = coef * xtab(i)^a; elseif ( a ~= 0 && b ~= 0 ) value = coef * xtab(i)^a * ytab(i)^b; end quad = quad + 0.5 * weight(i) * value; end exact = 1.0; err = abs ( exact - quad ); fprintf ( 1, ' %4d %14f %14f\n', rule, quad, err ); end end end return end