function  contour_result_on_mesh (i_p, Type, row) %CL
Pts_wide=2;
%b function  contour_result_on_mesh (i_p, Type)
% 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 ;
pre_p = 1;
pre_r = 0 ; %CL
if ( nargin == 0 )
  i_p = 0 ; Type = 0 ; row = 1 ;
elseif ( nargin == 1 )
  Type = 0 ; row = 1 ;
elseif ( nargin == 2 )
  row = 1 ; pre_r = 2 ; %CL
elseif ( nargin == 3 )
  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  msh_typ_nodes.tmp ; % nod_per_el nodes per element
nt = size (msh_typ_nodes,1);           % number of elements
nt = nt - less ; %CL
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
Col_1 = pre_e+2 ; %CL
Col_2 = nod_per_el+pre_e+1 ; %CL

  el_type (nt) = 1       ;       % pre-allocate array el_type
  Max_type = 1 ; %CL
  if ( nargin < 3 ) % mixed types possible %CL
%   get the element types %CL
    el_type = msh_typ_nodes (:, pre_e+1) ; %CL
    Max_type = max (el_type) ; %CL
    fprintf (' Maximum element type number is %g \n', Max_type) %CL
  end % if type %CL

%   count various types
  Type_count (1:Max_type) = 0 ;
  if ( Max_type > 1 )
    for j = 1:nt
      k = el_type (j) ;
      Type_count ( k ) = Type_count ( k ) + 1 ;
    end % for j
    for j = 1:Max_type
      fprintf ('Type %g has %g elements \n', j, Type_count (j))
    end % for j
  end % if Max_type

load node_results.tmp
nr = size (node_results, 1) ; 
max_p = size (node_results, 2) - pre_r; % number of columns %CL
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

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

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

%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_N, L_N] = min (f) ;
[V_X, L_X] = max (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');
%b F = TriScatteredInterp (x, y, f) ; % new method

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)
Pts_wide=3 ;
c = contour(X,Y,F,'LineWidth',Pts_wide);
clabel (c);
%b hold on 
%b if ( 0==0 )  %b
  %b error 'manual stop' %b
%b end % if manual %b

%    Loop over all elements
  % Type = 0 ; % debugging
  Old_type = 0;
  for it = 1:nt  ;
   This_type = el_type (it) ;
   if ( This_type == Type | Type == 0 ) % then plot it

%      Extract corner connectivity
    %b t_nodes = msh_typ_nodes (it, (pre_e+2):(nod_per_el+pre_e+1));
    t_nodes = msh_typ_nodes ((it+less), Col_1:Col_2 ) ; %CL

%      Get the type, nodes per type and loop array
     if ( Old_type ~= This_type  ) % get type
       Old_type = This_type    ;
       nod_per_t = nod_per_el ;
       [Low_node, Low_I] = min (t_nodes) ;
       if ( Low_node <= 0 )
         nod_per_t = Low_I - 1 ;
       end % if shorter topology list
       type_plus = nod_per_t + 1 ;

%                set constants
       clear loop;
       [loop] = get_El_Loop(nod_per_t) ;
     end % get type

%      Extract corner coordinates
    t_x  = x (t_nodes (1:nod_per_t)) ; % x at those nodes, only
    t_y  = y (t_nodes (1:nod_per_t)) ; % y at those nodes, only

%      Plot this polygon
    c_x  = t_x (loop(1:type_plus)); % x for nod_per_el line polygon
    c_y  = t_y (loop(1:type_plus)); % y for nod_per_el line polygon
    plot (c_x, c_y)           % plot nod_per_el lines
   end % if plot This_type
  end % for over all elements, nt

if ( i_p >= 1 )
  title(['FEA Component\_', int2str(i_p), ...
         '  (', num2str(V_N), ' to ', num2str(V_X),')'])
else % i_p = 0, get root mean sq
  title(['FEA RMS Value',...
         '  (', num2str(V_N), ' to ', num2str(V_X),')'])
end % if get RMS value
xlabel (['X: ', int2str(nt),' Elements with ', ...
        int2str(nod_per_el), ' nodes'])
ylabel (['Y on ', int2str(nt),' Elements']);

%       -depsc -tiff            % for an eps version
  %bprint ('-dpsc', ['contour_result_', int2str(i_p), '_on_mesh'])
  %bv_text = ['Created contour_result_', int2str(i_p),'_on_mesh.ps'] ;
  %bfprintf (1,'%s', v_text) ; fprintf (1, ' \n' )
%b end % of contour_result_on_mesh