这些代码，请您抽空注解一下或者部分注解一下。分析一下程序的结构
本帖最后由 417332551 于 2011-7-30 14:14 编辑TOPLSM 199-line version
港大的那个网站上的代码。。
您可以重新发个帖子，注解一下 这些代码。然后我把这个帖子处理掉。
%================================================%
% 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
% Department of Automation and Computer-Aided Engineering,
% The Chinese University of Hong Kong
% Email: yuwang@acae.cuhk.edu.hk
% The code can be downloaded from the webpage:
% http://www2.acae.cuhk.edu.hk/~cmdl/download.htm
%================================================%
function TOPLSM(DomainWidth, DomainHight, ENPR, ENPC, LV, LCur,FEAInterval, PlotInterval, TotalItNum)
% Step 1: Data initialization
EW = DomainWidth / ENPR; %The width of each finite element.
EH = DomainHight / ENPC; % The hight of each finite element.
M = [ ENPC + 1 , ENPR + 1 ];
[ x ,y ] = meshgrid( EW * [ -0.5 : ENPR + 0.5 ] , EH * [ -0.5 : ENPC + 0.5 ]);
[ FENd.x, FENd.y, FirstNdPCol ] = MakeNodes(ENPR,ENPC,EW,EH);
Ele.NdID = MakeElements( ENPR, ENPC, FirstNdPCol );
LSgrid.x = x(：）; LSgrid.y = y(：）; % The coordinates of each Level Set grid
for i = 1 : length(Ele.NdID(:,1))
Ele.LSgridID(i) = find((LSgrid.x - FENd.x(Ele.NdID(i,1)) - EW/2).^2 +...
(LSgrid.y - FENd.y(Ele.NdID(i,1)) - EH/2).^2 <= 100*eps);
end;
cx = DomainWidth / 200 * [ 33.33100166.670 66.67133.3320033.33100166.670 66.67133.3320033.33100166.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.')).'; % The level set function value on each Level Set grid
LSgrid.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');
F = sparse( 2 * (ENPR + 1 ) * (ENPC + 1), 1 );
F( 2 * (ENPR + 1 ) * (ENPC + 1) - ENPC ) = -1;
ItNum = 1;
while ItNum <= TotalItNum
%Step 2: Finite element analysis
disp(' '); disp( ['Finite Element Analysis No. ', num2str(ItNum), ' starts ... '] );
FEAStartTime = clock;
[ FENd.u, FENd.v , MeanCompliance] = FEA(1, 1e-3, F, ENPR, ENPC, EW, EH, FENd, Ele, 0.3);
LSgrid.Curv = CalcCurvature(reshape(LSgrid.Phi,M + 1), EW, EH);
%Step 3: Shape sensitivity analysis and normal velocity field calculation
LSgrid.Beta = zeros(size(LSgrid.x));
for e = 1: ENPR * ENPC
LSgrid.Beta(Ele.LSgridID(e)) = SensiAnalysis( 1, 1e-3, Ele.NdID, FENd.u , FENd.v, EW, EH,...
0.3, e, LV, LCur, LSgrid.Phi(Ele.LSgridID(e)), LSgrid.Curv(Ele.LSgridID(e)));
end;
LSgrid.Vn = LSgrid.Beta/ max(abs(LSgrid.Beta));
disp(['FEA costs ',num2str( etime(clock, FEAStartTime)), ' seconds']);
%Step 4: Level set surface update and reinitialization
LSgrid.Phi = LevelSetEvolve( reshape(LSgrid.Phi,M + 1), reshape(LSgrid.Vn,M + 1), EW, EH, FEAInterval);
if (ItNum == 1)||mod(ItNum, 5) ==0
LSgrid.Phi = Reinitialize(reshape(LSgrid.Phi,M + 1), EW, EH, 20);
end;
FENd.Phi = griddata( LSgrid.x, LSgrid.y, LSgrid.Phi, FENd.x, FENd.y, 'cubic');
% Step 5: Results visualization
if (ItNum == 1)||(mod(ItNum,PlotInterval) == 0)
subplot(1,2,1);
contourf( reshape(FENd.x , M), reshape(FENd.y , M), reshape(FENd.Phi , M), );
axis equal;grid on;
subplot(1,2,2); contourf( x, y, reshape(-LSgrid.Phi , M + 1), );alpha(0.05);hold on;
h3=surface(x, y, reshape(-LSgrid.Phi , M + 1));view();axis equal;grid on;
set(h3,'FaceLighting','phong','FaceColor','interp', 'AmbientStrength',0.6); light('Position',,'Style','infinite');
pause( 0.5 );
end;
%Step 6 : Go to step 2.
ItNum = ItNum + 1;
end;
function =AddBoundCondition(K , EleNumPerCol)
for i = 1 : EleNumPerCol + 1
K( 2 * i-1 : 2 * i,：） = 0;
K( : , 2 * i-1: 2 * i) = 0;
K(2 * i - 1, 2 * i - 1) = 1;
K(2 * i , 2 * i) = 1;
end;
function =Assemble(K,ke,elementsNodeID,EleID)
m = elementsNodeID(EleID,：）;
for i = 1 : length(m)
for j = 1 : length(m)
K(2*m(i)-1 : 2*m(i), 2*m(j)-1: 2*m(j)) = K(2 * m(i) - 1 : 2 * m(i) , 2 * m(j) - 1: 2*m(j) )+ ke( 2*i-1: 2*i , 2*j-1:2*j);
end;
end;
function = BasicKe(E,nu, a, b)
k = [-1/6/a/b*(nu*a^2-2*b^2-a^2), 1/8*nu+1/8, -1/12/a/b*(nu*a^2+4*b^2-a^2), 3/8*nu-1/8, ...
1/12/a/b*(nu*a^2-2*b^2-a^2), -1/8*nu-1/8, 1/6/a/b*(nu*a^2+b^2-a^2), -3/8*nu+1/8];
KE = E/(1-nu^2) * [ k(1) k(2) k(3) k(4) k(5) k(6) k(7) k(8)
k(2) k(1) k(8) k(7) k(6) k(5) k(4) k(3)
k(3) k(8) k(1) k(6) k(7) k(4) k(5) k(2)
k(4) k(7) k(6) k(1) k(8) k(3) k(2) k(5)
k(5) k(6) k(7) k(8) k(1) k(2) k(3) k(4)
k(6) k(5) k(4) k(3) k(2) k(1) k(8) k(7)
k(7) k(4) k(5) k(2) k(3) k(8) k(1) k(6)
k(8) k(3) k(2) k(5) k(4) k(7) k(6) k(1)];
function [ Curvature ] = CalcCurvature(Phi, dx, dy)
Matrices = Matrices4Diff( Phi );
Phix = (Matrices.i_plus_1 - Matrices.i_minus_1)/(2 * dx);
Phiy = (Matrices.j_plus_1 - Matrices.j_minus_1)/(2 * dy);
= UpwindDiff(Phiy , dx, 'x');
Phixy = (Phiyx_Bk + Phiyx_Fw)/2;
Phixx =(Matrices.i_plus_1 - 2*Phi + Matrices.i_minus_1)/dx^2;
Phiyy =(Matrices.j_plus_1 - 2*Phi + Matrices.j_minus_1)/dy^2;
Curvature = (Phixx .* Phiy.^2 - 2 * Phix .* Phiy .* Phixy + Phiyy .* Phix.^2) ./ ((Phix.^2 + Phiy.^2).^1.5 + 100*eps);
Curvature = Curvature(：）;
function = EleStiffMatrices(EW, EH,E1,E0,NU,Phi,EleNds,i)
if min(Phi(EleNds(i,：）)) > 0
ke = BasicKe(E1, NU, EW, EH);
elseif max(Phi(EleNds(i,：）)) < 0
ke = BasicKe(E0, NU, EW, EH);
else
[ s , t ] = meshgrid([-1 : 0.1 : 1],[-1 : 0.1 : 1]);
tmpPhi = (1 - s(：）).*(1 - t(：）)/4 * Phi(EleNds(i,1)) + (1 + s(：）).*(1 - t())/4 * Phi(EleNds(i,2))...
+ (1 + s(：）).*(1 + t(：）)/4 * Phi(EleNds(i,3)) + (1-s(：）).*(1 + t(：）)/4 * Phi(EleNds(i,4));
ke = length(find( tmpPhi >= 0 ))/length(s(：）) * BasicKe(E1,NU,EW, EH);
end;
function [ u, v , MeanCompliance] = FEA(E1, E0, F, ENPR, ENPC, EW, EH, FENds, Ele, NU)
K = sparse(2 * (ENPR + 1 )*(ENPC + 1), 2 * (ENPR + 1 )*(ENPC + 1) );
for i=1:ENPR * ENPC
ke = EleStiffMatrices(EW, EH, E1,E0, NU, FENds.Phi, Ele.NdID, i);
K = Assemble(K, ke, Ele.NdID, i);
end;
K = AddBoundCondition(K , ENPC);
U = K \ F;
MeanCompliance = F( 2 * (ENPR + 1 )*(ENPC + 1) - ENPC ) *U( 2 * (ENPR + 1 )*(ENPC + 1) - ENPC );
for i = 1: 0.5 * length(U)
u(i,1) = U(2 * i - 1);
v( i ,1) = U(2 * i );
end;
function = LevelSetEvolve(Phi0, Vn, dx, dy, LoopNum)
DetT = 0.3 * min(dx,dy)/max(abs(Vn(：）));
for i = 1 : LoopNum
= UpwindDiff(Phi0 , dx , 'x');
= 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(：）;
function =MakeElements(EleNumPerRow,EleNumPerCol,FirstNdPCol)
EleNds = zeros(EleNumPerRow * EleNumPerCol , 4);
for i=1:EleNumPerRow
EleNds( , 4) = .';
end;
EleNds(:,1)=EleNds(:,4)- 1;
EleNds(:,2)=EleNds(:,1)+ EleNumPerCol + 1;
EleNds(:,3)=EleNds(:,2)+ 1;
function = MakeNodes(EleNumPerRow,EleNumPerCol,EleWidth,EleHight)
[ x , y ]= meshgrid( EleWidth * [ 0 : EleNumPerRow ], EleHight * );
FirstNdPCol = find( y(：） == max(y(：）));
NodesX = x(：）; NodesY = y(：）;
function = Matrices4Diff( Phi )
Matrices.i_minus_1 = zeros(size(Phi));
Matrices.i_plus_1 = Matrices.i_minus_1;
Matrices.j_minus_1 = Matrices.i_minus_1;
Matrices.j_plus_1 = Matrices.i_minus_1;
Matrices.i_minus_1(:, 1) = Phi(:, end);
Matrices.i_minus_1(:, 2:end) = Phi(:,1:end-1);
Matrices.i_plus_1(:, end) = Phi(:,1);
Matrices.i_plus_1(:, 1:end-1) = Phi(:,2:end);
Matrices.j_minus_1(1, ) = Phi(end, ：）;
Matrices.j_minus_1(2:end , ：） = Phi(1:end-1,：）;
Matrices.j_plus_1(end,：） = Phi(1,）;
Matrices.j_plus_1(1:end-1, ：） = Phi(2:end, ：）;
function = Reinitialize(Phi0, dx, dy, LoopNum)
for i = 1 : LoopNum + 1
= UpwindDiff(Phi0 , dx , 'x');
= UpwindDiff(Phi0 , dy , 'y');
S = Phi0 ./ (sqrt(Phi0.^2 + (((Dx_L + Dx_R)/2).^2 + ((Dy_L + Dy_R)/2).^2) * dx^2) + eps);
DetT = 0.5 * min(dx,dy)/max(abs(S(：）));
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(S, 0) .* Grad_Plus + min(S, 0) .* Grad_Minus) - S);
end;
SignDistPhi = Phi0(：）;
function [ Beta ] = SensiAnalysis( E1, E0, EleNds, u, v, EW , EH,NU, EleID , L4Vol , L4Curv , Phi , curvature)
ae = [ u(EleNds(EleID,1)); v(EleNds(EleID,1)); u(EleNds(EleID,2)); v(EleNds(EleID,2));
u(EleNds(EleID,3)); v(EleNds(EleID,3)); u(EleNds(EleID,4)); v(EleNds(EleID,4))];
B = 1/2*[ -1/EW, 0, 1/EW, 0, 1/EW, 0, -1/EW, 0;
0, -1/EH, 0, -1/EH, 0, 1/EH, 0, 1/EH;
-1/EH, -1/EW, -1/EH, 1/EW, 1/EH, 1/EW, 1/EH, -1/EW ];
if Phi> 0.75 * EW
E = E1;
elseif Phi< -0.75 * EW
E = E0;
else
E = E1 * 0.75 * ( 1.0 - 1e-3 ) * ( Phi / (0.75 * EW) - (Phi / (0.75 * EW))^3 / 3.0 ) + 0.5 * ( 1 + 1e-3 );
end;
D = E / (1 - NU^2) * [ 1, NU , 0; NU , 1 , 0; 0 , 0, (1 - NU)/2 ];
Beta = L4Vol - (B * ae).' * D * (B * ae) - L4Curv * curvature;
function = UpwindDiff(Phi , dx , strDirection)
Matrices = Matrices4Diff( Phi );
if strDirection == 'x'
BackDiff = (Phi - Matrices.i_minus_1)/dx;
FawdDiff = (Matrices.i_plus_1 - Phi)/dx;
elseif strDirection == 'y'
BackDiff = (Phi - Matrices.j_minus_1)/dx;
FawdDiff = (Matrices.j_plus_1 - Phi)/dx;
end;
%=================================================================%
% The program, TOP_LevelSet, is developed for the implimentation oftopology optimization
% of structures with Level Set Methods.
% First Version : July 10, 2004
% Second Version: September 27, 2004
% Last Modification:October 31, 2005, optimize the code.
%
% Developed by: Michael Yu WANG, Shikui CHEN and Qi XIA
% Department of Automation and Computer-Aided Engineering,
% The Chinese University of Hong Kong
% Email: yuwang@acae.cuhk.edu.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.
% ************************************Disclaimer********************************************%
% The authors reserve all rights for the programs. The programs may be distributed and
% used for academic and educational purposes. The authors do not guarantee that the
% code is free from errors, and they shall not be liable in any event caused by the use
% of the programs.
%=================================================================% 等我回校看有没有时间，这个程序其实是在99行经典程序的基础上写的，大致的结构差不多啊 回复 2# yygao56
我回校了。
貌似这个程序有点长啊~~~
程序中也有了部分英文注释
你看你具体那一部分看不懂，我解释一下~~ 我问个简单的问题，用OptiStruct做优化，设计变量可以是原始模型的任意设计参数嘛，比如某条边的长度，孔的直径大小等。。。。:) 回复 3# yygao56
版主终于回学校了，想死我们了:lol
页:
[1]