% This function forms global stiffness matrix (plane stress/strain)
%
% Output: K - global stiffness matrix
%
% Input: IEN - element topology: numbering of control points
% P - coordinates of control points
% K - empty global stiffness matrix
% E - constitutive matrix
% t - thickness
% nx - number of elements in xi direction
% ny - number of elements in eta direction
% ncp - number of control points
% poly - polynomial order
function K=FormK_(IEN,P,K,E,t,nx,ny,ncp,poly)
[G,W]=Gauss(poly); % Call Gauss points and weights
e=1;
xi_v=(0:nx); % Vectors defining parameter space
eta_v=(0:ny);
for s=1:ny
for r=1:nx
IEN_e=nonzeros(IEN(e,:))'; % Element topology of current el.
eDof=[IEN_e IEN_e+ncp]; % eDof: first x, then y
k=0;
for g=1:size(G,1) % For each Gauss point
xi_n=G(g,1); % Gauss coord. reference element
eta_n=G(g,2); % Gauss coord. reference element
% Gauss coordinates in parameter space
xi=xi_v(r)+(xi_n+1)*(xi_v(r+1)-xi_v(r))/2;
eta=eta_v(s)+(eta_n+1)*(eta_v(s+1)-eta_v(s))/2;
[Nxi,dxi,Neta,deta]=BasisFunc(xi,eta,nx,ny,poly);
dxi=dxi./2; % Mapping of derivatives from
deta=deta./2; % parameter space to reference el.
Nxi=Nxi(r:r+poly); % Collect basis functions and
dxi=dxi(r:r+poly); % derivatives which support the el.
Neta=Neta(s:s+poly);
deta=deta(s:s+poly);
[J,dxy]=Jacobian_(Nxi,dxi,Neta,deta,P,IEN_e,poly);
B=zeros(3,2*length(IEN_e)); % B matrix
B(1,1:length(IEN_e))=dxy(1,:);
B(2,length(IEN_e)+1:2*length(IEN_e))=dxy(2,:);
B(3,1:length(IEN_e))=dxy(2,:);
B(3,length(IEN_e)+1:2*length(IEN_e))=dxy(1,:);
k=k+B'*E*B*t*det(J)*W(g); % Numerical integration of k
end
K(eDof,eDof)=K(eDof,eDof)+k;
e=e+1;
end
end
end
评论0