function  contour_result (i_p, row) %CL
Pts_wide=2;
%b function  contour_result (i_p)
% 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
% pre_e  = Element items before connectivity list
% pre_p  = Nodal items before coordinates
% pre_r  = Nodal items before results %CL
% row    = row number where cord & connectivity start, 1 or 2 %CL

% clear
pre_e  = 0 ; % items before type, connectivity list
pre_p = 1;
pre_r = 0 ; %CL
if ( nargin == 0 )
  i_p = 0 ; row = 1 ; %CL
elseif ( nargin == 1 )
  row = 1 ; %CL
elseif ( nargin == 2 )
  pre_r = 2 ; %CL
end % if no arguments
less = row - 1 ; %CL

load msh_bc_xyz.tmp;
np = size (msh_bc_xyz, 1);
np = np - less ;   %CL
fprintf ('Read %g mesh coordinate pairs \n', np)
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 nodal solution values \n', nr)
fprintf ('     with %g components each \n', max_p)
if ( i_p > max_p )
  error ('i_p > available data')
end % if error

  load  msh_typ_nodes.tmp ; % nod_per_el nodes per element
  nt = size (msh_typ_nodes,1) ; % number of elements in mesh 
  nod_per_el = size(msh_typ_nodes,2) - pre_e -1; %nodes per el

% fprintf ('Begin component %g value carpet plots: \n', i_p)

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

%b x = msh_bc_xyz (:, 2); y = msh_bc_xyz (:, 3);
x = msh_bc_xyz (row:(np+less), (pre_p+1)) ; % extract x  %CL
y = msh_bc_xyz (row:(np+less), (pre_p+2)) ; % extract y %CL
 
% f = node_results (:, i_p);
if ( i_p >= 1 )
  f = node_results(:, pre_r+i_p) ; %CL
else % i_p = 0, get root mean sq
  for k = 1:np
    f (k) = sqrt(sum(node_results (k, pre_r+1:pre_r+max_p).^2)) ; %CL
  end % for k
end % if get RMS value
V_X = max(f); V_N=min(f);
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');
% F = griddata (x, y, f, X, Y, 'linear') ;
F = griddata (x, y, f, X, Y, 'nearest') ;
% F = griddata (x, y, f, X, Y, 'v4') ;

clf
hold on
xmax = max (x)      ; xmin = min (x) ;
ymax = max (y)      ; ymin = min (y) ;
%b ymin=-0.19
%b ymax = 2.19
axis ([xmin, xmax, ymin, ymax])   % set axes
% axis('equal')
grid
% meshc (X, Y, F)
c = contour(X,Y,F,'LineWidth',Pts_wide);
% V=[0.25, 0.5, 0.75, 1., 1.25, 1.5, 1.75]
% V=[1., 2., 3., 4., 5.]
% V=[0.5,1.,1.5,2.,2.5,3.,3.5,4.,4.5,5.,5.5] ;
% V=[0.2,0.4,0.6,0.8,1.0] ;
  % clabel (c,V);
  clabel (c);
if ( i_p >= 1 )
  title(['FEA Solution Component\_', int2str(i_p), ...
         ' at ',int2str(np),' Nodes']);
else % i_p = 0, get root mean sq
  title(['FEA Solution RMS Value at ', ...
         int2str(np),' Nodes']);
end % if get RMS value
% xlabel (['X: ', int2str(np),' Nodes'])
xlabel (['X Coordinate for ', int2str(nt),' Elements with ', ...
        int2str(nod_per_el), ' nodes'])
ylabel (['Y Coordinate   (Contour min : max  ', num2str(V_N), ' : ', ...
         num2str(V_X), ') ']);
%       -depsc -tiff            % for an eps version
  % print ('-dpsc', ['contour_result_', int2str(i_p)])
  % v_text = ['Created contour_result_', int2str(i_p),'.ps'] ;
  % fprintf (1,'%s', v_text) ; fprintf (1, ' \n' )
% end of contour_result