bilapmorley.edp

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

load "Morley"

int n=100,nn=n+10;
real[int] xx(nn),yy(nn);

mesh Th=square(40,40);  // mesh definition of $\Omega$
fespace Vh(Th,P2Morley);      // finite element space
// correct june 21 2007 (Thank to  Hui Zhang [e-mail] who find the error)
macro bilaplacien(u,v) ( dxx(u)*dxx(v)+dyy(u)*dyy(v)+2.*dxy(u)*dxy(v)) // fin macro 
real f=1;
Vh [u,ux,uy],[v,vx,vy];

solve bilap([u,ux,uy],[v,vx,vy]) =
    int2d(Th)(  bilaplacien(u,v) )
   - int2d(Th)(f*v)
   + on(1,2,3,4,u=0,ux=0,uy=0) 
; 
   
plot(u,cmm="u", wait=1,fill=1);
plot(ux,wait=1,cmm="u_x");
plot(uy,wait=1,cmm="u_y");


for (int i=0;i<=n;i++)
 {
   xx[i]=real(i)/n;
   yy[i]=u(0.5,real(i)/n); // value of uh at point (0.5, i/10.) 
 }
 plot([xx(0:n),yy(0:n)],wait=1);