matlab实现水平集方法

【实例简介】
学习水平集方法的基础知识，用matlab实现水平集，与matlab99行拓扑优化类似。入门级程序代码
Discrete level-set topology optimization code written in matlab 455 the position of the center of element e, then the dis- Using these shape sensitivities, the normal velocity u relized level-set function y satisfies within element e at iteration k of the algorithm is o if x = l L (ce)1z0 If xe=0 ue keue-d I 89 (V(x)-Vr).(10) The discrete level-set function can then be updated to Based on Burger et al.(2004), the forcing term g find a new structure by solving (3)numerically. y is should be taken as initialized as a signed distance function and an upwind finite difference scheme is used so that the evolution 8=-sign (TL equation can be accurately solved. In addition, the time step for the finite difference scheme is chosen to satisfy where STL is the topological sensitivity of the the cfl stability condition Lagrangian L. For compliance minimization nucleat ing solid areas within the void regions of the design is pointless because such solid regions will not take any △t Inax v load. Therefore holes should only be nucleated within the solid structure and where h is the minimum distance between adjacent mum is taken over all gridpoints (e.g, Sethian 1999 8=oTL if v<0 gridpoints in the spatial discretization and the maxi 0ify≥0. (12) Osher and上 eakin2002) As stated above, the two scalar fields u and g are The topological sensitivity of the compliance objective typically chosen based on the shape and topological function in two dimensions with traction-free boundar sensitivities of the optimization objective, respectively. conditions on the nucleated hole and the unit ball as the In order to satisfy the volume constraint. here they are model hole is(e.g. Allaire et al. 2005) chosen using the shape and topological sensitivities of the e lagrangian T(+2 2(λ+) (4uu kee+(-u)u? (kTr)eue) L=c(x)+×(V(x)-V)+ k V(x)-Vr」, (13) (6) where up(kTr)eue is the finite element approximation to the product tr(o tr(e) where o is the stress tensor where nk and A are parameters which change with and e is the strain tensor. In (13), n and u are the each iteration k of the optimization algorithm. They are Lame constants for the solid material. The topological updated using the scheme sensitivity of the volume v(x) when the model hole is he unit ball k+1 (V(x)-Vreg) aA,(7) (14) where C E(0, 1)is a fixed parameter. This implements the augmented Lagrangian multiplier method for con- Putting these results into the definition of the strained optimization, as used for topology optimiza- Lagrangian and using(12)gives the source term g lion with the level-set method by Luo et al.(2008a) Those readers interested in more discussion of the The normal velocity v is chosen as a descent direc- algorithm are referred to the publications referenced tion for the Lagrangian L using its shape derivative above. Also, many more details are revealed by the (e.g., Allaire et al. 2004; Wang et al. 2003). In the case Matlab implementation and are discussed in the follow- of traction-tree boundary conditions on the moving ing section boundary the shape sensitivity of the compliance ob jective c(x)is the negative of the strain energy density (e.g, Allaire et al. 2004) 3 Explanation of Matlab implementation Sc le=-u,kele. 8Q (8) Following Sigmund(2001), several assumptions are made to simplify the matlab code. In particular, the The shape sensitivity of the volume v(x)is optimization domain is assumed to be rectangular and split into square finite clements. There arc nelx (9) elements along the horizontal direction and nely 456 V.J. Challis entirely solid (line 4). The level-set function lsf is ini lialized from the structure using the reinit function (line 5)which is described in detail in Section 3.2. The shape and topological sensitivity arrays shapesens and topsens are initialized as zero arrays(line 6) Lastly, the stiffness matrix and other material para mctcrs arc sct with a call to materialInfo (linc 7) Fig. I Left: load and supports for the simple bridge prob- before iteration of the algorithm begins 3,2,2) Iterations are counted with iterNum and continue for a maximum of 200 iterations(line 9). Inside the iteration loop, the finite element analysis routine FE elements along the vertical direction. The total num- is called (line 11)to allow calculation of the shape ber of elements is N=nelx x nely nd topological sensitivities for the current structure The code is envoked with the call (lines 12-21). The node numbering and method of finding Ue(=ue)in lines 12-16 are drawn directly from top levelset(nelx, nely, volReq, Sigmund(2001). Line 17 calculates the shape sensitivity stepLength, numReinit, topWeight)i of the compliance from( 8).Lines 18-19 calculate the where volReg is the required solid fraction(equivalent Lopological sensitivity of the compliance from(13).The to N), and the remaining three inputs are parameters compliance objective value is calculated from the shape of the algorithm. The parameter stepLength specifies sensitivity and stored in the vector objective (line the time interval over which the level-set function is 23), and the current solid volume fraction is calculatcd evolved in a design update. measured relative to the as volCurr (line 24). The objective and volume in CFL time step. The parameter numReinit determines formation for the current structure is displayed to the how often the level-set function is reinitialized to a scrccn and the current structurc is displayed as a figurc signed distance function: level-set reinitialization is per- (lines 25-27), similar to Sigmund(2001) formed at the end of each numreinit-th iteration A convergence check (lines 28-32) may terminate of the algorithm. The last parameter, topweight, is the algorithm before 200 iterations are completed. The the weighting o of the forcing term in the evolution convergence check is not performed for the first five it- equation erations of the algorithm(line 29). After these first five The boundary conditions and loads are specified iterations the optimization terminates if the volume is within the Matlab code. The version of the code given within 0.005 of the required value volreg (ine 29)and here performs the optimization for a simple bridge. a also the previous five objective function values are all suggested first call of the program is within a 1% tolerance of the current objective value tcp_-eve1set(50,30,0.3,3,2,2); If the algorithm does not terminate, a design update The boundary conditions and optimization result for is performed. Prior to the design update, the parame this call to the code are given in Fig. 1 ters for imposing the volume constraint with the aug Details of the code are explained in the following mented Lagrangian method are updated (lines 33-38) subsections. Wherever possible the code is similar to The matlab variable la stores the current value of k that of Sigmund(2001) La stores the current value of A, and alpha stores the parameter a. The sensitivities of the additional terms in 3.1 Main program--lines 1-48 the lagrangian in(6) are added to the sensitivities of the objective function (lines 39-41)and then the design The main program includes initialization. an iteration update is performed with a call to updateStep(line loop to perform the optimization and a check for 43). The reinitialization of the level-set function is then convergence carried out, but only for iteration numbers that are a The initialization phase defines the array struc of multiple of the user-specified parameter numReinit size nely x nelx that stores the current values xe (line 44-47). The optimization loop then repeats defining the structure. The initial structure is set to be 2The factor max(struc(ely, elx),00001)gives a very Names in typewriter font correspond to Matlab variables or small stiffness to the void phase which improves the robustness functions of the algorithm Discrete level-set topology optimization code written in matlab 457 3. 2 Reinitialization of the leyel-set function-lines ment function for those readers without the toolbox 49-53 is given in Appendix b 2. Set the sensitivities to zero for those elements The reinitialization of the level-set function to a signed which must not change phase during the evolution distance function is important to ensure accuracy in In particular, any load-bearing elements must re solving the evolution equation (e.g., Sethian 1999; main solid (lines 59-61). If desired, passive ele- Osher and Fedkiw 2002). However, to allow nuclcation ments can also be specificd to remain void b of holes under evolution of (3), the level-set function setting their sensitivities to zero should not be reinitialized too often(Burger et al. 3. Evolve the level-set function to find the new struc 2004). This is controlled with the user-specified pa ture using the function evolve(lines 62-63). The rameter numreinit. When required, reinitialization negative of the shape sensitivity shapesens is is achieved in the Matlab implementation using the passed into the evolve function as the normal ve function reinit, which accepts the current structure locity u for the boundary of the design The expres struc and returns the level-set function lsf(=v) for sion topsens *(lsf(2: end-1, 2: end-1)<0) that structure (line 50) in line 63 calculates the forcing term g from The level-set array Isf has size (nely +2)x the topological sensitivity of the Lagrangian (nelx+2)-this includes a border of void elements topSens as in(12). The user-specified parameter around the structure. the border is necessary to ensure topWeight is passed into evolve as the parame the boundary of the structure is accurately represented ter The evolve function is described nexu by the level-set function. struc is therefore enlarged to strucfull which has a border of void clements 3. 4 evolution of the lcvel-sct function -lincs 64-86 (line 51) prior to calculating lsf For each element, the level-set function is calculated The function evolve solves the evolution equation 3) as the minimum euclidean distance to an clement of to give a ncw structurc. The function takcs in the cur he opposite phase, as measured from the center of the rent level-set function lsf(=4), the normal velocity element. For solid elements the sign is chosen to be (v), the forcing term g(=8), the user-specified negative, whereas for void elements the sign is chosen parameter stepLength, and the parameter w(= o) to be positive, matching with the definition of the level- This function uses an upwind finite difference set function(4) scheme(e.g, Sethian 1999, page 65) to solve the evolu Matlab's bwdist function (part of the Image tion equation(3). The time-step is chosen to be a small Processing Toolbox) makes the reinitialization ex- portion(. 1) of the value given by the CFl condition tremely simple. Given an array, bwdist returns an in (5)(lines 69-70), and the number of time-steps is Tray of the same size where each entry is the distance chosen based on the user specified value stepLength to the nearest non-zero entry of the array. With care- (lines 71-72) ful use of signs and logical expressions, two calls to The finite difference scheme is fairly simple to imple bwdist generates the level-set lunction(line 53) ment using Matlabs circshift function. The velocity For readers without access to Matlab's Image v and source term g are extended with a border of zeros Processing Toolbox, the author has written a simple on lines 66-68 to match the level-set function which version of bwdist The code is given in Appendix b includes a border of void elements around the edge of the domain. The new structure is found from the sign of 3.3 Design update-lines 54-63 the updated level-set function(lines 85-86) The function updatestep implements a design update 3.5 Finite element implementation--lines 87-106 using the shape and topological sensitivity information The implementation involves the following steps The function Fe which implements the finite element method is almost identical to the implmentation by 1. Smooth the sensitivities (lines 56-58) to avoid Sigmund(2001), to which the reader is referred for very sharp changes between neighbouring elements more information (such as across the boundary of the structure). The only differences occur in the application of the Including this heuristic step is not necessary but loads and supports, which are specific to the bridge tends to give better results. The Matlab function problem (lines 100-102), and in the assembly of the padarray which is uscd in this smoothing step is global stiffness matrix bccausc there is no nced to pc- part of the Image Processing Toolbox. A replace- nalize intermediate densities. The factor of max(struc 458 V.J. Challis (elx, ely),00001)in the assembly of the globa stiffness matrix(line 97)is to ensure it is non-Singular, allowing solution of the system with Matlabs matrix divide operato 3.6 Element stiffness matrix-lines 107-129 The materialInfo function returns the element stiff- Fig. 2 Left: the load and supports for the half-wheel ness matrix KE (ke), the element matrix for calcu- problem, and right: the result of top levelset lating the trace of the stress tensor multiplied by the (601, 3, 2, 4)i with the code altered to optimize this half-wheel trace of the strain tensor kTr(=(kTr)e), and the lame constants lambda(=A)and mu(= u(line 108). The calculation of the element stiffness matrix is the same Note that the load for the half-wheel is the same as the s the implementation given by Sigmund(2001), but is bridge example so does not need to be changed rearranged to allow re-use of the code for calculating With these new lines, the half-wheel optimization the matrix ktr can be called with(for example) op1 evelset(60,30,0.3,3,2,4) 4 Simple extensions The final structure is displayed in Fig. 2 Next, we consider a simple cantilever probl It is easy to alter the code to consider different bound- picted in Fig 3, matching Sigmund (2001). For this ex ary conditions and loads, multiple load cases, keeping ample, the code must be changed in two places because certain areas void or solid, and different initial struc- both the boundary conditions and the placement of tures. Examples of the first two extensions are given the load have changed. The load bearing elements for h the canitlever must be set to remain solid during the optimization. This is done by setting the sensitivities 4.1 Other boundary conditions shapeSens and topsens of the elements to zero within the updatestep function Specifically, lines 59 It is straightforward to change the code to solve dif- 6l are replaced with ferent optimization problems. Two examples are pre 59 Load bearing elements must remain sented here solid Cantilever First, we change the boundary condition at the lower shape Sens(end, end)=0; right-hand corner of the bridge example of Fig. 1 to topSens(end, end)=0 consider the half-wheel design problem shown in Fig. 2 To restrict only the vertical motion of the gridpo int at In addition, the applied force and supports must be set the lower right-hand corner of the design domain line within the finite element implementation FE, replacing 102 in the code is replaced with: lines 100-102 102 fixeddofs=[(2*(ne1y+1)-1):(2★(ne1y+1)) 100 &Define loads and supports - Cantilever 2★(ne1x+1)★(ne1y+1)]; F((neLx+1)*(nely+1)*2,1)=1 fixeddofs =1: 2*(nely+1) The comment in line 100 should also be changed to reflect the new boundary conditions and be replaced with 100%Define loads and supports -Half-wheel: It is possible to ascribe zero stiff ness to the void phase, and then solve the equation KU=F via a conjugate gradient algorithm As noted earlier, a small nonzero stiffness for the void phase improves the robustness of the algorithm Readers should be careful if copying and pasting code snippets rom this article into Matlab. The author has found that spaces Fig. 3 Left: the load and supports for the cantilever prob are sometimes added during this process to give incorrect. Matlab lem, and right: the result of top levelset (32, 20,0.4, code 2,3, 2)i with the code altered to optimize this cantilever Discrete level-set topology optimization code written in matlab 459 With these new lines, the cantilever optimization can The two additional loads must be specined arter line be called with(for example) 101as top leve1sct(32,20,0.4,2,3;2); 10laF(2大( round(ne1x/2)+1)*(ne1y+1)-1,2) The final structure is displayed in Fig. 3 F(2*( round(ne1x/2)+1)*(nely+1)-1,3) C.5 4.2 Multiple load cases Figure 4 is the result of this multi-load version of the Allowing for multiple load cases is also not difficult, and code called with the following the code can be extended for this purpose in a very similar fashion to the Matlab SIMP implementation top_1 verset(60,30,0.3,4,2,3) Sigmund 2001). Here horizontal loads are added to the bridge problem to improve the stability of the bridge The bridge optimized with the three loads applied sep The new loads are depicted in Fig. 4. The new objective arately is obviously more stable under horizontal load is the sum of the compliance from the three load cases than the one optimized with the three loads applied The following line is added after line ll to set together as in Fig.1 the shape and topological sensitivities to zero before adding in the sensitivites from each load case 5 Discussion 1la shapeSens(: )=0; topSens(: )=0; To accumulate shapesens and topsens from the This discussion addresscs the various paramcters of the three load cases, lines 16-19 are replaced with the Matlab implementation of the level-set algorithm following loop also describes some additional extensions that might bc uscful for tackling morc complicated topology op 16oYi=1:3 timization problems Ue=U(L2*n1-1;2*n1;2*n2-1;2*n2;2*n2 +1;2*n2+2;2*n1+1;2*n1+2],i a table of suggested parameter values for the shape Sens(ely, elx)= shapeSens(ely, elx six parameters which are specified in the call to max(struc(ely, elx),00001)*Ue * KE* Ue; top_levelset are given in Table 1. The upper and topSens(ely, elx)= topSens(ely, elx)+ lower bounds for the parameters have been determined struc(ely, elx )*pi/2* by experimenting with the Matlab code on various (lambda+2*mu)/mu/(lambda+mu)*(4*mu*Ue simple compliance minimization problems. Choosing KE*Ue+(lambda-mu)*Ue*KTr*Ue different parameters within the suggested ranges can result in different topology optimization results. How- Finally, changes are required in the finite element sec- ever, using values within the suggested ranges tend tion of the code to account for now having three load to give sensible optimization results that have objec cases. line 91 becomes Live values close to optimal objective value. Excep- 91F= sparse(2*(ne1y+1)*(ne1x+1),3);= tions can occur, and in particular the two parameters zeros(2*(nely+1)*(nelx+1),3) numReinit and topWeight both affect the ability of the algorithm to nucleate holes in the design so need to be chosen carefully. Specifically, choosing both parameters at the lower end of the suggested ranges (numReinit=2 and topweight= 1)can mean that new holes cannot be nucleated during the optimization and give poor optimized designs Other parameters are hard-coded into the Matlab code. It is not necessary to complicate the program call y including them because the hard-coded values have been suitable for all of the optimization problems tested by the author. The meanings of these parameters along ig. 4 Left: the loads and supports for the multiple load bridge problem-the three loads are applied separately as indicated by with any relevant comments are given in Table 2 the numbers 1, 2 and 3, and right: the result of top levelset As well as altering parameters to see how this affects (6030,0.3,42,)i, with the code changed to consider the outcome of the optimization, readers may wish to this multi-load bridge problem implement new features into the Matlab code. Based 460 V.J. Challis 2 on the author's experience, there are three extensions which particularly come to mind E月9 Improving robustness of the algorithm by checking the change in the objective at each iteration. E巨§ Extending the code in order to impose multiple constraints Utilizing symmetry to speed up the finite element 5三g月 当82分目器罗 calculations during the optimization process a brief description of how to implement these features follows Checking the change in the objective prior to ac cepting an iteration of the optimization algorithm is straight-forward to implement but damages the simple 3月95≌E ①0 teration loop of the presented Matlab code so has not been presented here by the author. If the objective 59四 increases too much after the design has been updated with a call to updatestep, the new structure can be re jected and the author would suggest an adaptive correc 8 tion of the parameter stepLength prior to a new call to updateStep. By reducing stepLength to a frac- tion of its previous value, the new call to updatestep 号 will result in less change in the design. Note that it is ncccssary for some incrcasc in the objcctive function to be acceptable, therefore the tolerance on this checking process should not be set too tightly. The reason for this is two-fold: some ability to go up-hill can help escap local minima, and, when an entirely solid initial design is used the optimization will inevitably go up-hill To tackle more complicated optimization problems 8 it is important to be able to deal with more than one 反司 constraint in the optimization problem. The augmented Lagrangian technique applied in the matlab code to im plement the volume constraint can be used with more caint. For example, the technique ha been used with two constraints in conjunction with the d of topology opt tion by luo et al 目 (2008b). The author and colleagues have developed an alternative approach, which was successful for imposing 523 an isotropy constraint when designing materials in both 3535>5 two and three dimensions The reader is referred to Challis et al. (2008a) 自 9日 In the matlab implementation of the simp method (Sigmund 2001), it is easy to utilize symmetry to reduce the size of the design domain and improve the speed of 5E目 the finite element calculations. This only requires im 2 plementation of the correct boundary conditions along the line of symmetry. To utilize symmetry with in the level-set approach, both the reinitialization of the level- set function and the finite difference evolution scheme must be modified. This is necessary so that the level sct function accurately reflects the symmetry of the optimization problem. The code changes can be made Discrete level-set topology optimization code written in matlab 461 Table 2 Description of parameters hard-coded into the matlab code Parameter Meaning and relevant comments ones (nely, nelx)on line4 Initial structure from which to begin the optimization 200 on line 9 Maximum number of iterations of the optimization algorithm 0.0001 on line 17 Stiff ness of the void phase compared to the solid phase. This can be zero but a small nonzero value improves robustness of the algorithm 5 on line 29 Number of iterations at the start of the optimization for which the convergence criteria are not checked. This is included so that the condition on line 30 is only checked when it is valid otherwise the array index end-5 does not make sense) 0.005 on line 29 Tolerance for satisfaction of the volume constraint s on line 30 Number of iterations which must have similar objective values for termination of the algorithm. If hanged, the user must also change the 5 on line 29 0. 01 on line 30 Relative tolerance on the objective values for termination of the algorithm -0.01, 1000 and 0.9 on line 35 Initial parameters for imposing the volume constraint with the augmented Lagrangian method. If the initial value of la is too small or the value of alpha is too small then the optimization can become unstable array on lines 57 and 58 Convolution filter to smooth the sensitivities, currently chosen as 1/6*[0 10;1 2 1;0 1 0 0.1 on line 70 Fraction of the cfl time step to use as a time step for solving the evolution equation. This must b less than or equal to l to satisfy the CFL condition. Also, to ensure the correct interpretation of the parameter steplength, this value must be the reciprocal of the parameter 10 on line 72. If that value is changed, this must be changed to match 10 on line 72 Number of time steps required to evolve the level-set function for a time equal to the cFl time step To ensure the correct interpretation of the parameter stepLength, this value must be the reciprocal of the parameter (. 1 on line 70. If this value is changed, that must be changed to match 0.0001 on line 97 the global stiffness matrix is not singular, and should match the parameter 0.0001 on line 1? So that Stiff ness of the void phase compared to the solid phase. This must be a positive nonzero value 1 on line 110 The Youngs modulus of the solid phase 0.3 on line 110 The Poissons ratio of the solid phase with care on a case-by-case basis. These changes are anonymous reviewers for their useful comments that resulted in also required when designing unit cells for material significant improvements to the manuscript. The author has re- microstructures, so that periodicity of the design is ceived an Australian Postgraduate Award and funding from the Australian Research Council (grant DP0878785). This financial taken into account support is gratefully acknowledged 6 Conclusion Appendix a--Matlab code This paper presents a simple Matlab implementation 198 TOPOLOGY OPTIMIZATION USING THE LEVEL of the level-set method for topology optimization. The SET METHOD, VIVIEN U. CHALLIS 2009 code can be used to demonstrate the capabilities of the 2 function [struc] top levelset(nelx algorithm as well as specific implementation details nely, volReq, stepLength, numReinit The code could easily be used alongside the simP tp啊 eight) Matlab code written by Sigmund(2001)for simple com 3号工门itia1 ization parisons of the two methods in an educational setting 4 struc ones(nely, nelx); Despite having more lines than the simp implementa- 5 [lsf] reinit(struc) tion the level-set Matlab code is still very compact, and 6 shapeSens zeros(nely, nelx)i topsens zeros(nely, nelx); neglecting comments and the finite element code is only 7 [KE, KTr, lambda, mu]= materialInfo()i 63 lines ions 吕Main1oop 9 for iterNum = 1: 200 10 FE-analysis, calculate sensitivities Acknowledgements The author would like to acknowledge 11 [U= FE(struc, KE Prof. M.P. Bendsge, Prof. J.K. Guest, Dr. A. P. Roberts, Prof 2 for ely = 1: nely O. Sigmund and Dr. A.H. Wilkins for their helpful discussions 13 for elx =1: nelx and support of this work. The author would also like to thank the n1=(ne1y+1)*(e_x-1)+e1y; 462 V.J. Challis n2=(ne1y+1)*e1x+e1y; 55 function [struc, lsf] updateStep(lsf Ue-U([2*n1-1;2*n1;2*n2-1;2*n2;2*n2 shapeSens, topSens, stepLength, topWeight) +1;2*n2+2;2*n1+1;2*n1+2],1) 56 Smooth the sensitivities shapeSens(ely, elx)=-max(struc(ely, 57 [shapeSens] conv2(padarray( shapeSens elx),0.0001)*Ue′★KE为Ue ,[1,1], replicate),1/6*[010;121; topsens(ely, elx)= struc(ely, elx)* 010],′va1ia′) pi/2*(lambda+2*mu)/mu/(lambda+mu)* 58 [topSens] conv2(padarray (copsens, [ 1] plicate ) 1/6 21;010] (4*m1*Ue′*KF*UJe+(1 ambda-m)*Ue’*KTr 1id′) Ue 59 9 Load bearing pixels must remain solid Bridge 21 end 60 shapeSens(end, [1, round(end/ 2): round(end 22 g store data, print plot information 2+1),end])=0 23 objective(iterNum) sum(shapeSens(: ))i 61 topSens(end, [1, round(end/ 2): round(end 24 volCurr sum(struc(: ))/(nelx*nely) /2+1),end])=0 25 disp([! It: num2str (iterNum 628 Design update via evolution sprintf(号10.4f′, objective 63 [struc, lsf]=evolve(-shapeSens, topSens iterUm)) *(1sf(2:end-1,2:end-1)<0},1sE sprintf(′36.3f’, volCurr)]) stepLength, topweight) 27 colormap (gray); images(-struc,[-1,0]); 64号岩---- EVOLUTION OF LEVEL- SET FUNCTION axis equal; axis tight; axis off; drawnow 65 function [struc, lsf] evolve(v,g, lst 28 Check tor convergence steplength, w) 29 if iterUm 5&&( abs(volCurr-volReq 66 Extend sensitivites using a zero border <0.005)SS 67 VFull= zeros(size(v)+2)i vFull(2: end all( abs(objective(end)-objective( 1,2:end-1) end-5: end-1))<0.01*abs(objective 68 gFull zeros(size(g)+2)i gFull(2: end end))) 1,2:end-1)=g return 698 Choose time step for evolution based on 32 end CFL value 8 Set augmented Lagrangian parameters 70 dt =0.1/max(abs(v(: 34 if iterUm 1 1o Evolve for total time stepLength CFL 1a=-0.01;La=1000;a1pha=0.9 value 36 else for i=1:(10*stepLength 1/La *(volCurr volReq)i La 73 g Calculate derivatives on the grid a⊥pha*La 74 dpx circshift(is,[0,-1])-lsf Isf circshift(lsf,[o 1 2 Include volume sensitivities 76 dpy circshift(lsf, [-1,0]) 40 shapeSens- shapeSens -la +1/La*( 7 dmy lsf- circshift(lsf,[1,0])i 78号 Update leve1 set func七 ion using 41 topSens topSens pi*(la-1/La*( upwind scheme volCurr-volReq)) 79 lsf=lsf- dt min(yFull,0 Design update 80 sqrt( min(dmx, 0).2+max(dpx,0).2+min( [struc, lsf] pdateStep(Isf, shapesens dmy,0).^2+max(dpy,0).2) topSens, stepLength, topWeight)i dt max(fUll,0).* 44 Reinitialize level-set function grt( max(drx, 0).2+min(dpx, 0).2+max( 45 if -mod(iterNum, numReinit) dmy, 0).2+min(dpy, 0).2) TIsf]= reinit(struc) W★dt*gFu11 47 end 48 end 85 New structure obtained from lsf 49- REINITIALIZATION OF LEVEL-SET 86 strucFull (lsf<o)i struc strucFull(2: FUNCTION -- end-1,2:end-1); 50 function [lsf]- reinit(struc) 87号吕 FINTTE ELEMENT ANALYSIS 51 strucFull= zeros(size(struc)+2) 88 function [U]=FE(struc, KE strucFull(2: end-1 2: end-1)= struc 89 [nely, nelx]=size(struc)i 52 8 Use "bwdist "(Image Processing Toolbox. oK= sparse(2*(nelx+1)*(nely+1),2*(nelx 53 lsf =(-strucFull)*(waist(strucFull 1)*(ne1y+1)); 0.5)-strucFull*(bwdist(strucFull-1) 91 F= sparse(2*(nely+1)*(nelx+1),1);U 0.5 zers(2*(ne1y+1)*(ne1x+1),1); DESIGN UPDATE 92 for elx =1:nelx 【实例截图】
【核心代码】