function  contour_pt_ave_error_est ()
Pts_wide=2;
% Copyright 2000, J.E. Akin.  All rights reserved.
% plot finite element ave node error est as a surface 
% 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

load msh_bc_xyz.tmp;
np = size (msh_bc_xyz, 1);

load pt_ave_error_est.tmp
nr = size (pt_ave_error_est, 1);
if ( nr == 0 )
  error ('Error missing file pt_ave_error_est.tmp')
end % if error
 
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	

%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);
 
f = 100. * pt_ave_error_est (:);
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(nt),' Elements (with ', ...
        int2str(nod_per_el), ' nodes)'])
ylabel (['Y: ', int2str(np),' Nodes'])
title (['FEA SCP Energy Norm ', ...
       'Error Estimate Averaged at Nodes, % * 100']);

%     -depsc -tiff            % for an eps version
%bprint -dpsc contour_pt_ave_error_est
hold off
%bfprintf ('Created file contour_pt_ave_error_est.ps \n')
% end of contour_pt_ave_error_est