function  contour_exact_error_pc (i_p)
Pts_wide=2;
% Copyright 2000, J.E. Akin.  All rights reserved.
% plot finite element result 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 )
  error ' no argument' % 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

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

load exact_node_solution.tmp
nrr = size (exact_node_solution, 1);
max_pp = size (exact_node_solution, 2) ; % number of columns
if ( nrr == 0 )
  error ('Error missing file exact_node_solution.tmp')
end % if error
fprintf ('Read %g exact nodal value sets \n', nrr)
fprintf ('     with %g components each \n', max_pp)

if ( nr ~= nrr | max_p ~= max_pp )
  error ('FEA and exact results not the same sizes')
end % if same sizes

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

%b fprintf ('Begin component %g exact contour: \n', i_p)

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

x = msh_bc_xyz (:, 2); y = msh_bc_xyz (:, 3);
 
node_error = node_results - exact_node_solution ;

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
for j=1:np
  if ( exact_node_solution (:, i_p) ~= 0. ) 
    f (j) = f (j) / exact_node_solution (j, i_p) * 100.;
  else 
    f (j) = 0. ;             
  end % if 
end % for

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
hold on
xmax = max (x)      ; xmin = min (x) ;
ymax = max (y)      ; ymin = min (y) ;
axis ([xmin, xmax, ymin, ymax])   % set axes
axis('equal')
grid
% meshc (X, Y, F)
c = contour(X,Y,F,'LineWidth',Pts_wide);
clabel (c);
xlabel (['X: ', int2str(np),' Nodes'])
ylabel ('Y');
if ( i_p >= 1 )
  title(['Smoothed Exact Error % in Solution Component\_', ...
        int2str(i_p) ]); 
else
   title(['Smoothed Exact Solution Error RMS Value'])
end % if

%       -depsc -tiff            % for an eps version
  %b print ('-dpsc', ['contour_exact_error_', int2str(i_p)])
  %b v_text = ['Created contour_exact_error_', int2str(i_p),'.ps'] ;
  %b fprintf (1,'%s', v_text) ; fprintf (1, ' \n' )
hold off
% end of contour_exact_error_pc