// // Copyright O. Pantz, S. Afaf, April 2014 // //==============Parameters Definition==========================================================================================// verbosity=0;int wait=0;string name="mast"; 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 18. === // int Nx=61, Ny=91; // Nb of nodes along x and y int NbUpWind=10; // Nb of iteration during HJ int kx=2.,ky=2.; // Nb of holes in the initial Shape along x and y real lag=15; // Lagrange multiplier for the weight real Lx=4, Ly=6; // Length/Height of the domain real epsilonV=1.e-4; // Regularization of the velocity int Niter=800; // 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,Sh; real deltax=Lx/(Nx-1),deltay=Ly/(Ny-1); {//initial mesh definition Sh=square(Nx-1,Ny-1,[Lx*x-Lx/2.,Ly*y-Ly/2.],flags=1); fespace Xh(Sh,P0); // for cut function of grip part that as to be removed from the mesh Xh ToRemove=0; for(int j=0;j<(Ny-1)*2/3;j++){ for(int i=0;i<(Nx-1)/4;i++) {ToRemove[](2*((Nx-1)*j+i))=1;ToRemove[](2*((Nx-1)*j+i)+1)=1;} for(int i=3*(Nx-1)/4;i