//////////////////////// PARAMETERS //////////////////////////// verbosity=0; macro g [0,-1] // Applied surface load real Y=15; // Young Modulus real nu=0.35; // Poisson Ratio real lambda=Y*nu/((1.+nu)*(1.-2.*nu)); // Lame Moduli real mu=Y/(2.*(1.+nu)); int Niter=50; // Number of iterations real penVolume=0.1; // Penalization of the volume real GeoPrec=0.02; // scale of the geometric mesh real step=8.; // initial descent step string name="bridge"; bool poursuit=false,adapt=false; // Restart at previous computed state ? / Adapt mesh wrt the displacement ? /////////////////////////////////////////// //// Definition of the initial mesh //// /////////////////////////////////////////// int Free,Neumann1,Neumann2,Dirichlet,Dirichlet0; // Labels for optimized boundaries, Neumann and Dirichlet condidtions {Free=1;Neumann1=2;Neumann2=3;Dirichlet=4;Dirichlet0=8;}; int FreezeX,FreezeY,FreezeXY; // Freeze degree of freedom of theta {FreezeX=5;FreezeY=6;FreezeXY=7;} mesh Sh; // The Total mesh real H=0.2; // height of the bridge real L=1.; // length of the bridge if(poursuit) {Sh=readmesh(name+"-Sh.msh");} else { // The definition of the initial mesh could be modified as will // // -- a bridge -- real support=0.1; // support of the bridge (Dirichlet) real Sneumann=0.1; // suport of applied loads real deltax=L/5.,deltay=H/5.; border Bottom1(t=0,support-support/10){x=t;y=0;label=Dirichlet;} border Bottom10(t=support-support/10,support){x=t;y=0;label=Dirichlet0;} border Bottom2(t=support,L/3.-Sneumann/2.){x=t;y=0;label=Free;} border Bottom3(t=L/3.-Sneumann/2.,L/3.+Sneumann/2.){x=t;y=0;label=Neumann1;} border Bottom4(t=L/3.+Sneumann/2.,2.*L/3.-Sneumann/2.){x=t;y=0;label=Free;} border Bottom5(t=2.*L/3.-Sneumann/2.,2.*L/3.+Sneumann/2.){x=t;y=0;label=Neumann2;} border Bottom6(t=2.*L/3.+Sneumann/2.,L-support){x=t;y=0;label=Free;} border Bottom70(t=L-support,L-support+support/10){x=t;y=0;label=Dirichlet0;} border Bottom7(t=L-support+support/10,L){x=t;y=0;label=Dirichlet;} border Arc(t=pi,0){x=cos(t)*L/2.+L/2.;y=sin(t)*H;label=Free;} real r=H/10.; border Hole1(t=0.,2.*pi){x=r*cos(t)+L/3;y=r*sin(t)+H/2.;label=Free;} border Hole2(t=0.,2.*pi){x=r*cos(t)+2.*L/3.;y=r*sin(t)+H/2.;label=Free;} border Hole3(t=0.,2.*pi){x=r*cos(t)+L/2.+L/3;y=r*sin(t)+H/2.;label=Free;} // -- BIG BOX containing the shape --// int n=10; border DLeft(t=2.*H,-deltay){x=-deltax;y=t;label=FreezeX;} border DBottom(t=-deltax,L+deltax){x=t;y=-deltay;label=FreezeY;} border DRight(t=-deltay,2.*H){x=L+deltax;y=t;label=FreezeX;} border DTop(t=L+deltax,-deltax){x=t;y=2.*H;label=FreezeY;} Sh=buildmesh( DLeft(H*n)+DBottom(L*n)+DRight(H*n)+DTop(L*n) +Bottom1(n*support)+Bottom10(1)+Bottom2(n*(L/3))+Bottom3(n*Sneumann)+Bottom4(n*(L/3))+Bottom5(n*Sneumann)+Bottom6(n*(L/3))+Bottom70(1)+Bottom7(n*support) +Arc(L*n)+Hole1(10)+Hole2(10));//+Hole3(10)); } mesh Th; // The mesh of the shape int ShapeLabel; // Region label of the shape fespace Zh(Sh,P0); Zh reg=region; // Region labels if(poursuit) { ifstream Label(name+"-label"); Label>>ShapeLabel; } else { real px=L/2.; real py=H/1000.; // a point of the shape ShapeLabel=reg(px,py); } Th=trunc(Sh,reg==ShapeLabel); plot(Th,wait=1); //////////////////////////////////////// // end of initial mesh definition // //////////////////////////////////////// ////////////////////// END OF PARAMETERS LIST /////////////////// fespace Xh(Sh,[P1,P1]); // finite element space for the convection field fespace Wh(Th,[P2,P2]); // finite element space for the displacement int quadraticOrder=3; // Some macro macro u1 [u11,u12] // displacement macro v [v1,v2] // test function for elasticity macro e(u) [dx(u[0]),(dx(u[1])+dy(u[0]))/sqrt(2.),dy(u[1])] // deformation tensor macro div(u) (dx(u[0])+dy(u[1])) // divergence operator macro theta [theta1,theta2] // descent direction macro ptheta [ptheta1,ptheta2] // previous descent direction macro dtheta [dtheta1,dtheta2] // test function fot descent direction macro Theta [Theta1,Theta2] // interpolated descent direction Wh u1,v; problem elasticity1(u1,v)= int2d(Th,qforder=quadraticOrder)(2.*mu*(e(u1)'*e(v))+lambda*div(u1)*div(v)) -int1d(Th,Neumann1,Neumann2)(g'*v) +on(Dirichlet,u11=0,u12=0); // we remove Dirichlet0 (for regularity problems) Xh theta,dtheta; problem Computetheta(theta,dtheta)= int2d(Sh,qforder=3)(2.*mu*(e(theta)'*e(dtheta))+lambda*div(theta)*div(dtheta)) +int1d(Sh,Free)(jump(reg)*(2.*mu*(e(u1)'*e(u1))+lambda*div(u1)^2-penVolume)*([N.x,N.y]'*dtheta)) +on(Dirichlet0,Neumann1,Neumann2,theta1=0,theta2=0) +on(FreezeX,theta1=0) +on(FreezeY,Dirichlet,theta2=0); varf dJ(theta,dtheta)= int1d(Sh,Free)(jump(reg)*(2.*mu*(e(u1)'*e(u1))+lambda*div(u1)^2-penVolume)*([N.x,N.y]'*dtheta)) +on(Dirichlet0,Neumann1,Neumann2,theta1=0,theta2=0) +on(FreezeX,theta1=0) +on(FreezeY,Dirichlet,theta2=0); real J; // Initialisation of the Loop for(int iter=0;iterminT0/100.) { Sh=movemesh(pSh,[X,Y]);moveOK=true; } else { step/=2; cout<< " Triangle reversed = step/2 ="<0) { step/=2.; cout<<"GO BACK !!!! step="<