Dear GetFEM++ users,
I am working on adaptive refinement of a linear elasticity problem using a posteriori error estimate.
At each adaptive refinement step, I need to compute some quantity of interest defined on a subdomain $\Omega$. Let's call this quantity of interest J(u) = \int_{\Omega} div(u) dx.
I use 'asm_generic' of GetFEM++ to compute this quantity. The python code looks like:
QoI = "(Div_u)*Test_t"
QoI_asm = gf.asm_generic(mim,1,QoI,OMEGA
,"u",False,mfu,md.variable('u')
,"t",True,mfef, np.zeros(mfef.nbdof()))
QoI_asm_elem = QoI_asm [ QoI_asm.size - mfef.nbdof() : QoI_asm.size ]
qoi = abs(np.sum(QoI_asm_elem))
with
mfef = gf.MeshFem(m,1)
mfef.set_fem(gf.Fem('FEM_PK_DISCONTINUOUS(2,{d})'.format(d=0)))
I compute then the relative error of the quantity of interest (J(u)-J(u_h))/J(u), with u is the solution computed from very fine mesh, u_h is the solution of the adaptive mesh.
What I got is that (J(u)-J(u_h))/J(u) does not converge well under mesh refinement. Secondly, J(u)-J(u_h) differs from J(u-u_h) for adaptive refinement case, which is not acceptable since J is linear. For the uniform refinement case, they are however identical.
I checked that the region of interest OMEGA is refined, and is updated correctly at each refinement step.
I do not know where the problem comes from. Do you think the generic assembly code I wrote is correct? Any hint would be very helpful and appreciated.
Best,
Phuoc