/* Variational formulations !! 
Here is beam bending for Lame's equations */

r0 := 1.0;
r1 := 2.0;
border(1,0,2,5) begin x:=0; y:=2-t end;
border(2,0,10,30) begin x:=t; y:=0 end;
border(3,0,2,5) begin x:=10; y:=t end;
border(2,0,10,30) begin x:=10-t; y:=2 end;

buildmesh(800);
savemesh('lamemesh');
E= 21.5e10;
sigma := 0.29;
mu =E/(2*(1+sigma));
lambda = E*sigma/((1+sigma)*(1-2*sigma));
nu = lambda+sigma;

varsolve(u,v;w,s)
begin
     onbdy(1) u=0;
     onbdy(1) v=0;
     onbdy(3) dnu(v)=1e9;

     
e11 = dx(u);
e22 = dy(v);
e12 = 0.5*(dx(v)+dy(u));
e21 = e12;
dive = e11 + e22;
s11w=2*(lambda*dive+2*mu*e11)*dx(w);
s22s=2*(lambda*dive+2*mu+e22)*dy(s);
s12s = 2*mu*e12*(dy(w)+dx(s));
s21w = s12s;
a = s11w+s22s+s12s+s21w +0.1*s;
end : intt[a];
plot(u); plot(v);

save('u.dta',u); save('v.dta',v);
x = x + 0.1 * u;
y = y + 0.1 * v;
f = 0;
plot(f); /* trick to plot the new mesh */