function  contour_pt_ave_ex_err_mesh ()
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);
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 pt_ave_ex_error.tmp
nr = size (pt_ave_ex_error, 1);
if ( nr == 0 )
  error ('Error missing file pt_ave_ex_error.tmp')
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);
 
f = 100. * pt_ave_ex_error (:);
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);

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

%    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

%      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

xlabel (['X: ', int2str(nt),' Elements (with ', ...
        int2str(nod_per_el), ' nodes)'])
ylabel (['Y: ', int2str(np),' Nodes'])
%xlabel ('X'); ylabel ('Y');
title (['Smoothed FEA Exact Energy Norm ', ...
       'Error Averaged at Nodes, % * 100']);

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