// // Copyright O. Pantz, S. Afaf, April 2014 // //==============Parameters Definition==========================================================================================// verbosity=0;int wait=0;string name="cantilever"; real E=1,nu=0.3; real mu=E/(2.*(1+nu)),lambda=E*nu/((1+nu)*(1-2.*nu)); // Lame coefficient real weak=0.001; // multiplicative coefficient for the "weak" material mimicking void // ==== Fig 2. === // int Nx=81, Ny=41; // Nb of nodes along x and y int NbUpWind=10; // Nb of iteration during HJ int kx=3,ky=4; // Nb of holes in the initial Shape along x and y real lag=100; // Lagrange multiplier for the weight /* // ==== Fig 3. ==== // int Nx=81, Ny=41; // Nb of nodes along x and y int NbUpWind=10; // Nb of iteration during HJ int kx=1,ky=4; // Nb of holes in the initial Shape along x and y real lag=100; // Lagrange multiplier for the weight // ==== Fig 11.a. ==== // int Nx=161, Ny=81; // Nb of nodes along x and y int NbUpWind=10; // Nb of iteration during HJ int kx=3,ky=4; // Nb of holes in the initial Shape along x and y real lag=150; // Lagrange multiplier for the weight // ==== Fig 11. ==== // int Nx=161, Ny=81; // Nb of nodes along x and y int NbUpWind=5; // Nb of iteration during HJ int kx=5,ky=6; // Nb of holes in the initial Shape along x and y real lag=150; // Lagrange multiplier for the weight */ real Lx=2, Ly=1; // Length/Height of the domain real epsilonV=1.e-4; // Regularization of the velocity int Niter=100; // Nb of iterations macro u [u1,u2]// macro v [v1,v2]// macro g [g1,g2]// macro grad(u) [dx(u),dy(u)] // macro div(u) (dx(u[0])+dy(u[1])) // macro e(u) [dx(u[0]),(dx(u[1])+dy(u[0]))/sqrt(2.),dy(u[1])] // // =============== Definition of the initial mesh ============================================================================// mesh Th;int Dirichlet; {//initial mesh definition Th=square(Nx-1,Ny-1,[Lx*x-Lx/2.,Ly*y-Ly/2.],flags=1); Dirichlet=4; }//end of initial mesh definition // =============== Initialisation of tje Level Set module ========================================================================= // load "../LevelSet/LevelSet" real deltax=Lx/(Nx-1),deltay=Ly/(Ny-1); include "../LevelSet/shift.edp"; //===================Ploting parameters=========================================================================================// real[int] viso=[-100,-1.e-4,1.e-4,100]; real[int] hsv2= [ 1.,0.,0, 1.,0.,0., 1.,0.,1, 1.,0.,1 ]; real[int] hsv= [ 1.,0.,0, 1.,0.,0., 1.,0.,1, 1.,0.,1, 1.,0.,1, 1.,0.,1 ]; //====================Space definition=========================================================================================// fespace Vh(Th,[P1,P1]),Wh(Th,P1); // Vh = displacement space; Wh = Level set space Vh g=[0,0]; g1[](Nx*Ny+Nx-1)=-1; // In case of P2 finite elements, this definition has to be changed Vh u; // Displacement Wh phi,V,dV; // Level Set / velocity //===================initial shape PHI0=======================================================================================// real r=0.25; phi=cos(pi*kx*x)*cos(pi*ky*y)-r^2; //initialization string input; if(name!="") { plot(phi,fill=1,viso=viso,hsv=hsv2,ps="eps/"+name+"0.eps",bb=[[-Lx/2,-Ly/2.],[Lx/2.,Ly/2.]],wait=1); input="eps/"+name+"0.eps"; } plot(phi,fill=1,wait=1,viso=viso,cmm="phi initial",hsv=hsv); //=====================problem OF displacement sigma==========================================================================// macro chi() ((phi<0)*(1.-weak)+weak) // macro X() ((phi<0)*1.) // varf a(u,v) = int2d(Th)(chi*(2.*mu*e(u)'*e(v)+lambda*div(u)*div(v))) + on(Dirichlet,u1=0,u2=0); //=====================problem oF V regularization (transport velocity)========================================================// problem regularizationV(V,dV)= int2d(Th)(X*V*dV+epsilonV*grad(V)'*grad(dV)) -int2d(Th)(X*(2.*mu*e(u)'*e(u)+lambda*div(u)^2-lag)*dV); int NbReg=5; // Nb of regularization iteration of function phi real dtReg=0.5*min(deltax,deltay); // Regularization step //==============================Displacement of the function LevelSet at velocity V============================================// u=[0,0]; for(int i=1;i