// // Copyright O. Pantz, S. Afaf, April 2014 // //==============Parameters Definition==========================================================================================// verbosity=0;int wait=0;string name="bridge"; 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 real lag=30; // Lagrange multiplier for the weight // === Fig 14 === // int Nx=81,Ny=49; // Nb of nodes along x and y //int Nx=161,Ny=97; // Nb of nodes along x and y int NbUpWind=30; // Nb of iterations for HJ int kx=6,ky=6; //hole in PHI0. real Lx=2, Ly=1.2; // Length/Height of th domain real epsilonV=1.e-4; // Regularization of the velocity int Niter=200; // 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; {//initial mesh definition Th=square(Nx-1,Ny-1,[Lx*(x-0.5),Ly*(y-0.5)]); }//end of initial mesh definition // =============== Initialisation of the Level Set module ========================================================================= // load "../LevelSet/LevelSet" // chargement du module Level set real deltax=Lx/(Nx-1),deltay=Ly/(Ny-1); // taille des mailles selon x et y 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];g2[](Nx)=-1; // g = Ponctual applied load Vh u; // Displacement Wh phi,V,dV; // Level Set / velocity //==================================== Definition of the Penalization for Dirichlet BC =========================================// int[int] I=[1,2*(Nx-1)+1]; real[int] C=[1.e30,1.e30]; matrix D=[I,I,C]; // Dirichlet conditions D.resize(Vh.ndof,Vh.ndof); //=================== Initial shape PHI0=======================================================================================// real r=0.25; phi=cos(pi*kx*x/Lx+pi)*cos(pi*ky*y/Ly)-r^2; 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.]]); input="eps/"+name+"0.eps"; } plot(phi,viso=viso,fill=1,wait=1,cmm="phi initial",hsv=hsv); //=================== Rigidity bilinear form ============================================================================// 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))); //=====================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============================================// for(int i=0;i