Codes

Feb. 25th, 2009 11:58 pm
nikavia: (Default)
[personal profile] nikavia
SOOO i did alot of program work today, but it is amazing. no really here is the beautiful code:



function nyssagilkey2(file_name);

% OPEN DATA FILE

fid = fopen(file_name, 'r');

% READ DATA - READ NODE DATA - First read the number of nodes in the truss

Number_nodes=fscanf(fid, '%d', 1);

% Read the coordinates for each node

for i=1:Number_nodes
Node=fscanf(fid, '%d', 1);
Nodea(i,1)=Node;
Coordinate(Node, 1)=fscanf(fid, '%g', 1);
Coordinate(Node, 2)=fscanf(fid, '%g', 1);
end

% READ DATA - READ ELEMENT DATA - Now read the number of elements and the element definition

Number_elements=fscanf(fid, '%d', 1);
M=zeros(2*Number_nodes, 2*Number_nodes);

for i=1:Number_elements
Element=fscanf(fid, '%d', 1);
Node_from=fscanf(fid, '%d', 1);
Node_to=fscanf(fid, '%d', 1);
Area=fscanf(fid, '%d', 1);
Young=fscanf(fid, '%d', 1);
Elementa(i,1)=Element;
DOFa(i,1)=Node_from;
DOFa(i,2)=Node_to;
Areamap(i,1)=Area;
Youngmap(i,1)=Young;
angles(i,1)=atand((Coordinate(Node_to,2)-Coordinate(Node_from,2))/(Coordinate(Node_to,1)-Coordinate(Node_from,1)));
ca= cosd(angles(i,1));
sa= sind(angles(i,1));
Length=sqrt((Coordinate(Node_to,1)-Coordinate(Node_from,1))^2+(Coordinate(Node_to,2)-Coordinate(Node_from,2))^2);
Lengthmap(i,1)=Length;

k=Area*Young/Length;

DOF1=2*(Node_from)-1;
DOF2=2*(Node_from);
DOF3=2*(Node_to)-1;
DOF4=2*(Node_to);

M(DOF1,DOF1)=M(DOF1,DOF1)+ca^2*k;
M(DOF1,DOF2)=M(DOF1,DOF2)+ca*sa*k;
M(DOF1,DOF3)=M(DOF1,DOF3)-ca^2*k;
M(DOF1,DOF4)=M(DOF1,DOF4)-ca*sa*k;

M(DOF2,DOF1)=M(DOF2,DOF1)+ca*sa*k;
M(DOF2,DOF2)=M(DOF2,DOF2)+sa^2*k;
M(DOF2,DOF3)=M(DOF2,DOF3)-ca*sa*k;
M(DOF2,DOF4)=M(DOF2,DOF4)-sa^2*k;

M(DOF3,DOF1)=M(DOF3,DOF1)-ca^2*k;
M(DOF3,DOF2)=M(DOF3,DOF2)-ca*sa*k;
M(DOF3,DOF3)=M(DOF3,DOF3)+ca^2*k;
M(DOF3,DOF4)=M(DOF3,DOF4)+ca*sa*k;

M(DOF4,DOF1)=M(DOF4,DOF1)-ca*sa*k;
M(DOF4,DOF2)=M(DOF4,DOF2)-sa^2*k;
M(DOF4,DOF3)=M(DOF4,DOF3)+ca*sa*k;
M(DOF4,DOF4)=M(DOF4,DOF4)+sa^2*k;

end

% READ DATA - READ REACTION DATA - Now read in the reactions

Number_reactions = fscanf(fid, '%d', 1);
Constraint = zeros(2*Number_nodes, 1);

if (2*Number_nodes ~= (Number_elements + Number_reactions))
error('Invalid number of nodes, elements, and reactions');
end

% CREATE CONSTRAINT MATRIX

for i=1:Number_reactions;
Reaction = fscanf(fid, '%d', 1);
Node=fscanf(fid, '%d', 1);
Direction = fscanf(fid, '%s', 1);
if ((Direction == 'y') || (Direction == 'Y'))
Constraint(2*Node) = 1;
elseif ((Direction == 'x')||(Direction == 'X'))
Constraint(2*Node-1)= 1;
else
error('Invalid direction for reaction')
end
end

External=zeros(2*Number_nodes,1);
Number_forces=fscanf(fid, '%d', 1);

% READ DATA - READ FORCE DATA - Now read in the forces

for i=1:Number_forces
Node =fscanf(fid, '%d', 1);
Force=fscanf(fid, '%g', 1);
Direction = fscanf(fid, '%g', 1);
External(2*Node-1)=External(2*Node-1)+Force*cos(Direction*(pi/180));
External(2*Node) =External(2*Node) +Force*sin(Direction*(pi/180));
end

New_row=1;

% REDUCE GLOBAL DOF MATRIX - Reduce the global DOF matrix to calculate Q

for i = 1:Number_nodes*2
if (Constraint(i)==0);
New_col=1;
for j=1:Number_nodes*2;
if (Constraint(j)==0);
MR(New_row, New_col)=M(i,j);
New_col=New_col+1;
end
end
New_row = New_row+1;
end
end

% REDUCE EXTERNAL FORCE MATRIX - Reduce external force to calculate Q

New_row1=1;

for i = 1:Number_nodes*2
if (Constraint(i)==0);
ER(New_row1,1)=External(i,1);
New_row1=New_row1+1;
end
end

% CALCULATE Q

Q=MR\ER;

% RECREATE Q MATRIX - Create Q matrix with constraints

Qmap=zeros(2*Number_nodes,1);

New_row2=1;

for i = 1:Number_nodes*2;
if (Constraint(i)==0);
Qmap(i,1)=Q(New_row2,1);
New_row2=New_row2+1;
end
end

% REDUCE GLOBAL DOF MATRIX - Reduce the global DOF matrix to calculate reactions

New_row3=1;
for i = 1:Number_nodes*2
if (Constraint(i)==1);
New_col2=1;
for j=1:Number_nodes*2;
MR2(New_row3, New_col2)=M(i,j);
New_col2=New_col2+1;
end
New_row3 = New_row3+1;
end
end

% CALCULATE REACTIONS

React=MR2*Qmap;

% RECREATE REACTION MATRIX - Create Reaction matrix with constraints

New_row4=1;
Reactmap=zeros(2*Number_elements,1);
for i = 1:Number_nodes*2;
if (Constraint(i)==1);
Reactmap(i,1)=React(New_row4,1);
New_row4=New_row4+1;
end
end

% CALCULATE STRESSES - Calculate the stresses in each element

cosa=cosd(angles);
sina=sind(angles);

for i=1:Number_elements
Ang=[-cosa(i) -sina(i) cosa(i) sina(i)];
Nodef=DOFa(i,1);
Nodet=DOFa(i,2);
Qsigma=[Qmap(2*Nodef-1,1);Qmap(2*Nodef,1);Qmap(2*Nodet-1,1);Qmap(2*Nodet,1)];
sigma(i,1)=Youngmap(i,1)/Lengthmap(i,1)*Ang*Qsigma;
end

% CALCULATE DISPLACEMENTS - Calculate the change in length

for i=1:Number_elements;
displace(i,1)=Length*sigma(i,1)/Young;
end

% REDO REACTION MATRIX - Break in X and Y for each element

react=zeros(2,Number_elements);
for i=1:Number_elements*2;
react(i)=Reactmap(i,1);
end
reacts=react';

% REDO DISPLACEMENT MATRIX - Break in X and Y for each element

Qrec=zeros(2,Number_elements);
for i=1:Number_elements*2;
Qrec(i)=Qmap(i,1);
end
Qlast=Qrec';

% PRINT ALL THE RESULTS

fprintf('\nNode X Y \n')
for i=1:Number_elements
fprintf('%d %3d %3d\n',Nodea(i,1), Coordinate(i,1), Coordinate(i,2))
end

fprintf('\nElement From To Area E \n')
for i=1:Number_elements
fprintf('%d %3d %3d %d %d\n',Elementa(i,1), DOFa(i,1), DOFa(i,2), Areamap(i,1), Youngmap(i,1))
end

fprintf('\nNodal Displacement \nNode DX DY \n')
for i=1:Number_elements
fprintf('%d %-5f %-5f\n',Nodea(i,1), Qlast(i,1), Qlast(i,2))
end

fprintf('\nElement From To Stress Length Change \n')
for i=1:Number_elements
fprintf('%d %3d %3d %5.0f %6.4f\n',Elementa(i,1), DOFa(i,1), DOFa(i,2), sigma(i,1), displace(i,1))
end

fprintf('\nReaction or Constraint Forces \nNode ForceX ForceY \n')
for i=1:Number_elements
fprintf('%d %5.f %5.f\n',Nodea(i,1), reacts(i,1), reacts(i,2))
end

end % End of program
Uh. Other than that, i got to see an electron microscope today. That was pretty neat. I kinda wish we got to keep the neat little pictures.
My Printer ran out of ink. That sucked too.
I also am hopefully going to be inducted into Pi Tau Sigma, the mechanical engineering society. Which is awesome. Went to the meeting, found out about some things. Go me.


Monday 2/23 Tuesday 2/24 Wednesday 2/25 Thursday 2/26 Friday 2/27
  Materials Exam  Materials Lab Due    
  ProE Due  ME 360 Due    
  Anniversary      
  Grandpa's Birthday      
Monday 3/2 Tuesday 3/3 Wednesday 3/4 Thursday 3/5 Friday 3/6
  ProE due ME314 Hwk Due Music Quiz  
  Scuba Tank Due      
  Vibes Due      
Laters
Nyssa

Currently Reading: Genius Squad
Books: 9== 9%
Words: 10000/10000=100%


This account has disabled anonymous posting.
If you don't have an account you can create one now.
HTML doesn't work in the subject.
More info about formatting

Profile

nikavia: (Default)
nikavia

September 2010

S M T W T F S
   1234
567891011
12131415161718
192021 22232425
2627282930  

Most Popular Tags

Style Credit

Expand Cut Tags

No cut tags
Page generated Jul. 14th, 2025 04:12 am
Powered by Dreamwidth Studios