function  contour_n_ave_flux_grad (i_p, i_d)
Pts_wide=2;
% Copyright 2000, J.E. Akin.  All rights reserved.
% plot node scp ave i_d gradient component of i_p flux component
% convert any type of mesh to a structured square mesh

% The flux gradients are stored as N_R_B flux components
% associated with N_SPACE grident components. Here N_SPACE = 2
% Poisson Eq: U,x,x U,y,x : U,y,x U,y,y :
% Plane stress: e_x,x e_y,x e_xy,x : e_x,y e_y,y, e_xy,y
% or U,x,x U,y,x (U,y,x + U,x,x) : U,x,y U,y,y (U,y,y + U,x,y)

% x,y   ==  nodal coordinates of the FEA
% f     ==  nodal solution of the FEA
% X,Y   ==  nodal coordinates of the structured square mesh
% F     ==  interpolated solution for the structured square mesh
% NX    ==  number of nodes in x-direction for structured mesh 
% NY    ==  number of nodes in y-direction for structured mesh

if ( nargin == 0 )
  i_p = 1 ; i_d = 1 ;
end % if no arguments
if ( nargin == 1 )
  i_d = 1 ;
end % if no arguments
format short

load msh_bc_xyz.tmp;
np = size (msh_bc_xyz, 1);
fprintf ('Read %g mesh coordinate pairs \n', np)
pre_p = 1;
ns = size (msh_bc_xyz,2) - pre_p ;   % space dimension
if ( ns < 2 )
  error ('This is not a 2D mesh')
end % if not 2D data

pre_e  = 0 ;
load  msh_typ_nodes.tmp ; % nod_per_el nodes per element
nt         = size (msh_typ_nodes,1);           % number of elements
if ( nt == 0 )
  error ('Error missing file msh_typ_nodes.tmp')
end % if error
nod_per_el = size (msh_typ_nodes,2) - pre_e -1 ; % nodes per element

load pt_ave_grad_flux.tmp; % x,y followed by flux list
ng = size (pt_ave_grad_flux, 1) ;     % count points
nf = size (pt_ave_grad_flux, 2) ;  % count components
if ( ng == 0 )
  error ('Error missing file pt_ave_grad_flux.tmp')
end % if error
if ( nf <= 0 )
  error ('Error, no fluxes in file pt_ave_grad_flux.tmp')
else
  fprintf ('%g flux components with 2 gradient components \n', nf/2)
  nf = nf/2  ;
end % if error

% NX, NY = 31 are default values, these can be altered by user
NX = 41; NY = 41;

x = msh_bc_xyz (:, 1+pre_p);
y = msh_bc_xyz (:, 2+pre_p);
 
if ( i_d > 0 )
  i_c = (i_d - 1)*nf + i_p ;
  f =  pt_ave_grad_flux (:, i_c) ;
else  % average cross derivatives
  i_c = (abs(i_d) - 1)*nf + i_p  ;
  f =  pt_ave_grad_flux (:, i_c) ;
  i_c = (i_p - 1)*nf + abs(i_d) ;
  f =  0.5 * (f + pt_ave_grad_flux (:, i_c)) ;
  % for j = 1:ng
    % f(j) = sqrt ( sum(pt_ave_grad_flux (j, 3:(nf+2)).^2) ) ;
  % end % for ng points
end % which component

xlin = linspace (min(x), max(x), NX);
ylin = linspace (min(y), max(y), NY);
[X, Y] = meshgrid (xlin, ylin);
F = griddata (x, y, f, X, Y, 'cubic');

clf
c = contour (X,Y,F,'LineWidth',Pts_wide);
clabel (c);
hold on
grid

xlabel (['X: ', int2str(nt),' Elements'])
ylabel (['Y: ', int2str(np),' Nodes'])
% xlabel ('X'); ylabel ('Y');
if ( i_d < 0 )
  title (['Matlab Smoothed FE Average (U,', int2str(i_p), ',', ...
         int2str(abs(i_d)), ' and U,', int2str(abs(i_d)), ',', ...
         int2str(i_p), ') at the Nodes' ]);
  %b title (['Matlab Smoothed FE Grad\_', int2str(i_d),  ...
          %b ' of Flux\_', int2str(i_p), ' at the Nodes' ]); 
else
  title (['Matlab Smoothed FE U,', int2str(i_p), ',', ...
         int2str(i_d), ' at the Nodes' ]);
 %title (['Matlab Smoothed FEA QP Flux RMS Value at ', ...
 %        int2str(ng), ' Gauss Points'])
end % which component
%       -depsc -tiff            % for an eps version
  %bprint ('-dpsc', ['contour_n_ave_flux_grad_', int2str(i_p)])
  hold off
  %bv_text = ['Created contour_n_ave_flux_grad_', int2str(i_p), '.ps'];
  %bfprintf (1,'%s', v_text) ; fprintf (1, ' \n' )

% end of contour_n_ave_flux_grad