[Getfem-users] neumann boundary condition

[Alice Nicolas]

Dear all,

I am trying try to compute Laplace(u)=u in 1D, with u(0)=0, du/dx(0)=1. 
The solution should be u=x^3/6+x. I fail in programming the Neumann 
boundary condition. It seems that it is not taken into account, and I 
get as a solution something that is really well fitted by -x^3/6+x/2. It 
seems that I do not use correctly the source term brick. May I ask your 
help to debug this stubborn mistake?

Here is the part of my code that seems to be wrong:
#I have first defined the mesh on [0,1] and assigned to 0 the boundary 
region "left"
# create a MeshFem for u and rhs fields of dimension 1 (i.e. a scalar field)
mfp   = gf.MeshFem(m, 1)
mfrhs = gf.MeshFem(m, 1)
mfp.set_fem(gf.Fem('FEM_PK(1,2)'))
mfrhs.set_fem(gf.Fem('FEM_PK(1,2)'))
mim = gf.MeshIm(m, gf.Integ('IM_GAUSS1D(4)'))

# interpolate the source term
F1 = mfrhs.eval('x[0]')
#interpolate the boundary conditions
#Neumann
F2=mfrhs.eval('1.')
#Dirichlet
F3=mfrhs.eval('0.')

# model
md = gf.Model('real')
md.add_fem_variable('p', mfp)
md.add_Laplacian_brick(mim, 'p')

# volumic source term
md.add_initialized_fem_data('VolumicData', mfrhs, F1)
md.add_source_term_brick(mim, 'p', 'VolumicData')
#Neumann condition on the left dp/dx=1
md.add_initialized_fem_data('NeumannData', mfrhs, F2)
md.add_source_term_brick(mim,'p','NeumannData',left)
# Dirichlet condition on the left.
md.add_initialized_fem_data("DirichletData", mfrhs, F3)
md.add_Dirichlet_condition_with_penalization(mim, 'p', 1e20, 
left,'DirichletData')

# assembly of the linear system and solve.
md.solve()
P = md.variable('p')
mfp.export_to_pos('p.pos',P,'')

Thank you very much for your help!
Best regards,
Alice Nicolas