Dear Freefem community.
I wonder if there is any way to implement the blue-green adaptive refinement strategy using Adaptmesh:
- Start with a coarse mesh Th
- Solve a discrete problem in the actual mesh Th
- Compute eta_K for each triangle in Th
- Use blue-green marking to refine each K’ in Th whose indicator eta_K’ satisfies eta_K’>=param*max{eta_K: K in Th}, where param is a parameter from 0 to 1.
- Define the resulting mesh as the actual mesh and go to step 2.
All the steps, except N°4, can be done in Freefem++. This is what I’m looking for. Any help is appreciated.
The following is a code inspired by the tutorials on the Poisson equation.
verbosity=0;
// Parameters
real hinit = 1; //initial mesh size
func f=(x-y);
// Mesh
border ba(t=0, 1.0){x=t; y=0; label=1;}
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));
// Fespace
fespace Vh(Th, P1); //for the mesh size and solution
Vh h = hinit; //the FE function for the mesh size
Vh u, v;
fespace Ph(Th, P0); //for the error indicator
//Build a mesh with the given mesh size hinit
Th = adaptmesh(Th, h, IsMetric=1, splitpbedge=1, nbvx=10000);
// plot(Th, wait=1);
// Problem
problem Poisson (u, v)
= int2d(Th, qforder=5)(
u*v*1.0e-10
+ dx(u)*dx(v)
+ dy(u)*dy(v)
)
- int2d(Th, qforder=5)(
f*v
)
;
varf indicator2 (unused, chiK)
= intalledges(Th)(
chiK*lenEdge*square(jump(N.x*dx(u) + N.y*dy(u)))
)
+ int2d(Th)(
chiK*square(hTriangle*(f + dxx(u) + dyy(u)))
)
;
// Mesh adaptation loop. Set to 0 for test.
for (int i = 0; i < 0; i++){
u = u;
// Solve
Poisson;
Ph etak;
etak[] = indicator2(0, Ph);
etak[] = sqrt(etak[]);
real maxetak=etak[].max;
real[int] mark(etak.n);
// The loop to mark elements
for (int j = 0; j < etak.n; j++){
if (etak[](j) >=0.5*maxetak){
mark[j]=1;
}
else{
mark[j]=0;
}
}
// Here we will refine the mesh using Adaptmesh and the mark array. How?
}