convects.edp

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

// file convects.edp
// Characteristics Galerkin
border C(t=0, 2*pi) { x=cos(t);  y=sin(t); };
mesh Th = buildmesh(C(100));
fespace Uh(Th,P1);
Uh cold, c = exp(-10*((x-0.3)^2 +(y-0.3)^2));

real dt = 0.17,t=0;
Uh u1 = y, u2 = -x;
for (int m=0; m<2*pi/dt ; m++) {
    t += dt;     cold=c;
    c=convect([u1,u2],-dt,cold);
    plot(c,cmm=" t="+t + ", min=" + c[].min + ", max=" +  c[].max);
}
// Now with Discontinuous Galerkin
fespace Vh(Th,P1dc);

Vh w, ccold, v1 = y, v2 = -x, cc = exp(-10*((x-0.3)^2 +(y-0.3)^2));
real u, al=0.5;  dt = 0.05;

macro n(u)(N.x*v1+N.y*v2) //

problem  Adual(cc,w) = int2d(Th)((cc/dt+(v1*dx(cc)+v2*dy(cc)))*w)
  + intalledges(Th)((1-nTonEdge)*w*(al*abs(n(u))-n(u)/2)*jump(cc))
//  - int1d(Th,C)((n(u)<0)*abs(n(u))*cc*w)  // unused because cc=0 on d(Omega)^-
  - int2d(Th)(ccold*w/dt);

for ( t=0; t< 2*pi ; t+=dt)
{
 ccold=cc; Adual;
 plot(cc,fill=1,cmm="t="+t + ", min=" + cc[].min + ", max=" +  cc[].max);
};
real [int] viso=[-0.1,0,0.5,0.1,0.5,0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,0.75,0.8,0.9,1];
plot(c,wait=1,fill=1,value=1,ps="convectCG.eps",viso=viso);
plot(cc,wait=1,fill=1,value=1,ps="convectDG.eps",viso=viso);

// the same DG very much faster
varf aadual(cc,w) = int2d(Th)((cc/dt+(v1*dx(cc)+v2*dy(cc)))*w)
        + intalledges(Th)((1-nTonEdge)*w*(al*abs(n(u))-n(u)/2)*jump(cc));
varf bbdual(ccold,w) =  - int2d(Th)(ccold*w/dt);
matrix  AA= aadual(Vh,Vh);
matrix BB = bbdual(Vh,Vh);
set (AA,init=t,solver=UMFPACK);
Vh rhs=0;
for ( t=0; t< 2*pi ; t+=dt)
{
  ccold=cc;
  rhs[] = BB* ccold[];
  cc[] = AA^-1*rhs[];
  plot(cc,fill=0,cmm="t="+t + ", min=" + cc[].min + ", max=" +  cc[].max);
};
plot(cc,wait=1,fill=1,value=1,ps="convectDG.eps",viso=viso);