Hello,
I have the following code written with the many help of Frédéric Hecht. Other than calculating the deformation, I need to calculate the internal forces that are applied to each node of the mesh.
How would you suggest to calculate the internal forces?
// Parameters
real E1 = 21.5;
real sigma = 0.29;
real gravity = -0.05;// 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;// 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 = Esigma/((1 + sigma)(1 - 2sigma));
//Not used cout << "Lambda = " << lamda << endl;
//Not used cout << "Mu = " << mu << endl;
cout << "Gravity = " << gravity << endl;
solve Elasticity ([uu, vv], [w, s], solver=sparsesolver)
= int2d(th)(
lambdadiv(w, s)div(uu, vv)
+ 2.mu(epsilon(w, s)’ * epsilon(uu, vv))
)
- int2d(th)(gravitys)
+ brhs
+ on(1, uu=0, vv=0)
;
cout << "Max displacement = " << uu.linfty << 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,Ph,“export_vh.txt”);
ffSaveData(E, “efile.txt”)