////////////////////////////////////// // // // Thickness or sizing optimization // // of an elastic plate // // // // Ecole Polytechnique, MAP 562 // // Copyright G. Allaire, 2004 // // modified by T. Wick // // MAP562: Winter 2017 // // // ////////////////////////////////////// // Program: // 1. Compliance minimization: alternate minimization in sigma and h // without using the adjoint explicitly thanks to // the self-adjointness of the equations // (see lecture notes and Homework 4) // 2. TODO: Displacement tracking // 2.0 Define geometry (done) // 2.1 Re-define cost functional: matching displacement // 2.2 Define adjoint problem // 2.3 Update functional derivative via // J'(h)(\delta h) = \int \delta h \sigma \nabla p // or simply in strong form: J'(h) = \sigma \nabla p // 2.4 Algorithm // h_{n+1} = Proj (h_n - \mu J'(h_n)) // 2.5 Take care of the two projections: // 1. hmax and hmin // 2. volume constraint // // Idea for the development: keep compliance // cost functional and change first only the algorithm // from the optimality criteria method to // a gradient method (thus step 2.3 - 2.5) // and check if the result is the same as before. // If yes, change now to the new cost functional. // Runtime parameters and definitions ofstream compli("plate.obj") ; string save="plate"; // Prefix of backup files int niter=30; // Number of iterations real lagrange=1.; // Lagrange multiplier for the volume constraint real lagmin, lagmax ; // Bounds for the Lagrange multiplier real J; // Cost functional, here the compliance real volume0; // Initial volume real volume,volume1; // Volume of the current shape string caption; // Caption for the graphics real E=100; // Young modulus real nu=0.3; // Poisson coefficient (between -1 and 1/2) real[int] vviso(21); for (int i=0;i<21;i++) vviso[i]=i*0.05 ; // Problem data (traction force on a section of the boundary) func g1=0; // Applied forces func g2=100; macro g [g1,g2]// ////////////////////////////////////// // Computation of Lame coefficients // ////////////////////////////////////// real lambda=E*nu/((1.+nu)*(1.-2.*nu)); real mu=E/(2.*(1.+nu)); //////////////////////////////////////////////////// // Lower and upper bounds for the plate thickness // //////////////////////////////////////////////////// real hmin=0.1; real hmax=1.0; real hmoy=0.5; func h0=hmoy ; /////////////////////// 1:Dirichlet boundary condition // Domain definition // 2:Neumann boundary or free boundary condition /////////////////////// 3:Non-homogeneous Neumann or applied load // 4:mixed conditions for // diff. comp. on top boundary mesh Th, Dh; border bd(t=0,1) { x=+1; y=t;label=2; }; // right boundary of the domain border bg1(t=0.8,0.7) { x=-1; y=t;label=2; }; // left boundary of the domain border bg2(t=0.7,0) { x=-1; y=t;label=2; }; // left boundary of the domain border bs(t=1,-0.5) { x=t; y=1;label=4; }; // top boundary of the domain border b1(t=-1,-0.9) { x=t; y=0;label=2; }; // lower boundary of the domain 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=2; }; border b4(t=0.1,0.9) { x=t; y=0;label=2; }; border b5(t=0.9,1) { x=t; y=0;label=3; }; // Corner top left boundary border corn1(t=1,0.8) { x=-0.5; y=t; label=2; }; border corn2(t=-0.5,-1) { x=t; y=0.8; label=2; }; // Interior borders border i1(t=0.8,0.7) { x=-0.5; y=t; }; border i2(t=-1.0,-0.5) { x=t; y=0.7; }; // Cylinder (homogeneous Dirichlet boundary conditions) border cc(t=0,2*pi) {x=0.5 + 0.1 * cos(t); y= 0.5 + 0.1* sin(t); label=1; } /////////////////////// // Building the mesh // /////////////////////// int n = 1 ; int m = 1 ; Th= buildmesh (bd(10*n)+bs(20*n)+corn1(10*n)+corn2(10*n)+bg1(m)+bg2(10*n)+b1(m)+b2(8*n)+b3(2*m)+b4(8*n)+b5(m)+i1(m)+i2(m)+cc(-50)); plot(Th,wait=0,ps=save+".mesh.eps"); // Obtaining material ids for the two domains int Omegabig = Th(-0.72,0.4).region; int Omegacorner = Th(-0.75,0.85).region; ///////////////////////////////////////////// // Defining macros for short hand notation // ///////////////////////////////////////////// macro u [u1,u2] // macro phi [phi1,phi2]// macro grad(u) [dx(u), dx(u)] // macro e(u) [dx(u[0]),dy(u[1]),(dx(u[1])+dy(u[0]))/sqrt(2.0)]// macro div(u) (dx(u[0])+dy(u[1]))// ///////////////////////////////////////////// // Definition of the finite element spaces // ///////////////////////////////////////////// // Vector-valued displacements fespace Vh2(Th,[P2,P2]); // The thickness h fespace Vh0(Th,P0); Vh2 u,phi; Vh0 h,hold,density; h = h0 ; /////////////////////// // Elasticity system // ////////////////////// // Define state equation problem elasticity(u,phi) = int2d(Th)(2.0 * h * mu * e(u)' * e(phi) + h * lambda * div(u) * div(phi)) -int1d(Th,3)(g' * phi) +on(4,u2=0) +on(1,u1=0,u2=0) ; //////////////////// // Initial volume // //////////////////// volume0=int2d(Th)(h); // Initial value of the cost functional elasticity; J=int1d(Th,3)(g' * u); compli << J << endl ; cout<<"Initialization. Compliance: "< volume0) { lagmax = lagrange ; while (volume > volume0) { lagrange = 2*lagrange ; h = sqrt(density*hold/lagrange) ; h = min(h,hmax) ; h = max(h,hmin) ; volume=int2d(Th)(h) ; }; lagmin = lagrange ; }; // if (volume1 == volume0) { lagmin = lagrange ; lagmax = lagrange ; }; /////////////////////////////////////////////////////// // Dichotomy on the value of the Lagrange multiplier // /////////////////////////////////////////////////////// int inddico=0; while ((abs(1.-volume/volume0)) > 0.000001) { lagrange = (lagmax+lagmin)*0.5 ; h = sqrt(density*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<<"number of iterations of the dichotomy: "<