// // Copyright O. Pantz, S. Afaf, April 2014 // //==============Parameters Definition==========================================================================================// verbosity=0;int wait=0;string name="Lbeam"; 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 int Nx=81, Ny=81; // Nb of nodes along x and y int NbUpWind=10; // Nb of iteration during HJ int kx=7,ky=8; // Nb of holes in the initial Shape along x and y real lag=100; // Lagrange multiplier for the weight real Lx=1, Ly=1; // Length/Height of the domain real epsilonV=1.e-4; // Regularization of the velocity int Niter=400; // 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); int Dirichlet=3; {//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=(Ny-1)/2;j<(Ny-1);j++){ for(int i=(Nx-1)/2;i