1// Mesh 2border ba(t=0, 1.0){x=t; y=0; label=1;} 3border bb(t=0, 0.5){x=1; y=t; label=2;} 4border bc(t=0, 0.5){x=1-t; y=0.5;label=3;} 5border bd(t=0.5, 1){x=0.5; y=t; label=4;} 6border be(t=0.5, 1){x=1-t; y=1; label=5;} 7border bf(t=0.0, 1){x=0; y=1-t;label=6;} 8mesh Th = buildmesh (ba(6) + bb(4) + bc(4) +bd(4) + be(4) + bf(6)); 9savemesh(Th, "th.msh"); 10 11// Fespace 12fespace Vh(Th, P2); 13fespace Nh(Th, P0); 14Vh u, v; 15Nh rho, logrho; 16 17// Variables 18real[int] viso(21); 19for (int i = 0; i < viso.n; i++) 20 viso[i]=10.^(+(i-16.)/2.); 21 22// Functions 23func f = (x-y); 24 25// Problem 26problem Probem1 (u, v, solver=CG, eps=1.0e-6) 27 = int2d(Th, qforder=5)( 28 u * v * 1.0e-10 29 + dx(u)*dx(v) 30 + dy(u)*dy(v) 31 ) 32 + int2d(Th, qforder=10)(-f*v) 33 ; 34 35varf indicator2 (uu, chiK) 36 = intalledges(Th)( 37 chiK*lenEdge*square(jump(N.x*dx(u) + N.y*dy(u))) 38 ) 39 + int2d(Th)( 40 chiK*square(hTriangle*(f+dxx(u)+dyy(u))) 41 ) 42 ; 43 44// Adaptation loop 45real error = 0.01; 46for (int i = 0; i < 4; i++) { 47 // Solve 48 Probem1; 49 cout << u[].min << " " << u[].max << endl; 50 51 // Plot 52 plot(u, wait=1); 53 54 // Indicator 55 rho[] = indicator2(0, Nh); 56 rho = sqrt(rho); 57 logrho = log10(rho); 58 cout << "rho min = " << rho[].min << ", rho max = " << rho[].max << endl; 59 60 // Plot 61 plot(rho, fill=1, wait=1, cmm="indicator density ", value=1, viso=viso, nbiso=viso.n); 62 plot(logrho, fill=1, wait=1, cmm="log 10 indicator density ", value=1, nbiso=10); 63 64 // Mesh adaptation 65 Th = adaptmesh(Th, [dx(u), dy(u)], err=error, anisomax=1); 66 plot(Th, wait=1); 67 68 // Interpolation 69 u = u; 70 rho = rho; 71 72 // Update error 73 error = error/2; 74} 75