function Efrac = frcenshl(Nodes,shlnodes,phi,interpopt,numdiv,mode_anal,nu,xi_int,yi_int,pxx,pxy,pyy)
%
%  Name:  frcenbar.m
%
%  Usage: 
%
%    Calculates the fractional strain energy for a shell element
%
%  Input:
%  Nodes = List of nodes in DIAMOND geometry format
%  shlnodes = 4-element vector containing indices of nodes for this shell
%	phi = Mode shape value at each end of the bar element, for each mode
%  interpopt = 'linear','spline', or 'cubic'
%  numdiv = number of subdivisions for the bar
%  mode_anal=vector of modes to use
%
%  Output:
%	Etot(i)=total strain energy in the element for mode i
% 	Efrac(i,j)=energy in the jth subdivision of the element for the ith mode
%

%  Version SWD970403
%
%  Author(s): Phillip Cornwell, cornwell@lanl.gov
%	modified by Scott W. Doebling, doebling@lanl.gov
%	
%  Institution: Los Alamos National Laboratory,
%               ESA-EA
%
%  Copyright (c) 1996-1997, The Regents of the University
%       of California. (See top level of code for
%       full copyright statement.)
%

%
% Define the parameters for the numerical integration
%

nn=3;  % Numerically integrate at 5 points across each element subdivision

%
% Parameters for Node 1 - Node 2 direction
%

%L12 = norm(Nodes(shlnodes(2),2:4) - Nodes(shlnodes(1),2:4));
%x12 = linspace(0,L12,(numdiv*nn));
%dx12 = x12(2) - x12(1);

%
% Node 1 - Node 4 direction
%

%L14 = norm(Nodes(shlnodes(4),2:4) - Nodes(shlnodes(1),2:4));
%x14 = linspace(0,L14,(numdiv*nn))';
%dx14 = x14(2) - x14(1);



% Initialize Efrac

%Efrac = zeros(length(mode_anal),numdiv,numdiv);

%
% find derivatives
%

%%
%%  Determine which DOF in Nodes_int are associated with this element
%%

   x1int = linspace(Nodes(shlnodes(1),2),Nodes(shlnodes(2),2),numdiv+1);
   y1int = linspace(Nodes(shlnodes(1),3),Nodes(shlnodes(2),3),numdiv+1);
   
   x2int = linspace(Nodes(shlnodes(2),2),Nodes(shlnodes(3),2),numdiv+1);
   y2int = linspace(Nodes(shlnodes(2),3),Nodes(shlnodes(3),3),numdiv+1);
   
   x3int = linspace(Nodes(shlnodes(3),2),Nodes(shlnodes(4),2),numdiv+1);
   y3int = linspace(Nodes(shlnodes(3),3),Nodes(shlnodes(4),3),numdiv+1);
   
   x4int = linspace(Nodes(shlnodes(4),2),Nodes(shlnodes(1),2),numdiv+1);
   y4int = linspace(Nodes(shlnodes(4),3),Nodes(shlnodes(1),3),numdiv+1);
   
   inodeint = [];
   
   for in = 1:numdiv,
      
         [in1i,in1j] = find(xi_int == x1int(in) & yi_int == y1int(in));
         [in2i,in2j] = find(xi_int == x2int(in) & xi_int == y2int(in));
         [in3i,in3j] = find(xi_int == x3int(in) & xi_int == y3int(in));
         [in4i,in4j] = find(xi_int == x4int(in) & xi_int == y4int(in));
         
         inodeint = [inodeint;[in1i,in1j];[in2i,in2j];[in3i,in3j];[in4i,in4j]];
         
   end

for imod = 1:length(mode_anal),
  
%   [x12i,x14i,phii]=griddata([0 L12 L12 0],...
%      [0 0 L14 L14],phi(:,mode_anal(imod)),x12,x14);
   
%   [px,py]=gradient(phii,dx12,dx14);
%   [pxx,pxy]=gradient(px,dx12,dx14);
%   [pyx,pyy]=gradient(py,dx12,dx14);

%
% Determine the integrand
%

fx=pxx(:,:,imod).^2+pyy(:,:,imod).^2+...
   2*nu*pxx(:,:,imod).*pyy(:,:,imod)+2*(1-nu)*pxy(:,:,imod).^2;

%
%  Do the integration
%

for i=1:numdiv
   for j=1:numdiv
      xvalues=x12((i-1)*nn+1:i*nn);
      yvalues=x14((j-1)*nn+1:j*nn);
      zvalues=fx((j-1)*nn+1:j*nn,(i-1)*nn+1:i*nn);
      temp = zeros(1,nn);
      for k=1:nn
         temp(k)=trapz(yvalues,zvalues(:,k));
      end
      Efrac(imod,i,j)=trapz(xvalues,temp);
      
   end
end


end

return
