function fem2d_pack_test06 ( ) %*****************************************************************************80 % %% TEST06 tests DIV_Q4. % % Licensing: % % This code is distributed under the GNU LGPL license. % % Modified: % % 14 February 2003 % % Author: % % John Burkardt % m = 21; n = 13; fprintf ( 1, '\n' ); fprintf ( 1, 'TEST06\n' ); fprintf ( 1, ' DIV_Q4 estimates divergence and vorticity\n' ); fprintf ( 1, ' using 4 node quadrilateral elements.\n' ); fprintf ( 1, '\n' ); fprintf ( 1, ' Original U, V data forms %d rows and ', m ); fprintf ( 1, ' %d columns.\n', n ); % % Set the nodes of the quadrilateral. % q(1:2,1:4) = [ ... 0.0, 0.0; ... 1.0, 0.0; ... 1.0, 1.0; ... 0.0, 1.0 ]'; ilo = 5; jlo = 5; ihi = 8; jhi = 8; % % Put dummy data into U and V at the data nodes. % fprintf ( 1, '\n' ); fprintf ( 1, ' I, J, X, Y, U(I,J), V(I,J)\n' ); fprintf ( 1, '\n' ); for i = 1 : m for j = 1 : n y(i,j) = ( ( m - i ) * q(2,1) ... + ( i - 1 ) * q(2,3) ) ... / ( m - 1 ); x(i,j) = ( ( n - j ) * q(1,1) ... + ( j - 1 ) * q(1,3) ) ... / ( n - 1 ); end end r8mat_print_some ( m, n, x, ilo, jlo, ihi, jhi, ' Some of X:' ); r8mat_print_some ( m, n, y, ilo, jlo, ihi, jhi, ' Some of Y:' ); for i = 1 : m for j = 1 : n u(i,j) = x(i,j) * y(i,j); v(i,j) = sin ( x(i,j) * x(i,j) + y(i,j) * y(i,j) ); end end r8mat_print_some ( m, n, u, ilo, jlo, ihi, jhi, ' Some of U:' ); r8mat_print_some ( m, n, v, ilo, jlo, ihi, jhi, ' Some of V:' ); % % Get DIV and VORT. % [ div, vort ] = div_q4 ( m, n, u, v, q ); % % Compare computed and known values at the centers of the elements. % for i = 1 : m-1 for j = 1 : n-1 ym(i,j) = ( ( 2 * m - 2 * i - 1 ) * q(2,1) ... + ( 2 * i - 1 ) * q(2,3) ) ... / ( 2 * m - 2 ); xm(i,j) = ( ( 2 * n - 2 * j - 1 ) * q(1,1) ... + ( 2 * j - 1 ) * q(1,3) ) ... / ( 2 * n - 2 ); end end r8mat_print_some ( m-1, n-1, xm, ilo, jlo, ihi, jhi, ' Some of XM:' ); r8mat_print_some ( m-1, n-1, ym, ilo, jlo, ihi, jhi, ' Some of YM:' ); fprintf ( 1, '\n' ); fprintf ( 1, ' I, J, DIV(I,J), Exact Divergence, Difference\n' ); fprintf ( 1, '\n' ); for i = ilo : ihi for j = jlo : jhi dudx = ym(i,j); dudy = xm(i,j); dvdx = 2.0 * xm(i,j) * cos ( xm(i,j) * xm(i,j) + ym(i,j) * ym(i,j) ); dvdy = 2.0 * ym(i,j) * cos ( xm(i,j) * xm(i,j) + ym(i,j) * ym(i,j) ); fprintf ( 1, '%4d %4d %7f %7f %12f %12f %12f\n', ... i, j, xm(i,j), ym(i,j), div(i,j), dudx + dvdy, ... div(i,j) - ( dudx + dvdy ) ); end fprintf ( 1, '\n' ); end fprintf ( 1, '\n' ); fprintf ( 1, ' I, J, VORT(I,J), Exact Vorticity, Difference\n' ); fprintf ( 1, '\n' ); for i = ilo : ihi for j = jlo : jhi dudx = ym(i,j); dudy = xm(i,j); dvdx = 2.0 * xm(i,j) * cos ( xm(i,j) * xm(i,j) + ym(i,j) * ym(i,j) ); dvdy = 2.0 * ym(i,j) * cos ( xm(i,j) * xm(i,j) + ym(i,j) * ym(i,j) ); fprintf ( 1, '%4d %4d %7f %7f %12f %12f %12f\n', ... i, j, xm(i,j), ym(i,j), vort(i,j), dvdx - dudy, ... vort(i,j) - ( dvdx - dudy ) ); end fprintf ( 1, '\n' ); end return end