Won't Plot in 3D

Question

Ariana Garcia 2019-4-20

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/457598-won-t-plot-in-3d

评论： Walter Roberson 2019-4-20

%-----------------------------------------------------------------
%.... Name:Ariana Garcia
%.... Date:April 21,2017
%.... Project:CP4
%.... Recitation Section:Mon
%-----------------------------------------------------------------
%-----------------------------------------------------------------
%.... 
%.... Truss Analysis
%.... 
%-----------------------------------------------------------------
%***** READ THIS *****
% 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.
clear; clc;
      
% Print Banner to command window
	fprintf('%s\n','*-----------------------------------------------*')
	fprintf('%s\n','| CEE210 Statics, Computing Project 4           |')
	fprintf('%s\n','| Original by Amie Baisley, 11.01.2013          |')
	fprintf('%s\n','| Modified by Chris A. Lawrence, 04.02.2016     |')
	fprintf('%s\n','*-----------------------------------------------*')
	fprintf('\n\n\n')
%% ***********************************************************************
% Create a separate m file that contains several different truss input
% sections, any one of which can be copied and pasted into the truss
% analysis program, without having to modify the rest of the code.
% The file should be named CEE210_NNN_Truss_Inputs.m where NNN is your
% initials.
% The input file should contain every truss that was used to test your
% truss analysis code and must be submitted with your project.
% ************************************************************************
%% **********************************************************************
% 3D Space Trust from CP3
% Units used: 
	LUnit = 'ft';	% Length = ft or m
	FUnit = 'lb';	% Force = lb or kip or N or kN
%.. Input problem data
    nNodes = 7;             % Number of nodes (pin joints) in the truss
    nMembers = 15;          % Number of members
    nReactions = 3;         % Number of force vector reactions
    nLoads = 1;             % 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 node coordinates
    Node(1,:) = [2,-3,0];
    Node(2,:) = [2,3,0];
    Node(3,:) = [-4,3,0];
    Node(4,:) = [1,-1.5,6];
    Node(5,:) = [1,1.5,6];
    Node(6,:) = [-2,1.5,6];
    Node(7,:) =[0,0,12];
    %.... Finish inputting nodes
    
%.... Input member connections
    Member(1,:) = [1,2];
    Member(2,:) = [2,3];
    Member(3,:) = [3,4];
    Member(4,:) = [3,5];
    Member(5,:) = [3,6];
    Member(6,:) = [1,3];
    Member(7,:) = [2,5];
    Member(8,:) = [2,4];
    Member(9,:) = [1,4];
    Member(10,:) = [4,7];
    Member(11,:) = [5,7];
    Member(12,:) = [6,7];
    Member(13,:) = [5,6];
    Member(14,:) = [4,5];
    Member(15,:) = [5,7];
    %.... Finish inputting member connections
    
    % 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,0];		% Pin at Node 1
    Reaction(2,:) = [2, 0,1,0];		% Roller at Node 5
    Reaction(3,:) = [3, 1,1,1];	
    
% 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,:) = [5, 0, 200, -500];		% 1000 lb load at Node 2
    
  
	
% End of Input Section for CP3 Space Trust.......
%*****************************************************************************************
%%  Compute the unit vector (direction) for every member
PosVec = zeros(nMembers, 3);		%Initialize Position Vector matrix
UnVec = zeros(nMembers, 3);		%Initialize Unit Vector matrix
OppUnVec = zeros(nMembers, 3);	%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.  *****
    PosVec(i,:) = Node(EN,:)- Node(SN,:);   % postion vector
        length=norm(PosVec(i,:));
        UnVec(i,:)=PosVec(i,:)./norm(PosVec(i,:));    
	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(3*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(3*i-1, j) = UnVec(j,1);
			C(3*i, j) = UnVec(j,2);
            C(3*i, j) = UnVec(j,3);
		elseif (EN == i)  % If member j ends at node i
			C(3*i-1, j) = OppUnVec(j,1);
			C(3*i, j) = OppUnVec(j,2);
            C(3*i, j) = OppUnVec(j,3);
		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);
    F(3*j) = Load(i,4);
    
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, 3*nNodes);
for (j = 1 : nReactions)
	for (i = 1 : nNodes)
		if (Reaction(j,1) == i)
			if (Reaction(j,2) == 1)
				Crest(nReaction,:) = C(3*i-1,:);
				nReactionEquation(3*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(3*i,:);
				nReactionEquation(3*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 : 3*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 = -(Cfree)\(Ffree)' ;     % finish this equation
% Since all values in the Crest and T matrices are now known, you can solve
% for the support reactions.
Reactions = -Crest*T ;    % 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(3*i-1) == 1)
		SupportReaction(i,1) = i;
		SupportReaction(i,3) = Reactions(nReaction);
		nReaction = nReaction+1;   
    end
    if (nReactionEquation(3*i-2) == 1)
		SupportReaction(i,1) = i;
		SupportReaction(i,2) = Reactions(nReaction);
		nReaction = nReaction+1;   
	end
	if (nReactionEquation(3*i) == 1)
		SupportReaction(i,1) = i;
		SupportReaction(i,4) = Reactions(nReaction);
		nReaction = nReaction+1;   
	end
end
%% Create output for command window 
fprintf('%s\n'     ,  '----------------------------------------')
fprintf('%s\n'     ,  '-------------    Truss    --------------')
fprintf('%s\n'     ,  '----------------------------------------')
fprintf('%s',' Node Coordinates (', LUnit, ')')
fprintf ('\n')
for i = 1 : nNodes
	fprintf('%s%4i%8.3f%8.3f%8.3f\n','   Node: ', i, Node(i,:)')
end % for i, nNodes
	fprintf ('\n')
	fprintf('%s',' Member Forces (', FUnit, ')')
	fprintf ('\n')
for i = 1 : nMembers
	fprintf('%s%4i%12.3f\n','   Member Force: ',i, T(i)) 
end % for i, nMembers
	fprintf ('\n')
for i = 1 : nReactions
	for j = 1 : nNodes
		if (Reaction(i,1) == j)
			fprintf('%s%s%s%4i%12.3f%12.3f%12.3f\n',...
				'   Support Reaction (',FUnit,') at Node ',...
				SupportReaction(j,:)) 
			fprintf ('\n')
		end % if Reaction
	end % for j, nNodes
end % for i, nReactions
fprintf('\n%s%8.3f\n','   Residual Error: ', Res)
fprintf('\n\n\n')
      
      
%% 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,')' )); 
zlabel(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 = [Node(SN,1); Node(EN,1)];
	Y = [Node(SN,2); Node(EN,2)];
    Z = [Node(SN,3); Node(EN,3)];
  if nDim== 3  
	p = plot3(X,Y,Z);
  end
  if nDim== 2  
	p = plot(X,Y);
  end
 
	set(p, 'Color', MColor, 'LineWidth', 6/nDim);
end % for m, nMembers
		
%.. Plot external reaction forces
Rlength = 0.1 * max(max(Node));  % 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,:) = Node(j,:);
					xR(2,:) = Node(j,:) - Rlength * e(k, 1:nDim);
                    xR(3,:) = Node(j,:) - Rlength * e(k, 1:nDim);
                    if nDim== 2
					p = plot(xR(:,1),xR(:,2));
                    end
                    if nDim== 3
                    	p = plot3(xR(:,1),xR(:,2),xR(3,:));
                    end
                 
					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(Node));  % 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,:) = Node(j,:);
				xL(2,:) = Node(j,:) + Llength * F * e(k, 1:nDim); 
                xL(3,:) = Node(j,:) + Llength * F * e(k, 1:nDim); 
				if nDim== 2
                p = plot(xL(:,1), xL(:,2));
                end
                if nDim== 3
                p = plot3(xL(:,1), xL(:,2),xL(3,:));
                end
                
			set(p,'Color','c','LineWidth',8/nDim);
			end % for k
		end % if Load == j
	end % for j, nNodes
end % for i, nLoads
      
%.. Plot nodes 
if nDim==2
plot(Node(:,1), Node(:,2), 'o', 'LineWidth', 2,...
	'MarkerFaceColor','w',...
	'MarkerEdgeColor','k',...
	'MarkerSize',Marker);
view(2);
end
if nDim == 3
plot3(Node(:,1), Node(:,2),Node(:,3), 'o', 'LineWidth', 2,...
	'MarkerFaceColor','w',...
	'MarkerEdgeColor','k',...
	'MarkerSize',Marker);
view(2);
end
		
% End of Program