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