Parker/Gambrel Truss Algorithm- undefined function or variable 'm'

I kept getting a 'preallocating for speed error' on line 73 for r(m).

I solved that issue by setting:

r(m) = zeros; %(before the for loop function)

NOW I'm getting an 'undefined function or variable 'm' error' on line 65.

NEED to make this code plot a Parker/Gambrel truss like in the picture:

what is wrong with my code?

clear, clf, clc

choice1=menu('Truss Program would like you to enter truss dimensions manually and define its dimensions?'...
,'Enter Dimensions Manually');

switch choice1
    case 1
        J=input('Please enter the number of joints:\n');
        M=(2*J)-3; %Joint-member relationship for a simple, planar truss.
        C=zeros(J,M); %Preallocates the connection matrix
        x=zeros(J,1); %Preallocates the joint position vectors
        y=zeros(J,1); %Prompts user for coordinates for joint positions through a loop

    for j=1:J
            promptx=sprintf('Please enter the x coordinate of joint %d in FEET:\n',j);
            x(j)=input(promptx);
            prompty=sprintf('Now enter the y coordinate of joint %d in FEET:\n',j);
            y(j)=input(prompty);
    end
    
    %Builds the connection matrix from a loop of queries asking the
    %user which members are connected to each joint.
        
    cnct{1,J}=[];  %Each cell in 'cnct' contains the indices of the corresponding row in connection matrix 'C' that are to be assigned the value of 1.
    for j=1:J
            str=sprintf('How many members are connected to joint %d?\n',j);
            choice2=menu(str,'1','2','3','4','5','6');
            cnct{j}=zeros(choice2,1);
            for i=1:choice2
                rem=choice2-(i-1);
                str2=sprintf('Which members? (%d members remaining for joint %d.)\n',rem,j);
                cnct{j}(i)=input(str2);
            end
            C(j,cnct{j})=1;
    end
end
     
%--------------------DEFINE SUPPORT AND LOAD MATRICES---------------------%

tstart=tic;
%Creates support force connection matrix
Sx=zeros(J,3);
Sy=zeros(J,3);
Sx(1,1)=1;
Sy(1,2)=1;
Sy(J,3)=1;

%Creates a load vector describing the location and strength of the
%load on the truss
L=zeros(2*J,1);
choice3=menu('Would you like to set a load at a joint?','Set a Load');
switch choice3 
    case 1
        Lj=input('What number joint will the load be placed at?\n');
        F=input('Please enter the value of the load in Newtons: \n');
        L(J+Lj)=F;
end

%---------------------CONSTRUCT EQUILIBRIUM EQUATIONS---------------------%

%Preallocation

deltax=zeros(M,1);
deltay=zeros(M,1);
ux=zeros(M,1);
uy=zeros(M,1);
Ax=zeros(J,M);
Ay=zeros(J,M);
r=zeros(1,M);

%Calculates coefficients for matrix A.

for m=1:M
    vec=find(C(:,m));
    deltax(m)=abs(x(vec(2))-x(vec(1)));
    deltay(m)=abs(y(vec(2))-y(vec(1)));
    r(m)= magnitude(deltax(m),deltay(m)); 
    ux(m)=deltax(m)/r(m);
    uy(m)=deltay(m)/r(m);
    if x(vec(1))<x(vec(2))
        Ax(vec(1),m)=ux(m);
        Ax(vec(2),m)=-ux(m);
    else
        Ax(vec(1),m)=-ux(m);
        Ax(vec(2),m)=ux(m);
    end
    if y(vec(1))<y(vec(2))
        Ay(vec(1),m)=uy(m);
        Ay(vec(2),m)=-uy(m);
    else
        Ay(vec(1),m)=-uy(m);
        Ay(vec(2),m)=uy(m);
    end
end

A=[Ax Sx;Ay Sy]; %Creates matrix A
T=A\L; %Solves equilibrium equations
    
%--------------------------ANALYSIS OF RESULTS----------------------------%

f=plottruss(x,y,C,critmember,Lj);
 set(f,'Visible','on','Name','NameTest')

[Daniel M]

You haven't defined m at that point yet. Nor have you defined r, so there is no need to index into it anyways. You also did not choose the size of the variable using inputs to the function zeros.

Also you shouldn't name a variable rem because it is the name of a native function rem().

[Ruben Mora]

Yes, these errors are as a result to the fact that I am very new to MATLAB and trying to complete a project for school. I updated the size of the variable using inputs:

r=zeros(1,M);

is that correct? b/c now I am getting:

'undefined function or variable 'magnitude' error'

how would I define m and r?

and not sure where you see that I used 'rem()'

Thanks for all your help.

Would you be able to help me code this correctly- I just need to make a graphic showing a plot of the truss with the point load forces at the nodes as seen in the inserted image

P.S. - Thanks for your quick response!

[Daniel M]

I can't run your code because I don't know what inputs to enter. Can you post something with the inputs already written in?

[Ruben Mora]

now I am getting this error with this code:

Unable to perform

assignment because the

indices on the left side

are not compatible with the

size of the right side.

Error in inputs (line 88)

r(m)=zeros(1,M);

here is the code w/inputs:

clear, clf, clc
choice1=menu('Truss Program would like you to select a truss type and define its dimensions?'...
    ,'Select Truss Type');
switch choice1
    case 1
        %Select truss geometry (i.e. warren truss, pratt truss, parker truss, etc.)
        type=input('Please type the name of the truss in single quotes:\n');
        J=input('Please enter the number of joints:\n');
        M=(2*J)-3;
        C=zeros(J,M); %Preallocates the connection matrix
        %Parker Truss
        if type=='P'
            a=input('Please enter the base dimension (cm):\n');
            b=a*sin(pi/3); %height of the truss
            bottom=(J/2)+.5; %number of joints on top and bottom
            top=(J/2)-.5;
            v=[1 1 0 1 1];
            m=2;
         for j=3:(J-2)
             C(j,m:(m+4))=v;
             m=m+2;
         end
         C(1,1:2)=[1 1];
         C(2,1:4)=[1 0 1 1];
         C((J-1),(M-3):M)=[1 1 0 1];
         C(J,(M-1):M)=[1 1];
        end
        
    %Creates x and y location vectors for the joints of the Parker
    %Truss
    x=zeros(1,J);
    y=zeros(1,J);
    
    n=0;
    for j=1:2:J
        x(j)=a*n;
        n=n+1;
    end
    
    n=0;
    for j=2:2:J-1
        x(j)=(a/2)+a*n;
        y(j)=b;
        n=n+1;
    end
end

%--------------------DEFINE SUPPORT AND LOAD MATRICES---------------------%

tstart=tic;
%Creates support force connection matrix
Sx=zeros(J,3);
Sy=zeros(J,3);
Sx(1,1)=1;
Sy(1,2)=1;
Sy(J,3)=1;

%Creates a load vector describing the location and strength of the
%load on the truss
L=zeros(2*J,1);
choice3=menu('Would you like to use the default truss loading or set a load at a joint?','Default Loading','Set a Load');
switch choice3
    case 1
        F=(500*9.81)/1000;
        Lj=input('What number joint will the load be placed at?\n');
        L(J+Lj)=F;
    case 2
        Lj=input('What number joint will the load be placed at?\n');
        F=input('Please enter the value of the load in Newtons: \n');
        L(J+Lj)=F;
end

%---------------------CONSTRUCT EQUILIBRIUM EQUATIONS---------------------%

%Preallocation
deltax=zeros(M,1);
deltay=zeros(M,1);
ux=zeros(M,1);
uy=zeros(M,1);
Ax=zeros(J,M);
Ay=zeros(J,M);
r(m)=zeros(1,M);

%Calculates coefficients for matrix A.
for m=1:M
    vec=find(C(:,m));
    deltax(m)=abs(x(vec(2))-x(vec(1)));
    deltay(m)=abs(y(vec(2))-y(vec(1)));
    r(m)=magnitude(deltax(m),deltay(m)); 
    ux(m)=deltax(m)/r(m);
    uy(m)=deltay(m)/r(m);
    if x(vec(1))<x(vec(2))
        Ax(vec(1),m)=ux(m);
        Ax(vec(2),m)=-ux(m);
    else
        Ax(vec(1),m)=-ux(m);
        Ax(vec(2),m)=ux(m);
    end
    if y(vec(1))<y(vec(2))
        Ay(vec(1),m)=uy(m);
        Ay(vec(2),m)=-uy(m);
    else
        Ay(vec(1),m)=-uy(m);
        Ay(vec(2),m)=uy(m);
    end
end

A=[Ax Sx;Ay Sy]; %Creates matrix A
T=A\L; %Solves equilibrium equations

%--------------------------ANALYSIS OF RESULTS----------------------------%

%Shows which members are in tension/compression

Tm=T(1:length(T)-3);
tension=find(Tm>0);
compression=find(Tm<0);
zeroforce=find(Tm==0);

%Graph the truss
 
 f=plottruss(x,y,C,critmember,Lj);
 set(f,'Visible','on','Name','NameTest')

[Daniel M]

Please upload an .m file that contains a test case already programmed in. I do not know what inputs you want, and it makes it more annoying to debug. Make it easy for people to help you.

[Ruben Mora]

% Load
F = 10; % ton

c = cosd(45);
s = sind(45);

% Coefficient matrix
A = [-c  0  0  0  c;
     -s  0  -1 0  -s;
      0  -1 0  1  0;
      0  0  1  0  0;
      0  0  0  -1 -c];
 
 % Load vector
 b = [0;
      0;
      0;
      F;
      0 ];
  
  % Solve linear system
  f = A\b;
  
  % GaussPivot: Gauss elimination pivoting
  %   x = GaussPivot(A,b): Gauss elimination with pivoting.
  % input: 
  %   A = coefficient matrix
  %   b = right hand side vector
  % output:
  %   x = solution vector 
  [m,n] = size(A);
  if m == n
  end
  nb = n + 1;
  Aug = [A b]; 
  % forward elimination
  for k = 1:minus(n,1)
      % partial pivoting
      [big,i] = max(abs(Aug(k:n,k)));
      ipr = i + k;
      if ipr  == k
          Aug([k,ipr],:) = Aug([ipr,k],:);
      end
       for i = k+1:n
           factor  =  Aug(i,k)/Aug(k,k);
           Aug(i,k:nb) = minus(Aug(i,k:nb) , factor*Aug(k,k:nb));
       end
       theta = angle(30);
      
  end
  % back substitution
  x = zeros(n,1);
  x(n) = Aug(n,nb)/Aug(n,n);
  for i = minus(n , 1):uminus(1):1
      x(i) = ((Aug(i,nb) - Aug(i,i+1:n)*x(i+1:n)))/Aug(i,i);
  end  
  % Post processing
  plotTruss(f);
  

  function plotTruss(f)
  % Scale range
  xmin = min(f);
  xmax = max(f);
  % Heatmap
  L = 10; %ft
  setupHeatmapPlot(xmin,xmax, length(f));
  plotMember(0, 0,  L, L, f(1),xmin, xmax );
  plotMember(0, 0,  L, 0, f(2),xmin, xmax );
  plotMember(L, 0,  L, L, f(3),xmin, xmax );
  plotMember(L, 0, 2*L,0, f(4),xmin, xmax );
  plotMember(L, L, 2*L,0, f(5),xmin, xmax );
  end
  % Plot structure member i-j colored by load Fij
  function plotMember(xi,yi,xj,yj,Fij,Fmin,Fmax)
  plot([xi xj], [yi yj],'-','LineWidth',5,'Color',scaleColor(Fij,Fmin,Fmax));
  plot([xi xj], [yi yj],'-','LineWidth',5,'Color',[0 0 0]);
  end
  
 % Return a scaled color of value v b/w range [vmin vmax]
 function c = scaleColor(v,vmin,vmax)
 m = colormap();
 i = linspace(vmin,vmax,length(m));
 c = interp1(i,m,v,'linear');
 end
 
 %Prepare figure to display heatmap
 function b = setupHeatmapPlot(vmin,vmax,n)
 colormap(hsv(n));
 b = colorbar();
 b.TickLabels = vmin + b.Ticks*(vmax - vmin);
 axis equal
 hold on
 grid on
 end

[Daniel M]

This seems to be completely different code. It also executes without error. Is your issue solved?

[Ruben Mora]

I need it to look like the one on the image. Can't quite figure it out..

Thank You very much-

I do appreciate your patience!!

[Daniel M]

If you have a specific question, I can try to help.

[Ruben Mora]

I made another code and it's giving an error. Could you run it and check it, please?

NodeCoordinates={{0,0,0},{10,5,0},{10,0,0},{20,8,0},{20,0,0},{30,9,0},...
    {30,0,0},{40,8,0},{40,0,0},{50,5,0},{50,0,0},{60,0,0}};

ElemNodes={{1,3},{3,5},{5,7},{7,9},{9,11},{11,12},...
    {1,2},{2,4},{4,6},{6,8},{8,10},{10,12},...
    {2,3},{4,5},{6,7},{8,9},{10,11},...
    {2,5},{4,7},{7,8},{9,10}};

PrintSpaceTrussNodeCoordinates(NodeCoordinates,'Node coordinates:',[1 0.000000 0.000000 0.000000; 2 10.000000 5.000000 0.000000; 3 10.000000 0.000000 0.000000; 4 20.000000 8.000000 0.000000; 5 20.000000 0.000000 0.000000; 6 30.000000 9.000000 0.000000; 7 30.000000 0.000000 0.000000; 8 40.000000 8.000000 0.000000; 9 40.000000 0.000000 0.000000; 10 50.000000 5.000000 0.000000; 11 50.000000 0.000000 0.000000; 12 60.000000 0.000000 0.000000
]);
numnod=Length(NodeCoordinates);  numele=Length(ElemNodes);
Em=1000; Abot=2; Atop=10; Abat=3; Adia=1;
ElemMaterials= Table(Em,{numele});
ElemFabrications={Abot,Abot,Abot,Abot,Abot,Abot,Atop,Atop,Atop,Atop,...
    Atop,Atop,Abat,Abat,Abat,Abat,Abat,Adia,Adia,Adia,Adia};
PrintSpaceTrussElementData(ElemNodes,ElemMaterials,ElemFabrications,...
    'Element data:',[1 {1, 3} 1000.00 2.00 2 {3, 5} 1000.00 2.00 3 {5, 7} 1000.00 2.00 4 {7, 9} 1000.00 2.00 5 {9, 11} 1000.00 2.00 6 {11, 12} 1000.00 2.00 7 {1, 2} 1000.00 10.00 8 {2, 4} 1000.00 10.00 9 {4, 6} 1000.00 10.00 10 {6, 8} 1000.00 10.00 11 {8, 10} 1000.00 10.00 12 {10, 12} 1000.00 10.00 13 {2, 3} 1000.00 3.00 14 {4, 5} 1000.00 3.00 15 {6, 7} 1000.00 3.00 16 {8, 9} 1000.00 3.00 17 {10, 11} 1000.00 3.00 18 {2, 5} 1000.00 1.00 19 {4, 7} 1000.00 1.00 20 {7, 8} 1000.00 1.00 21 {9, 10} 1000.00 1.00
]);
ProcessOptions= {True};

view ={ {0,1,0},{0,0,1} };
labels={{True,0.06,-1.5,1.5},{True,0.12},{"Times",11,"Roman"}};
PlotSpaceTrussElementsAndNodes(NodeCoordinates,ElemNodes,...
    'bridge mesh',{view,-1,labels}) ;

NodeDOFTags=  Table({0,0,1},{numnod});
NodeDOFValues=Table({0,0,0},{numnod});
NodeDOFValues(3)=({0,-10,0}); 
NodeDOFValues(5)= {0,-10,0};
NodeDOFValues(7)={0,-16,0};
NodeDOFValues(9)={0,-10,0};
NodeDOFValues(11)={0,-10,0};
NodeDOFTags(1)= {1,1,1};     %(* fixed node 1 *)
NodeDOFTags(numnod)={0,1,1}; %(* hroller @ node 12 *)
PrintSpaceTrussFreedomActivity(NodeDOFTags,NodeDOFValues,...
    'DOF Activity:',[1 1 1 1 0 0 0 2 0 0 1 0 0 0 3 0 0 1 0 -10 0 4 0 0 1 0 0 0 5 0 0 1 0 -10 0 6 0 0 1 0 0 0 7 0 0 1 0 -16 0 8 0 0 1 0 0 0 9 0 0 1 0 -10 0 10 0 0 1 0 0 0 11 0 0 1 0 -10 0 12 0 1 1 0 0 0
]);

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

NodeDisplacements,NodeForces,ElemForces,ElemStresses = SpaceTrussSolution( NodeCoordinates,ElemNodes,ElemMaterials,ElemFabrications,NodeDOFTags,NodeDOFValues,ProcessOptions );

PrintSpaceTrussNodeDisplacements(NodeDisplacements,...
    'Computed node displacements:',[1 0.000000 0.000000 0.000000 2 0.809536 -1.775600 0.000000 3 0.280000 -1.792260 0.000000 4 0.899001 -2.291930 0.000000 5 0.560000 -2.316600 0.000000 6 0.847500 -2.385940 0.000000 7 0.847500 -2.421940 0.000000 8 0.795999 -2.291930 0.000000 9 1.135000 -2.316600 0.000000 10 0.885464 -1.775600 0.000000 11 1.415000 -1.792260 0.000000 12 1.695000 0.000000 0.000000
]);
PrintSpaceTrussNodeForces(NodeForces,...
    'Node forces including reactions:',[1 0.0000 28.0000 0.0000 2 0.0000 0.0000 0.0000 3 0.0000 -10.0000 0.0000 4 0.0000 0.0000 0.0000 5 0.0000 -10.0000 0.0000 6 0.0000 0.0000 0.0000 7 0.0000 -16.0000 0.0000 8 0.0000 0.0000 0.0000 9 0.0000 -10.0000 0.0000 10 0.0000 0.0000 0.0000 11 0.0000 -10.0000 0.0000 12 0.0000 28.0000 0.0000]);
PrintSpaceTrussElemForcesAndStresses(ElemForces,ElemStresses,...
    'Int Forces and Stresses:',[1 56.0000 28.0000 2 56.0000 28.0000 3 57.5000 28.7500 4 57.5000 28.7500 5 56.0000 28.0000 6 56.0000 28.0000 7 -62.6100 -6.2610 8 -60.0300 -6.0030 9 -60.3000 -6.0300 10 -60.3000 -6.0300 11 -60.0300 -6.0030 12 -62.6100 -6.2610 13 10.0000 3.3330 14 9.2500 3.0830 15 12.0000 4.0000 16 9.2500 3.0830 17 10.0000 3.3330 18 1.6770 1.6770 19 3.2020 3.2020 20 3.2020 3.2020 21 1.6770 1.6770
]);

view={{0,1,0},{0,0,1}};
box={{0,-4,0},{60,-4,0},{60,10,0},{0,10,0}};
PlotSpaceTrussDeformedShape(NodeCoordinates,ElemNodes,NodeDisplacements,...
    {0,1},box,'deformed shape (unit magnif)',{view,-1,{'black','blue'}});
PlotSpaceTrussStresses(NodeCoordinates,ElemNodes,ElemStresses,1,box,...
    "axial stresses in truss members",{view,0,labels});

[Daniel M]

You didn't post the functions that you're calling. So I get an error on line 7 for an undefined function PrintSpaceTrussNodeCoordinates.