diff --git a/CEopt-1.0/TrussConst.m b/CEopt-1.0/TrussConst.m new file mode 100644 index 0000000..5d6524f --- /dev/null +++ b/CEopt-1.0/TrussConst.m @@ -0,0 +1,78 @@ +% ----------------------------------------------------------------- +% Truss10Const.m +% ----------------------------------------------------------------- +% programmer: Marcos Vinicius Issa +% marcos.issa@uerj.br +% +% Originally programmed in: Apr 04, 2024 +% Last updated in: Jul 31, 2024 +% ----------------------------------------------------------------- +% Constraints for truss optimization problem. +% ----------------------------------------------------------------- +function [G,H] = TrussConst(A,MyTruss) + + % truss structure parameters + rho = MyTruss.rho; + E = MyTruss.E; + AddedMass = MyTruss.AddedMass; + NODES = MyTruss.NODES; + ELEM = MyTruss.ELEM; + Nnodes = MyTruss.Nnodes; + Nelem = MyTruss.Nelem; + Ndofs = MyTruss.Ndofs; + + % initialize inequatily constraint + G = zeros(1,3); + + % initialize equatily constraint + H = zeros(1,Ndofs); + + % preallocate memory for global matrices + K = zeros(Ndofs,Ndofs); + M = zeros(Ndofs,Ndofs); + + % assembly global matrices + for e = 1:Nelem + % distance between NODESs + dx = NODES(ELEM(e,2),1) - NODES(ELEM(e,1),1); + dy = NODES(ELEM(e,2),2) - NODES(ELEM(e,1),2); + l = sqrt(dx^2+dy^2); + % strain matrix + c = dx/l; + s = dy/l; + B = 1/l*[-c -s c s]; + % element DoFs + eDof = [2*ELEM(e,1)-1, 2*ELEM(e,1),... + 2*ELEM(e,2)-1, 2*ELEM(e,2)]; + % local stiffness and mass matrices + Ke = B'*E*B*A(e)*l; + Me = rho*l*A(e)*[2 0 1 0; + 0 2 0 1; + 1 0 2 0; + 0 1 0 2]/6; + % assembly global stiffness and mass matrices + K(eDof,eDof) = K(eDof,eDof) + Ke; + M(eDof,eDof) = M(eDof,eDof) + Me; + end + + % update mass matrix with the added mass + M = M + AddedMass*eye(Ndofs,Ndofs); + + % fixed and free nodes + FIX = union(2*[5 6]-1,2*[5 6]); + FREE = setdiff(1:2*Nnodes,FIX); + + % solve the generalized eigenvalue problem + [V,D] = eig(K(FREE,FREE),M(FREE,FREE)); + + % sort frequencies + [omega2,Idx] = sort(diag(D)); + + % equilibrium constraints + H = 0; + + % frequency constraints + G(1) = 1 - sqrt(omega2(1))/(2*pi*7); + G(2) = 1 - sqrt(omega2(2))/(2*pi*15); + G(3) = 1 - sqrt(omega2(3))/(2*pi*20); +% ----------------------------------------------------------------- \ No newline at end of file diff --git a/CEopt-1.0/TrussMass.m b/CEopt-1.0/TrussMass.m new file mode 100644 index 0000000..ad01574 --- /dev/null +++ b/CEopt-1.0/TrussMass.m @@ -0,0 +1,29 @@ +% ----------------------------------------------------------------- +% Truss10Mass.m +% ----------------------------------------------------------------- +% programmer: Marcos Vinicius Issa +% marcos.issa@uerj.br +% +% Originally programmed in: Apr 04, 2024 +% Last updated in: Jul 31, 2024 +% ----------------------------------------------------------------- +% Objective function for truss optimization problem. +% ----------------------------------------------------------------- +function M = TrussMass(A,MyTruss) + + % truss structure parameters + rho = MyTruss.rho; + NODES = MyTruss.NODES; + ELEM = MyTruss.ELEM; + Nelem = MyTruss.Nelem; + + % compute mass + M = 0.0; + for e = 1:Nelem + dx = NODES(ELEM(e,2),1) - NODES(ELEM(e,1),1); + dy = NODES(ELEM(e,2),2) - NODES(ELEM(e,1),2); + l = sqrt(dx^2+dy^2); + M = M + A(e)*l*rho; + end +end +% ----------------------------------------------------------------- \ No newline at end of file diff --git a/CEopt-1.0/TrussPlot.m b/CEopt-1.0/TrussPlot.m new file mode 100644 index 0000000..abec337 --- /dev/null +++ b/CEopt-1.0/TrussPlot.m @@ -0,0 +1,126 @@ +% ----------------------------------------------------------------- +% TrussPlot.m +% ----------------------------------------------------------------- +% programmer: Marcos Vinicius Issa +% marcos.issa@uerj.br +% +% Originally programmed in: Apr 04, 2024 +% Last updated in: Jul 31, 2024 +% ----------------------------------------------------------------- +% Plot a truss structure given the elements cross-section area. +% ----------------------------------------------------------------- +function TrussPlot(x,MyTitle) + + % inch to meter convertion factor + Inch2Meter = 0.0254; + + % truss geometric dimensions (m) + l1 = 360*Inch2Meter; % 1st length + l2 = 360*Inch2Meter; % 2nd length + h = 360*Inch2Meter; % heigth + + % mesh points in global coordinates (Nnodes x 2) + XY_Mesh = [ (l1+l2) h; + (l1+l2) 0; + l1 h; + l1 0; + 0 l1; + 0 0]; + + % global coordinates of element nodes (Nelem x Nlnodes) + GlobalNodeID = [3 5; + 1 3; + 4 6; + 2 4; + 3 4; + 1 2; + 4 5; + 3 6; + 2 3; + 1 4]; + + grayColor = [.7 .7 .7]; + + figure('DefaultAxesFontSize',10) + clf + hold on + + % Support 1 + xl1 = [-1.0 -0.7] ; yl1 = [-1.0 1.0]; + [X1,Y1] = hatch_coordinates(xl1, yl1, 0.15) ; + plot(X1,Y1,'red','linewidth',1.0); + a11 = [0; -0.7; -0.7]; + a12 = [0; 1.0; -1.0]; + patch(a11,a12,grayColor,'EdgeColor','none'); + plot([a11(2) a11(3)],[a12(2) a12(3)],'red','linewidth',2.5); + + % Support 2 + xl2 = [-1.0 -0.7] ; yl2 = [h-1.0 h+1.0]; + [X2,Y2] = hatch_coordinates(xl2, yl2, 0.15) ; + plot(X2,Y2,'red','linewidth',1.0); + a21 = [0; -0.7; -0.7]; + a22 = [h; h+1; h-1]; + patch(a21,a22,grayColor,'EdgeColor','none'); + plot([a21(2) a21(3)],[a22(2) a22(3)],'red','linewidth',2.5); + + for i = 1:length(x) + patch('Faces' ,GlobalNodeID(i,:),... + 'Vertices' ,XY_Mesh,... + 'EdgeColor' ,grayColor,... + 'FaceColor' ,'blue', ... + 'LineWidth' ,x(i)); + + patch('Faces' ,GlobalNodeID(i,:),... + 'Vertices' ,XY_Mesh,... + 'EdgeColor' ,'none',... + 'FaceColor' ,'none', ... + 'MarkerEdgeColor','blue',... + 'Marker' ,'o',... + 'MarkerFaceColor','white',... + 'MarkerSize' ,10,... + 'LineWidth' ,3); + end + + + set(gca,'xtick',[],'ytick',[]) + set(gca,'XColor', 'none','YColor','none') + title(MyTitle,'FontSize',18) + pause(1) +end +%------------------------------------------------------------ + +%------------------------------------------------------------ +function [X,Y] = hatch_coordinates(xlim, ylim, xstep, ystep, merge) + %// function [X,Y] = hatch_coordinates(xlim,ylim,xstep,ystep,merge) + %// + %// Return coordinates for plotting a hatch pattern + %// The angle of the lines can be adjusted by varying the + % ratio xstep/ystep + %// xlim and ylim are vectors with two elements, + % where the first element needs to be smaller than the second. + + % // set default options + if nargin < 3 ; xstep = 1 ; end + if nargin < 4 ; ystep = xstep ; end + if nargin < 5 ; merge = true ; end + + % // define base grid + xpos = xlim(1):xstep:xlim(2) ; nx = numel(xpos) ; + ypos = ylim(1):ystep:ylim(2) ; ny = numel(ypos) ; + + % // Create the coordinates + nanline = NaN*ones(1,nx+ny-3) ; + X = [ [ xpos(1)*ones(1,ny-2) xpos(1:end-1) ] ; ... + [ xpos(2:end) xpos(end)*ones(1,ny-2) ] ; ... + nanline ] ; + Y = [ [ypos(end-1:-1:1) ylim(1)*ones(1,nx-2) ] ; ... + [ypos(end)*ones(1,nx-1) ypos(end-1:-1:2)] ; ... + nanline ] ; + + % // merge if asked too + if merge + X = X(:) ; + Y = Y(:) ; + end +end +%------------------------------------------------------------ \ No newline at end of file