水平集拓扑优化经典程序.doc

function TOPLSM(DomainWidth, DomainHight, EleNumPerRow, EleNumPerCol, LV, LCur,FEAInterval, PlotInterval, TotalItNum) %=% TOPLSM, a 199-line Matlab program, is developed and presented here for the% mean compliance optimization of structures in 2D, with the classical level set method. % Developed by: Michael Yu WANG, Shikui CHEN and Qi XIA% First Version : July 10, 2004% Second Version: September 27, 2004% Last Modification:October 31, 2005, optimize the code.% The code can be downloaded from the webpage:% .hk/cmdl/download.htm% Department of Automation and Computer-Aided Engineering, % The Chinese University of Hong Kong% Email: .hk%Main references:% (1.)M.Y. Wang, X. M. Wang, and D. M. Guo,A level set method for structural topology optimization,% Computer Methods in Applied Mechanics and Engineering, 192(1-2), 227-246, January 2003%(2.) M. Y. Wang and X. M. Wang, PDE-driven level sets, shape sensitivity, and curvature flow for structural topology optimization,% CMES: Computer Modeling in Engineering & Sciences, 6(4), 373-395, October 2004.%(3.) G. Allaire, F. Jouve, A.-M. Toader, Structural optimization using sensitivity analysis and a level-set method , % J. Comp. Phys. Vol 194/1, pp.363-393 ,2004. %Parameters:% DomainWidth : the width of the design domain;% DomainHight : the hight of the design domain;% EleNumPerRow : the number of finite elements in horizontal direction;% EleNumPerCol : the number of finite elements in vertical direction;% LV : Lagrange multiplier for volume constraint;% LCur : Lagrange multiplier for perimeter constraint whose shape sensitivity is curvature;% FEAInterval : parameters to specify the frequency of finite element% analysis;% PlotInterval : parameters to specify the frequency of plotting;% TotalItNum : total iteration number.%=% Step 1: Data initializationEW = DomainWidth / EleNumPerRow; % The width of each finite element.EH = DomainHight / EleNumPerCol; % The hight of each finite element.M = EleNumPerCol + 1 , EleNumPerRow + 1 ; % the number of nodes in each dimension x , y = meshgrid( EW * -0.5 : EleNumPerRow + 0.5 , EH * -0.5 : EleNumPerCol + 0.5 ); FENd.x, FENd.y, FirstNodePerCol = MakeNodes(EleNumPerRow,EleNumPerCol,EW,EH); % get the coordinates of the finite element nodesEle.NodesID = MakeElements( EleNumPerRow, EleNumPerCol, FirstNodePerCol ); LSgrid.x = x(:); LSgrid.y = y(:); % The coordinates of each Level Set grid for i = 1 : length(Ele.NodesID(:,1) Ele.LSgridID(i) = find(LSgrid.x - FENd.x(Ele.NodesID(i,1) - EW/2).2 +. % get the ID of the level set grid that lies in the middle of a finite element (LSgrid.y - FENd.y(Ele.NodesID(i,1) - EH/2).2 = 100*eps); end;cx = DomainWidth / 200 * 33.33 100 166.67 0 66.67 133.33 200 33.33 100 166.67 0 66.67 133.33 200 33.33 100 166.67;cy = DomainHight / 100 * 0 0 0 25 25 25 25 50 50 50 75 75 75 75 100 100 100; for i = 1 : length( cx ) tmpPhi( :, i ) = - sqrt ( ( LSgrid . x - cx ( i ) ) .2 + ( LSgrid . y - cy ( i ) ) .2 ) + DomainHight/10; end;LSgrid.Phi = - (max(tmpPhi.).; % define the initial level set functionLSgrid.Phi(LSgrid.x - min(LSgrid.x) .* (LSgrid.x - max(LSgrid.x) .* (LSgrid.y - max(LSgrid.y) .* (LSgrid.y - min(LSgrid.y) = 100*eps) = -1e-6;FENd.Phi = griddata( LSgrid.x, LSgrid.y, LSgrid.Phi, FENd.x, FENd.y, cubic); % project Phi values from the level set function to the finite element nodesF = sparse( 2 * (EleNumPerRow + 1 ) * (EleNumPerCol + 1), 1 ); F( 2 * (EleNumPerRow + 1 ) * (EleNumPerCol + 1) - EleNumPerCol ) = -1; % construct force vectorItNum = 1;while ItNum 0 % the case that the element is inside the boudary ke = BasicKe(E1, NU, EW, EH);elseif max(Phi(EleNodeID(i,:) = 0 )/length(s(:);ke = AreaRatio * BasicKe(E1,NU,EW, EH);end;% THE END OF THIS FUNCTIONfunction u, v , MeanCompliance = FEA(E1, E0, F, EleNumPerRow, EleNumPerCol, EW, EH, FENodes, Ele, NU)%=% function u, v , MeanCompliance = FEA(E1, E0, F, EleNumPerRow, EleNumPerCol, EW, EH, FENodes, Ele, NU)% returns the 2-dimensional displacement field and the mean compliance;% E1 : Youngs modulus of elastic material;% E0 : Youngs modulus of void material;% F : the force vector;% EW: geometric width of a finite element;% EH : geometric hight of a finite element;% FENodes : the structure storing finite element nodes;% Ele : the structure storing finite element;% NU: Poisson ratio, assume that the two materials have the same NU;% %=%K = sparse(2 * (EleNumPerRow + 1 )*(EleNumPerCol + 1), 2 * (EleNumPerRow + 1 )*(EleNumPerCol + 1) ); for i=1:EleNumPerRow * EleNumPerCol ke = EleStiffMatrix(EW, EH, E1,E0, NU, FENodes.Phi, Ele.NodesID, i); K = Assemble(K, ke, Ele.NodesID, i); end;K = AddBoundCondition(K , EleNumPerCol);U = K F;tmp = 2 * (EleNumPerRow + 1 )*(EleNumPerCol + 1) - EleNumPerCol;MeanCompliance = F( tmp ) *U( tmp );for i = 1: 0.5 * length(U)u(i,1) = U(2 * i - 1); v( i ,1) = U(2 * i );end;% THE END OF THIS FUNCTIONfunction Phi1 = LevelSetEvolve(Phi0, Vn, dx, dy, LoopNum)%=% function Phi1 = LevelSetEvolve(Phi0, Vn, dx, dy, LoopNum)% updates the level set surface using a first order space convex.% Phi0: is an m-by-n matrix. Its the level set surface before evolution.% Vn: is an m-by-n matrix.Its the normal velocity field.% Phi1: is an m*n-by-1 vector. Its the level set surface after evolution.% %=%DetT = 0.5 * min(dx,dy)/max(abs(Vn(:);for i = 1 : LoopNumDx_L, Dx_R = UpwindDiff(Phi0 , dx , x);Dy_L, Dy_R = UpwindDiff(Phi0 , dy , y);Grad_Plus = (max(Dx_L,0).2 + (min(Dx_R , 0).2 + (max(Dy_L,0).2 + (min(Dy_R,0).2 ).0.5;Grad_Minus = (min(Dx_L,0).2 + (max(Dx_R , 0).2 + (min(Dy_L,0).2 + (max(Dy_R,0).2 ).0.5;Phi0 = Phi0 - DetT .* (max(Vn, 0) .* Grad_Plus + min(Vn, 0) .* Grad_Minus);end;Phi1 = Phi0(:);% THE END OF THIS FUNCTIONfunction EleNodesID=MakeElements(EleNumPerRow,EleNumPerCol,FirstNodePerCol)%=% function EleNodesID=MakeElements(EleNumPerRow,EleNumPerCol,FirstNodePerCol) % is used to produce finite elements.% EleNodesID(1:NumEle, 1:4): stands for node ID that make up each element;% O-O% |4 3 |% | |% | |% |1 2 |% O-O% EleNumPerRow: number of element per row% EleNumPerCol: number of element per column% FirstNodePerCol: the serial number of the first node in each column% %=%EleNodesID = zeros(EleNumPerRow * EleNumPerCol , 4);for i=1:EleNumPerRow EleNodesID(i * EleNumPerCol : -1:(i-1) * EleNumPerCol + 1 , 4) = FirstNodePerCol(i): -1:FirstNodePerCol(i) - EleNumPerCol + 1.;end;EleNodesID(:,1)=EleNodesID(:,4)- 1;EleNodesID(:,2)=EleNodesID(:,1)+ EleNumPerCol + 1;EleNodesID(:,3)=EleNodesID(:,2)+ 1;% THE END OF THIS FUNCTIONfunction NodesX, NodesY, FirstNodePerCol = MakeNodes(EleNumPerRow,EleNumPerCol,EleWidth,EleHight)%=% function nodesPosition,FirstNodePerCol=MakeNodes(EleNumPerRow,EleNumPerCol,EleWidth,EleHight)% is used to make nodesPosition% NodesX: an n-by-1 vector discribing the x coordinates of FE nodes% NodesY: an n-by-1 vector discribing the x coordinates of FE nodes% EleNumPerRow: the number of elements in each row;% EleNumPerCol: the number of elements in each col% EleWidth: is the width of an element. Here EleWidth=1;% EleHight: is the hight of an element. Here EleHight=1;% %=% x , y = meshgrid( EleWidth * 0 : EleNumPerRow , EleHight * 0 : EleNumPerCol);FirstNodePerCol = find( y(:) = max(y(:);NodesX = x(:); NodesY = y(:);% THE END OF THIS FUNCTIONfunction Matrix = Matrix4diff( Phi )%=% function Matrix = Matrix4diff( Phi ) produces a structure used for% upwind finite diffence.% Phi: is an m-by-n matrix;% Matrix: a structure which includes 4 matrixes used to calculate finite% difference.%=%Matrix.i_minus_1 = zeros(size(Phi);Matrix.i_plus_1 = Matrix.i_minus_1; Matrix.j_minus_1 = Matrix.i_minus_1; Matrix.j_plus_1 = Matrix.i_minus_1;Matrix.i_minus_1(:, 1) = Phi(:, end); Matrix.i_minus_1(:, 2:end) = Phi(:,1:end-1);Matrix.i_plus_1(:, end) = Phi(:,1); Matrix.i_plus_1(:, 1:end-1) = Phi(:,2:end);Matrix.j_minus_1(1, :) = Phi(end, :); Matrix.j_minus_1(2:end , :) = Phi(1:end-1,:);Matrix.j_plus_1(end,:) = Phi(1,:); Matrix.j_plus_1(1:end-1, :) = Phi(2:end, :);% THE END OF THIS FUNCTIONfunction SignDistPhi = Reinitialize(Phi0, dx, dy, LoopNum)%=%function SignDistPhi = Reinitialize(Phi0, dx, dy, LoopNum) is used to% regulize the level set function to be a signed distance function. We % adopt the PDE-based method proposed by Peng, Merriman, and Osher% etc.,where% Phi0: is an m-by-n matrix. Its the level set surface before reinitialization.% dx : the interval between two adjacent grids in axis X.% dy : the interval between two adjacent grids in axis Y.% SignDistPhi : is an m*n-by-1 vector. Its the level set surface after reinitialization.% %=%for i = 1 : LoopNum + 1 Dx_L, Dx_R = UpwindDiff(Phi0 ,