///////////////////////////////////////////////////// // Méthode SIMP pour la minimisation // // de la compliance d'un pont // // // // Ecole Polytechnique, MAP 562 // // Copyright G. Allaire, 2004 // ///////////////////////////////////////////////////// ofstream fichobj("pont.obj") ; string sauve="pont"; // nom du fichier de sauvegarde int niternp=10; // nombre minimum d'it?rations non p?nalis?es int niter=100; // nombre d'it?rations total int n=3; // Taille du maillage real compliance,compliance0,compliance1; // Compliance string legende; // L?gende pour les sorties graphiques real E=1; // Module de Young real nu=0.3; // Coefficient de Poisson real lagrange=1; // Multiplicateur de Lagrange pour le volume real lagmin, lagmax; // Encadrement du multiplicateur de Lagrange real masse0; // Masse initiale real masse,masse1; // Masse de la forme courante int inddico ; // Indice dans la recherche du multiplicateur de Lagrange real thetamin=0.01; // Densit? minimale de mati?re real thetamoy=0.3; // Densit? moyenne de mati?re real pi=4*atan(1); // Valeur de pi=3.14159... real p=1; // Exposant de p?nalisation real gris,gris0; // Indicateur de p?nalisation int ii0=0; // Indice de mise à jour de l'exposant de p?nalisation func f1=0; // Forces appliqu?es func f2=-1; // D?finition des isovaleurs de trac? 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)); ///////////////////////////// 1:Condition de Dirichlet // Définition du domaine // 2:Condition Libre ///////////////////////////// 3:Condition de Neuman non nulle mesh Th; border bd(t=0,1.2) { x=+1; y=t;label=2; }; // bord droit de la forme border bg(t=1.2,0) { x=-1; y=t;label=2; }; // bord gauche de la forme border bs(t=1,-1) { x=t; y=1.2;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 // ////////////////////////////// Th= buildmesh (bd(12*n)+bs(20*n)+bg(12*n)+b1(n)+b2(8*n)+b3(2*n)+b4(8*n)+b5(n)); plot(Th,wait=1,ps=sauve+".maillage.eps"); savemesh(Th,sauve+".msh"); fespace Vh0(Th,P0); fespace Vh2(Th,[P2,P2]); Vh2 [u,v],[w,s]; Vh0 coef,theta,eu,ev,euv,intermed; // theta = densit? de mati?re // coef = theta^p func theta0=thetamoy; theta=theta0; coef=theta^p; // Masse initiale real volume=int2d(Th)(1.); masse0 = int2d(Th)(theta); masse0 = masse0/volume ; /////////////////////////// // Définition du système // /////////////////////////// problem elasticite([u,v],[w,s]) = int2d(Th)( 2.*coef*mu*(dx(u)*dx(w)+dy(v)*dy(s)+((dx(v)+dy(u))*(dx(s)+dy(w)))/2.) +coef*lambda*(dx(u)+dy(v))*(dx(w)+dy(s)) ) -int1d(Th,3)(f1*w+f2*s) +on(1,v=0) ; ////////////////////////////////// // Boucle d'optimisation // ////////////////////////////////// int iter; for (iter=1;iter< niter+1;iter=iter+1) { elasticite; compliance1 = compliance ; compliance = int1d(Th,3)(f1*u+f2*v); cout <<"Iteration " < masse0) { lagmax = lagrange ; while (masse > masse0) { lagrange = 2*lagrange ; theta = (p*intermed*coef/lagrange)^(1/(p+1)); theta = min(1,theta); theta = max(theta,thetamin); masse = int2d(Th)(theta)/volume ; }; lagmin = lagrange ; }; // if (masse1 == masse0) { lagmin = lagrange ; lagmax = lagrange ; }; // dichotomie sur le multiplicateur de Lagrange inddico=0; while ((abs(1.-masse/masse0)) > 0.000001) { lagrange = (lagmax+lagmin)*0.5 ; theta = (p*intermed*coef/lagrange)^(1/(p+1)); theta = min(1,theta); theta = max(theta,thetamin); masse = int2d(Th)(theta)/volume ; inddico=inddico+1; if (masse < masse0) {lagmin = lagrange ;} ; if (masse > masse0) {lagmax = lagrange ;} ; }; cout<<"Nombre d'iterations de Lagrange "<= niternp) { if ( (abs(compliance-compliance1) < 0.01*compliance0)&(ii0 > 4) ) {gris=int2d(Th)((theta-thetamin)*(1.-theta))*4/volume ; p=p*(1+0.3^(1.+gris*0.5)) ; p=min(p,4.); ii0 = 0 ; cout << "gris = "<