//////////////////////////////////////////////// // Grip Optimization // // by the boundary variation method // // Copyright O.Pantz, G. Allaire (2005) // //////////////////////////////////////////////// // Optimization Parameters int niter=60; //Number of iterations real errelas=0.005; //Error for the resolution of the variational problems real metricmsh=10.; //Metric used for the coarse mesh: Bigger=finer mesh real initialstep=1.; //Initial step real maxstep=3.; //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="grip"; //Name under which the results are stored int nsave=1; //Number of iterations to be saved real volume0=-1.; //=0 Constant Volume during the iterations, >0 Target Volume, <0 No Volume's constraint real penVolume=0.05; //Volume penalization func cutVolume=(x<5.0); // construction of the initial mesh mesh Sh,Th,Shprev,Thprev; //Boundaries with (non null)-Neumann and Dirichlet conditions are not optimized int neumann=1; int dirichlet=2; int teeth=4; //Boundaries with free conditions are optimized int free=3; border a(t=3,-3){x=0;y=t;label=free;}; //Left Boundary (Optimized Boundary) border b1(t=0,1){x=t;y=-3;label=neumann;}; //Low Boundary (with Neumann of Free Conditions) border c1(t=0,1){x=1+5*t;y=-3+2.5*t;label=free;}; border dd1(t=6,5){x=t;y=-0.5;label=teeth;}; //Teeth of the grip border e(t=-0.5,0.5){x=5;y=t;label=free;}; border dd2(t=5,6){x=t;y=0.5;label=teeth;}; border c2(t=1,0){x=1+5*t;y=3-2.5*t;label=free;}; //Upper Boundary (with Neumann of Free Conditions) border b2(t=1,0){x=t;y=3;label=neumann;}; //A Hole with Dirichlet conditions (not Optimized) border GammaD(t=0,2*pi){x=cos(t)*0.2+3;y=sin(t)*0.2;label=dirichlet;}; //Optimized Holes border Hole1(t=0,2*pi){x=cos(t)*0.5+2.;y=sin(t)*0.5+1.5;label=free;}; border Hole2(t=0,2*pi){x=cos(t)*0.5+2.;y=sin(t)*0.5-1.5;label=free;}; Sh= buildmesh (a(15)+b1(5)+c1(10)+dd1(4)+e(4)+dd2(4)+c2(10)+b2(5)+GammaD(10)+Hole1(-10)+Hole2(-10)); Th=Sh; Thprev=Th; // cut function used to regularize the gradient shape near singularities func cut=1.-(x<6)*(x>5)*(y>-0.5)*(y<0.5); // Material parameters real E=15.; //Young Modulus (always positif) real nu=0.4; //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=1.; func g2=0.; //Target displacement func f1target=0.; func f2target=-y*(x>3.); //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 WSh [u1target,u2target]; WSh [p1,p2]; //Adjoint //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 Pression; //Pression on the griped body real Cost; //Value of the cost function 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 grapical 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)+on(teeth,u2=0); /////////////////////// // Adjoint Problem // /////////////////////// problem adjoint([p1,p2],[v1,v2]) = int2d(Sh)( 2.*mu*(dx(p1)*dx(v1)+dy(p2)*dy(v2)+((dx(p2)+dy(p1))*(dx(v2)+dy(v1)))/2.) +lambda*(dx(p1)+dy(p2))*(dx(v1)+dy(v2))) -int2d(Sh)( 2.*mu*(dx(u1target)*dx(v1)+dy(u2target)*dy(v2)+((dx(u2target)+dy(u1target))*(dx(v2)+dy(v1)))/2.) +lambda*(dx(u1target)+dy(u2target))*(dx(v1)+dy(v2))) +on(dirichlet,p1=0,p2=0)+on(teeth,p2=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)*Pression/(compliance*compliance) +(2.*mu*((dx(u1target)-dx(p1))*dx(u1)+(dy(u2target)-dy(p2))*dy(u2) +((dx(u2target)-dx(p2)+dy(u1target)-dy(p1))*(dx(u2)+dy(u1)))/2.) +lambda*(dx(u1target)-dx(p1)+dy(u2target)-dy(p2))*(dx(u1)+dy(u2)))/compliance+penVolume*cutVolume // /////////////////////// // 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,teeth,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 [u1target,u2target]=[f1target,f2target]; adjoint; // Solving Adjoint Equation compliance=int1d(Sh,neumann)(g1*u1+g2*u2); // Initial Compliance Pression=int2d(Sh)( 2.*mu*(dx(u1target)*dx(u1)+dy(u2target)*dy(u2)+((dx(u2target)+dy(u1target))*(dx(u2)+dy(u1)))/2.) +lambda*(dx(u1target)+dy(u2target))*(dx(u1)+dy(u2))); Cost=Pression/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 [u1target,u2target]=[f1target,f2target]; adjoint; // Solving Adjoint Equation compliance=int1d(Sh,neumann)(g1*u1+g2*u2); //Initial compliance Pression=int2d(Sh)( 2.*mu*(dx(u1target)*dx(u1)+dy(u2target)*dy(u2)+((dx(u2target)+dy(u1target))*(dx(u2)+dy(u1)))/2.) +lambda*(dx(u1target)+dy(u2target))*(dx(u1)+dy(u2))); Cost=Pression/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. "<