mesh Th;int Dirichlet; real L=20,H=1; { int dn=5; Th=square(L*dn,H*dn,[L*x,H*y]); Dirichlet=4; } real sqrt2=sqrt(2.); macro e(u) [dx(u[0]),dy(u[1]),(dx(u[1])+dy(u[0]))/sqrt2]// macro div(u) (dx(u[0])+dy(u[1]))// fespace Wh(Th,[P1,P1]); macro u [u1,u2]// macro up [up1,up2]// macro uP [uP1,uP2]// macro v [v1,v2]// Wh u,v,up,uP; real t=0,dt=10.,T=2000; real mu=1,lambda=1,rho=1; Wh [U1,U2]=[-x,0]; // initial displacement Wh [V1,V2]=[0,0.1*sin(2.*x/L)]; // initial velocity //Wh [V1,V2]=[0,0]; // initial velocity real dt2=square(dt); problem elasticity(u,v,init=t)=int2d(Th)(rho*u'*v/dt2+2.*mu*e(u)'*e(v)+lambda*div(u)*div(v)) +int2d(Th)(rho*(-2*up+uP)'*v/dt2) +on(Dirichlet,u1=0,u2=0); // initialisation uP=[U1,U2]; up=[U1+dt*V1,U2+dt*V2]; real exa=0.2; int iter=0; for(t=0;t