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.

% 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;
E     = propstruct.youngmodulus;
nu    = propstruct.poissonratio;
bodyf = propstruct.bodyforce;

% 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

% 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
    [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
U = unrotateDOF(numnp,U,nangle);                   % return to (x,y) components

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

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