function StaticElast(fighandle)
% StaticElast(FIGHANDLE)
% This function accepts the handle FIGHANDLE to the main GUI figure and
% runs the 2-D static elasticity finite element problem defined there.
% 3-noded elements, force or displacement boundary conditions, body forces,
% uniform temperature change of the body, and plane stress or plane strain
% are all supported.

%%%%%%%%%%%%%%%%%%%%
%%%%% DO NOT CHANGE vvvvvvvvvv
%%%%%%%%%%%%%%%%%%%%

% get handles struct
handles=guidata(fighandle);

% define solution method as static
handles.solutionexists = 'Displ. done.';


% unpack necessary variables
meshstruct = handles.mesh;
meshdata=meshstruct.size;
numels = meshdata(1);
numnp = meshdata(2);
nnpe = meshdata(3);
numeq = meshdata(4);
ndof = meshdata(5);
x=meshstruct.xcoords;
y=meshstruct.ycoords;
np=meshstruct.connect;
edges = meshstruct.edges;

npts = nnpe;

bcinfostruct = handles.bcsi;
bccode = bcinfostruct.bccodes;
nangle = bcinfostruct.bcnangles;
boundt = bcinfostruct.bcboundts;
boundn = bcinfostruct.bcboundns;

propstruct = handles.prop;
ipsos = propstruct.twoDflag;
dT    = propstruct.deltatemp;
E     = propstruct.youngmodulus;
nu    = propstruct.poissonratio;
bodyf = propstruct.bodyforce;
alpha = propstruct.thermalcoeff;
if isfield(propstruct,'anisomods')
    anisomods = propstruct.anisomods;
end
%%%%%%%%%%%%%%%%%%%%
%%%%% DO NOT CHANGE ^^^^^^^
%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%
%%%%% BEGIN SPARSE 
%%%%%%%%%%%%%%%%%%%%
% initialize the problem
Fsr=zeros(numels*2*nnpe,1); Fsv=zeros(numels*2*nnpe,1); 
Ksr=zeros(numels*4*nnpe*nnpe,1); Ksc=zeros(numels*4*nnpe*nnpe,1); 
Ksv=zeros(numels*4*nnpe*nnpe,1);
fe     = zeros(1,ndof);    % elemental force vector
fer    = zeros(1,ndof);    % rotated elemental force vector
ke     = zeros(ndof);      % elemental stiffness matrix
ker    = zeros(ndof);      % rotated elemental stiffness matrix
elarea = 0.0  ;            % elemental area
Nint   = zeros(2,6);       % Shape function array
B      = zeros(3,6);       % shape function derivative array
r=zeros(3);

% h is a correction factor to produce the right shear strain energy
% term in the matrix multiplication.
h=[ 1 0 0; 0 1 0; 0 0 2];
if E<0
    r = matlb(ipsos,nu,anisomods); % define material stiffness matrix
else
    r = matlb(ipsos,nu,E); % define material stiffness matrix
end

% now loop over the elements
fc=1; kc=1;
for iele =1:numels

    [elarea,Nint,B] = sf_arrs(iele,nnpe,np,x,y);  % find shape fcn, s.f. derivative arrays
    ke = elarea*B'*h*r*B;                         % form elemental stiffness matrix
    fe = Nint'*bodyf;                             % form elemental force vector
    fet= elarea*alpha*B'*r*dT*[1;1;0];            % form thermal pseudo force
    fe = fe+fet;                                  % add pseudo force to fe
    [ker,fer] = rotateDOF(iele,nnpe,np,nangle,ke,fe);   % move (x,y) comp's to (norm,tang) comp's
    % perform stiffness matrix and force vector assembly
    for  i=1:nnpe
        i2=2*np(iele,i); % find d.o.f. corresponding to y- or norm-component
        i1=i2-1;         % find d.o.f. corresponding to x- or tang-component
        Fsr(fc)=i1; Fsv(fc)=fer(2*i-1); fc=fc+1;
        Fsr(fc)=i2; Fsv(fc)=fer(2*i);   fc=fc+1;
        for j=1:nnpe
            j2=2*np(iele,j);
            j1=j2-1;
            Ksr(kc)=i1; Ksc(kc)=j1; Ksv(kc)=ker(2*i-1,2*j-1); kc=kc+1;
            Ksr(kc)=i1; Ksc(kc)=j2; Ksv(kc)=ker(2*i-1,  2*j); kc=kc+1;
            Ksr(kc)=i2; Ksc(kc)=j2; Ksv(kc)=ker(2*i  ,  2*j); kc=kc+1;
            Ksr(kc)=i2; Ksc(kc)=j1; Ksv(kc)=ker(2*i  ,2*j-1); kc=kc+1;
        end
    end
    % assembly complete
end
Fs=sparse(Fsr,1,Fsv);
Ks=sparse(Ksr,Ksc,Ksv);
[Ks,Fs] = applybc(numnp,bccode,boundn,boundt,Ks,Fs);   % apply boundary conditions

Us=full(Ks\Fs);
%%%%%%%%%%%%%%%%%%%%
%%%%% END SPARSE 
%%%%%%%%%%%%%%%%%%%%
% % % % % % % %%%%%%%%%%%%%%%%%%%%
% % % % % % % %%%%% BEGIN ORIGINAL METHOD
% % % % % % % %%%%%%%%%%%%%%%%%%%%
% % % % % % % % initialize the problem
% % % % % % % K      = zeros(numeq);     % global stiffness matrix
% % % % % % % U      = zeros(1,numeq);   % global nodal displacements
% % % % % % % F      = zeros(1,numeq);   % global force vector
% % % % % % % fe     = zeros(1,ndof);    % elemental force vector
% % % % % % % fer    = zeros(1,ndof);    % rotated elemental force vector
% % % % % % % ke     = zeros(ndof);      % elemental stiffness matrix
% % % % % % % ker    = zeros(ndof);      % rotated elemental stiffness matrix
% % % % % % % elarea = 0.0  ;            % elemental area
% % % % % % % Nint   = zeros(2,6);       % Shape function array
% % % % % % % B      = zeros(3,6);       % shape function derivative array
% % % % % % % r=zeros(3);

% % % % % % % % h is a correction factor to produce the right shear strain energy
% % % % % % % % term in the matrix multiplication.
% % % % % % % h=[ 1 0 0; 0 1 0; 0 0 2];
% % % % % % % r = matlb(ipsos,nu,E); % define material stiffness matrix
% % % % % % % junk='regular method'
% % % % % % % tic
% % % % % % % % now loop over the elements
% % % % % % % for iele =1:numels
% % % % % % % 
% % % % % % %     [elarea,Nint,B] = sf_arrs(iele,nnpe,np,x,y);  % find shape fcn, s.f. derivative arrays
% % % % % % %     ke = elarea*B'*h*r*B;                         % form elemental stiffness matrix
% % % % % % %     fe = Nint'*bodyf;                             % form elemental force vector
% % % % % % %     fet= elarea*alpha*B'*r*dT*[1;1;0];            % form thermal pseudo force
% % % % % % %     fe = fe+fet;                                  % add pseudo force to fe
% % % % % % %     [ker,fer] = rotateDOF(iele,nnpe,np,nangle,ke,fe);   % move (x,y) comp's to (norm,tang) comp's
% % % % % % %     % perform stiffness matrix and force vector assembly
% % % % % % %     for  i=1:nnpe
% % % % % % %         i2=2*np(iele,i); % find d.o.f. corresponding to y- or norm-component
% % % % % % %         i1=i2-1;         % find d.o.f. corresponding to x- or tang-component
% % % % % % %         F(i1)=F(i1)+fer(2*i-1);
% % % % % % %         F(i2)=F(i2)+fer(2*i);
% % % % % % %         for j=1:nnpe
% % % % % % %             j2=2*np(iele,j);
% % % % % % %             j1=j2-1;
% % % % % % %             K(i1,j1)    =   K(i1,j1)  +  ker(2*i-1,2*j-1);
% % % % % % %             K(i1,j2)    =   K(i1,j2)  +  ker(2*i-1,  2*j);
% % % % % % %             K(i2,j2)    =   K(i2,j2)  +  ker(2*i  ,  2*j);
% % % % % % %             K(i2,j1)    =   K(i2,j1)  +  ker(2*i  ,2*j-1);
% % % % % % %         end
% % % % % % %     end
% % % % % % %     % assembly complete
% % % % % % % end
% % % % % % % [K,F] = applybc(numnp,bccode,boundn,boundt,K,F);   % apply boundary conditions
% % % % % % % U=F/K;                                             % solve the system of equations
% % % % % % % toc
% % % % % % % %%%%%%%%%%%%%%%%%%%%
% % % % % % % %%%%% END ORIGINAL METHOD
% % % % % % % %%%%%%%%%%%%%%%%%%%%
U = unrotateDOF(numnp,Us,nangle)';                   % return to (x,y) components

% create the output structure
displstruct = struct('displ', U);
handles.displ = displstruct;

% update the handles structure
guidata(fighandle,handles);
