function  contour_test (i_p)
% Copyright 2000, J.E. Akin.  All rights reserved.
% plot finite element exact as a surface for 1 <= i_p<= n_g_fre
% convert any type of mesh to a structured square mesh

% 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

% clear
if ( nargin == 0 )
  i_p = 0 ;
end % if no arguments

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 node_results.tmp
nr = size (node_results, 1);
  if ( nr == 0 )
error ('Error missing file node_results.tmp')
end % if error
max_p = size (node_results, 2) ; % number of columns
fprintf ('Read %g fea value sets \n', nr)
fprintf ('     with %g components each \n', max_p)

load exact_node_solution.tmp
nr = size (exact_node_solution, 1);
max_p = size (exact_node_solution, 2) ; % number of columns
if ( nr == 0 )
  error ('Error missing file exact_node_solution.tmp')
end % if error
fprintf ('Read %g nodal solution values \n', nr)
fprintf ('     with %g components each \n', max_p)

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

if ( i_p > max_p )
  error ('i_p > available data')
end % if error

fprintf ('Begin component %g error contour plots: \n', i_p)

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

% msh_bc_xyz has: pre_p items then: x, y
  x = msh_bc_xyz (1:np, (pre_p+1)) ; % extract x column
  y = msh_bc_xyz (1:np, (pre_p+2)) ; % extract y column 

  node_error = node_results ; %b- exact_node_solution ;
 
% f = node_error (:, i_p);
if ( i_p >= 1 )
  f = node_error(:, i_p) ;
else % i_p = 0, get root mean sq
  for k = 1:np
    f (k) = sqrt ( sum (node_error (k, 1:max_p).^2)) ;
  end % for k
end % if get RMS value
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
%    Initialize plots
xmax = max (x)      ; xmin = min (x) ;
ymax = max (y)      ; ymin = min (y) ;
axis ([xmin, xmax, ymin, ymax])   % set axes
axis ('equal')                    % true shape style
hold on
grid
% meshc (X, Y, F)
c = contour(X,Y,F);
fprintf ('Select contours to label with mouse using space bar to \n') 
fprintf ('enter contours and return to terminate: \n')
clabel (c, 'manual') ;

% draw the mesh
%                set constants
  [loop] = get_El_Loop (nod_per_el) ;

for i = 1:np;
  p_text = sprintf ('  %g', i);
  text (x(i), y(i), p_text) 
end % for all points

%    Loop over all elements
  for it = 1:nt  ;

%      Extract corner connectivity
    t_nodes = msh_typ_nodes (it, (pre_e+2):(nod_per_el+pre_e+1))%;

%      Extract corner coordinates
    t_x  = x (t_nodes) %;         % x at those nodes, only
    t_y  = y (t_nodes) %;         % y at those nodes, only

% est results from grad
    r4 = 0.25*( (t_x(4) - t_x(1)) * scp_node_ave_fluxes (t_nodes(1),1) ...
              + (t_y(4) - t_y(1)) * scp_node_ave_fluxes (t_nodes(1),2) ...
              + (t_x(4) - t_x(2)) * scp_node_ave_fluxes (t_nodes(2),1) ...
              + (t_y(4) - t_y(2)) * scp_node_ave_fluxes (t_nodes(2),1) )
%      Plot this polygon
    c_x  = t_x (loop) ;       % x for nod_per_el line polygon
    c_y  = t_y (loop) ;       % y for nod_per_el line polygon
    plot (c_x, c_y)           % plot nod_per_el lines
  end % for over all elements

if ( i_p >= 1 )
  title(['Matlab Smoothed Exact Error in Component\_', int2str(i_p), ...
         ':  ', int2str(nt), ' Elements, ', int2str(np),' Nodes']);
else % i_p = 0, get root mean sq
  title(['Matlab Smoothed Exact Error in RMS Value:  ', ...
        int2str(nt), ' Elements, ', int2str(np),' Nodes' ]);
end % if get RMS value
xlabel (['X: ', int2str(nt),' Elements'])
ylabel (['Y: ', int2str(np),' Nodes'])

%       -depsc -tiff            % for an eps version
  print ('-dpsc', ['contour_test_', int2str(i_p), '_on_mesh'])
  v_text = ['Created contour_test_', int2str(i_p),'_on_mesh.ps'] ;
  fprintf (1,'%s', v_text) ; fprintf (1, ' \n' )
% end of contour_test