1load "Element_HCT"
2load "splitmesh3"
3load "Element_P3dc"
4
5macro DD(f,hx,hy) ( (f(x1+hx,y1+hy)-f(x1-hx,y1-hy))/(2*(hx+hy))) //
6macro dn(f) ( N.x*dx(f)+N.y*dy(f)) //
7macro dnn(f) ( N.x*f#2+N.y*f#3) //
8mesh Th=square(1,1,flags=3);//,[10*(x+y/3),10*(y-x/3)]);
9// Th=trunc(Th,x<0.5);
10mesh Th3=splitmesh3(trunc(Th,1,split=10));
11for(int i=0; i<3;++i)
12cout <<i << " "  <<" "<<  Th[0][i].x << " " << Th[0][i].y  << endl;
13real x1=0.7,y1=0.9, h=1e-6;
14int it1=Th(x1,y1).nuTriangle;
15
16func ff =  2*x*x + 3*y*y + 4.5*y*x + 5*x + 6*y + 7;
17func ffx = 4*x + 4.5*y + 5;
18func ffy = 6*y + 4.5*x + 6;
19func ffxx = 4.;
20func ffyy = 6.;
21func ffxy =4.5;
22
23fespace Vh(Th,HCT);
24fespace Wh(Th3,P3dc);
25fespace Eh(Th,P0edge);
26
27Eh  edges;
28
29Vh [a1,a2,a3],[b1,b2,b3],[c1,c2,c3];
30
31
32[a1,a2,a3] = [ ff,ffx,ffy];
33cout << a1(0,0) << endl;
34if(1)
35for(int k=0; k<12; ++k)
36{
37    a1[]=0;
38    a1[][Vh(0,k)]=1;
39    Wh aa=a1,dxaa=dx(aa),dyaa=dy(aa);
40    real err=int2d(Th)( square(aa-a1));
41    real errx=int2d(Th)( square(dx(aa)-a2));
42    real erry=int2d(Th)( square(dy(aa)-a3));
43    real errxx=int2d(Th)( square(dx(dxaa)-dxx(a1)));
44    cout << int2d(Th)( square(dxaa));
45    cout << int2d(Th)( square(dxx(a1)));
46  assert(errxx<1e-9);
47
48    real erryy=int2d(Th)( square(dy(dyaa)-dyy(a1)));
49    real errxy=int2d(Th)( square(dy(dxaa)-dxy(a1)));
50    cout << " err= " <<k << ":" << err << " /" << errx << " "<< erry << " / " << errxx << " "<< errxy << " "
51    << erryy << endl;
52    assert(int2d(Th)( square(dy(a2)-dyx(a1))) < 1e-9);
53    assert(int2d(Th)( square(dy(a3)-dyy(a1))) < 1e-9);
54    assert(int2d(Th)( square(dx(a2)-dxx(a1))) < 1e-9);
55    assert(int2d(Th)( square(dx(a3)-dxy(a1))) < 1e-9);
56    assert(errx<1e-9);
57    assert(erry<1e-9);
58   assert(errxx<1e-9);
59   assert(errxy<1e-9);
60   assert(erryy<1e-9);
61   [b1,b2,b3]=[a1,a2,a3];
62   cout << " a2, a3 :" << a2(0.5,0.5) << " " << a3(0.5,0.5) <<endl;
63   cout << k <<  " b1 = " << b1[] << endl;
64   b1[]-=a1[];
65   assert(b1[].linfty < 1e-6);
66}
67 cout << a1(0.6,0.2) << " == " << ff(0.6,0.2) << endl;
68 cout << a2(0.6,0.2) << " == " << ffx(0.6,0.2) << endl;
69  cout << a3(0.6,0.2) << " == " << ffy(0.6,0.2) << endl;
70  cout << dxx(a1)(0.6,0.2) << " xx== " << ffxx(0.6,0.2) << endl;
71  cout << dyy(a1)(0.6,0.2) << " yy== " << ffyy(0.6,0.2) << endl;
72  cout << dxy(a1)(0.6,0.2) << " xy== " << ffxy(0.6,0.2) << endl;
73  cout << dyx(a1)(0.6,0.2) << " yx== " << ffxy(0.6,0.2) << endl;
74 cout << a1(0.2,0.6) << " == " << ff(0.2,0.6) << endl;
75 cout << a2(0.2,0.6) << " == " << ffx(0.2,0.6) << endl;
76  cout << a3(0.2,0.6) << " == " << ffy(0.2,0.6) << endl;
77 cout << " 00 = " << int2d(Th)(square(a1-ff)) << endl;
78 cout << " 00 = " << int2d(Th)(square(a2-ffx)) << endl;
79 cout << " 00 = " << int2d(Th)(square(a3-ffy)) << endl;
80plot(a1,wait=1);
81
82//Th=square(1,1,[10*(x+y/3),10*(y-x/3)]);
83
84varf vFlux([a],[e]) = intalledges(Th)( dn(a1)*e*(jump(real(nuTriangle))<= 0));
85varf vMean([a],[e]) = intalledges(Th)( (a1)*e*(jump(real(nuTriangle))<= 0)/lenEdge);
86
87
88for (int i=0;i<Vh.ndofK;++i)
89	cout << i << " -> " << Vh(0,i) << endl;
90for (int i=0;i<Vh.ndofK;++i)
91{
92  cout << " ***  node " << i << " of Traingle " << it1 << endl;
93  a1[]=0;
94  int j=Vh(it1,i);
95  a1[][j]=1;
96  for (int k=0; k< 3; ++k)
97    cout <<" v at " << Th[it1][k] <<  " == " << a1(Th[it1][k].x, Th[it1][k].y) << " " << Th[it1][k].x<< " " <<  Th[it1][k].y <<endl;
98  edges[]=vFlux(0,Eh);
99  cout << "Flux  edges = " << edges[] << endl;
100  edges[]=vMean(0,Eh);
101  cout << " Mean   edges = " << edges[] << endl;
102
103  plot(a1, wait=1,cmm="w_"+i);
104  [b1,b2,b3]=[a1,a2,a3];
105
106  plot(a1,b1,cmm="w"+i, wait=1);
107
108  c1[] = a1[] - b1[];
109  cout << " int_1 " <<  int1d(Th,1) (dn(a1)) << endl;
110  cout << " int_3 " <<  int1d(Th,3) (dn(a1)) << endl;
111  cout << " int_2 " <<  int1d(Th,2) (dn(a1)) << endl;
112  cout << " int_4 " <<  int1d(Th,4) (dn(a1)) << endl;
113
114  cout << " int_1 " <<  int1d(Th,1) (dnn(a)) << endl;
115  cout << " int_3 " <<  int1d(Th,3) (dnn(a)) << endl;
116  cout << " int_2 " <<  int1d(Th,2) (dnn(a)) << endl;
117  cout << " int_4 " <<  int1d(Th,4) (dnn(a)) << endl;
118
119  cout << " ---------" << i << " " << c1[].max << " " << c1[].min << endl;
120  cout << " a = " << a1[] <<endl;
121  cout << " b = " << b1[] <<endl;
122
123  assert(c1[].max < 1e-5 && c1[].min > -1e-9);
124
125  cout << " dx(a1)(x1,y1) = " << dx(a1)(x1,y1) << " == " << DD(a1,h,0) << " == " << a2(x1,y1) << endl;
126  cout << " dy(a1)(x1,y1) = " << dy(a1)(x1,y1) << " == " << DD(a1,0,h) << " == " << a3(x1,y1) << endl;
127
128  cout << " dx(a2)(x1,y1) = " << dx(a2)(x1,y1) << " == " << DD(a2,h,0)  << " == " << dxx(a1)(x1,y1) << endl;
129  cout << " dy(a2)(x1,y1) = " << dy(a2)(x1,y1) << " == " << DD(a2,0,h)  << " == " << dxy(a1)(x1,y1) << endl;
130  cout << " dx(a3)(x1,y1) = " << dx(a3)(x1,y1) << " == " << DD(a3,h,0)  << " == " << dxy(a1)(x1,y1) << endl;
131  cout << " dy(a3)(x1,y1) = " << dy(a3)(x1,y1) << " == " << DD(a3,0,h)  << " == " << dyy(a1)(x1,y1) << endl;
132
133  assert( abs(dx(a1)(x1,y1)-DD(a1,h,0) ) < 1e-4);
134  assert( abs((a2)(x1,y1)-DD(a1,h,0) ) < 1e-4);
135  assert( abs((a3)(x1,y1)-DD(a1,0,h) ) < 1e-4);
136  assert( abs(dx(a2)(x1,y1)-DD(a2,h,0) ) < 1e-4);
137  assert( abs(dy(a1)(x1,y1)-DD(a1,0,h) ) < 1e-4);
138  assert( abs(dy(a2)(x1,y1)-DD(a2,0,h) ) < 1e-4);
139
140
141}
142
143