/////////////////////////////////////// // // // Optimisation de la forme // // d'un radiateur de refroidissement // // engl. Optimal radiator // // // fonction cout = compliance = // // énergie de dissipation thermique // // // // Ecole Polytechnique, MAP 562 // // Copyright G. Allaire, 2004 // // // // modified by T. Wick // // MAP 562 // // Winter 2016/2017 // // // // // /////////////////////////////////////// /////////////////////////////////////// // This example demonstrates the // results of lecture 8 (chapter7ctd.pdf), // page 39 - 41 /////////////////////////////////////// // As usually we first define // global variables, material properties // and problem data (forces) ////////////////////////////////////// int niter=100; // Number of total iterations int npen= 30; // Starting iterations with penalization int n=3; // Size of the mesh real lagrange=1.; // Lagrange multiplier for the volume real lagmin, lagmax ; // Lower and upper bounds of the Lagrange multiplier int inddico ; // Counter of iterations of the dichonomy problem real compliance; // Compliance real volume0; // Initial volume real volume,volume1; // Current volume // Some output commands ofstream objectif("radiateur.obj") ; string sauve="radiateur"; // security file string legende; // Légende pour les sorties graphiques real[int] vviso(21); for (int i=0;i<21;i++) vviso[i]=i*0.05 ; // Definition of pi and applied forces on the top boundary real pi = 4*atan(1) ; func g=1; // Temperature flux ///////////////////////////////// // Defining the two phases // and the characteristic function theta ///////////////////////////////// real alpha=0.01; real beta=1.0; func theta0=(x>0.1)*(x<0.9)*(y<0.9) ; //////////////////////////////// 1:Dirichlet condition // Definition of the domain // 2:Neumann condition //////////////////////////////// 3:Non-homogen. Neumann (temp. flux) mesh Th ; border bd(t=0,1) { x=1; y=t;label=2; }; // Right boundary border bg(t=1,0) { x=0; y=t;label=2; }; // Left boundary border bs(t=1,0) { x=t; y=1;label=3; }; // Top boundary border b1(t=0,0.1) { x=t; y=0;label=1; }; // Bottom boundary border b2(t=0.1,0.9) { x=t; y=0;label=2; }; // Inner boundary border b3(t=0.9,1) { x=t; y=0;label=1; }; // Inner boundary ////////////////////////////// // Construction of the mesh // ////////////////////////////// Th= buildmesh (bd(10*n)+bs(10*n)+bg(10*n)+b1(1*n)+b2(8*n)+b3(1*n)); plot(Th,wait=0,ps=sauve+".maillage.eps"); savemesh(Th,sauve+".msh"); ////////////////////////////////// // Defining the finite elements // ////////////////////////////////// fespace Vh0(Th,P0); fespace Vh2(Th,P2); // Solution variable: T // Test function: v Vh2 T,v; // P0-space (constants) for // the characteristic functions // and the density (similar to geometric optimization) Vh0 theta,density, thetb; //////////////////////////////////// // Proportions of the two phases // //////////////////////////////////// theta = theta0 ; thetb = 1-theta ; /////////////////////////////// // Implementation of the PDE // /////////////////////////////// problem conduction(T,v) = int2d(Th)( (alpha*theta+beta*(1-theta))*(dx(T)*dx(v)+dy(T)*dy(v)) ) -int1d(Th,3)(g*v) +on(1,T=0) ; //////////////////////////////// // In the following, we perform // initialization of the // volume and compliance //////////////////////////////// // Initial volume volume0=int2d(Th)(1-theta0); // Solve once the PDE in order // to get the first time a temperature. // This temperature is then used // to initialize the compliance. conduction; // Compute now the compliance compliance = int1d(Th,3)(g*T); objectif << compliance << endl ; cout<<"Initialization. Compliance = "< npen) { theta = ( 1 - cos(pi*theta) )/2 ; }; // Similar to geometric optimization, // we calculate the Lagrange multiplier // and make sure that it stays within // its bounds. volume=int2d(Th)(1-theta); volume1 = volume ; if (volume1 < volume0) { lagmin = lagrange ; while (volume < volume0) { lagrange = lagrange/2 ; theta = ( beta - sqrt(density*(beta-alpha)/(lagrange)) )/(beta-alpha); theta = min(theta,1) ; theta = max(theta,0) ; if (iter > npen) { theta = ( 1 - cos(pi*theta) )/2 ; }; volume=int2d(Th)(1-theta) ; }; lagmax = lagrange ; }; if (volume1 > volume0) { lagmax = lagrange ; while (volume > volume0) { lagrange = 2*lagrange ; theta = ( beta - sqrt(density*(beta-alpha)/(lagrange)) )/(beta-alpha); theta = min(theta,1) ; theta = max(theta,0) ; if (iter > npen) { theta = ( 1 - cos(pi*theta) )/2 ; }; volume=int2d(Th)(1-theta) ; }; lagmin = lagrange ; }; if (volume1 == volume0) { lagmin = lagrange ; lagmax = lagrange ; }; // Solving the dichonomy problem // to compute a Lagrange multiplier // in order to enforce the volume // constraint (see chapter5.pdf, page 32) inddico=0; while ((abs(1.-volume/volume0)) > 0.000001) { lagrange = (lagmax+lagmin)*0.5 ; theta = ( beta - sqrt(density*(beta-alpha)/(lagrange)) )/(beta-alpha); theta = min(theta,1) ; theta = max(theta,0) ; if (iter > npen) { theta = ( 1 - cos(pi*theta) )/2 ; }; volume=int2d(Th)(1-theta) ; inddico=inddico+1; if (volume < volume0) {lagmin = lagrange ;} ; if (volume > volume0) {lagmax = lagrange ;} ; }; //cout<<"Number of iteration of the dichonomy problem: "<