//////////////////////////////////// // // // Optimisation de l'épaisseur // // d'une plaque élastique // // // // Ecole Polytechnique, MAP 562 // // Copyright G. Allaire, 2004 // // Mise à jour O. Pantz 2013 // // // // Testé sous FreeFem++ v3.20 // // // //////////////////////////////////// int wait=0; // =0 pas de pause entre chaque affichage; =1 pause à chaque sortie graphique string sauve="plaque"; // Nom du fichier de sauvegarde (si vide, pas de sauvegarde effectué) real coef=0.2; // Coefficient d'exagération de la déformation 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;func g2=-100; // Forces appliquées int niter=30; // Nombre d'itérations int n=2,m=3; // Paramètres controlant la Taille du maillage real hmin=0.1,hmax=1.0; // Epaisseur minimale et maximale de la plaque func h0=0.5; // Epaisseur initiale verbosity=0; // Silence !!! ofstream compli(sauve+".obj") ; ///////////////////////////////////// // Calcul des coefficients de Lamé // ///////////////////////////////////// real lambda=E*nu/((1.+nu)*(1.-2.*nu)); real mu=E/(2.*(1.+nu)); //////////////////////////////////////////// // Paramètres pour les sorties graphiques // //////////////////////////////////////////// real[int] vviso(21); for (int i=0;i<21;i++) vviso[i]=i*0.05*(hmax-hmin)+hmin ; real[int] colors=[0,0,1,0,0,0]; ///////////////////////////// 1:Condition de Dirichlet // Définition du domaine // 2:Condition Libre ///////////////////////////// 3:Condition de Neumann non nulle mesh Th, Dh; // Maillage de référence, Maillage déformé int Dirichlet=1,Libre=2,Neumann=3; // Labels pour les conditions aux limites {// Construction du maillage // border bd(t=0,1) { x=+1; y=t;label=Libre; }; // bord droit de la forme border bg(t=1,0) { x=-1; y=t;label=Libre; }; // bord gauche de la forme border bs(t=1,-1) { x=t; y=1;label=Libre; }; // bord supérieur de la forme border b1(t=-1,-0.9) { x=t; y=0;label=Dirichlet; }; // bord inférieur de la forme border b2(t=-0.9,-0.1) { x=t; y=0;label=Libre; }; border b3(t=-0.1,0.1) { x=t; y=0;label=Neumann; }; border b4(t=0.1,0.9) { x=t; y=0;label=Libre; }; border b5(t=0.9,1) { x=t; y=0;label=Dirichlet; }; Th= buildmesh (bd(10*n)+bs(20*n)+bg(10*n)+b1(m)+b2(8*n)+b3(2*m)+b4(8*n)+b5(m)); }// Construction du maillage // //if (sauve!="") savemesh(Th,sauve+".msh"); plot(Th,wait=wait); // Opérateur différentiels // real sqrt2=sqrt(2.); macro e(u) [dx(u[0]),(dy(u[0])+dx(u[1]))/sqrt2,dy(u[1])] // macro div(u) (dx(u[0])+dy(u[1]))// // Espaces et variables // macro u [u1,u2]// macro v [v1,v2]// macro g [g1,g2]// fespace Vh0(Th,P0),Vh2(Th,[P2,P2]); //Définition des espaces P0 et P2*P2 Vh0 h=h0; // Initialisation de l'épaissuer Vh2 u,v; // Déplacement et fonction test associée /////////////////////////// // Définition du système // /////////////////////////// problem elasticite(u,v) = int2d(Th)(2.*h*mu*e(u)'*e(v)+h*lambda*div(u)*div(v))-int1d(Th,Neumann)(g'*v) +on(Dirichlet,u1=0,u2=0); ////////////////////////////// // Calcul du volume initial // ////////////////////////////// real volume,volume0=int2d(Th)(h); // Volume courant et Volume initial ///////////////////////// // Compliance initiale // ///////////////////////// elasticite; real compliance=int1d(Th,Neumann)(g'*u);// Compliance de la plaque 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 // ////////////////////////////////////////////////// int inddico=0; while ((abs(1.-volume/volume0)) > 0.000001) {// Boucle Dichotomique 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 (volumevolume0) lagmax = lagrange; }// Boucle Dichotomique //cout<<"nbre iterations dichotomie: "<