//////////////////////////////////////////////// // Compliance Optimization // // by the boundary variation method // // Copyright O.Pantz, G. Allaire (2005) // //////////////////////////////////////////////// // Optimization Parameters int niter=120; //Number of iterations real errelas=0.01; //Error for the resolution of the variational problems real metricmsh=5.; //Metric used for the coarse mesh: Bigger=finer mesh real initialstep=0.2; //Initial step real maxstep=-1.; //Maximum step (<0 => no limit) real lagrangestep=3.; //Updating step of the Lagrange multiplier associated to the Volume's constraint bool tocontinue=false; //false=start a new computation, true=continue a previous computation string backupname="cantilever"; //Name under which the results are stored int nsave=1; //Number of iterations to be saved real volume0=95.; //=0 Constant Volume during the iterations, >0 Target Volume, <0 No Volume's constraint // construction of the initial mesh mesh Sh,Th,Shprev,Thprev; //Definition of the boundary labels int neumann=1; int dirichlet=2; int free=3; // right boundary (with neumann conditions) border a(t=-1,1) { x=20; y=t;label=neumann;}; // left boundary (with free or dirichlet conditions) border c1(t=4,2) { x=0; y=t;label=dirichlet; }; border c2(t=2,-2) { x=0; y=t;label=free; }; border c3(t=-2,-4) { x=0; y=t;label=dirichlet; }; // upper and lower boundary (with free conditions) border b(t=1,0) { x=20.*t; y=4.-3.*t;label=free; }; border d(t=0,1) { x=20.*t; y=-4.+3.*t;label=free; }; // Holes border Hole1(t=0,2*pi) {x=0.2*cos(t)+6;y=0.2*sin(t)+2;label=free;}; border Hole2(t=0,2*pi) {x=0.2*cos(t)+12;y=0.2*sin(t)-1;label=free;}; Sh= buildmesh (b(30)+c3(10)+c2(10)+c1(10)+d(30)+a(20)+Hole1(-11)+Hole2(-11)); Th=Sh; Thprev=Th; // cut function used to regularize the gradient shape near singularities func cut=1-(x>19.5)*(abs(y)<1.5)-(x<0.5)*(y>1.5)*(y<4.5)-(x<0.5)*(y<-1.5)*(y>-4.5); // Material parameters real E=15.; //Young Modulus (always positif) real nu=0.35; //Poisson ratio (between -1 and 1/2) real lambda=E*nu/((1.+nu)*(1.-2.*nu)); //Lame Moduli real mu=E/(2.*(1.+nu)); //Applied Loads func g1=0.; func g2=-1.; //Definition of the variables and intialization //Finite elements spaces on the finer mesh Sh fespace Wh(Sh,P1); fespace WSh(Sh,[P2,P2]); //Finite elements spaces on the coarse mesh Th fespace WTh(Th,[P2,P2]); Wh gradient; //Gradient WSh [d1,d2]; //extension field WSh [theta1,theta2]; //test functions for the extension field WSh [u1,u2],[v1,v2]; //Displacement and Test functions //Interpolated Finite elements on the coarse mesh WTh [D1,D2],[Dprev1,Dprev2]; WTh [U1,U2],[Uprev1,Uprev2]; int nsaved=0; //Initialization counter of the back-up copies int nsave0=0; //Current number of back-up copies real compliance; //Compliance real volume; //Volume of the current shape real perimeter; //Perimeter of the optimized boundary string caption; //Caption for graphic outputs real GradNorm; //Gradient Norm real GradProduct; //Scalar product between the current gradient and the previous one bool rollback=false; //Did we roll back ? real actualstep; //Actual step size real lagrange; //Lagrange Multiplier associated to the Volume constraint real[int] isoval(20); //isovalue used for graphic outputs int constrained=1; if (volume0<0) {constrained=0;}; ///////////////////////// // Elasticity Problem // ///////////////////////// problem elasticity([u1,u2],[v1,v2]) = int2d(Sh)( 2.*mu*(dx(u1)*dx(v1)+dy(u2)*dy(v2)+((dx(u2)+dy(u1))*(dx(v2)+dy(v1)))/2.) +lambda*(dx(u1)+dy(u2))*(dx(v1)+dy(v2))) -int1d(Sh,neumann)(g1*v1+g2*v2) +on(dirichlet,u1=0,u2=0); //////////////////////////////////////// // Expression of the shape's Gradient // //////////////////////////////////////// macro gradientexp() -2.*mu*(dx(u1)^2+dy(u2)^2+((dx(u2)+dy(u1))^2)/2.)-lambda*(dx(u1)+dy(u2))^2 // /////////////////////// // Extension Problem // /////////////////////// // H1 scalar product between vector-valued functions macro prodscal(t1,t2,p1,p2) dx(t1)*dx(p1)+dy(t1)*dy(p1)+dx(t2)*dx(p2)+dy(t2)*dy(p2)+t1*p1+t2*p2 // problem extension([d1,d2],[theta1,theta2]) = int2d(Sh)(prodscal(d1,d2,theta1,theta2)) +int1d(Sh,free) ((theta1*N.x+theta2*N.y)*(gradientexp + lagrange*constrained)*cut) +on(dirichlet,neumann,d1=0,d2=0); ////////////////////////////// // Initialization // ////////////////////////////// int iter,iter0; string backupnameNb; if(!tocontinue){ iter0=1; if (volume0==0) {volume0=int1d(Sh)(x*N.x+y*N.y)/2;} // Computation of the initial Volume elasticity; // Solving the State Equation compliance=int1d(Sh,neumann)(g1*u1+g2*u2); // Initial Compliance gradient=-(gradientexp); // Computation of the shape derivative //Evaluation of the Lagrande Multiplier associated to the Volume constraint lagrange=int1d(Sh,free)(gradient*cut)/int1d(Sh,free)(cut); gradient=(gradient-lagrange*constrained)*cut; //Mesh Adaptation Th=Sh; [U1,U2]=[u1,u2]; //interpolation of u on the coarse mesh Th [Uprev1,Uprev2]=[0,0]; Uprev1[]=U1[]; Uprev2[]=U2[]; Sh=adaptmesh (Th,Uprev1,Uprev2,err=errelas); actualstep=initialstep; //Initialization of the initial step //We save the intial shape if (nsave!=0) { //isovalue for graphic output real gmax=gradient[].max; real gmin=gradient[].min; for (int j=0;j>volume0;} else {real next; f >>next;}; f >>lagrange; //Lagrange Multiplier f >>actualstep; //Actual step size f >>rollback; //Did we roll back ? f >>iter0; //iteration de départ f >>nsave0; //Number of previous back-up copies f >>isoval; f >>Dprev1[];f>>Dprev2[];//Previous extenstion fields nsaved=nsave0; nsave0=nsave0+1;}; ///////////////////////////////////// // Initialization of the Algorithm // ///////////////////////////////////// elasticity; //Solving the state equation compliance=int1d(Sh,neumann)(g1*u1+g2*u2); //Initial compliance gradient=-(gradientexp); //Computation of the shape derivative //Evaluation of the Lagrange Multiplier associated to the Volume constraint lagrange=int1d(Sh,free)(gradient*cut)/int1d(Sh,free)(cut); gradient=(gradient-lagrange*constrained)*cut; /////////////////////////////// // Optimization Loop // /////////////////////////////// for (iter=iter0;iter< niter+iter0;iter=iter+1){ cout <<"Iteration " <0)&(GradProduct<(GradNorm/2.))){ actualstep=actualstep*2*GradNorm/(2*GradNorm-GradProduct); cout<<"We increase the step size"<(GradNorm/2.)){ actualstep=actualstep*4./3.; cout<<"We increase the step size"<0.) actualstep=min(actualstep,maxstep); Thprev=Th; //We keep a copy of the meshs (in case of backup at the next loop's step) Shprev=Sh; ////////////////////////////////// // We move the coarse mesh Th // ////////////////////////////////// if (!rollback){ real aa,minaire=checkmovemesh (Th,[x,y])/10000.; while (minaire > (aa=checkmovemesh(Th,[x+actualstep*D1,y+actualstep*D2])) ){ cout << " problem reversed triangle => actualstep= actualstep/2. "<