/////////////////////////////////////////////////////// // Level set method for the geometric optimization // // of a minimal compliance cantilever // // // // Ecole Polytechnique, MAP 562 // // Copyright G. Allaire, C. Dousset 2007 // /////////////////////////////////////////////////////// // La fonction ligne de niveaux est appelée phi. La // forme est caractérisée par phi<0. La fonction phi // est convectée par un schéma caractéristique implicite // (cf. J. Strain, JCP 151, 1999). // The level set function is denoted by phi. The shape // is characterized by phi<0. The function phi is advected // by a caracteristic implicit scheme (cf. J. Strain, // JCP 151, 1999). /////////////////////////////////////////////////////// int n=2; //mesh size parameter real mu=8; //Lame coefficient real lambda=1; //Lame coefficient real g1=0; //horizontal component of the applied force on the boundary real g2=-5; //vertical component of the applied force on the boundary real lag=0.5; //Lagrange multiplier for the weight real T=.1; //descent step or time step for the advection real TT=0.01; //time step for the re-initialization real weak=0.0001; //multiplicative coefficient for the "weak" material mimicking void real eps1=0.001; //small parameter avoiding division by zero in the normal computation real eps2=0.001; //regularization parameter for the velocity real tol=0.0001; //tolerance for the increase of the objective function real compliance; real weight; real objective; real objectiveinitial; real compliancem; real weightm; real objectivem; int kx=4; //number of holes in the horizontal direction of the initial shape int ky=3; //number of holes in the vertical direction of the initial shape int N=30; //number of iterations for the gradient algorithm int Ninit=5; //number of iterations for the re-initialization int iter; int iterinit; string caption; string save="cantilever-ls"; ofstream object("cantilever-ls.obj"); //Definition of the boundary and of the mesh // right boundary of the mesh border bd1(t=-0.5,-0.05) { x=1; y=t; label=1; }; border bd2(t=-0.05,0.05) { x=1; y=t; label=3; }; border bd3(t=0.05,0.5) { x=1; y=t; label=1; }; //left boundary of the mesh border bg(t=0.5,-0.5) { x=-1; y=t; label=2; }; // upper boundary of the mesh border bs(t=1,-1) { x=t; y=0.5; label=1; }; // lower boundary of the mesh border bi(t=-1,1) { x=t; y=-0.5; label=1; }; // label meanings: 1=Neumann, 3=applied force, 2=Dirichlet mesh Th= buildmesh (bd1(9*n)+bd2(2*n)+bd3(9*n)+ bg(20*n)+bs(40*n)+bi(40*n)); //Definition of finite element spaces fespace Vh(Th,[P1,P1]); fespace Vh1(Th,P1); Vh [u,v],[w,s],[um,vm]; Vh1 N1, N2, phi, V, phim, M1, M2, S, phiinit, nabla, Xapprox, X, Vreg, vv, Xapproxm, Xm; // Definition of the initial shape (perforated by holes according to kx and ky) func phi0=-0.1+sin(pi*kx*x)*sin(pi*ky*(y-0.5)); phi=phi0; plot(phi,wait=0,fill=1,value=1,ps=save+".init.eps"); // Re-initialization of phi such that it becomes the signed distance // function to the shape boundary. Implicit linearized solving of : // (d phi)/(d t) + sign(phi) ( |gradphi)| - 1 ) = 0 Vh1 h1=hTriangle; real h=h1[].max; //maximal size of a mesh triangle cout<<"maximal size of a mesh triangle h="<objective*(1+tol)) { T=T/2; cout <<"descent step rejected"<