function  smooth_exact_solution (i_p)
% Copyright 2000 J.E. Akin. All rights reserved.
% plot finite element nodal exact values 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

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_pp)
    error ('i_p > available data')
  end % if error

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

x = msh_bc_xyz (:, 2); y = msh_bc_xyz (:, 3);

if ( i_p >= 1 )
  f = exact_node_solution (:, i_p) ;
else % i_p = 0, get root mean sq
  for k = 1:np
    f (k) = sqrt ( sum (exact_node_solution (k, 1:max_pp).^2)) ;
  end % for k
end % if get RMS value

%    Cite max, min values
[V_X, L_X] = max (f) ;
[V_N, L_N] = min (f) ;
fprintf ('Max value is %g at node %g \n', V_X, L_X)
fprintf ('Min value is %g at node %g \n', V_N, L_N)
 
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
meshc (X, Y, F)
hold on
if ( i_p >= 1 )
  zlabel (['Component ', int2str(i_p), ' (max = ', ...
  num2str(V_X), ', min = ', num2str(V_N), ')'])
  title(['Matlab Smoothed Exact Value of Solution Component\_', ...
  int2str(i_p),': ',int2str(np),' Nodes'])
  %int2str(i_p),': ',int2str(nt),' Elements, ',int2str(np),' Nodes'])
          % int2str(i_p) ]); 
else % i_p = 0, get root mean sq
  zlabel (['RMS Value (max = ', ...
  num2str(V_X), ', min = ', num2str(V_N), ')'])
  title(['Matlab Smoothed RMS Value of the Exact Solution: ', ...
          int2str(np),' Nodes'])
          %int2str(nt),' Elements, ', int2str(np),' Nodes'])
  %title(['Matlab Smoothed RMS Value of the Exact Solution'])
end % if get RMS value
xlabel ('X'); ylabel ('Y');

% label max min points
  v_text = sprintf ('----min') ;
  text  (x(L_N), y(L_N), V_N, v_text) ;
  v_text = sprintf ('----max') ;
  text  (x(L_X), y(L_X), V_X, [v_text]) ;

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

hold off
% end of smooth_exact_solution