//////////////////////////////////// // // // Optimisation de l'épaisseur // // d'une plaque élastique // // // // Ecole Polytechnique, MAP 562 // // Copyright G. Allaire, 2004 // // // //////////////////////////////////// ofstream compli("plaque.obj") ; string sauve="plaque"; // Nom du fichier de sauvegarde int niter=30; // Nombre d'itérations int n=2; // Taille du maillage real lagrange=1.; // Multiplicateur de Lagrange pour le volume real lagmin, lagmax ; // Encadrement du multiplicateur de Lagrange int inddico ; real compliance; // Compliance 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; // Legende 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 nulle mesh Th, Dh; 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; }; ////////////////////////////// // Construction du maillage // ////////////////////////////// n = 2 ; int m = 3 ; Th= buildmesh (bd(10*n)+bs(20*n)+bg(10*n)+b1(m)+b2(8*n)+b3(2*m)+b4(8*n)+b5(m)); plot(Th,wait=0,ps=sauve+".maillage.eps"); savemesh(Th,sauve+".msh"); 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]; Vh0 h,hold,densite; h = h0 ; /////////////////////////// // 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) ; ////////////////////////////// // Calcul du volume initial // ////////////////////////////// volume0=int2d(Th)(h); ///////////////////////// // Compliance initiale // ///////////////////////// elasticite; compliance=int1d(Th,3)(g1*u+g2*v); compli << compliance << endl ; cout<<"Initialisation. Compliance: "< volume0) { lagmax = lagrange ; while (volume > volume0) { lagrange = 2*lagrange ; h = sqrt(densite*hold/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 = sqrt(densite*hold/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: "<