An issue for h-adaptation in FreeFem++

[zhaog6]

Hello,

I’d like to implement h-adaptation on structured grids. I find the mesh quality doesn’t seem to be good when using splitmesh to refine mesh locally, because splitmesh(...) is not based on the partition principle of the longest edge of the grid elements. I wonder if there is any way to improve it.
For example:

macro grad(u) [dx(u), dy(u)] //

//----------------------+
//    Exact Solution    |
//----------------------+
func f = 1;
func g = 0;

//---------------+
//     Mesh      |
//---------------+
int nn = 1;
int[int] labs = [1, 2, 3, 4];
border bd1(t = 0, 1) {x = t; y = 0; label=1;}
border bd2(t = 0, 1) {x = 1; y = t/2.; label=1;}
border bd3(t = 0, 1) {x = 1-t/2.; y = 0.5; label=1;}
border bd4(t = 0, 1) {x = 0.5; y = t/2.+0.5; label=1;}
border bd5(t = 0, 1) {x = 0.5-t/2.; y = 1; label=1;}
border bd6(t = 0, 1) {x = 0; y = 1-t; label=1;}
mesh Th = buildmesh(bd1(2*nn) + bd2(nn) + bd3(nn) + bd4(nn) + bd5(nn) + bd6(2*nn));
/*
for (int it = 0; it < 10; ++it)
{
	Th = adaptmesh(Th, 1./nn, IsMetric=1, nbvx=1e+5);
}
*/
//mesh Th = square(nn, nn);
plot(Th, cmm="Initial Mesh", wait=1);

//----------------------------+
//    Finite Element Space    |
//----------------------------+
fespace Vh(Th, P1), Ph(Th, P0);
Vh uh;
Vh xx = x, yy = y;
Ph etaK;

//----------------------------------------------------+
//    Solve Poisson Equations by h-adaptive Methods   |
//----------------------------------------------------+
varf vA(u, v) = int2d(Th) (grad(u)'*grad(v)) + on(labs, u=g);
varf vl(unused, v) = int2d(Th) (f * v) + on(labs, unused=g);
varf vetaK(unused, chiK) = intalledges(Th) ( chiK * lenEdge * square(jump(grad(uh)'*[N.x, N.y])) )
	+ int2d(Th) ( chiK * square(hTriangle * (f + dxx(uh) + dyy(uh))) );
matrix A;
for (int it = 0; it < 100; ++it)
{
	etaK[] = vetaK(0, Ph);
	Th = splitmesh(Th, 1+(etaK>1.25e-4/2^it));
	uh = uh;
	etaK = etaK;
	A = vA(Vh, Vh);
	real[int] b = vl(0, Vh);
	uh[] = A^-1 * b;
	plot(Th, cmm="Refine Mesh", wait=1);
}