//////////////////////////////////// // // // Optimisation d'une console // // (cantilever en anglais) // // // // Ecole Polytechnique, MAP 562 // // Copyright G. Allaire, O. Pantz // // 2004 // // // //////////////////////////////////// ofstream compli("cantilever.obj") ; ofstream grad("cantilever.grad") ; string sauve="cantilever"; // Nom du fichier de sauvegarde int niter=20; // Nombre d'itérations real errelas=0.05; // Niveau d'erreur souhaité pour l'adaptation de maillage real pasinit=1.; // Pas de descente initial real pasef; // Pas de descente effectif real muv; // Multiplicateur de Lagrange associé au volume real pasmu=2; // Pas de mise à jour du multiplicateur de Lagrange int nsaved=0; // Initialisation du compteur des sauvegardes real compliance; // Compliance real volume0; // Volume initial real volume; // Volume de la forme déformée string legende; // Légende pour les sorties graphiques real norme; // Norme du gradient real produit; // Produit scalaire entre le gradient actuel et le gradient précédent bool arriere=false; // Paramètre de controle du retour en arrière real E=15; // Module de Young (toujours positif) real nu=0.35; // Coefficient de Poisson (toujours compris entre -1 et 1/2) real pi=4*atan(1) ; func g1=0; // Forces appliquées func g2=-1; ///////////////////////////////////// // Calcul des coefficients de Lamé // ///////////////////////////////////// real lambda=E*nu/((1.+nu)*(1.-2.*nu)); real mu=E/(2.*(1.+nu)); ////////////////////////////////////////////////// // Fonction définissant la zone à optimiser // // (=1 sur la frontière optimisable, =0 sinon) // ////////////////////////////////////////////////// func cut=1-(x>8.5)*(abs(y)<1.5)-(x<0.5)*(y>1.5)*(y<4.5)-(x<0.5)*(y<-1.5)*(y>-4.5); // définition de la frontière à optimiser func real d1(real & tt) { return 9.0*tt;; }; func real d2(real & tt) {return -4.+3.*tt;}; func real b1(real & tt) {real r=1.-tt; return d1(r); }; func real b2(real & tt) {real r=1.-tt; return -d2(r);}; ///////////////////////////// 1:Condition de Dirichlet // Définition du domaine // 2:Condition Libre ///////////////////////////// 3:Condition de Neumann non nulle mesh Sh; border a(t=-1,1) { x=9; y=t;label=3; }; // bord droit de la forme border c1(t=4,2) { x=0; y=t;label=1; }; // bord gauche de la forme border c2(t=2,-2) { x=0; y=t;label=2; }; border c3(t=2,4) { x=0; y=-t;label=1; }; border b(t=0,1) { x=b1(t); y=b2(t);label=2; };// partie supérieure de la forme border d(t=0,1) { x=d1(t); y=d2(t);label=2; };// partie inférieure de la forme real ray=0.5 ; border z1(t=0,2*pi) {x=7.5+ray*cos(t);y=0.+ray*sin(t);label=2;}; border z2(t=0,2*pi) {x=6.+ray*cos(t);y=0.8+ray*sin(t);label=2;}; border z3(t=0,2*pi) {x=6.+ray*cos(t);y=-0.8+ray*sin(t);label=2;}; border z4(t=0,2*pi) {x=4.5+ray*cos(t);y=0.+ray*sin(t);label=2;}; border z5(t=0,2*pi) {x=3.+ray*cos(t);y=-1.+ray*sin(t);label=2;}; border z6(t=0,2*pi) {x=3.+ray*cos(t);y=1.+ray*sin(t);label=2;}; border z7(t=0,2*pi) {x=1.5+ray*cos(t);y=0.+ray*sin(t);label=2;}; ///////////////////////////////////// // Construction des deux maillages // // (Sh = fin, Th = grossier) // ///////////////////////////////////// //Sh= buildmesh (b(30)+c3(10)+c2(20)+c1(10)+d(30)+a(20)+z1(-10)+z2(-10)+z3(-10)+z4(-10)+z5(-10)+z6(-10)); Sh= buildmesh (b(30)+c3(10)+c2(10)+c1(10)+d(30)+a(20)+z1(-10)+z2(-10)+z3(-10)+z4(-10)+z5(-10)+z6(-10)+z7(-10)); //Sh= buildmesh (b(30)+c3(10)+c2(10)+c1(10)+d(30)+a(20)+z1(-10)+z2(-10)+z3(-10)+z4(-10)); //Sh= buildmesh (b(30)+c3(10)+c2(10)+c1(10)+d(30)+a(20)); mesh Th=Sh; fespace Vh1(Sh,P1); //Définition de l'espace P1 fespace Vh2(Sh,[P2,P2]); //Définition de l'espace P2*P2 fespace Wh2(Th,[P2,P2]); //Définition de l'espace P2*P2 Vh2 [u,v],[w,s],[eu,ev]; Vh1 gradient, energ; Wh2 [ui,vi],[eui,evi],[ui2,vi2],[eui2,evi2]; /////////////////////////// // Résolution du système // /////////////////////////// problem elasticite([u,v],[w,s]) = int2d(Sh)( 2.*mu*(dx(u)*dx(w)+dy(v)*dy(s)+((dx(v)+dy(u))*(dx(s)+dy(w)))/2.) +lambda*(dx(u)+dy(v))*(dx(w)+dy(s)) ) -int1d(Sh,3)(g1*w+g2*s) +on(1,u=0,v=0) ; /////////////////////////////////////////////// // Problème d'extension du gradient de forme // // (qui doit être étendu à l'intérieur) // /////////////////////////////////////////////// problem extension([eu,ev],[w,s]) = int2d(Sh)( dx(eu)*dx(w)+dy(eu)*dy(w)+dx(ev)*dx(s)+dy(ev)*dy(s)+eu*w+ev*s ) -int1d(Sh,2) ((w*N.x+s*N.y)*(2*mu*(dx(u)^2+dy(v)^2+((dx(v)+dy(u))^2)/2.) +lambda*(dx(u)+dy(v))^2 - muv)*cut) +on(1,3,eu=0,ev=0) ; ////////////////////////////// // Calcul du volume initial // // (qui sera volume cible) // ////////////////////////////// volume0=int1d(Sh)(x*N.x+y*N.y)/2; cout <<"Volume initial = "<0), on augmente le pas de descente. // Si le gradient de la fonction cout n'est pas dans la même direction que le // gradient de l'itération précédente (produit<0), on fait marche arrière et // on diminue le pas de descente (on évite ainsi les oscillations). if ((produit<0)&(!arriere)) {pasef=pasef*norme/(norme-produit)/4.; cout<<"******* Pas de descente trop important: Marche arriere !!! *******"<0)&(produit<(norme/2.))) {pasef=pasef*2*norme/(2*norme-produit); cout<<"On augmente le pas de descente"<(norme/2.)) {pasef=pasef*4./3.; cout<<"On augmente le pas de descente"< (aa=checkmovemesh(Th,[x+pasef*eui,y+pasef*evi])) ) { cout << "******* Probleme de triangle retourne: on diminue le pas de descente, " << minaire << " > " << aa<<" *******"<< endl; pasef= pasef/2; } Th = movemesh (Th,[x+pasef*eui,y+pasef*evi]); cout<<"Pas de descente effectif = "<