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
53
54
55
56
57
58
59
60
61
|
load "Element_P4dc"
macro DD(f,hx,hy) ( (f(x1+hx,y1+hy)-f(x1-hx,y1-hy))/(2*(hx+hy))) //
macro DD2(f,hx,hy) ( (-2*f(x1,y1)+f(x1+hx,y1+hy)+f(x1-hx,y1-hy))/(square(hx+hy))) //
macro dn(f) ( N.x*dx(f)+N.y*dy(f)) //
macro dnn(f) ( N.x*f#2+N.y*f#3) //
mesh Th=square(1,1,[10*(x+y/3),10*(y-x/3)]);
real x1=0.7,y1=0.9, h=1e-4;
int it1=Th(x1,y1).nuTriangle;
fespace Vh(Th,P4dc);
fespace Eh(Th,P0edge);
Eh edges;
Vh a1,b1,c1;
varf vFlux([a],[e]) = intalledges(Th)( dn(a1)*e*(jump(real(nuTriangle))<= 0));
varf vMean([a],[e]) = intalledges(Th)( (a1)*e*(jump(real(nuTriangle))<= 0)/lenEdge);
for (int i=0;i<Vh.ndofK;++i)
cout << i << " -> " << Vh(0,i) << endl;
for (int i=0;i<Vh.ndofK;++i)
{
cout << " *** node " << i << " of Traingle " << it1 << endl;
a1[]=0;
int j=Vh(it1,i);
a1[][j]=1;
edges[]=vFlux(0,Eh);
cout << "Flux edges = " << edges[] << endl;
edges[]=vMean(0,Eh);
cout << " Mean edges = " << edges[] << endl;
plot(a1, wait=1,cmm="w_"+i);
b1=a1;
plot(a1,b1,cmm="w"+i, wait=1);
c1[] = a1[] - b1[];
cout << " ---------" << i << " " << c1[].max << " " << c1[].min << endl;
cout << " a = " << a1[] <<endl;
cout << " b = " << b1[] <<endl;
assert(c1[].max < 1e-5 && c1[].min > -1e-9);
cout << " dx(a1)(x1,y1) = " << dx(a1)(x1,y1) << " == " << DD(a1,h,0) << endl;
cout << " dy(a1)(x1,y1) = " << dy(a1)(x1,y1) << " == " << DD(a1,0,h) << endl;
cout << " dxx(a1)(x1,y1) = " << dxx(a1)(x1,y1) << " == " << DD2(a1,h,0) << endl;
cout << " dyy(a1)(x1,y1) = " << dyy(a1)(x1,y1) << " == " << DD2(a1,0,h) << endl;
assert( abs(dx(a1)(x1,y1)-DD(a1,h,0) ) < 1e-4);
assert( abs(dxx(a1)(x1,y1)-DD2(a1,h,0) ) < 1e-4);
assert( abs(dy(a1)(x1,y1)-DD(a1,0,h) ) < 1e-4);
assert( abs(dyy(a1)(x1,y1)-DD2(a1,0,h) ) < 1e-4);
}
|
