-
Notifications
You must be signed in to change notification settings - Fork 197
/
Copy pathexample28.edp
28 lines (22 loc) · 972 Bytes
/
example28.edp
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
// example28.edp, from file chapt3/convects.edp
// With Discontinuous Galerkin
border C(t=0, 2*pi) { x=cos(t); y=sin(t); };
mesh Th = buildmesh(C(100));
fespace Vh(Th,P1dc);
Vh w, ccold, v1 = y, v2 = -x, cc = exp(-10*( (x-0.3)^2 + (y-0.3)^2) );
real u, alpha=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*( alpha*abs(n) - n/2 )*jump(cc) )
// - int1d(Th,C)( (n(u)<0)*abs(n(u) )*cc*w) // unused: cc=0 on d(Omega)^-
- int2d(Th)( ccold*w/dt );
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];
for (real t=0; t< 2*pi ; t+=dt){
ccold=cc;
Adual;
plot(cc, fill=true, dim=3, viso=viso,
cmm="t="+t + ", min=" + cc[].min + ", max=" + cc[].max);
};
plot(cc,wait=true,fill=false,value=true,viso=viso,
cmm="example28, min=" + cc[].min + ", max=" + cc[].max);