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