Re: [Getfem-users] Generic assembly of GetFEM++

[Yves Renard]

Dear Phuoc,

If you just want to compute J(u) = \int_{\Omega} div(u) dx, then I would say 
that the more straigthforward computation is

gf.asm_generic(mim,0,"Div_u",OMEGA
                                 ,"u",False,mfu,md.variable('u')
                                 ,"t",True,mfef, np.zeros(mfef.nbdof()))

However, if mfef is a Lagrange finite element, what you wrote will also do the 
job, may be except that it sums all the components of course.

Best regards,

Yves.

----- Original Message -----
From: "Huu Phuoc BUI" <address@hidden>
To: address@hidden
Sent: Wednesday, April 18, 2018 1:03:59 PM
Subject: [Getfem-users] Generic assembly of GetFEM++

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