/////////////////////////////////////// // // // Optimization of the cross-section // // of an elastic bar under torsion // // // // Ecole Polytechnique, MAP 562 // // Copyright G. Allaire, 2004 // // // /////////////////////////////////////// ofstream compli("torsion.obj") ; string save="torsion"; // Prefix of backup files int niter=30; // Number of iterations int n=5; // Size of the mesh real lagrange=1.; // Lagrange multiplier for the volume constraint real lagmin, lagmax ; // Bounds for the Lagrange multiplier int inddico ; real compliance; // Compliance real volume0; // Initial volume real volume,volume1; // Volume of the current shape string caption; // Caption for the graphics func g=1; // Applied force (torsion momentum) real[int] vviso(21); for (int i=0;i<21;i++) vviso[i]=i*0.05 ; //////////////////////////////////////////// // Inverse shear moduli of the two phases // //////////////////////////////////////////// real alpha=0.5; real beta=1.; /////////////////////////////////////////////////////// // Initialization of the proportion of phase alpha // /////////////////////////////////////////////////////// real thetamoy=0.5; func theta0=thetamoy ; ////////////////////////////////// // Domain definition // // Dirichlet boundary condition // ////////////////////////////////// mesh Th ; border bd(t=0,1) { x=1; y=t;label=1; }; // right boundary of the domain border bg(t=1,0) { x=0; y=t;label=1; }; // left boundary of the domain border bs(t=1,0) { x=t; y=1;label=1; }; // upper boundary of the domain border bi(t=0,1) { x=t; y=0;label=1; }; // lower boundary of the domain /////////////////////// // Building the mesh // /////////////////////// Th= buildmesh (bd(10*n)+bs(10*n)+bg(10*n)+bi(10*n)); plot(Th,wait=0,ps=save+".mesh.eps"); ///////////////////////////////////////////// // Definition of the finite element spaces // ///////////////////////////////////////////// fespace Vh0(Th,P0); fespace Vh1(Th,P1); Vh1 u,v; Vh0 theta,density; //////////////////////////////// // Proportion of phase alpha // //////////////////////////////// theta = theta0 ; /////////////////////////////// // Torsion elasticity system // /////////////////////////////// problem torsion(u,v) = int2d(Th)( alpha*beta*(dx(u)*dx(v)+dy(u)*dy(v))/(beta*theta+alpha*(1-theta)) ) -int2d(Th)(g*v) +on(1,u=0) ; //////////////////// // Initial volume // //////////////////// volume0=int2d(Th)(theta); ///////////////////////// // Initial compliance // ///////////////////////// torsion; compliance=int2d(Th)(g*u); compli << -compliance << endl ; cout<<"Initialization. Compliance: "< volume0) { lagmax = lagrange ; while (volume > volume0) { lagrange = 2*lagrange ; theta = ( sqrt(density*alpha*beta*(beta-alpha)/lagrange) - alpha )/(beta-alpha); theta = min(theta,1) ; theta = max(theta,0) ; volume=int2d(Th)(theta) ; }; lagmin = lagrange ; }; // if (volume1 == volume0) { lagmin = lagrange ; lagmax = lagrange ; }; /////////////////////////////////////////////////////// // Dichotomy on the value of the Lagrange multiplier // /////////////////////////////////////////////////////// inddico=0; while ((abs(1.-volume/volume0)) > 0.000001) { lagrange = (lagmax+lagmin)*0.5 ; theta = ( sqrt(density*alpha*beta*(beta-alpha)/lagrange) - alpha )/(beta-alpha); theta = min(theta,1) ; theta = max(theta,0) ; volume=int2d(Th)(theta) ; inddico=inddico+1; if (volume < volume0) {lagmin = lagrange ;} ; if (volume > volume0) {lagmax = lagrange ;} ; }; //cout<<"number of iterations of the dichotomy: "<