function sedamshl
%
% Function: sedamshl.m
% 
% Usage: sedamshl
% 
% Description:
%    Computes and plots the strain energy
%    damage indicators for shell geometry
%

% 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.                                 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

global HANDLES GEOM DAMAGEID

numdiv = str2num(get(HANDLES.hsemb(10),'string'));
nu = str2num(get(HANDLES.hsemb(17),'string'));
iopt = get(HANDLES.hsemb(12),'value');

if iopt==1,
	interpopt='spline';
elseif iopt==2,
	interpopt='cubic';	
end,

%  Check each SEMquad to make sure it has valid curvature


%     First, make a vector of coordinates in x and y directions

nquads = length(GEOM.SEMquads);

for iquad = 1:nquads,
   
   thisquad = GEOM.SEMquads{iquad}(2:end);
   
   nnodes = length(thisquad); % Nodes which define this quad
   
   rnod = zeros(nnodes,1);
   
   for inod = 1:nnodes,
      
      rnod(inod) = find(GEOM.Nodes(:,1) == thisquad(inod));
      
   end
   
   xvec = zeros(nnodes,nquads);
   yvec = zeros(nnodes,nquads);
   
   for i = 1:nnodes,
      
      xvec(i,iquad) = GEOM.Nodes(rnod(i),2);
      yvec(i,iquad) = GEOM.Nodes(rnod(i),3);
      
   end
   
end
   

for i = 1:nquads,
   
   %     Next, find the number of unique coordinates in each direction
   
   x_unique = unique(xvec(:,i));
   y_unique = unique(yvec(:,i));
   
   %     If either number is less than 3, we have invalid curvature
   
   if  length(x_unique) < 3,
      disp(['WARNING: SEMquad #',num2str(i),' does not have at least three points in the x direction(invalid curvature)'])
   end
   
   if  length(y_unique) < 3,
      disp(['WARNING: SEMquad #',num2str(i),' does not have at least three points in the y direction(invalid curvature)'])
   end
   
end



%
%  Compute the damage index for each subdivision of each element
%

[DImat,DIloc] = dishls(GEOM.Nodes,GEOM.SEMquads,DAMAGEID.undammode,DAMAGEID.dammode,DAMAGEID.respDOF,interpopt,numdiv,DAMAGEID.mode_anal,nu);

%
%  Set plot parameters
%

ngrad = numdiv;

inodal = 0;  % DI values are elemental

colormap cool

ptdamshl(GEOM.Nodes,GEOM.Shells,DImat,DAMAGEID.haxes,inodal,ngrad,DIloc);

%  Plot undeformed geometry
hold on

updtgeom(DAMAGEID.haxes,GEOM.Nodes)

hold off

damscale(DAMAGEID.hlegax,DImat);

axes(DAMAGEID.haxes)
DAMAGEID.damtitstr = damtitle(['Strain Energy Method - Shell Elements'],DAMAGEID.undampath,DAMAGEID.undamfile,...
   DAMAGEID.dampath,DAMAGEID.damfile,['Number of Divs:',num2str(DAMAGEID.numdiv),', ','Modes: ',num2str(DAMAGEID.mode_anal)]);
title(['Strain Energy Method - Shell Elements'])

rotate3d on