You must modify this code so that it works for any truss that the user

wishes to analyze, whether 2-D or 3-D.

Therefore, after your code is modified, all trusses should be assumed as

3-D Space Trusses.

For 2-D planar trusses, the z-coordinate of each joint should be zero

and the e3 magnitude of all external forces and reactions should be zero.

.M file given that must be modify

% Units used:

LUnit = 'ft'; % Length = ft or m

FUnit = 'lb'; % Force = lb or kip or N or kN

%.. Input problem data

nNodes = 8; % Number of nodes (pin joints) in the truss

nMembers = 11; % Number of members

nReactions = 2; % Number of force vector reactions

nLoads = 3; % Number of externally applied loads

nDim = 3; % Dimensions 2=2D, 3=3D

% *** Once your code is modified and tested, nDim must always equal 3 and

% the variable nDim cannot be used in your code.

nDOF = nDim * nNodes; % Number of degrees of freedom

tol = 1.e-8; % Remainder tolerance on norm of residual force

% Input coordinates of each node

x(1,:) = [0,0,0]; % A

x(2,:) = [3,3,3]; % B

x(3,:) = [6,6,6]; % C

x(4,:) = [9,9,9]; % D

x(5,:) = [18, 0]; % G

x(6,:) = [9, 0,0]; % H

x(7,:) = [12, 6,3]; % E

x(8,:) = [15, 3,3]; % F

% Input member connections

%.. First input is the start node; Second input is the end node

Member(1,:) = [1,2]; % AB

Member(2,:) = [1,6]; % AH

Member(3,:) = [2,6]; % BH

Member(4,:) = [2,3]; % BC

Member(5,:) = [3,6]; % CH

Member(6,:) = [3,4]; % CD

Member(7,:) = [4,6]; % DH

Member(8,:) = [4,5]; % DG

Member(9,:) = [5,6]; % GH

Member(10,:) = [7,6]; % EH

Member(11,:) = [8,6]; % FH

% Input support reactions

%.. First input is the node number

%.. Second input is resistance in the e1 direction

%.. Third input is resistance in the e2 direction

%.. Fourth input is resistance in the e3 direction (zero for 2-D trusses)

%*** 1 for resistance in that direction, 0 for no resistance ***

Reaction(1,:) = [1, 1, 1]; % Pin at Node 1

Reaction(2,:) = [5, 0, 1]; % Roller at Node 5

% Input externallly applied loads

%.. First input is the node number, remaining inputs are the magnitude of

%.. the force in the e1, e2 and e3 directions, respectively.

Load(1,:) = [2, 500, -866.025]; % 1000 lb load at Node 2

Load(2,:) = [3, 866.025, -500]; % 1000 lb load at Node 3

Load(3,:) = [4, 1000, 0]; % 1000 lb load at Node 4

% End of Input Section for Assessment 04 16Sp.(2-D).......

% ************************************************************************

%% Compute the unit vector (direction) for every member

PosVec = zeros(nMembers, 2); %Initialize Position Vector matrix

UnVec = zeros(nMembers, 2); %Initialize Unit Vector matrix

OppUnVec = zeros(nMembers, 2); %Initialize opposite direction matrix

for i = 1 : nMembers

SN = Member(i,1);

EN = Member(i,2);

% ***** This section is where your CP3 code should go. *****

%% ---------------------------------------------------------------

%....

%.... PART 1

%.... PLANAR Truss

%....

%-----------------------------------------------------------------

%.... problem data

nNodesP = 12; %number of nodes

nMembersP = 21; %number of members

NodePoly = zeros(nNodesP, 3); %Initialize memory storage

MemberPoly = zeros(nMembersP, 2);

%.... node coordinates

NodePoly(1,:) = [0,0,0]; %Node 1

NodePoly(2,:) = [2,0,0]; %NOde 2

NodePoly(3,:) = [4,0,0]; %Node 3

NodePoly(4,:) = [6,0,0]; %Node 4

NodePoly(5,:) = [8,0,0]; %Node 5

NodePoly(6,:) = [10,0,0]; %Node 6

NodePoly(7,:) = [12,0,0]; %Node 7

NodePoly(8,:) = [10,3/2,0]; %Node 8

NodePoly(9,:) = [8,3,0]; %Node 9

NodePoly(10,:) = [6,3,0]; %Node 10

NodePoly(11,:) = [4,3,0]; %Node 11

NodePoly(12,:) = [2,3/2,0]; %Node 12

%.... member connections

MemberPoly(1,:) = [1,2];

MemberPoly(2,:) = [2,3];

MemberPoly(3,:) = [3,4];

MemberPoly(4,:) = [4,5];

MemberPoly(5,:) = [5,6];

MemberPoly(6,:) = [6,7];

MemberPoly(7,:) = [8,7];

MemberPoly(8,:) = [8,6];

MemberPoly(9,:) = [9,8];

MemberPoly(10,:) = [9,6];

MemberPoly(11,:) = [9,5];

MemberPoly(12,:) = [10,9];

MemberPoly(13,:) = [10,5];

MemberPoly(14,:) = [10,4];

MemberPoly(15,:) = [10,3];

MemberPoly(16,:) = [11,10];

MemberPoly(17,:) = [11,3];

MemberPoly(18,:) = [11,2];

MemberPoly(19,:) = [12,11];

MemberPoly(20,:) = [12,2];

MemberPoly(21,:) = [1,12];

%% .. Use a for loop to compute the position vector between connected

PosVecP = zeros(nMembersP, 3); %Initialize memory storage

UnitVecP = zeros(nMembersP, 3); %unit vector

OppPosVecP = zeros(nMembersP, 3);

for (j = 1 : nMembersP)

SN = MemberPoly(j,1); %SN = start node

EN = MemberPoly(j,2); %EN = end node

PosVec(j,:) =NodePoly(EN,:)-NodePoly(SN,:);

UnitVecP(j,:)=PosVecP(j,:)/norm(PosVecP(j,:));

PosVecP(j,:)=-PosVecP(j,:)/norm(PosVecP(j,:));

end % for j (members)

%% .Create output for command window

fprintf('%s ' , '----------------------------------------')

fprintf('%s ' , '-------------Planar Truss---------------')

fprintf('%s ' , '----------------------------------------')

for (i = 1 : nNodesP)

fprintf('%s%8i%8.3f%8.3f%8.3f ',' Node: ',i, NodePoly(i,:)')

end

fprintf(' ')

for (i = 1 : nMembersP)

fprintf('%s%8i%8.3f%8.3f%8.3f ',' Unit Vector: ',i, UnitVecP(i,:)')

end

fprintf(' ')

%% .Create the plot

fig1 = figure(1); clf; grid on; axis equal; hold on;

xlabel('X'); ylabel('Y'); title('Planar Truss');

for m = 1 : nMembersP

SN = MemberPoly(m,1);

EN = MemberPoly(m,2);

X = [NodePoly(SN,1); NodePoly(EN,1)];

Y = [NodePoly(SN,2); NodePoly(EN,2)];

Z = [NodePoly(SN,3); NodePoly(EN,3)];

p = plot(X,Y);

set(p, 'Color', 'k', 'LineWidth', 2);

end % for m, nMembersP

scatter(NodePoly(:,1),NodePoly(:,2),'fill','red');

%% ---------------------------------------------------------------

%....

%.... PART 1b

%.... SPACE Truss

%....

%-----------------------------------------------------------------

%....truss data

nNodesS = 7; %number of nodes

nMembersS = 15; %number of members

NodeS = zeros(nNodesS, 3); %Initialize memory storage

nMemberS = zeros(nMembersS, 2);

%.... node coordinates

NodeS(1,:) = [2,-3,0]; %Node 1

NodeS(2,:) = [2,3,0]; %Node 2

NodeS(3,:) = [-4,3,0]; %Node 3

NodeS(4,:) = [1,-3/2,6]; %Node 4

NodeS(5,:) = [1,3/2,6]; %Node 5

NodeS(6,:) = [-2,3/2,6]; %Node 6

NodeS(7,:) = [0,0,12]; %Node 7

%.... members coordinates

MemberS(1,:) = [1,2];

MemberS(2,:) = [2,3];

MemberS(3,:) = [2,4];

MemberS(4,:) = [2,5];

MemberS(5,:) = [1,3];

MemberS(6,:) = [1,4];

MemberS(7,:) = [3,6];

MemberS(8,:) = [3,5];

MemberS(9,:) = [3,4];

MemberS(10,:) = [5,4];

MemberS(11,:) = [5,6];

MemberS(12,:) = [6,4];

MemberS(13,:) = [4,7];

MemberS(14,:) = [5,7];

MemberS(15,:) = [6,7];

%% .. Use a for loop to compute the position vector between connected

PosVecPoly = zeros(nMembersS, 3); %Initialize memory storage

UnitVecPoly = zeros(nMembersS, 3); %unit vector

OppPosVecPoly = zeros(nMembersS, 3); %negative position vector

for (j=1:nMemberS)

SN=MemberS(j,1); %start node

EN=MemberS(j,2); %end node

PosVecPoly(j,:)=NodeS(EN,:)-NodeS(SN,:); %compute position vector

UnitVecPoly(j,:)=PosVecPoly(j,:)/norm(PosVecPoly(j,:)); %compute unit vector

OppPosVecPoly(j,:)=-PosVecPoly(j,:)/norm(PosVecPoly(j,:)); %compue negative position vector

end

%% ..Create output for command window

fprintf('%s ' , '---------------------------------------')

fprintf('%s ' , '-------------Space Truss---------------')

fprintf('%s ' , '---------------------------------------')

for (i = 1 : nNodesS)

fprintf('%s%8i%8.3f%8.3f%8.3f ',' Node: ',i, NodeS(i,:)')

end

fprintf(' ')

for (i = 1 : nMembersS)

fprintf('%s%8i%8.3f%8.3f%8.3f ',' Unit Vector: ',i, UnitVecPoly(i,:)')

end

fprintf(' ')

%% .Create the plot

fig2 = figure(2); clf; grid on; axis equal; hold on;

xlabel('X'); ylabel('Y'); zlabel('Z');title('Space Truss');

for m = 1 : nMembersS

SN = MemberS(m,1);

EN = MemberS(m,2);

X = [NodeS(SN,1); NodeS(EN,1)];

Y = [NodeS(SN,2); NodeS(EN,2)];

Z = [NodeS(SN,3); NodeS(EN,3)];

p = plot3(X,Y,Z);

set(p, 'Color', 'k', 'LineWidth', 2);

end % for m, nMemberS

scatter3(NodeS(:,1),NodeS(:,2),NodeS(:,3),'fill','blue');

%% ---------------------------------------------------------------

%....

%.... PART 2

%.... Polygon

%....

%-----------------------------------------------------------------

%.... Input problem data (radius, # of points, # of members, theta)

r = 10; %Radius

nNodesPoly = 6;

nMembersPoly = nNodesPoly - 1;

theta = pi/nMembersPoly; %a semi-circular arch

%% .. Use a for loop to compute the x, y and z coordinates of each point and

for (i=1:nNodesPoly)

x(i)=[r*cos((i-1)*theta)]; %x as a function of r and theta

y(i)=[r*sin((i-1)*theta)]; %y as a function of r and theta

z(i)=[0]; %z as a function of r and theta

end

NodePoly=[x;y;z]';

%.... Use a for loop to compute the members or connections and the unit

%.... vector for each member.

PosVecPoly = zeros(nMembersPoly, 3); %Initialize memory storage

UnitVecPoly = zeros(nMembersPoly, 3);

OppPosVecPoly = zeros(nMembersPoly, 3);

for j=1:nMembersPoly

MemberPoly(j,:) = [j,j+1];

SN=MemberPoly(j,1); %Start node

EN=MemberPoly(j,2); %End node

PosVecPoly(j,:)=NodePoly(EN,:)-NodePoly(SN,:);

UnitVecPoly(j,:)=PosVecPoly(j,:)/norm(PosVecPoly(j,:));

OppPosVecPoly(j,:)=-PosVecPoly(j,:)/norm(PosVecPoly(j,:));

end

%% .Create output for command window

fprintf('%s ' , '---------------------------------------')

fprintf('%s ' , '--------------Polygon------------------')

fprintf('%s ' , '---------------------------------------')

for (i = 1 : nNodesPoly)

fprintf('%s%8i%8.3f%8.3f%8.3f ',' Node: ', NodePoly(i,:)')

end

fprintf(' ')

for i = 1 : nMembersPoly

fprintf('%s%8i%8.3f%8.3f%8.3f ',' Unit Vector: ',i, UnitVecPoly(i,:)')

end

fprintf(' ')

%% ..Create the plot

fig3 = figure(3); clf; grid on; axis equal; hold on;

xlabel('X'); ylabel('Y'); title('Polygon Truss');

p = plot(NodePoly(:,1), NodePoly(:,2));

set(p,'Color','red','LineWidth',2);

scatter(NodePoly(:,1), NodePoly(:,2),'fill','green');

% End of Program

OppUnVec(i,:) = -PosVec(i,:) / norm(PosVec(i,:));

end

%% The coefficient matrix, C, contains the unit vectors (direction vectors)

% for each member force associated with each node. The unit vectors are

% then separated into e1, e2, e3 components and placed in separate rows.

% Therefore, each node has 2 rows (for 2-D) or 3 rows (for 3-D) associated

% with it; one for each of the e1, e2, e3 components.

% This corresponds with the Degrees of Freedom, nDOF, which is 2*nNodes

% for 2-D or 3*nNodes for 3-D.

% So the coefficient matrix has 2 or 3 rows for each node and one column

% for each member force.

C = zeros(2*nNodes, nMembers); %Initialize coefficient matrix

% Loop through all nodes, create the vector force equation for that node

% and store it in the proper row of the coefficient matrix, C.

for i = 1 : nNodes

for j = 1 : nMembers

SN = Member(j,1);

EN = Member(j,2);

if (SN == i) % If member j starts at node i

C(2*i-1, j) = UnVec(j,1);

C(2*i, j) = UnVec(j,2);

elseif (EN == i) % If member j ends at node i

C(2*i-1, j) = OppUnVec(j,1);

C(2*i, j) = OppUnVec(j,2);

end % if SN

end % for nMembers

end % for nNodes

%% The External Force matrix, F, contains the magnitude of all externally

% applied loads, (input as e1, e2, e3 components), stored in the

% proper rows to corespond with the node it is applied at.

% Therefore, the F matrix has 2 or 3 rows for each node and one column.

F = zeros(3*nNodes, 1); %Initialize External Force matrix

% Loop through all externally applied loads and place each component of the

% load (e1, e2, e3) in the proper rows in the external load matrix, F.

for i = 1 : nLoads

j = Load(i,1); % Which node, j is the load i on

F(3*j-1) = Load(i,2);

F(3*j) = Load(i,3);

end % for i,nLoads

%% The Restraint matrix, Crest, contains the unit vectors (directions)

% for each member force associated with those nodes that are restrained by

% external supports.

% The Crest matrix has one row for each reaction component and one column

% for each member.

%.... Loop through all nodes and determine if there is a reaction at that

%.... node and seperate the reaction node equations from the force load

%.... equations.

nReaction = 1; %Used to count the number of reaction equations

%.... The nReactionEquation vector determines if a row should be put in the

%.... Crest matrix or the Cfree matrix

nReactionEquation = zeros(1, 2*nNodes);

for (j = 1 : nReactions)

for (i = 1 : nNodes)

if (Reaction(j,1) == i)

if (Reaction(j,2) == 1)

Crest(nReaction,:) = C(2*i-1,:);

nReactionEquation(2*i-1) = 1;

% If nReactionEquation is 1 it should be in the Crest

% matrix and not the Cfree matrix.

nReaction = nReaction+1;

end

if (Reaction(j,3) == 1)

Crest(nReaction,:) = C(2*i,:);

nReactionEquation(2*i) = 1;

% If nReactionEquation is 1 it should be in the Crest

% matrix and not the Cfree matrix.

nReaction = nReaction+1;

end

end

end

end

nForce = 1; %Used to count the number of Cfree equations

for (i = 1 : 2*nNodes)

if (nReactionEquation(i) == 1)

else

Cfree(nForce,:) = C(i,:);

Ffree(nForce) = F(i);

nForce = nForce + 1;

end

end

%% Solve first for the member forces and then solve for the support

% reactions.

% All of the values in the Cfree and Ffree matrices are known so you can

% now solve for the T matrix, which is the tension force in each member.

T = % finish this equation

% Since all values in the Crest and T matrices are now known, you can solve

% for the support reactions.

Reactions = % finish this equation

% Calculate the residual to find the amount of error and ensure that

% equilibrium was satisfied.

Residual1 = (Cfree * T + Ffree');

Residual2 = (Crest * T + Reactions);

Res = norm(Residual1) + norm(Residual2);

% Create the support reaction matrix.

nReaction = 1;

for (i = 1 : nNodes)

if (nReactionEquation(2*i-1) == 1)

SupportReaction(i,1) = i;

SupportReaction(i,2) = Reactions(nReaction);

nReaction = nReaction+1;

end

if (nReactionEquation(2*i) == 1)

SupportReaction(i,1) = i;

SupportReaction(i,3) = Reactions(nReaction);

nReaction = nReaction+1;

end

end

%% Create output for command window

fprintf('%s ' , '----------------------------------------')

fprintf('%s ' , '------------- Truss --------------')

fprintf('%s ' , '----------------------------------------')

fprintf('%s',' Node Coordinates (', LUnit, ')')

fprintf (' ')

for i = 1 : nNodes

fprintf('%s%4i%8.3f%8.3f%8.3f ',' Node: ', i, x(i,:)')

end % for i, nNodes

fprintf (' ')

fprintf('%s',' Member Forces (', FUnit, ')')

fprintf (' ')

for i = 1 : nMembers

fprintf('%s%4i%12.3f ',' Member Force: ',i, T(i))

end % for i, nMembers

fprintf (' ')

for i = 1 : nReactions

for j = 1 : nNodes

if (Reaction(i,1) == j)

fprintf('%s%s%s%4i%12.3f%12.3f%12.3f ',...

' Support Reaction (',FUnit,') at Node ',...

SupportReaction(j,:))

fprintf (' ')

end % if Reaction

end % for j, nNodes

end % for i, nReactions

fprintf(' %s%8.3f ',' Residual Error: ', Res)

fprintf(' ')

%% Create the plot............................................

% Once your code is modified the plot will be 3-D and so you may need the

% "Rotate 3D" command to view the truss figure properly.

%.. Plot members with color indicating tension or compression

Marker = ceil(25 / sqrt(nNodes));

fig1 = figure(1); clf; grid on; axis equal;

hold on;

xlabel(cat(2, 'X (',LUnit,')' ));

ylabel(cat(2, 'Y (',LUnit,')' ));

title('Truss Analysis');

small = max(T)*tol;

for m = 1 : nMembers

if (T(m) < -small)

MColor = 'b'; % Color compression members blue

elseif (T(m) > small)

MColor = 'r'; % Color tension members red

else

MColor = 'k'; % Color "zero force members" black

end % if T

SN = Member(m,1);

EN = Member(m,2);

X = [x(SN,1); x(EN,1)];

Y = [x(SN,2); x(EN,2)];

p = plot(X,Y);

set(p, 'Color', MColor, 'LineWidth', 6/nDim);

end % for m, nMembers

%.. Plot external reaction forces

Rlength = 0.1 * max(max(x)); % establish size of line

e = [ 1, 0, 0 ; 0, 1, 0 ; 0, 0, 1];

for i = 1 : nReactions

for j = 1 : nNodes

if (Reaction(i,1) == j)

for k = 1 : 2

if(Reaction(i, k+1) == 1)

xR(1,:) = x(j,:);

xR(2,:) = x(j,:) - Rlength * e(k, 1:nDim);

p = plot(xR(:,1),xR(:,2));

set(p,'Color','g','LineWidth',8/nDim);

end % if Reaction == 1

end % for k

end % if Reaction == j

end % for j, nNodes

end % for i, nReactions

%.. Plot loads

Llength = 0.1 * max(max(x)); % establish size of line

Fmax = max(max(abs(Load))); % establish maximum force level

e = [ 1, 0, 0 ; 0, 1, 0 ; 0, 0, 1];

for i = 1 : nLoads

for j = 1 : nNodes

if (Load(i,1) == j)

for k = 1 : 2

F = Load(i, k+1) / Fmax;

xL(1,:) = x(j,:);

xL(2,:) = x(j,:) + Llength * F * e(k, 1:nDim);

p = plot(xL(:,1), xL(:,2));

set(p,'Color','c','LineWidth',8/nDim);

end % for k

end % if Load == j

end % for j, nNodes

end % for i, nLoads

%.. Plot nodes

plot(x(:,1), x(:,2), 'o', 'LineWidth', 2,...

'MarkerFaceColor','w',...

'MarkerEdgeColor','k',...

'MarkerSize',Marker);

view(2);

% End of Program