function [elarea,Nint,B_sym,B_skw] = sf_arrs(iele,nnpe,np,x,y)
% define the shape function and shape function derivative arrays

% find the nodes belonging to this element
I=np(iele,:);

% subtract x- and y-components of node 1
x1 = 0.0;% = x(I(1))-x(I(1)); 
x2 = x(I(2))-x(I(1));
x3 = x(I(3))-x(I(1));
y1 = 0.0;% = y(I(1))-y(I(1));
y2 = y(I(2))-y(I(1));
y3 = y(I(3))-y(I(1));

%calculate elemental area
elarea = (x2*y3 - x3*y2)/2;% = Area of the element
nfac   = elarea/3;         % = integration of s.f. over element
fac    = elarea*2;         % = factor used in finding s.f. derivatives
%if any area is negative, return an error and stop running
if elarea<=0 
    error('Negative element area!');
end
    
% Nint is the matrix used to form fe,
% = array of shape function integrals
Nint=nfac*[1 0 1 0 1 0;
           0 1 0 1 0 1];

%calculate s.f. derivatives (function listed below)
for ifor=1:nnpe
	dndx(ifor) = sfdx(y2,y3,fac,ifor);
	dndy(ifor) = sfdy(x2,x3,fac,ifor);
end

% the U vector is arranged x1 y1 x2 y2 x3 y3 etc.
% B_sym matrix is built accordingly.
% B_sym matrix multiplies displacements to get strains.
B_sym=[     dndx(1)      0  dndx(2)      0  dndx(3)      0 ;
             0  dndy(1)      0  dndy(2)      0  dndy(3);
   0.5*[dndy(1) dndx(1) dndy(2) dndx(2) dndy(3) dndx(3)]];

% B_skw matrix multiplies displacements to get rotation.
B_skw = [ 0 0 0 0 0 0 ; 
         0 0 0 0 0 0 ; 
   0.5*[dndy(1) -dndx(1) dndy(2) -dndx(2) dndy(3) -dndx(3)]];
 
function dsfdx = sfdx(y2,y3,fac,i)
%dndx function statements
switch i
 case 1
  dsfdx = -(y3-y2)/fac;
 case 2
  dsfdx = y3/fac;
 case 3
  dsfdx = -y2/fac;
end
function dsfdy = sfdy(x2,x3,fac,i)
%  dndy function statement
switch i
 case 1
  dsfdy = (x3-x2)/fac;
 case 2
  dsfdy = -x3/fac;
 case 3
  dsfdy = x2/fac;
end
