Mercurial > forge
changeset 8460:464dc6d279e4 octave-forge
Initial join of ocframe with mechanics
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/main/mechanics/inst/ocframe/MSNForces.m Fri Sep 30 15:19:38 2011 +0000 @@ -0,0 +1,205 @@ +## Copyright (C) 2010 Johan Beke +## +## This software is free software; you can redistribute it and/or modify it +## under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 3 of the License, or (at +## your option) any later version. +## +## This software is distributed in the hope that it will be useful, but +## WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +## General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with this software; see the file COPYING. If not, see +## <http://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## @deftypefn {Function File} [@var{x}, @var{M}, @var{S}, @var{N}] = MSNForces(@var{joints}, @var{members}, @var{dist}, @var{point}, @var{MemF}, @var{membernum}, @var{divisions}) +## +## This function returns the internal forces of a member for each position x. The member +## is divided in 20 subelements if the argument is not given. The used sign convention is displayed in the help file. +## +## Input parameters are similar as with SolveFrame and PlotFrame with extra arguments: +## membernum = Number of the member to calculate +## divisions = Number of divisions for the member +## +## @end deftypefn + +function [x,M,S,N]=MSNForces(joints,members,dist,point,MemF,membernum,divisions) + if (nargin < 6) + usage("[x,M,S,N]=MSNForces(joints,members,dist,point,MemF,membernum,divisions) use the help command for more information"); + end + if (nargin<7) + divisions=20; + end + Near=members(membernum,1); + Far=members(membernum,2); + L=sqrt((joints(Near,1)-joints(Far,1))^2+(joints(Near,2)-joints(Far,2))^2); + x=0:L/divisions:L; + M=zeros(1,columns(x)); + S=zeros(1,columns(x)); + N=zeros(1,columns(x)); + FxN_end=MemF(membernum,1); + FyN_end=MemF(membernum,2); + MzN_end=MemF(membernum,3); + + %Moment forces + for i=1:columns(x) + %due to end forces + M(i)+=-MzN_end+FyN_end*x(i); + %due to dist + for j=1:rows(dist) + if (dist(j,1)==membernum) + FxN=dist(j,2); + FyN=dist(j,3); + FxF=dist(j,4); + FyF=dist(j,5); + a=dist(j,6); + b=dist(j,7); + local=dist(j,8); + if (local==0) + %convert FxN,... to local axes + temp=[FxN,FyN,0,FxF,FyF,0]; + temp=(TransformationMatrix(joints(Near,1),joints(Near,2),joints(Far,1),joints(Far,2))')^-1*temp'; + FxN=temp(1); + FyN=temp(2); + FxF=temp(4); + FyF=temp(5); + end + if (x(i)>L-b) + M(i)+= (FyN+FyF)/2*(L-a-b)* (x(i)-((FyN+2*FyF)*L+(2*a-b)*FyN+(a-2*b)*FyF)/(3*FyN+3*FyF)); + end + if (x(i)>a && x(i)<=L-b) + %% bug :S zie bug.m + %M(i)+=(((3*FyN*x(i)+3*a*FyN)*L+(2*FyF-2*FyN)*x(i)^2+((-3*b-2*a)*FyN-a*FyF)*x(i)+(-3*a*b-2*a^2)*FyN-a^2*FyF)*((2*FyN*x(i)*L+(FyF-FyN)*x(i)^2+(-2*b*FyN-2*a*FyF)*x(i))/(2*L-2*b-2*a)-(2*a*FyN*L+(-2*a*b-a^2)*FyN-a^2*FyF)/(2*L-2*b-2*a)))/(6*FyN*L+(3*FyF-3*FyN)*x(i)+(-6*b-3*a)*FyN-3*a*FyF); + M(i)+=((3*FyN*x(i)^2-6*a*FyN*x(i)+3*a^2*FyN)*L+(FyF-FyN)*x(i)^3+(-3*b*FyN-3*a*FyF)*x(i)^2+((6*a*b+3*a^2)*FyN+3*a^2*FyF)*x(i)+(-3*a^2*b-2*a^3)*FyN-a^3*FyF)/(6*L-6*b-6*a); + end + end + end + %due to point load on member + for j=1:rows(point) + if (point(j,1)==membernum) + Fx=point(j,2); + Fy=point(j,3); + a=point(j,4); + local=point(j,5); + if (local==0) + %convert to global axes + temp=[Fx,Fy,0,0,0,0]; + temp=(TransformationMatrix(joints(Near,1),joints(Near,2),joints(Far,1),joints(Far,2))')^-1*temp'; + Fx=temp(1); + Fy=temp(2); + end + if (x(i)>a) + M(i)+=(x(i)-a)*Fy; + end + end + end + + end + M=M.*-1; %sign convention changed during development + + %Shear forces + for i=1:columns(x) + %due to end forces + S(i)+=FyN_end; + %due to dist + for j=1:rows(dist) + if (dist(j,1)==membernum) + FxN=dist(j,2); + FyN=dist(j,3); + FxF=dist(j,4); + FyF=dist(j,5); + a=dist(j,6); + b=dist(j,7); + local=dist(j,8); + if (local==0) + %convert FxN,... to local axes + temp=[FxN,FyN,0,FxF,FyF,0]; + temp=(TransformationMatrix(joints(Near,1),joints(Near,2),joints(Far,1),joints(Far,2))')^-1*temp'; + FxN=temp(1); + FyN=temp(2); + FxF=temp(4); + FyF=temp(5); + end + if (x(i)>L-b) + S(i)+=((FyN+FyF)*L^2+((-2*b-2*a)*FyF-2*b*FyN)*L+b^2*FyN+(b^2+2*a*b)*FyF)/(2*L-2*b-2*a)-(2*a*FyN*L+(-2*a*b-a^2)*FyN-a^2*FyF)/(2*L-2*b-2*a); + end + if (x(i)>a && x(i)<=L-b) + S(i)+=(2*FyN*x(i)*L+(FyF-FyN)*x(i)^2+(-2*b*FyN-2*a*FyF)*x(i))/(2*L-2*b-2*a)-(2*a*FyN*L+(-2*a*b-a^2)*FyN-a^2*FyF)/(2*L-2*b-2*a); + end + end + end + %due to point load on member + for j=1:rows(point) + if (point(j,1)==membernum) + Fx=point(j,2); + Fy=point(j,3); + a=point(j,4); + local=point(j,5); + if (local==0) + %convert to global axes + temp=[Fx,Fy,0,0,0,0]; + temp=(TransformationMatrix(joints(Near,1),joints(Near,2),joints(Far,1),joints(Far,2))')^-1*temp'; + Fx=temp(1); + Fy=temp(2); + end + if (x(i)>a) + S(i)+=Fy; + end + end + end + end + %Normal forces + for i=1:columns(x) + %due to end forces + N(i)+=-FxN_end; + %due to dist + for j=1:rows(dist) + if (dist(j,1)==membernum) + FxN=dist(j,2); + FyN=dist(j,3); + FxF=dist(j,4); + FyF=dist(j,5); + a=dist(j,6); + b=dist(j,7); + local=dist(j,8); + if (local==0) + %convert FxN,... to local axes + temp=[FxN,FyN,0,FxF,FyF,0]; + temp=(TransformationMatrix(joints(Near,1),joints(Near,2),joints(Far,1),joints(Far,2))')^-1*temp'; + FxN=temp(1); + FyN=temp(2); + FxF=temp(4); + FyF=temp(5); + end + if (x(i)>L-b) + N(i)-=((FxN+FxF)*L^2+((-2*b-2*a)*FxF-2*b*FxN)*L+b^2*FxN+(b^2+2*a*b)*FxF)/(2*L-2*b-2*a)-(2*a*FxN*L+(-2*a*b-a^2)*FxN-a^2*FxF)/(2*L-2*b-2*a); + end + if (x(i)>a && x(i)<=L-b) + N(i)-=(2*FxN*x(i)*L+(FxF-FxN)*x(i)^2+(-2*b*FxN-2*a*FxF)*x(i))/(2*L-2*b-2*a)-(2*a*FxN*L+(-2*a*b-a^2)*FxN-a^2*FxF)/(2*L-2*b-2*a); + end + end + end + %due to point load on member + for j=1:rows(point) + if (point(j,1)==membernum) + Fx=point(j,2); + Fy=point(j,3); + a=point(j,4); + local=point(j,5); + if (local==0) + %convert to global axes + temp=[Fx,Fy,0,0,0,0]; + temp=(TransformationMatrix(joints(Near,1),joints(Near,2),joints(Far,1),joints(Far,2))')^-1*temp'; + Fx=temp(1); + Fy=temp(2); + end + if (x(i)>a) + N(i)-=Fx; + end + end + end + end +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/main/mechanics/inst/ocframe/PlotDiagrams.m Fri Sep 30 15:19:38 2011 +0000 @@ -0,0 +1,198 @@ +## Copyright (C) 2010 Johan Beke +## +## This software is free software; you can redistribute it and/or modify it +## under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 3 of the License, or (at +## your option) any later version. +## +## This software is distributed in the hope that it will be useful, but +## WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +## General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with this software; see the file COPYING. If not, see +## <http://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## @deftypefn {Function File} PlotDiagrams(@var{joints}, @var{members}, @var{dist}, @var{point}, @var{MemF}, @var{diagram}, @var{divisions}, @var{scale}) +## +## This function plots the internal forces for all members. The force to be plotted can be selected with @var{diagram} +## which will be "M", "S" or "N" for the moment, shear or normal forces. +## +## Input parameters are similar as with SolveFrame and PlotFrame. +## @end deftypefn + +function PlotDiagrams(joints,members,dist,point,MemF,diagram,divisions,scale) + %during development some errors occured on windows: http://wiki.octave.org/wiki.pl?OctaveForWindows see oct2mat and plot package + + %% scaling goes wrong when there al members have end moment 0 + if (nargin < 6) + usage("PlotDiagrams(joints,members,dist,point,MemF,diagram,scale) use the help command for more information"); + end + if (nargin<7) + divisions=20; + end + if (diagram=="M") + if (nargin<8) + scale=2/max(max(abs(MemF(:,[3,6])))) + end + %plot frame first + newplot(); + clf(); + hold off; + set (gca(), "dataaspectratio", [1,1,1]); + hold on; + + for i=1:rows(members) + Near=members(i,1); %near joint + Far=members(i,2); %far joint + %make line from near to far joint + line([joints(Near,1),joints(Far,1)],[joints(Near,2),joints(Far,2)],"color","blue"); + + end + %plot coord numbers + for i=1:rows(joints) + text(joints(i,1),joints(i,2),int2str(i)); %'color','blue' + hold on; + end + %plot member numbers + for i=1:rows(members) + text((joints(members(i,1),1)+joints(members(i,2),1))/2,(joints(members(i,1),2)+joints(members(i,2),2))/2,int2str(i),'color','blue'); + end + + for i=1:rows(members) + [x,M]=MSNForces(joints,members,dist,point,MemF,i,divisions); + %finding maximum and minimum + [M_max,x_max]=max(M); + [M_min,x_min]=min(M); + %scale example scale factor: 2/max moment + M=M.*scale; % moment curve consistent with deformation + Near=members(i,1); + Far=members(i,2); + + %transformation to real member position! + L=sqrt((joints(Near,1)-joints(Far,1))^2+(joints(Near,2)-joints(Far,2))^2); + l=(joints(Far,1)-joints(Near,1))/L; + m=(joints(Far,2)-joints(Near,2))/L; + T=[l,m;-m,l]'; + for i=1:columns(x) + temp=T*[x(i);M(i)]; %rotation + x(i)=temp(1)+joints(Near,1); %displacement x and y + M(i)=temp(2)+joints(Near,2); + end + plot(x,M,'r-'); + %plot max and min + text(x(x_max),M(x_max),num2str(M_max,'%g'),'color','red'); + text(x(x_min),M(x_min),num2str(M_min,'%g'),'color','red'); + end + end + if (diagram=="S") + if (nargin<8) + scale=2/max(max(abs(MemF(:,[2,5])))); + end + %plot frame first + newplot(); + clf(); + hold off; + set (gca(), "dataaspectratio", [1,1,1]); + hold on; + + for i=1:rows(members) + Near=members(i,1); %near joint + Far=members(i,2); %far joint + %make line from near to far joint + line([joints(Near,1),joints(Far,1)],[joints(Near,2),joints(Far,2)],"color","blue"); + + end + %plot coord numbers + for i=1:rows(joints) + text(joints(i,1),joints(i,2),int2str(i)); %'color','blue' + hold on; + end + %plot member numbers + for i=1:rows(members) + text((joints(members(i,1),1)+joints(members(i,2),1))/2,(joints(members(i,1),2)+joints(members(i,2),2))/2,int2str(i),'color','blue'); + end + + for i=1:rows(members) + [x,M,S]=MSNForces(joints,members,dist,point,MemF,i,divisions); + %finding maximum and minimum + [S_max,x_max]=max(S); + [S_min,x_min]=min(S); + %scale example scale factor: 2/max moment + S=S.*scale; + Near=members(i,1); + Far=members(i,2); + + %transformation to real member position! + L=sqrt((joints(Near,1)-joints(Far,1))^2+(joints(Near,2)-joints(Far,2))^2); + l=(joints(Far,1)-joints(Near,1))/L; + m=(joints(Far,2)-joints(Near,2))/L; + T=[l,m;-m,l]'; + for i=1:columns(x) + temp=T*[x(i);S(i)]; + x(i)=temp(1)+joints(Near,1); + S(i)=temp(2)+joints(Near,2); + end + plot(x,S,'r-'); + %plot max and min + text(x(x_max),S(x_max),num2str(S_max,'%g'),'color','red'); + text(x(x_min),S(x_min),num2str(S_min,'%g'),'color','red'); + end + end + if (diagram=="N") + if (nargin<8) + scale=2/max(max(abs(MemF(:,[1,4])))); + end + %plot frame first + newplot(); + clf(); + hold off; + set (gca(), "dataaspectratio", [1,1,1]); + hold on; + + for i=1:rows(members) + Near=members(i,1); %near joint + Far=members(i,2); %far joint + %make line from near to far joint + line([joints(Near,1),joints(Far,1)],[joints(Near,2),joints(Far,2)],"color","blue"); + + end + %plot coord numbers + for i=1:rows(joints) + text(joints(i,1),joints(i,2),int2str(i)); %'color','blue' + hold on; + end + %plot member numbers + for i=1:rows(members) + text((joints(members(i,1),1)+joints(members(i,2),1))/2,(joints(members(i,1),2)+joints(members(i,2),2))/2,int2str(i),'color','blue'); + end + + for i=1:rows(members) + [x,M,S,N]=MSNForces(joints,members,dist,point,MemF,i,divisions); + %finding maximum and minimum + [N_max,x_max]=max(N); + [N_min,x_min]=min(N); + %scale example scale factor: 2/max moment + N=N.*scale; + Near=members(i,1); + Far=members(i,2); + + %transformation to real member position! + L=sqrt((joints(Near,1)-joints(Far,1))^2+(joints(Near,2)-joints(Far,2))^2); + l=(joints(Far,1)-joints(Near,1))/L; + m=(joints(Far,2)-joints(Near,2))/L; + T=[l,m;-m,l]'; + for i=1:columns(x) + temp=T*[x(i);N(i)]; + x(i)=temp(1)+joints(Near,1); + N(i)=temp(2)+joints(Near,2); + end + plot(x,N,'r-'); + %plot max and min + text(x(x_max),N(x_max),num2str(N_max,'%g'),'color','red'); + text(x(x_min),N(x_min),num2str(N_min,'%g'),'color','red'); + end + end +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/main/mechanics/inst/ocframe/PlotFrame.m Fri Sep 30 15:19:38 2011 +0000 @@ -0,0 +1,78 @@ +## Copyright (C) 2010 Johan Beke +## +## This software is free software; you can redistribute it and/or modify it +## under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 3 of the License, or (at +## your option) any later version. +## +## This software is distributed in the hope that it will be useful, but +## WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +## General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with this software; see the file COPYING. If not, see +## <http://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## @deftypefn {Function File} {} PlotFrame(@var{joints}, @var{members}, @var{D}, @var{factor}) +## +## +## Plots a 2D frame (with displacements if needed) using +## the following input parameters: +## +## joints = [x , y, constraints ; ...] +## +## constraints=[x , y, rotation] free=0, supported=1 +## +## members = [nodeN, nodeF, E, I, A; ...] +## +## Optional arguments: +## +## D = [x,y,rotation;...] Displacements as returned by SolveFrame +## +## factor= Scaling factor for the discplacements (default: 10) +## +## @end deftypefn + +function PlotFrame(joints,members,D,factor) + if (nargin < 2) + usage("PlotFrame(joints,members,D,factor) D and factor ar optional. Use the help command for more information"); + end + %during development some errors occured on windows: http://wiki.octave.org/wiki.pl?OctaveForWindows see oct2mat and plot package + newplot(); + clf(); + set (gca(), "dataaspectratio", [1,1,1]); + for i=1:rows(members) + N=members(i,1); %near joint + F=members(i,2); %far joint + %make line from near to far joint + line([joints(N,1),joints(F,1)],[joints(N,2),joints(F,2)],"color","blue"); + hold on; + if (nargin>=3) + %plot also displacements + if (nargin==3) + % default factor + factor=10; + end + cN=node2code(N); + cF=node2code(F); + line([joints(N,1)+D(N,1)*factor,joints(F,1)+D(F,1)*factor],[joints(N,2)+D(N,2)*factor,joints(F,2)+D(F,2)*factor],"color","red"); + hold on; + end + end + %one meter extra as a border + set (gca (), "xlim", [min(joints(:,1))-1,max(joints(:,1))+1]); + set (gca (), "ylim", [min(joints(:,2))-1,max(joints(:,2))+1]); + + %plot coord numbers + for i=1:rows(joints) + text(joints(i,1),joints(i,2),int2str(i)); %'color','blue' + hold on; + end + %plot member numbers + for i=1:rows(members) + text((joints(members(i,1),1)+joints(members(i,2),1))/2,(joints(members(i,1),2)+joints(members(i,2),2))/2,int2str(i),'color','blue'); + end + hold off; +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/main/mechanics/inst/ocframe/SolveFrame.m Fri Sep 30 15:19:38 2011 +0000 @@ -0,0 +1,173 @@ +## Copyright (C) 2010 Johan Beke +## +## This software is free software; you can redistribute it and/or modify it +## under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 3 of the License, or (at +## your option) any later version. +## +## This software is distributed in the hope that it will be useful, but +## WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +## General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with this software; see the file COPYING. If not, see +## <http://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## @deftypefn {Function File} {[@var{Reactions}, @var{Displacements}, @var{MemF}]} = +## SolveFrame(@var{joints}, @var{members}, @var{nodeloads}, @var{dist}, @var{point}) +## +## +## Solves a 2D frame with the matrix displacement method for +## the following input parameters: +## +## joints = [x , y, constraints ; ...] +## +## constraints=[x , y, rotation] free=0, supported=1 +## +## members = [nodeN, nodeF, E, I, A; ...] +## +## nodeloads = [node, Fx, Fy, Mz; ...] +## +## loads on members: +## +## dist = [membernum,FxN,FyN,FxF,FyF,a,b,local ; ...] for distributed loads +## where FxN and FyN are the loads on distance a from the near node +## (same with far node and distance b) +## local=1 if loads are on local axis, 0 if global +## +## point = [membernum,Fx,Fy,a,local; ...] +## where Fx and Fy are the loads on distance a from the node near +## local=1 if loads are on local axis, 0 if global +## +## Output is formated as follows (rownumber corresponds to +## node or member number): +## +## Reactions = [Fx,Fy,Mz;...] where NaN it was a non supported dof +## +## Displacements = [x,y,rotation;...] +## +## MemF = [FxN, FyN, MzN, FxF, FyF, MzF; ...] +## @end deftypefn + +function [Reactions,Displacements,MemF]=SolveFrame(joints,members,nodeloads,dist,point) + if (nargin < 5) + usage("[Reactions,Displacements,MemF]=SolveFrame(joints,members,nodeloads,dist,point) Use the help command for more information"); + end + % Elements Axis + % y| + % | + % |___________________________ x + % Near far + %joints: [x, y, constraints; ...] 1= fixed + %members [nodeN,nodeF,E,I,A; ...] + %nodeloads [node,Fx,Fy,Mz; ...] + + P=D=zeros(rows(joints)*3,1); + %add nodal loads to P matrix + for load=1:rows(nodeloads) + c=node2code(nodeloads(load,1)); + for i=0:2 + P(c(1+i))=P(c(1+i))+nodeloads(load,2+i); + end + end + free=[]; %contains unconstrainted codenumbers + supported=[]; %contains constrainted codenumbers + for node=1:rows(joints) + c=node2code(node); + for i=3:5 + if (joints(node,i)==0) + free=[free,c(i-2)]; + else + supported=[supported,c(i-2)]; + end + end + end + + %% global equation + %% { P_f } [ K_ff K_fs ] { Delta_f } { P_F_f } + %% { } = [ ] . { } + { } + %% { P_s } [ K_sf K_ss ] { Delta_s } { P_F_s } + %% Solutions: + %% Delta_f = K_ff^-1 (P_f - P_F_f) + %% P_s = K_sf * Delta_f + P_F_s + + %stiffness matrix + tic + [ K_ff, K_sf ] = GlobalStiffnessMatrixRegrouped (joints, members,free,supported); %K_ss, K_fs are not used and if not calculated script is faster + toc + + %nodal forces and equivalent nodal ones + [P_F,MemFEF]=EquivalentJointLoads(joints,members,dist,point); + + + %reduced matrices + P_f=P(free,:); + P_s=P(supported,:); + + P_F_f=P_F(free,:); + P_F_s=P_F(supported,:); + + + %solution: find unknown displacements + %can be quicker for large systems with sparse matrix ! + + tic + [D_f, flag] = pcg(K_ff,P_f-P_F_f,1e-6,1000); + if (flag==1) + %max iteration + printf('max iteration'); + try + %A X = B => solve with cholesky => X = R\(R'\B) + r = chol (K_ff); + D_f=r\(r'\(P_f-P_F_f)); + catch + error("System is unstable because K_ff is singular. Please check the support constraints!"); + end_try_catch + end + if (flag==3) + %K_ff was found not positive defnite + error("System is unstable because K_ff is singular. Please check the support constraints!"); + end + toc + + + D(free,1)=D_f; + + %solution: find unknown (reaction) forces + P_s=K_sf*D_f+P_F_s; + P(supported,1)=P_s; + + + %find unknown member forces + MemF=[]; + for member=1:rows(members) + N=members(member,1); + F=members(member,2); + xN=joints(N,1); + yN=joints(N,2); + xF=joints(F,1); + yF=joints(F,2); + T=TransformationMatrix(xN,yN,xF,yF); + k=MemberStiffnessMatrix(sqrt((xN-xF)**2+(yN-yF)**2),members(member,3),members(member,4),members(member,5)); + c=[node2code(N),node2code(F)]; + MemF=[MemF;(k*T*D(c'))']; + end + MemF=MemF+MemFEF;%+fixed end forces + + %format codenumbers to real output + Displacements=[]; + Reactions=[]; + for i=0:rows(joints)-1 + Displacements=[Displacements;D(1+3*i:3+3*i)']; + Reactions=[Reactions;P(1+3*i:3+3*i)']; + end + for i=1:rows(Reactions) + for j=1:columns(Reactions) + if (joints(i,j+2)==0) + Reactions(i,j)=NaN; + end + end + end +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/main/mechanics/inst/ocframe/SolveFrameCases.m Fri Sep 30 15:19:38 2011 +0000 @@ -0,0 +1,158 @@ +## Copyright (C) 2010 Johan Beke +## +## This software is free software; you can redistribute it and/or modify it +## under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 3 of the License, or (at +## your option) any later version. +## +## This software is distributed in the hope that it will be useful, but +## WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +## General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with this software; see the file COPYING. If not, see +## <http://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## @deftypefn {Function File} {[@var{results}] =} SolveFrameCases (@var{joints}, @var{members}, @var{loadcases}) +## +## +## Solves a 2D frame with the matrix displacement method for +## the following input parameters: +## +## joints = [x , y, constraints ; ...] +## +## constraints=[x , y, rotation] free=0, supported=1 +## +## members = [nodeN, nodeF, E, I, A; ...] +## +## loadcases is a struct array with for each loadcase the fields +## +## - nodeloads = [node, Fx, Fy, Mz; ...] +## +## - dist = [membernum,FxN,FyN,FxF,FyF,a,b,local ; ...] +## +## - point = [membernum,Fx,Fy,a,local; ...] +## +## input is as for the function SolveFrame. +## +## Output is a struct array with the fields: Displacements, Reactions and MemF +## +## (output formated as for the function SolveFrame.) +## @end deftypefn + +function [results]=SolveFrameCases(joints,members,loadcases) + if (nargin != 3) + usage("[Reactions,Displacements,MemF]=SolveFrame(joints,members,loadcases) Use the help command for more information"); + end + + free=[]; %contains unconstrainted codenumbers + supported=[]; %contains constrainted codenumbers + for node=1:rows(joints) + c=node2code(node); + for i=3:5 + if (joints(node,i)==0) + free=[free,c(i-2)]; + else + supported=[supported,c(i-2)]; + end + end + end + + %stiffness matrix + [ K_ff, K_sf ] = GlobalStiffnessMatrixRegrouped (joints, members,free,supported); %K_ss, K_fs are not used and if not calculated script is faster + + %from now nodeloads,dist,point must depend on case: implement with ? !!!!!!!! + for lc=1:length(loadcases) + + %convert loadcases to other variables + nodeloads=loadcases(lc).nodeloads; + dist=loadcases(lc).dist; + point=loadcases(lc).point; + + P=D=zeros(rows(joints)*3,1); + %add nodal loads to P matrix + for load=1:rows(nodeloads) + c=node2code(nodeloads(load,1)); + for i=0:2 + P(c(1+i))=P(c(1+i))+nodeloads(load,2+i); + end + end + + %nodal forces and equivalent nodal ones + [P_F,MemFEF]=EquivalentJointLoads(joints,members,dist,point); + + + %reduced matrices + P_f=P(free,:); + P_s=P(supported,:); + + P_F_f=P_F(free,:); + P_F_s=P_F(supported,:); + + + %solution: find unknown displacements + %can be quicker for large systems with sparse matrix ! + + [D_f, flag] = pcg(K_ff,P_f-P_F_f,1e-6,1000); + if (flag==1) + %max iteration + printf('max iteration'); + try + %A X = B => solve with cholesky => X = R\(R'\B) + r = chol (K_ff); + D_f=r\(r'\(P_f-P_F_f)); + catch + error("System is unstable because K_ff is singular. Please check the support constraints!"); + end_try_catch + end + if (flag==3) + %K_ff was found not positive defnite + error("System is unstable because K_ff is singular. Please check the support constraints!"); + end + + + + D(free,1)=D_f; + + %solution: find unknown (reaction) forces + P_s=K_sf*D_f+P_F_s; + P(supported,1)=P_s; + + + %find unknown member forces + MemF=[]; + for member=1:rows(members) + N=members(member,1); + F=members(member,2); + xN=joints(N,1); + yN=joints(N,2); + xF=joints(F,1); + yF=joints(F,2); + T=TransformationMatrix(xN,yN,xF,yF); + k=MemberStiffnessMatrix(sqrt((xN-xF)**2+(yN-yF)**2),members(member,3),members(member,4),members(member,5)); + c=[node2code(N),node2code(F)]; + MemF=[MemF;(k*T*D(c'))']; + end + MemF=MemF+MemFEF;%+fixed end forces + + %format codenumbers to real output + Displacements=[]; + Reactions=[]; + for i=0:rows(joints)-1 + Displacements=[Displacements;D(1+3*i:3+3*i)']; + Reactions=[Reactions;P(1+3*i:3+3*i)']; + end + for i=1:rows(Reactions) + for j=1:columns(Reactions) + if (joints(i,j+2)==0) + Reactions(i,j)=NaN; + end + end + end + results(lc).Displacements=Displacements; + results(lc).Reactions=Reactions; + results(lc).MemF=MemF; + end +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/main/mechanics/inst/ocframe/ocframe_ex1.m Fri Sep 30 15:19:38 2011 +0000 @@ -0,0 +1,24 @@ +## -*- texinfo -*- +## @deftypefn {Function File} {} ocframe_ex1() +## Example of a planar frame. +## +## @end deftypefn +function [P,D,MemF]=ocframe_ex1() + joints=[0,0,1,1,1; + 4,4,0,0,0; + 8,4,1,1,0]; + + EIA=[210e9,23130*(10^-2)^4,84.5*(10^-2)^2];%IPE400 + + members=[1,2,EIA;2,3,EIA]; + + nodeloads=[]; + + dist=[1,0,-2e3,0,-2e3,0,0,0;2,0,-4e3,0,-4e3,0,0,1]; + point=[];%1,0,-3e4,3,1 + + [P,D,MemF]=SolveFrame(joints,members,nodeloads,dist,point); + PlotFrame(joints,members,D,10); + %plot moment diagram + PlotDiagrams(joints,members,dist,point,MemF,"S"); +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/main/mechanics/inst/ocframe/ocframe_ex2.m Fri Sep 30 15:19:38 2011 +0000 @@ -0,0 +1,36 @@ +## -*- texinfo -*- +## @deftypefn {Function File} {} ocframe_ex2() +## Example of a beam. +## +## @end deftypefn +function [P,D,MemF]=ocframe_ex2() + joints=[0,0,1,1,0; + 4,0,0,1,0; + 8,0,0,1,0]; + %Each 4 meter a node => 0,0 and 4,0 and 8,0 are the coordinates. + %The first node is a hinge and thus supported in the x and y direction => 1,1,0 for the constraints + %The following nodes are just rollers and thus supported in the y direction => 0,1,0 for the constraints + + EIA=[200e6,200e6*(10^-3)^4,6000*(10^-3)^2]; + %EIA as a single vector to be used afterwards + + members=[1,2,EIA;2,3,EIA]; + %2 members, connection node 1 to node 2 and node 2 to node 3 + + nodeloads=[]; + %there aren't nodal nodes in this example + + dist=[1,0,-4e3,0,-4e3,0,0,1; + 2,0,-4e3,0,-4e3,0,0,1]; + %both members have a distributed load which takes the full length of the member. Notice the - sign caused + %by the axes conventions + %as we are working with newtons and meters the load is -4e3 N and not -4 kN + + + [P,D,MemF]=SolveFrame(joints,members,nodeloads,dist,[]) + %solve the frame with + % P: reactions + % D: displacements + % MemF: the member end forces + +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/main/mechanics/inst/ocframe/ocframe_ex3.m Fri Sep 30 15:19:38 2011 +0000 @@ -0,0 +1,46 @@ +## -*- texinfo -*- +## @deftypefn {Function File} {} ocframe_ex3() +## Example of a planar frame. +## +## @end deftypefn +function [P,D,MemF]=ocframe_ex3() + UNP120=[210e9,363.990e-8,17.000e-4]; + RHS=[210e9,28.900e-8,8.730e-4]; + joints=[0,0,1,1,0; + 1.25,0,0,0,0; + 1.575,0,0,0,0; + 2.5,0,0,0,0; + 3.425,0,0,0,0; + 3.75,0,0,0,0; + 5.,0,0,1,0; + 0,1.5,0,0,0; + 1.25,1.5,0,0,0; + 2.5,1.5,0,0,0; + 3.75,1.5,0,0,0; + 5,1.5,0,0,0]; + members=[1,2,UNP120; + 2,3,UNP120; + 3,4,UNP120; + 4,5,UNP120; + 5,6,UNP120; + 6,7,UNP120; + 8,9,UNP120; + 9,10,UNP120; + 10,11,UNP120; + 11,12,UNP120; + 1,8,RHS; + 2,9,RHS; + 4,10,RHS; + 6,11,RHS; + 7,12,RHS]; + + nodeloads=[3,0,-30e3,0; + 5,0,-30e3,0; + 12,30e3,0,0]; + + + tic + [P,D,MemF]=SolveFrame(joints,members,nodeloads,[],[]); + toc + PlotFrame(joints,members,D,1); +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/main/mechanics/inst/ocframe/ocframe_exLC.m Fri Sep 30 15:19:38 2011 +0000 @@ -0,0 +1,68 @@ +## -*- texinfo -*- +## @deftypefn {Function File} {} ocframe_exLC() +## Example of a beam with generation of eurocode ULS +## load cases +## +## @end deftypefn + +function ocframe_exLC() + + joints=[0,0,1,1,0; + 4,0,0,1,0; + 8,0,0,1,0; + 12,0,0,1,0; + 16,0,0,1,0]; + + EIA=[200e6,200e6*(10^-3)^4,6000*(10^-3)^2]; + members=[1,2,EIA; + 2,3,EIA; + 3,4,EIA; + 4,5,EIA]; + + %there is a variable load 4kN/m to be combined with 2kN/m permanent + %eurocode ULS + %Permanent actions: if favourable 1 if not 1,35 + %Variable actions: if favourable 0 if not 1,5 + %Category A: domestic phi0 = 0,7 phi1 = 0,5 phi2 = 0,3 + % each case consists of permanent load * 1,35 plus one of the variable cases + % which makes 16 cases (try for each variable load 0 and 1,5 as the factor) + loadcases={}; + for i=0:15 + lc=toascii(dec2bin(i,4)).-48; + %no point loads and nodal loads + loadcases(i+1).nodeloads=[]; + loadcases(i+1).point=[]; + %dist load depends on case + dist=[]; + for j=1:4 + %add permanent + dist=[dist;[j,0,-2e3*1.35,0,-2e3*1.35,0,0,1]]; + if (lc(j)==1) + dist=[dist;[j,0,-4e3*0.7*1.5,0,-4e3*0.7*1.5,0,0,1]]; + end + loadcases(i+1).dist=dist; + end + end + + [results]=SolveFrameCases(joints,members,loadcases); + + %moment diagram envelope + moments=[]; + for lc=1:16 + x=[]; + M=[]; + for i=1:4 + [u,Mo]=MSNForces(joints,members,loadcases(lc).dist,loadcases(lc).point,results(lc).MemF,i,40); + x=[x,u.+(4*(i-1))]; + M=[M,Mo]; + end + moments=[moments;M]; + end + + plot(x,max(moments),x,min(moments)) + +end + + + + \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/main/mechanics/inst/ocframe/ocframe_railwaybridge.m Fri Sep 30 15:19:38 2011 +0000 @@ -0,0 +1,46 @@ +## -*- texinfo -*- +## @deftypefn {Function File} {} ocframe_railwaybridge() +## Example taken from a real railwaybridge. +## +## @end deftypefn +function [P,D,MemF]=ocframe_railwaybridge() + joints=[0:5.3636:59;zeros(4,12)]'; + joints(1,3:4)=[1,1]; + joints(12,3:4)=[1,1]; + + temp=[]; + %parabola with height = 8 m => f(x)=(32*x)/59-(32*x^2)/3481 + for i=2:11 + temp=[temp;joints(i,1),(32*joints(i,1))/59-(32*joints(i,1)^2)/3481,0,0,0]; + end + + joints=[joints;temp]; + + %EIA + beam1=[200e9,135.7*(0.1)^4,3.204*(0.1)^2]; + beam2=[200e9,80.5*(0.1)^4,2.011*(0.1)^2]; + members=[]; + for i=1:11 + members=[members;i,i+1,beam1]; + end + members=[members;1,13,beam1]; + for i=13:21 + members=[members;i,i+1,beam1]; + end + members=[members;22,12,beam1]; + + for i=2:11 + members=[members;i,i+11,beam2]; + end + + %own weight of beams neglected + + %some forces + nodeloads=[6,0,-50e3,0;7,0,-50e3,0]; + + tic + [P,D,MemF]=SolveFrame(joints,members,nodeloads,[],[]); + toc + + %PlotFrame(joints,members,D,100); +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/main/mechanics/inst/ocframe/private/EquivalentJointLoads.m Fri Sep 30 15:19:38 2011 +0000 @@ -0,0 +1,86 @@ +## Copyright (C) 2010 Johan Beke +## +## This software is free software; you can redistribute it and/or modify it +## under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 3 of the License, or (at +## your option) any later version. +## +## This software is distributed in the hope that it will be useful, but +## WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +## General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with this software; see the file COPYING. If not, see +## <http://www.gnu.org/licenses/>. + +function [P_F,MemFEF]=EquivalentJointLoads(joints,members,dist,point) + % joints: [x, y, constraints; ...] 1= fixed + % members [nodeN,nodeF,E,I,A; ...] + % dist containts distributed loads on members + % [membernum,FxN,FyN,FxF,FyF,a,b,local ; ...] + % point contains the point loads on members + % [membernum,Fx,Fy,a,local; ...] + % local=1 if in local coords + + P_F=zeros(rows(joints)*3,1); + MemFEF=zeros(rows(members),6); + %distributed loads + for i=1:rows(dist) + membernum=dist(i,1); + FxN=dist(i,2); + FyN=dist(i,3); + FxF=dist(i,4); + FyF=dist(i,5); + a=dist(i,6); + b=dist(i,7); + local=dist(i,8); + N=members(membernum,1); + F=members(membernum,2); + if (local==0) + %convert FxN,... to local axes + temp=[FxN,FyN,0,FxF,FyF,0]; + temp=(TransformationMatrix(joints(N,1),joints(N,2),joints(F,1),joints(F,2))')^-1*temp'; + FxN=temp(1); + FyN=temp(2); + FxF=temp(4); + FyF=temp(5); + end + [F1,F2,F3,F4,F5,F6]=FixedEndForcesDist(sqrt((joints(N,1)-joints(F,1))^2+(joints(F,2)-joints(N,2))^2),FxN,FyN,FxF,FyF,a,b); + F_local=[F1,F2,F3,F4,F5,F6]; + F_global=((TransformationMatrix(joints(N,1),joints(N,2),joints(F,1),joints(F,2))')*F_local')'; + + c=[node2code(N),node2code(F)]; + for j=1:6 + P_F(c(j))+=F_global(j); + MemFEF(membernum,j)+=F_local(j); + end + end + for i=1:rows(point) + %[membernum,Fx,Fy,a,local; ...] + membernum=point(i,1); + Fx=point(i,2); + Fy=point(i,3); + a=point(i,4); + local=point(i,5); + N=members(membernum,1); + F=members(membernum,2); + if (local==0) + %convert to global axes + temp=[Fx,Fy,0,0,0,0]; + temp=(TransformationMatrix(joints(N,1),joints(N,2),joints(F,1),joints(F,2))')^-1*temp'; + Fx=temp(1); + Fy=temp(2); + end + %FixedEndForcesPnt(l,a,Fx,Fy) + [F1,F2,F3,F4,F5,F6]=FixedEndForcesPnt(sqrt((joints(N,1)-joints(F,1))^2+(joints(F,2)-joints(N,2))^2),a,Fx,Fy); + F_local=[F1,F2,F3,F4,F5,F6]; + F_global=((TransformationMatrix(joints(N,1),joints(N,2),joints(F,1),joints(F,2))')*F_local')'; + + c=[node2code(N),node2code(F)]; + for j=1:6 + P_F(c(j))+=F_global(j); + MemFEF(membernum,j)+=F_local(j); + end + end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/main/mechanics/inst/ocframe/private/FixedEndForcesDist.m Fri Sep 30 15:19:38 2011 +0000 @@ -0,0 +1,39 @@ +## Copyright (C) 2010 Johan Beke +## +## This software is free software; you can redistribute it and/or modify it +## under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 3 of the License, or (at +## your option) any later version. +## +## This software is distributed in the hope that it will be useful, but +## WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +## General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with this software; see the file COPYING. If not, see +## <http://www.gnu.org/licenses/>. + +function [F1,F2,F3,F4,F5,F6]=FixedEndForcesDist(l,FxN,FyN,FxF,FyF,a,b) + % distributed load in local coords. + % a: distance from Near + % b: distance from Far + % FxN, FyN forces on start of load in local coordinates + % FxF, FyF forces on end of load in local coordinates + + % formulas found by integrating the formulas for a point load with maxima maxima.sourceforge.net + F1=-((2*FxN+FxF)*l^2+((-b-4*a)*FxN+(b-2*a)*FxF)*l+(-b^2+a*b+2*a^2)*FxN+(-2*b^2-a*b+a^2)*FxF)/(6*l); + F2=-1*(((7*FyN+3*FyF)*l^4+((-3*b-13*a)*FyN+(3*b-7*a)*FyF)*l^3+((-3*b^2+4*a*b-3*a^2)*FyN+(3*b^2-4*a*b+3*a^2)*FyF)*l^2+ +((-3*b^3+a*b^2+a^2*b+17*a^3)*FyN+(-17*b^3-a*b^2-a^2*b+3*a^3)*FyF)*l+(2*b^4-2*a*b^3+2*a^2*b^2-2*a^3*b-8*a^4)*FyN+(8*b^4+2*a*b^3-2*a^2*b^2+2*a^3*b-2*a^4)* +FyF)/(20*l^3)); + F3=-1*(((3*FyN+2*FyF)*l^4+((3*a-2*b)*FyN+(2*b-3*a)*FyF)*l^3+((-2*b^2+a*b-27*a^2)*FyN+(2*b^2-a*b-3*a^2)*FyF)*l^2+ +((-2*b^3-a*b^2+4*a^2*b+33*a^3)*FyN+(-18*b^3+a*b^2-4*a^2*b+7*a^3)*FyF)*l+(3*b^4-3*a*b^3+3*a^2*b^2-3*a^3*b-12*a^4)*FyN+ +(12*b^4+3*a*b^3-3*a^2*b^2+3*a^3*b-3*a^4)*FyF)/(60*l^2)); + F4=-((FxN+2*FxF)*l^2+((a-2*b)*FxN+(-4*b-a)*FxF)*l+(b^2-a*b-2*a^2)*FxN+(2*b^2+a*b-a^2)*FxF)/(6*l); + F5=-1*(((3*FyN+7*FyF)*l^4+((3*a-7*b)*FyN+(-13*b-3*a)*FyF)*l^3+((3*b^2-4*a*b+3*a^2)*FyN+(-3*b^2+4*a*b-3*a^2)*FyF)*l^2+ +((3*b^3-a*b^2-a^2*b-17*a^3)*FyN+(17*b^3+a*b^2+a^2*b-3*a^3)*FyF)*l+(-2*b^4+2*a*b^3-2*a^2*b^2+2*a^3*b+8*a^4)*FyN+(-8*b^4-2*a*b^3+2*a^2*b^2-2*a^3*b+2*a^4)* +FyF)/(20*l^3)); + F6=-1*(-((2*FyN+3*FyF)*l^4+((2*a-3*b)*FyN+(3*b-2*a)*FyF)*l^3+((-3*b^2-a*b+2*a^2)*FyN+(-27*b^2+a*b-2*a^2)*FyF)*l^2+ +((7*b^3-4*a*b^2+a^2*b-18*a^3)*FyN+(33*b^3+4*a*b^2-a^2*b-2*a^3)*FyF)*l+(-3*b^4+3*a*b^3-3*a^2*b^2+3*a^3*b+12*a^4)*FyN+ +(-12*b^4-3*a*b^3+3*a^2*b^2-3*a^3*b+3*a^4)*FyF)/(60*l^2)); +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/main/mechanics/inst/ocframe/private/FixedEndForcesPnt.m Fri Sep 30 15:19:38 2011 +0000 @@ -0,0 +1,33 @@ +## Copyright (C) 2010 Johan Beke +## +## This software is free software; you can redistribute it and/or modify it +## under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 3 of the License, or (at +## your option) any later version. +## +## This software is distributed in the hope that it will be useful, but +## WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +## General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with this software; see the file COPYING. If not, see +## <http://www.gnu.org/licenses/>. + +function [F1,F2,F3,F4,F5,F6]=FixedEndForcesPnt(l,a,Fx,Fy) + % point load in local coords. + % a: distance from Near + % l: length + % FxN, FyN forces on start of load in local coordinates + % FxF, FyF forces on end of load in local coordinates + + % b: distance from Far + b=l-a; + + F1=-b*Fx/l; + F2=-Fy*(b/l-a^2*b/l^3+a*b^2/l^3); + F3=-Fy*a*b^2/l^2; + F4=-a*Fx/l; + F5=-Fy*(a/l+a^2*b/l^3-a*b^2/l^3); + F6=Fy*a^2*b/l^2; +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/main/mechanics/inst/ocframe/private/GlobalStiffnessMatrix.m Fri Sep 30 15:19:38 2011 +0000 @@ -0,0 +1,85 @@ +## Copyright (C) 2010 Johan Beke +## +## This software is free software; you can redistribute it and/or modify it +## under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 3 of the License, or (at +## your option) any later version. +## +## This software is distributed in the hope that it will be useful, but +## WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +## General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with this software; see the file COPYING. If not, see +## <http://www.gnu.org/licenses/>. + +function K=GlobalStiffnessMatrix(joints,members) + #joints: [x, y, constraints; ...] 1= fixed + #members [nodeN,nodeF,E,I,A; ...] + K=sparse(rows(joints)*3,rows(joints)*3); + for member=1:rows(members) + N=members(member,1); + F=members(member,2); + xN=joints(N,1); + yN=joints(N,2); + xF=joints(F,1); + yF=joints(F,2); + %T=TransformationMatrix(xN,yN,xF,yF); + %k=transpose(T)*MemberStiffnessMatrix(sqrt((xN-xF)**2+(yN-yF)**2),members(member,3),members(member,4),members(member,5))*T; + c=[node2code(N),node2code(F)]; + + %TransformationMatrix stuff + L=sqrt((xN-xF)^2+(yN-yF)^2); + l=(xF-xN)/L; + m=(yF-yN)/L; + %MemberStiffnessMatrix stuff + E=members(member,3); + I=members(member,4); + A=members(member,5); + + %for i=1:rows(k) + % for j=1:columns(k) + % K(c(i),c(j))=K(c(i),c(j))+k(i,j); + % end + %end + + %fill Global Stiffness Matrix + K(c(1),c(1))+=(l^2*A*E)/L+(12*m^2*E*I)/L^3; + K(c(1),c(2))+=(l*m*A*E)/L-(12*l*m*E*I)/L^3; + K(c(1),c(3))+=-(6*m*E*I)/L^2; + K(c(1),c(4))+=-(l^2*A*E)/L-(12*m^2*E*I)/L^3; + K(c(1),c(5))+=(12*l*m*E*I)/L^3-(l*m*A*E)/L; + K(c(1),c(6))+=-(6*m*E*I)/L^2; + K(c(2),c(1))+=(l*m*A*E)/L-(12*l*m*E*I)/L^3; + K(c(2),c(2))+=(m^2*A*E)/L+(12*l^2*E*I)/L^3; + K(c(2),c(3))+=(6*l*E*I)/L^2; + K(c(2),c(4))+=(12*l*m*E*I)/L^3-(l*m*A*E)/L; + K(c(2),c(5))+=-(m^2*A*E)/L-(12*l^2*E*I)/L^3; + K(c(2),c(6))+=(6*l*E*I)/L^2; + K(c(3),c(1))+=-(6*m*E*I)/L^2; + K(c(3),c(2))+=(6*l*E*I)/L^2; + K(c(3),c(3))+=(4*E*I)/L; + K(c(3),c(4))+=(6*m*E*I)/L^2; + K(c(3),c(5))+=-(6*l*E*I)/L^2; + K(c(3),c(6))+=(2*E*I)/L; + K(c(4),c(1))+=-(l^2*A*E)/L-(12*m^2*E*I)/L^3; + K(c(4),c(2))+=(12*l*m*E*I)/L^3-(l*m*A*E)/L; + K(c(4),c(3))+=(6*m*E*I)/L^2; + K(c(4),c(4))+=(l^2*A*E)/L+(12*m^2*E*I)/L^3; + K(c(4),c(5))+=(l*m*A*E)/L-(12*l*m*E*I)/L^3; + K(c(4),c(6))+=(6*m*E*I)/L^2; + K(c(5),c(1))+=(12*l*m*E*I)/L^3-(l*m*A*E)/L; + K(c(5),c(2))+=-(m^2*A*E)/L-(12*l^2*E*I)/L^3; + K(c(5),c(3))+=-(6*l*E*I)/L^2; + K(c(5),c(4))+=(l*m*A*E)/L-(12*l*m*E*I)/L^3; + K(c(5),c(5))+=(m^2*A*E)/L+(12*l^2*E*I)/L^3; + K(c(5),c(6))+=-(6*l*E*I)/L^2; + K(c(6),c(1))+=-(6*m*E*I)/L^2; + K(c(6),c(2))+=(6*l*E*I)/L^2; + K(c(6),c(3))+=(2*E*I)/L; + K(c(6),c(4))+=(6*m*E*I)/L^2; + K(c(6),c(5))+=-(6*l*E*I)/L^2; + K(c(6),c(6))+=(4*E*I)/L; + end +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/main/mechanics/inst/ocframe/private/GlobalStiffnessMatrixModified.m Fri Sep 30 15:19:38 2011 +0000 @@ -0,0 +1,48 @@ +## Copyright (C) 2010 Johan Beke +## +## This software is free software; you can redistribute it and/or modify it +## under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 3 of the License, or (at +## your option) any later version. +## +## This software is distributed in the hope that it will be useful, but +## WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +## General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with this software; see the file COPYING. If not, see +## <http://www.gnu.org/licenses/>. + +function [ Kmod, Kcor ] = GlobalStiffnessMatrixModified(joints,members,supported) + %% \cite{Felippa2010} : "To apply support conditions without rearranging the equations we clear (set to zero) rows and columns + %% corresponding to prescribed zero displacements as well as the corresponding force components, and place + %% ones on the diagonal to maintain non-singularity. The resulting system is called the modified set of master + %% stiffness equations. " + + %% joints: [x, y, constraints; ...] 1= fixed + %% members [nodeN, nodeF, E, I, A; ...] + %% supported: vector with the supported (fixed) code numbers + + %% Kmod + Kcor = K + + %% Kcor is the correction matrix to make K + + Kmod=GlobalStiffnessMatrix(joints,members); + Kcor=sparse(rows(joints)*3,rows(joints)*3); + + %% Creat Kcor + for i=1:columns(supported) + Kcor(:,i)=Kmod(:,i); + Kcor(i,:)=Kmod(i,:); + end + %% delete rows and columns. This might be not the fastest method. + for i=1:columns(supported) + Kmod(:,i)=zeros(rows(joints)*3,1); + Kmod(i,:)=zeros(1,rows(joints)*3); + + %set 1 and -1 in Kmod and Kcor + Kmod(i,i)+=1; + Kcor(i,i)-=1; + end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/main/mechanics/inst/ocframe/private/GlobalStiffnessMatrixRegrouped.m Fri Sep 30 15:19:38 2011 +0000 @@ -0,0 +1,54 @@ +## Copyright (C) 2010 Johan Beke +## +## This software is free software; you can redistribute it and/or modify it +## under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 3 of the License, or (at +## your option) any later version. +## +## This software is distributed in the hope that it will be useful, but +## WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +## General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with this software; see the file COPYING. If not, see +## <http://www.gnu.org/licenses/>. + +function [ K_ff, K_sf, K_ss, K_fs ] = GlobalStiffnessMatrixRegrouped (joints, members,free,supported) + %% Returns the components of the global stiffness matrix regrouped as in the following equation: + %% { P_f } [ K_ff K_fs ] { Delta_f } + %% { } = [ ] . { } + %% { P_s } [ K_sf K_ss ] { Delta_s } + + %% joints: [x, y, constraints; ...] 1= fixed + %% members [nodeN, nodeF, E, I, A; ...] + %% free: vector with the free code numbers + %% supported: vector with the supported (fixed) code numbers + + %zeros(row,col) + + %number of free en supported joints + + %!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + %TODO: for large systems this function is still slow + %!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + K=GlobalStiffnessMatrix(joints,members); + + %K_ff + K_ff=K(free,free); + + %K_sf + K_sf=K(supported,free); + + + %K_ss + if (nargout>=3) + K_ss=K(supported,supported); + end + + %K_fs + if (nargout==4) + K_fs=K(free,supported); + end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/main/mechanics/inst/ocframe/private/GlobalStiffnessMatrixRegrouped.m.old Fri Sep 30 15:19:38 2011 +0000 @@ -0,0 +1,98 @@ +## Copyright (C) 2010 Johan Beke +## +## This software is free software; you can redistribute it and/or modify it +## under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 3 of the License, or (at +## your option) any later version. +## +## This software is distributed in the hope that it will be useful, but +## WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +## General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with this software; see the file COPYING. If not, see +## <http://www.gnu.org/licenses/>. + +function [ K_ff, K_sf, K_ss, K_fs ] = GlobalStiffnessMatrixRegrouped (joints, members,free,supported) + %% Returns the components of the global stiffness matrix regrouped as in the following equation: + %% { P_f } [ K_ff K_fs ] { Delta_f } + %% { } = [ ] . { } + %% { P_s } [ K_sf K_ss ] { Delta_s } + + %% joints: [x, y, constraints; ...] 1= fixed + %% members [nodeN, nodeF, E, I, A; ...] + %% free: vector with the free code numbers + %% supported: vector with the supported (fixed) code numbers + + %zeros(row,col) + + %number of free en supported joints + + %!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + %TODO: for large systems this function is still slow + % due to the indexof! + %!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + f=columns(free); + s=columns(supported); + + K_ff=spalloc(f,f); + K_ss=spalloc(s,s); + K_sf=spalloc(s,f); + K_fs=spalloc(f,s); + for member=1:rows(members) + N=members(member,1); + F=members(member,2); + xN=joints(N,1); + yN=joints(N,2); + xF=joints(F,1); + yF=joints(F,2); + T=TransformationMatrix(xN,yN,xF,yF); + k=transpose(T)*MemberStiffnessMatrix(sqrt((xN-xF)**2+(yN-yF)**2),members(member,3),members(member,4),members(member,5))*T;%global coords. + c=[node2code(N),node2code(F)];%code numbers of the near and far nodes of the current member + + for i=1:rows(k) + for j=1:columns(k) + if (k(i,j)!=0) + c_row=c(i); + c_col=c(j); + + %K_ff + row=indexof(free,c_row); + col=indexof(free,c_col); + if (row!=-1 && col!=-1) + K_ff(row,col)=K_ff(row,col)+k(i,j); + end + + %K_sf + row=indexof(supported,c_row); + col=indexof(free,c_col); + if (row!=-1 && col!=-1) + K_sf(row,col)=K_sf(row,col)+k(i,j); + end + + + %K_ss + if (nargout>=3) + row=indexof(supported,c_row); + col=indexof(supported,c_col); + if (row!=-1 && col!=-1) + K_ss(row,col)=K_ss(row,col)+k(i,j); + end + end + + %K_fs + if (nargout==4) + row=indexof(free,c_row); + col=indexof(supported,c_col); + if (row!=-1 && col!=-1) + K_fs(row,col)=K_fs(row,col)+k(i,j); + end + end + end + end + end + + end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/main/mechanics/inst/ocframe/private/MemberStiffnessMatrix.m Fri Sep 30 15:19:38 2011 +0000 @@ -0,0 +1,24 @@ +## Copyright (C) 2010 Johan Beke +## +## This software is free software; you can redistribute it and/or modify it +## under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 3 of the License, or (at +## your option) any later version. +## +## This software is distributed in the hope that it will be useful, but +## WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +## General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with this software; see the file COPYING. If not, see +## <http://www.gnu.org/licenses/>. + +function k=MemberStiffnessMatrix(L,E,I,A) + k=[A*E/L,0,0,-A*E/L,0,0; + 0,12*E*I/L^3,6*E*I/L^2,0,-12*E*I/L^3,6*E*I/L^2; + 0,6*E*I/L^2,4*E*I/L,0,-6*E*I/L^2,2*E*I/L; + -A*E/L,0,0,A*E/L,0,0; + 0,-12*E*I/L^3,-6*E*I/L^2,0,12*E*I/L^3,-6*E*I/L^2; + 0,6*E*I/L^2,2*E*I/L,0,-6*E*I/L^2,4*E*I/L]; +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/main/mechanics/inst/ocframe/private/TransformationMatrix.m Fri Sep 30 15:19:38 2011 +0000 @@ -0,0 +1,29 @@ +## Copyright (C) 2010 Johan Beke +## +## This software is free software; you can redistribute it and/or modify it +## under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 3 of the License, or (at +## your option) any later version. +## +## This software is distributed in the hope that it will be useful, but +## WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +## General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with this software; see the file COPYING. If not, see +## <http://www.gnu.org/licenses/>. + +function T=TransformationMatrix(x1,y1,x2,y2) + #Displacement transformation matrix = T + #Force transformation = T transposed + L=sqrt((x1-x2)^2+(y1-y2)^2); + l=(x2-x1)/L; + m=(y2-y1)/L; + T=[l,m,0,0,0,0; + -m,l,0,0,0,0; + 0,0,1,0,0,0; + 0,0,0,l,m,0; + 0,0,0,-m,l,0; + 0,0,0,0,0,1]; +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/main/mechanics/inst/ocframe/private/indexof.m Fri Sep 30 15:19:38 2011 +0000 @@ -0,0 +1,25 @@ +## Copyright (C) 2010 Johan Beke +## +## This software is free software; you can redistribute it and/or modify it +## under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 3 of the License, or (at +## your option) any later version. +## +## This software is distributed in the hope that it will be useful, but +## WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +## General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with this software; see the file COPYING. If not, see +## <http://www.gnu.org/licenses/>. + +function [index] = indexof(vector,number) + index=-1; + for i=1:columns(vector) + if (vector(i)==number) + index=i; + break + end + end +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/main/mechanics/inst/ocframe/private/matrix2tex.m Fri Sep 30 15:19:38 2011 +0000 @@ -0,0 +1,33 @@ +## Copyright (C) 2010 Johan Beke +## +## This software is free software; you can redistribute it and/or modify it +## under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 3 of the License, or (at +## your option) any later version. +## +## This software is distributed in the hope that it will be useful, but +## WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +## General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with this software; see the file COPYING. If not, see +## <http://www.gnu.org/licenses/>. + +function matrix2tex(m) + printf("[tex]\\left(\\begin{array}{") + for c=1:columns(m) + printf("c") + end + printf("}\n") + for r=1:rows(m) + for c=1:columns(m) + if (c!=columns(m)) + printf(" %3.3e &",m(r,c)) + else + printf(" %3.3e \\\\\n",m(r,c)) + end + end + end + printf("\\end{array}\\right)[/tex]\n\n") +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/main/mechanics/inst/ocframe/private/node2code.m Fri Sep 30 15:19:38 2011 +0000 @@ -0,0 +1,19 @@ +## Copyright (C) 2010 Johan Beke +## +## This software is free software; you can redistribute it and/or modify it +## under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 3 of the License, or (at +## your option) any later version. +## +## This software is distributed in the hope that it will be useful, but +## WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +## General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with this software; see the file COPYING. If not, see +## <http://www.gnu.org/licenses/>. + +function c=node2code(node) + c=(node-1)*3+1:(node-1)*3+3; +end \ No newline at end of file