// simulation of wafer heating by a laser./**************************************************************From:  Fabian DortuMaterials and Components Analysis group (MCA)SPDT divisionIMEC, Kapeldreef 75, B-3001 Leuven, Belgiumtel: +32/16 28 8774e-mail: Fabian.Dortu@imec.be <mailto:Fabian.Dortu@imec.be> **************************************************************/// The poisson equation (here the heat equation) is solved in cylindrical coordinates://   $ \Delta u = 1/r \partial_r(r \partial_r u) + \partial_{zz} u +  1/r^2 \partial_{\theta^2} u $// so the variationnal formulation is after integer on the 3d cylinder // rotation invariant given $\partial_{\theta} u =0$//  $ - \Delta u = f $//  on the 1/2 plan $ \theta =0$ //  $ \int r ( \partial_r u \partial_r v + \partial_z u \partial_z v ) = \int r f v  - \int_{\Gamma} r \partial_n u v $  //  the radial coordinate 'r' is called 'x'.//  the depth coordinate 'z' is called 'y'.// The wafer stuck is a rectangle with upper left corner at (0,0)//  and lower right corner at (radius,-thickness)// The laser beam it the surface from top to bottom at the center (0,0). // The units assumed for the distance is the micrometer (um)//***********************//*** User parameters ***//***********************real radius=1000;   // the wafer radius (um)real thick=600;     // the wafer thickness (um)real thC=1.5e-4;    // the thermal conductivity (W/K/um)real beamr=2;     // the laser beam radius (in the sens of a gaussian)real refl=0.55;     // reflection coefficient Air/Siliconreal Eg=1.12;       // Band gap of Silicon (eV)real En=1.49;       // Energy of laser (eV) - must be greater than Eg.real I0=1e-2;       // Laser beam intensity (W/um^2)real alpha=0.2;     // absorption coefficient (um^-1)func f=I0 * (1-refl) * (En-Eg)/En * alpha * exp(alpha*y) * exp(-x*x/beamr/beamr);                    // the heat generation function: a laser beam of radius 'beamr'                    // lateral shape is gaussian.                    // in depth shape is decreasing exponential because of absorption.real g=300;          // temperature at right and bottom surfaces (in Kelvin)//********************************//*** Geometry/Mesh definition ***//********************************//--- less basic mesh ---border bottom1(t=0,radius/5) {x=t; y=-thick; label=1;}border bottom2(t=radius/5,radius) {x=t; y=-thick; label=2;}border right1(t=-thick,0) {x=radius; y=t; label=3;};  // 'right' actually means external surfaceborder top1(t=radius,radius/5*3) {x=t; y=0; label=4;};border top2(t=radius/5*3,radius/5) {x=t; y=0; label=5;}border top3(t=radius/5,0) {x=t; y=0; label=6;};border left1(t=0,-thick/5) {x=0; y=t; label=7;}; // 'left' actually means center of cylinderborder left2(t=-thick/5,-thick) {x=0; y=t; label=8;}; // 'left' actually means center of cylindermesh Th = buildmesh ( bottom1(20) + bottom2(10) + right1(20) + top1(10) + top2(20) + top3(50) + left1(50) + left2(25) );plot(Th,wait=1,fill=1,ps="wafer_mesh2.eps");//savemesh(Th,"wafer_stuck.msh");//**************************//*** Problem definition ***//**************************fespace Vh(Th,P1);Vh u,v,zero;u=0;     // The temperature variable.v=0;     // The weight function.zero=0;  // used to set initial conditionint i=0; // variable used for mesh reconstructionreal error=0.1, coef=0.1^(1./5.);//--- Variational Form ----problem Problem1(u,v,solver=CG,init=i,eps=1.0e-8) =int2d(Th) ( x*dx(u)*dx(v) + x*dy(u)*dy(v) )  + int2d(Th) ( -v*x*f/thC )    // the source term = laser heating   + on( 1,2,3,u=g ) ;             // fixed temperature at bottom and right surface.//*************//*** Solve ***//*************real cpu=clock();cout << "***USER*** " << "Begin solve/adapt iterations" << endl;for (i=0;i<10;i++){  cout << "***USER*** " << "Iteration number: " << i << endl;  real d = clock();  Problem1;  plot(u,wait=1);  Th=adaptmesh(Th,u,inquire=1,err=error);  cout << " CPU = " << clock()-d << endl;  error = error * coef;  } ;plot(Th,wait=1,bb=[[0,0],[-10,-10]],nbiso=20,ps="adaptedmesh.eps"); // plot the last adapted meshplot(u,wait=1,bb=[[0,0],[-10,-10]],nbiso=20,ps="temperature.eps"); // plot the solution//Plot a cut section at r=0real[int] xx(10),yy(10);for (i=0;i<10;i++){  x=0.; y=i/10.; // this line is used by the next one   xx[i]=i;  yy[i]=u; // value of u at point (0. , i/10.)}plot([xx,yy],ps="likegnu.eps",wait=true);cout << " CPU = " << clock()-cpu << endl;