1load "Morley"
2macro DD(f,hx,hy) ( (f(x1+hx,y1+hy)-f(x1-hx,y1-hy))/(2*(hx+hy))) //
3macro dn(f) ( N.x*dx(f)+N.y*dy(f)) //
4macro dnn(f) ( N.x*f#2+N.y*f#3) //
5mesh Th=square(1,1,[10*(x+y/3),10*(y-x/3)]);
6
7real x1=0.7,y1=0.9, h=1e-6;
8int it1=Th(x1,y1).nuTriangle;
9
10func ff =  2*x*x + 3*y*y + 4.5*y*x + 5*x + 6*y + 7;
11func ffx = 4*x + 4.5*y + 5;
12func ffy = 6*y + 4.5*x + 6;
13func ffxx = 4.;
14func ffyy = 6.;
15func ffxy =4.5;
16
17fespace Vh(Th,P2Morley);
18fespace Eh(Th,P0edge);
19
20Eh  edges;
21
22Vh [a1,a2,a3],[b1,b2,b3],[c1,c2,c3];
23
24[a1,a2,a3] = [ ff,ffx,ffy];
25
26 cout << a1(0.6,0.2) << " == " << ff(0.6,0.2) << endl;
27 cout << a2(0.6,0.2) << " == " << ffx(0.6,0.2) << endl;
28  cout << a3(0.6,0.2) << " == " << ffy(0.6,0.2) << endl;
29  cout << dxx(a1)(0.6,0.2) << " xx== " << ffxx(0.6,0.2) << endl;
30  cout << dyy(a1)(0.6,0.2) << " yy== " << ffyy(0.6,0.2) << endl;
31  cout << dxy(a1)(0.6,0.2) << " xy== " << ffxy(0.6,0.2) << endl;
32  cout << dyx(a1)(0.6,0.2) << " yx== " << ffxy(0.6,0.2) << endl;
33 cout << a1(0.2,0.6) << " == " << ff(0.2,0.6) << endl;
34 cout << a2(0.2,0.6) << " == " << ffx(0.2,0.6) << endl;
35  cout << a3(0.2,0.6) << " == " << ffy(0.2,0.6) << endl;
36 cout << " 00 = " << int2d(Th)(square(a1-ff)) << endl;
37 cout << " 00 = " << int2d(Th)(square(a2-ffx)) << endl;
38 cout << " 00 = " << int2d(Th)(square(a3-ffy)) << endl;
39plot(a1,wait=1);
40
41Th=square(1,1,[10*(x+y/3),10*(y-x/3)]);
42
43varf vFlux([a],[e]) = intalledges(Th)( dn(a1)*e*(jump(real(nuTriangle))<= 0));
44varf vMean([a],[e]) = intalledges(Th)( (a1)*e*(jump(real(nuTriangle))<= 0)/lenEdge);
45
46
47for (int i=0;i<Vh.ndofK;++i)
48	cout << i << " -> " << Vh(0,i) << endl;
49for (int i=0;i<Vh.ndofK;++i)
50{
51  cout << " ***  node " << i << " of Traingle " << it1 << endl;
52  a1[]=0;
53  int j=Vh(it1,i);
54  a1[][j]=1;
55  edges[]=vFlux(0,Eh);
56  cout << "Flux  edges = " << edges[] << endl;
57  edges[]=vMean(0,Eh);
58  cout << " Mean   edges = " << edges[] << endl;
59
60  plot(a1, wait=1,cmm="w_"+i);
61  [b1,b2,b3]=[a1,a2,a3];
62
63  plot(a1,b1,cmm="w"+i, wait=1);
64
65  c1[] = a1[] - b1[];
66  cout << " int_1 " <<  int1d(Th,1) (dn(a1)) << endl;
67  cout << " int_3 " <<  int1d(Th,3) (dn(a1)) << endl;
68  cout << " int_2 " <<  int1d(Th,2) (dn(a1)) << endl;
69  cout << " int_4 " <<  int1d(Th,4) (dn(a1)) << endl;
70
71  cout << " int_1 " <<  int1d(Th,1) (dnn(a)) << endl;
72  cout << " int_3 " <<  int1d(Th,3) (dnn(a)) << endl;
73  cout << " int_2 " <<  int1d(Th,2) (dnn(a)) << endl;
74  cout << " int_4 " <<  int1d(Th,4) (dnn(a)) << endl;
75
76  cout << " ---------" << i << " " << c1[].max << " " << c1[].min << endl;
77  cout << " a = " << a1[] <<endl;
78  cout << " b = " << b1[] <<endl;
79
80  assert(c1[].max < 1e-5 && c1[].min > -1e-9);
81
82  cout << " dx(a1)(x1,y1) = " << dx(a1)(x1,y1) << " == " << DD(a1,h,0) << " == " << a2(x1,y1) << endl;
83  cout << " dy(a1)(x1,y1) = " << dy(a1)(x1,y1) << " == " << DD(a1,0,h) << " == " << a3(x1,y1) << endl;
84
85  cout << " dx(a2)(x1,y1) = " << dx(a2)(x1,y1) << " == " << DD(a2,h,0)  << " == " << dxx(a1) << endl;
86  cout << " dy(a2)(x1,y1) = " << dy(a2)(x1,y1) << " == " << DD(a2,0,h)  << " == " << dxy(a1) << endl;
87  cout << " dx(a3)(x1,y1) = " << dx(a3)(x1,y1) << " == " << DD(a3,h,0)  << " == " << dxy(a1) << endl;
88  cout << " dy(a3)(x1,y1) = " << dy(a3)(x1,y1) << " == " << DD(a3,0,h)  << " == " << dyy(a1) << endl;
89
90  assert( abs(dx(a1)(x1,y1)-DD(a1,h,0) ) < 1e-4);
91  assert( abs((a2)(x1,y1)-DD(a1,h,0) ) < 1e-4);
92  assert( abs((a3)(x1,y1)-DD(a1,0,h) ) < 1e-4);
93  assert( abs(dx(a2)(x1,y1)-DD(a2,h,0) ) < 1e-4);
94  assert( abs(dy(a1)(x1,y1)-DD(a1,0,h) ) < 1e-4);
95  assert( abs(dy(a2)(x1,y1)-DD(a2,0,h) ) < 1e-4);
96
97
98}
99
100