function haxes = ptdamshl(Nodes,Shells,DImat,haxes,inodal,ngrad,DIloc)
%
%  Name: ptdamshl
% 
%  Usage: ptdamshl(Nodes,Bars,DImat,haxes,inodal,ngrad,naxaz,naxel)
%
% 
%  Description:
%     Plot damage indicator DI for shell geometry
%     For nodal DI - use color gradient with ngrad increments
%     For elemental DI - use solid color for each element subdivision
%
%  Inputs:
%
%	Nodes   = List of nodes
%  Shells  = List of shell connectivity
%  DImat   = Matrix of DI values -- corresponding to rows of either Nodes or Bars
%             For elemental DI's: DImat second dim. corresponds to values of subelements 
%             in 1-2 direction, DImat third dim. corresponds to values of subelements in
%             1-4 direction)
%
%  haxes   = Handle of axes for plot
%  inodal  = 0, DI are elemental values
%          = 1, DI are nodal values
%  ngrad   = Number of gradient steps (in each direction)for each shell element 
%            (or each elemental subdivision, if applicable)
%  naxaz, naxel = 3-D viewing parameters

%  Version SWD970916
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This matlab source code was originally     %
% developed as part of "DIAMOND" at          %
% Los Alamos National Laboratory. It may     %
% be copied, modified, and distributed in    %
% any form, provided:                        %
%  a) This notice accompanies the files and  %
%     appears near the top of all source     %
%     code files.                            %
%  b) No payment or commercial services are  %
%     received in exchange for the code.     %
%                                            %
% Original copyright is reserved by the      %
% Regents of the University of California,   %
% in addition to Scott W. Doebling, Phillip  %
% J. Cornwell, Erik G. Straser, and Charles  %
% R. Farrar.                                 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


if ~exist('haxes')
	haxes = axes;
elseif ~ishandle(haxes),
	haxes = axes;
end

axes(haxes)

cla

%  Draw the undeformed geometry

updtgeom(haxes,Nodes)

%
%  Plot the damage indicator
%

%  Define the colors from the colormap

map = colormap;

%  Compute full range of DI

dimin = min(DImat(:));
dimax = max(DImat(:));

distep = (dimax - dimin) / (size(map,1)-1);

%  Map the colors to the DI values

%% ===================================
%%
%%    NODAL DAMAGE INDICATOR VALUES
%%
%% ===================================

if inodal == 1,
   
%  Loop over each bar element

   for ishl = 1:size(Shells,1),
   
      %  Select Nodes associated with that element

      n1 = Shells(ishl,2);
      n2 = Shells(ishl,3);
      n3 = Shells(ishl,4);
      n4 = Shells(ishl,5);

      %  Find locations of end-nodes in Node list

      in1 = find(Nodes(:,1) == n1);
      in2 = find(Nodes(:,1) == n2);
      in3 = find(Nodes(:,1) == n3);
      in4 = find(Nodes(:,1) == n4);

      %  Determine the coordinates and DI values at the nodes

      loc1 = Nodes(in1,2:4);
      loc2 = Nodes(in2,2:4);
      loc3 = Nodes(in3,2:4);
      loc4 = Nodes(in4,2:4);

      DI1  = DImat(in1);
      DI2  = DImat(in2);
      DI3  = DImat(in3);
      DI4  = DImat(in4);
   
      %  Define an interpolation vector along edges 14 and 23
      
      edge14 = linspace(DI1,DI4,ngrad+1);
      edge23 = linspace(DI2,DI3,ngrad+1);
      
      %  Define a matrix with all of the DI values
      
      DIsubmat = zeros(ngrad+1,ngrad+1);
      
      for i = 1:(ngrad+1),
         DIsubmat(:,i) = linspace(edge14(i),edge23(i),ngrad+1)';
      end   
      
      %  Compute the geometric and color increments
      
      geoinc12 = (loc2 - loc1) / ngrad;
      geoinc23 = (loc3 - loc2) / ngrad;
   
      diinc12 = (DI2 - DI1) / ngrad;
      diinc13 = (DI3 - DI1) / ngrad;
      diinc14 = (DI4 - DI1) / ngrad;
   
      di1old = DI1;
      di2old = DI2;
      di3old = DI3;
      di4old = DI4;
   
      diavg = (DI1 + DI2 + DI3 + DI4)/4;
   
      p1 = loc1;
      p2 = loc2;
      p3 = loc3;
      p4 = loc4;
   
      %  Loop over the number of gradient increments in each direction
   
      for igrad14 = 1:ngrad,

         for igrad12 = 1:ngrad,
         
            %  Determine the coordinates of the corners of the patch

            p1 = loc1 + (igrad12-1) * geoinc12 + (igrad14-1) * geoinc23;
            p2 = loc1 + igrad12 * geoinc12 + (igrad14-1) * geoinc23;
            p3 = loc1 + igrad12 * geoinc12 + igrad14 * geoinc23;
            p4 = loc1 + (igrad12-1) * geoinc12 + igrad14 * geoinc23;
            
            di1 = DIsubmat(igrad12,igrad14);
            di2 = DIsubmat(igrad12+1,igrad14);
            di3 = DIsubmat(igrad12+1,igrad14+1);
            di4 = DIsubmat(igrad12,igrad14+1);
            
            diavg = (di1+di2+di3+di4)/4;
      
            icolor = floor((diavg - dimin) / distep) + 1;

            segcolor = map(icolor,:);
      
            %  Draw the line segment with the appropriate color
  
            h = patch([p1(1) p2(1) p3(1) p4(1)], [p1(2) p2(2) p3(2) p4(2)],...
                [p1(3) p2(3) p3(3) p4(3)]);
      
            set(h,'facecolor',segcolor,'linestyle','none')
             
         end

      end

   end

%% ===================================
%%
%%  ELEMENTAL DAMAGE INDICATOR VALUES
%%
%% ===================================

elseif inodal == 0,
   
   %  Loop over each bar element

   for iquad = 1:size(DImat,1),
      
      for jseg = 1:size(DImat,2),
         
         for iseg = 1:size(DImat,3),
            
            %  Determine the x and y coordinates of each segment
            
            xloc = DIloc{iquad}.xloc{jseg,iseg};
            yloc = DIloc{iquad}.yloc{jseg,iseg};
            
            %   Determine the color of the segment
            
            icolor = floor((DImat(iquad,jseg,iseg) - dimin) / distep) + 1;
            
            segcolor = map(icolor,:);
            
            %   Draw a patch of the appropriate color
            
            h = patch([xloc(1) xloc(2) xloc(2) xloc(1)], [yloc(1) yloc(1) yloc(2) yloc(2)],...
               [0 0 0 0]);
            
            set(h,'facecolor',segcolor,'linestyle','none')
            
      
         end
         
      end
      
   end
   
end

set(haxes,'dataaspectratio',[1 1 1])
rotate3d on

%  Define axis labels

xlabel('X'),ylabel('Y'),zlabel('Z')

return