-
Notifications
You must be signed in to change notification settings - Fork 0
/
Ke_truss.m
60 lines (52 loc) · 1.59 KB
/
Ke_truss.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Ke_truss.m - November, 21 2002 %
% author: Jean H. Prevost %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Ke,ke,Qe]=Ke_truss(E,A,xn,ien,nen,ndf,nsd)
%------------------------------------------------------------------------
% Purpose:
% compute truss element stiffness
%
% Synopsis:
% [Ke,ke,Qe]=Ke_truss(E,A,xn,ien,nen,ndf,nsd)
%
% Variable Description:
% nen - number of nodes per element
% ndf - number of equations per node
% nsd - number of spatial dimensions
%------------------------------------------------------------------------
node1=ien(1); node2=ien(2);
% compute the length for each element
v=xn(:,node2)-xn(:,node1); L=norm(v);
% local stiffness
ke=zeros(2,2);
ke(:,:)=(E*A/L)*[[1 -1]; [-1 1]];
% Rotation matrix to global coordinate system
if ( nsd == 1) % 1D case
if (xn(1,node1) < xn(1,node2))Qe=1;
else Qe=-1;
end;
else
Qe=zeros(nen,ndf*nen);
v=v/norm(v); % normalize vector
if (nsd ==2) % 2D case
Qe(1,1)=v(1);
Qe(1,2)=v(2);
Qe(2,3)=v(1);
Qe(2,4)=v(2);
elseif (nsd ==3) % 3D case
Qe(1,1)=v(1);
Qe(1,2)=v(2);
Qe(1,3)=v(3);
Qe(2,4)=v(1);
Qe(2,5)=v(2);
Qe(2,6)=v(3);
end;
end;
% Global Stiffness Ke(i,j)
Ke=zeros(ndf*nen,ndf*nen);
if (nsd >1)
Ke(:,:)=Qe(:,:)'*ke(:,:)*Qe(:,:);
else
Ke(:,:)=Qe*ke(:,:);
end;