load "Element_P4dc" load "Element_P4" load "Element_P3pnc" func real cc(real aa) {real a=aa; if(abs(a)<1e-10) a=0; return a;} int[int] ne1=[1,2,0]; int[int] ne2=[2,0,1]; // the ref triangle int[int] ll=[2,0,0,1]; mesh Th=square(1,1,flags=2,label=ll); Th = trunc(Th,x<0.5,label=0); //Th = movemesh(Th,[x+0.1*y,y-0.2*x]); //Th = movemesh(Th,[x*2,y*2]); mesh Thg = trunc(Th,1,split=2,label=-1); plot(Th, wait=1); int it0=0; fespace Lh(Th,P1); Lh[int] l(3); l[0][][0]=1; l[1][][1]=1; l[2][][2]=1; fespace Wh(Th,P4dc); fespace Vh(Th,P3pnc); Wh[int] mn(12); int k=0; real cc7 = 3*4*5*6*7; func bk = (l[0] - l[1]) * (l[1] - l[2]) * (l[2] - l[0]); func l0 =l[0]; func l1 =l[1]; func l2 =l[2]; mn[k++]= l0 * l0 * l0; mn[k++]= l1 * l1 * l1; mn[k++]= l2 * l2 * l2; //3 mn[k++]= l0 * l0 * l1; mn[k++]= l0 * l0 * l2; mn[k++]= l1 * l1 * l0; mn[k++]= l1 * l1 * l2; mn[k++]= l2 * l2 * l0; mn[k++]= l2 * l2 * l1; //6 mn[k++]= l0 * l1 * l2; //2 mn[k++]= bk * l0; mn[k++]= bk * l1 ; // P4 element /* l0 * l0 * l0, l1 * l1 * l1, l2 * l2 * l2, //3 l0 * l0 * l1, l0 * l0 * l2, l1 * l1 * l0, l1 * l1 * l2, l2 * l2 * l0, l2 * l2 * l1, //6 l0 * l1 * l2, //2 bk * l0, bk * l1 // P4 element */ /* for(int i=0; i<3;++i) { int i1 = (i+1)%3, i2=(i+3)%3; mn[k++]= l[i]; mn[k++]= l[i1]*l[i2]; mn[k++]= l[i]*l[i]*l[i]; } mn[k++]= l[0]*l[1]*l[2]; mn[k++]= l[0]*bk; mn[k++]= l[1]*bk; */ Vh u,v; real[int,int] CC(12,12),C1(12,12); for (int j=0;j " << i << " : " ; u=0; u[][i]=1; for(int k=0; k<3;k++) { int i1= (k+1)%3, i2=(k+2)%3; if ( i2 < i1) swap(i1,i2); cout << " " << cc(int1d(Th,k,qforder=9)(u*l[i1]/lenEdge)) << " " << cc(int1d(Th,k,qforder=9)(u*l[i2]/lenEdge)) << " " << cc(int1d(Th,k,qforder=9)(u*l[i2]*l[i1]/lenEdge)) ; } cout << " " << cc(int2d(Th,qforder=9)( l[0]*u/area)) << " " << cc(int2d(Th,qforder=9)( l[1]*u/area)) <<" " << cc(int2d(Th,qforder=9)( l[2]*u/area)) << endl; //plot(u,wait=1); v=u; for (int j=0;j 1e-6) cout << i << " " << e << " " << e1 << " " << e2 << " :: " << dxu []. linfty << " " << dyu []. linfty << endl; err+= ee; } cout << " err=" << err << endl; ; assert( err < 1e-6);