I have an issue while using the elasticity module. When I print the displacements ux and uy (called uu and vv on my code) they turn out to be the identical even though the results while ploting the displacement say otherwise. My code follows bellow
real E1 = 21.5;
real sigma = 0.29;
real gravity = 0;
// Mesh
int bottombeam = 2;
border a(t=2, 0){x=0; y=t ;label=1;} // left beam
border b(t=0, 10){x=t; y=0 ;label=bottombeam;} // bottom of beam
border c(t=0, 2){x=10; y=t ;label=1;} // rigth beam
border d(t=0, 10){x=10-t; y=2; label=3;} // top beam
real xp = 5,yp=1;
real[int,int] P=[[xp,yp]]; // add point (xp,yp) in a mesh ..
mesh th = buildmesh(b(20) + c(5) + d(20) + a(5),points=P,fixeborder=1);
plot(th,wait=1);
fespace Ph(th,P0); // constant par triangle
Ph E=E1;
// Fespace
fespace Vh(th, [P1, P1]);
Vh [uu, vv], [w, s];
// find midlle (xp,yp) dof
w[]=0:w[].n-1; // put dof number in array
int dof1=w(xp,yp)+0.5;
int dof2=s(xp,yp)+0.5;
assert( abs(dof1-w(xp,yp)) + abs(dof2-s(xp,yp)) < 1e-6);
cout << " dof "<< dof1 << " "<< dof2 << endl;
real[int] brhs(Vh.ndof); // the add rhs ..
brhs = 0;
//brhs[dof2]= 10;
brhs[dof1]= 10;
// Macro
real sqrt2 = sqrt(2.);
macro epsilon(u1, u2) [dx(u1), dy(u2), (dy(u1) + dx(u2))/sqrt2] // EOM
macro div(u, v) (dx(u) + dy(v)) // EOM
// Problem (with solve)
func mu = E/(2*(1 + sigma));
func lambda = E*sigma/((1 + sigma)*(1 - 2*sigma));
cout << "Gravity = " << gravity << endl;
solve Elasticity ([uu, vv], [w, s], solver=sparsesolver)
= int2d(th)(
lambda*div(w, s)*div(uu, vv)
+ 2.*mu*(epsilon(w, s)' * epsilon(uu, vv))
)
- int2d(th)(gravity*s)
+ brhs
+ on(1, uu=0, vv=0)
;
cout << "Max displacement = " << uu[].linfty << endl;
cout << "Max displacement = " << vv[].linfty << endl;
//Stiffness Matrix
varf bin(uu,vv) = int2d(th)(dx(uu) * dx(vv) + dy(uu) * dy(vv));
matrix K = bin(Vh, Vh);
real[int] U(K.n);
for(int i=0; i < K.n / 2; i++){
U(2*i) = uu[](i);
U(2*i+1) = vv[](i);
}
real[int] F(K.n);
F = K * U;
real[int] f1(K.n / 2);
real[int] f2(K.n / 2);
for(int i=0; i < K.n / 2; i++){
f1(i) = F(2*i);
f2(i) = F(2*i+1);
}
{ofstream fout("U.txt");
fout << U << endl;
}
{ofstream fout("f1.txt");
fout << f1 << endl;
}
{ofstream fout("f2.txt");
fout << f2 << endl;
}
cout << "Max f = " << f1.linfty << endl;
{ofstream fout("Stiffness_matrix.txt");
fout << K << endl;
}
//Plot
plot([uu, vv], wait=1);
//plot([uu, vv], wait=1, bb=[[-0.5, 2.5], [2.5, -0.5]]);
// Move mesh
mesh th1 = movemesh(th, [x+uu, y+vv]);
plot(th1, wait=1);
include "ffmatlib-Edited.idp"
savemesh(th,"export_mesh.msh");
ffSaveVh(th,Vh,"export_vh.txt");
ffSaveData(E, "efile.txt")
ffSaveData(uu, "uufile.txt")
ffSaveData(vv, "vvfile.txt")