1// Beam under it own weight
2
3// Parameters
4real E = 21.5;
5real sigma = 0.29;
6real gravity = -0.05;
7
8// Mesh
9int bottombeam = 2;
10border a(t=2, 0){x=0; y=t ;label=1;}	// left beam
11border b(t=0, 10){x=t; y=0 ;label=bottombeam;}	// bottom of beam
12border c(t=0, 2){x=10; y=t ;label=1;}	// rigth beam
13border d(t=0, 10){x=10-t; y=2; label=3;}	// top beam
14mesh th = buildmesh(b(20) + c(5) + d(20) + a(5));
15
16// Fespace
17fespace Vh(th, [P1, P1]);
18Vh [uu, vv], [w, s];
19
20// Macro
21real sqrt2 = sqrt(2.);
22macro epsilon(u1, u2) [dx(u1), dy(u2), (dy(u1) + dx(u2))/sqrt2] // EOM
23macro div(u, v) (dx(u) + dy(v)) // EOM
24
25// Problem (with solve)
26real mu = E/(2*(1 + sigma));
27real lambda = E*sigma/((1 + sigma)*(1 - 2*sigma));
28cout << "Lambda = " << lambda << endl;
29cout << "Mu = " << mu << endl;
30cout << "Gravity = " << gravity << endl;
31solve Elasticity ([uu, vv], [w, s], solver=sparsesolver)
32	= int2d(th)(
33		  lambda*div(w, s)*div(uu, vv)
34		+ 2.*mu*(epsilon(w, s)' * epsilon(uu, vv))
35	)
36	- int2d(th)(gravity*s)
37	+ on(1, uu=0, vv=0)
38	;
39cout << "Max displacement = " << uu[].linfty << endl;
40
41// Plot
42plot([uu, vv], wait=1);
43plot([uu, vv], wait=1, bb=[[-0.5, 2.5], [2.5, -0.5]]);
44
45// Move mesh
46mesh th1 = movemesh(th, [x+uu, y+vv]);
47plot(th1, wait=1);
48