//////////////////////////////////// // // // Optimisation de l'épaisseur // // d'une plaque élastique // // avec état adjoint // // // //////////////////////////////////// ofstream objectif1("objectif1.txt") ; ofstream gradient1("gradient1.txt") ; string sauve="plaqadj"; // Nom du fichier de sauvegarde int niter=40; // Nombre d'itérations int n=3; // Taille du maillage real lagrange=0.01; // Multiplicateur de Lagrange pour le volume real pas=0.1 ; // Pas de descente real lagmin, lagmax ; int inddico ; real objectif; // Fonction cout real volume0; // Volume initial real volume,volume1; // Volume de la forme courante real exa=1; // Coefficient d'exagération de la déformation string legende; // Légende pour les sorties graphiques real E=100; //Module de Young (toujours positif) real nu=0.3; //Coefficient de Poisson (toujours compris entre -1 et 1/2) func g1=0; //Forces appliquées func g2=-100; real[int] vviso(21); for (int i=0;i<21;i++) vviso[i]=i*0.05 ; ///////////////////////////////////// // Calcul des coefficients de Lamé // ///////////////////////////////////// real lambda=E*nu/((1.+nu)*(1.-2.*nu)); real mu=E/(2.*(1.+nu)); ///////////////////////////////////////// // Epaisseur mini et maxi de la plaque // ///////////////////////////////////////// real hmin=0.1; real hmax=1.0; real hmoy=0.5; func h0=hmoy ; ///////////////////////////// 1:Condition de Dirichlet // Définition du domaine // 2:Condition Libre ///////////////////////////// 3:Condition de Neuman non nulles mesh Th; border bd(t=0,1) { x=+1; y=t;label=2; }; // bord droit de la forme border bg(t=1,0) { x=-1; y=t;label=2; }; // bord gauche de la forme border bs(t=1,-1) { x=t; y=1;label=2; }; // bord supérieur de la forme border b1(t=-1,-0.9) { x=t; y=0;label=1; }; // bord inférieur de la forme border b2(t=-0.9,-0.1) { x=t; y=0;label=2; }; border b3(t=-0.1,0.1) { x=t; y=0;label=3; }; border b4(t=0.1,0.9) { x=t; y=0;label=2; }; border b5(t=0.9,1) { x=t; y=0;label=1; }; real center=0.9 ; border a5(t=0,2*pi) { x=0.-0.1*sin(t);y=center+0.1*cos(t);label=5; }; ////////////////////////////// // Construction du maillage // ////////////////////////////// int m = int(n*1.5) ; Th= buildmesh (bd(10*n)+bs(20*n)+bg(10*n)+b1(m)+b2(8*n)+b3(2*m)+b4(8*n)+b5(m)+a5(5*m)); plot(Th,wait=0); fespace Vh0(Th,P0); //Définition de l'espace P0 fespace Vh2(Th,[P2,P2]); //Définition de l'espace P2*P2 Vh2 [u,v],[w,s],[p,q]; Vh0 h,hold,gradient; h = h0 ; /////////////////////// // Déplacement cible // /////////////////////// real u0x=0.; real cx=0.; real u0y=1.; //real cy=1.; func u0=u0x; func v0=u0y; Vh0 reg = region ; int indiceinclusion = reg (0.,center) ; Vh0 cy=1.*(region==indiceinclusion); /////////////////////////// // Définition du système // /////////////////////////// problem elasticite([u,v],[w,s]) = int2d(Th)( 2.*h*mu*(dx(u)*dx(w)+dy(v)*dy(s)+((dx(v)+dy(u))*(dx(s)+dy(w)))/2.) +h*lambda*(dx(u)+dy(v))*(dx(w)+dy(s)) ) -int1d(Th,3)(g1*w+g2*s) +on(1,u=0,v=0) ; ///////////////////////////// // Definition de l'adjoint // ///////////////////////////// problem adjoint([p,q],[w,s]) = int2d(Th)( 2.*h*mu*(dx(p)*dx(w)+dy(q)*dy(s)+((dx(q)+dy(p))*(dx(s)+dy(w)))/2.) +h*lambda*(dx(p)+dy(q))*(dx(w)+dy(s)) ) +int2d(Th)( 2*(cx*(u-u0)*w+cy*(v-v0)*s) ) +on(1,p=0,q=0) ; ////////////////////////////// // Calcul du volume initial // ////////////////////////////// volume0=int2d(Th)(h); //////////////////////////// // Fonction cout initiale // //////////////////////////// elasticite; objectif=int2d(Th)( cx*(u-u0)*(u-u0)+cy*(v-v0)*(v-v0) ); objectif1 << objectif << endl ; cout<<"Initialisation. Cout: "< volume0) { lagmax = lagrange ; while (volume > volume0) { lagrange = lagrange + 0.01 ; h = hold - pas*gradient - lagrange ; h = min(h,hmax) ; h = max(h,hmin) ; volume=int2d(Th)(h) ; }; lagmin = lagrange ; }; // if (volume1 == volume0) { lagmin = lagrange ; lagmax = lagrange ; }; ////////////////////////////////////////////////// // Dichotomie sur le multiplicateur de lagrange // ////////////////////////////////////////////////// inddico=0; while ((abs(1.-volume/volume0)) > 0.000001) { lagrange = (lagmax+lagmin)*0.5 ; h = hold - pas*gradient - lagrange ; h = min(h,hmax) ; h = max(h,hmin) ; volume=int2d(Th)(h) ; inddico=inddico+1; if (volume < volume0) {lagmin = lagrange ;} ; if (volume > volume0) {lagmax = lagrange ;} ; }; //cout<<"nbre iterations dichotomie: "<