membrane.edp

来自「FreeFem++可以生成高质量的有限元网格。可以用于流体力学」· EDP 代码 · 共 29 行

// file membrane.edp
real theta=4.*pi/3.;
real a=2.,b=1.; // the length of the semimajor axis and  semiminor axis
func z=x;

border Gamma1(t=0,theta)    { x = a * cos(t); y = b*sin(t); }
border Gamma2(t=theta,2*pi) { x = a * cos(t); y = b*sin(t); }
mesh Th=buildmesh(Gamma1(100)+Gamma2(50));   // construction of mesh

fespace Vh(Th,P2); // P2 conforming triangular FEM
Vh phi,w, f=1;

solve Laplace(phi,w)=int2d(Th)(dx(phi)*dx(w) + dy(phi)*dy(w))
                - int2d(Th)(f*w) + on(Gamma1,phi=z);  // resolution of laplace equation
plot(phi,wait=true, ps="membrane.eps"); //Plot Th and v
plot(Th,wait=true, ps="membraneTh.eps"); //Plot Th and v

// to build a gnuplot data file
{ ofstream ff("graph.txt");
   for (int i=0;i<Th.nt;i++)
   { for (int j=0; j <3; j++)
       ff<<Th[i][j].x  << "    "<< Th[i][j].y<< "  "<<phi[][Vh(i,j)]<<endl;
    ff<<Th[i][0].x  << "    "<< Th[i][0].y<< "  "<<phi[][Vh(i,0)]<<endl
      <<endl<<endl;
   }
}

savemesh(Th,"Th.msh");