/usr/share/doc/freefem++/examples/examples++-chapt3/convects.edp is in freefem++-doc 3.19.1-1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 | // 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()(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)-n/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)-n/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);
|