border ba(t=0,1.0){x=t; y=0; label=1;};// comment border bb(t=0,0.5){x=1; y=t; label=2;}; border bc(t=0,0.5){x=1-t; y=0.5;label=3;}; border bd(t=0.5,1){x=0.5; y=t; label=4;}; border be(t=0.5,1){x=1-t; y=1; label=5;}; border bf(t=0.0,1){x=0; y=1-t;label=6;}; mesh Th = buildmesh (ba(6) + bb(4) + bc(4) +bd(4) + be(4) + bf(6)); int n6= 600,n3=n6/2; mesh Tf = buildmesh (ba(n6) + bb(n3) + bc(n3) +bd(n3) + be(n3) + bf(n6)); fespace Vf(Tf,P1); savemesh(Th,"th.msh"); fespace Vh(Th,P1); fespace Nh(Th,P0); Vh u,v; Vf uf,vf; real error=0.01; func f=(x-y); problem Probem1(u,v,solver=CG,eps=1.0e-6) = int2d(Th,qforder=2)( u*v*1.0e-10+ dx(u)*dx(v) + dy(u)*dy(v)) + int2d(Th,qforder=2)( -f*v); solve Probemf(uf,vf,solver=CG,eps=1.0e-6) = int2d(Tf,qforder=2)( uf*vf*1.0e-10+ dx(uf)*dx(vf) + dy(uf)*dy(vf)) + int2d(Tf,qforder=2)( -f*vf); int i; Nh eta; varf indicator2(uu,eta) = intalledges(Th)(eta*lenEdge*square(jump(N.x*dx(u)+N.y*dy(u)))) +int2d(Th)(eta*square(f*hTriangle)); for (i=0;i< 4;i++) { Probem1; real errL2 = sqrt(int2d(Tf)(sqr(uf-u))); real errH1 = sqrt(int2d(Tf)(sqr(dx(uf)-dx(u))+sqr(dy(uf)-dy(u)))); cout << u[].min << " " << u[].max << " Err= " << errL2 << " " << errH1 << endl; plot(u,wait=1); cout << " indicator2 " << endl; eta[] = indicator2(0,Nh); eta=sqrt(eta); cout << eta[].min << " " << eta[].max << endl; plot(eta,fill=1,wait=1); Th=adaptmesh(Th,u,err=error,anisomax=1); plot(Th,wait=1); u=u; eta=eta; error = error/2; } ;