lame.edp

FreeFem++可以生成高质量的有限元网格。可以用于流体力学

// file lame.edd
mesh Th=square(10,10,[20*x,2*y-1]);
fespace Vh(Th,P2);
Vh u,v,uu,vv;
real sqrt2=sqrt(2.);
macro epsilon(u1,u2)  [dx(u1),dy(u2),(dy(u1)+dx(u2))/sqrt2] // EOM
// remark the $1/\sqrt2$ in term (dy(u1)+dx(u2)) is because 
// we want: 
//  epsilon(u1,u2)'* epsilon(v1,v2) == $ \varepsilon(\bm{u}): varepsilon(\bm{v})$
macro div(u,v) ( dx(u)+dy(v) ) // EOM


real E = 21e5, nu = 0.28, mu= E/(2*(1+nu)); 
real lambda = E*nu/((1+nu)*(1-2*nu)), f = -1; //  

solve lame([u,v],[uu,vv])= int2d(Th)(  
        lambda*div(u,v)*div(uu,vv)
        +2.*mu*( epsilon(u,v)'*epsilon(uu,vv) ) )	
        - int2d(Th)(f*vv) 
        + on(4,u=0,v=0);
real coef=100;
plot([u,v],wait=1,ps="lamevect.eps",coef=coef);

mesh th1 = movemesh(Th, [x+u*coef, y+v*coef]);
 plot(th1,wait=1,ps="lamedeform.eps");
real dxmin  = u[].min; 
real dymin  = v[].min; 

cout << " - dep.  max   x = "<< dxmin<< " y=" << dymin << endl; 
cout << "   dep.  (20,0)  = " << u(20,0) << " " << v(20,0) << endl;