[Getfem-users] Neumann BC in Python interface

[Arkadiusz Witek]

Dear Getfrem-Users,

I'm trying to translate transient problem form the heat_equation.cc file  to Python interface. Everything works fine until setting Neumann boundary conditions. The code breaks without any comment while solving first iteration. Could you please help me with this problem?
Underneath I paste the code up to the line where the error occurs. Maybe there is same error in stating FEM problem because I'm not very familiar with FEM details. I run the code with IDLE under UBUNTU.

import getfem as gf
import numpy as np
#
m=gf.Mesh('cartesian',np.arange(0,5.,0.1),np.arange(0,5.,0.1))
mf=gf.MeshFem(m,1)
mf.set_fem(gf.Fem('FEM_QK(2,1)'))
mim=gf.MeshIm(m,gf.Integ('IM_EXACT_PARALLELEPIPED(2)'))
#
top=101
down=102
flst=m.outer_faces()
fnor=m.normal_of_faces(flst)
ttop=abs(fnor[1,:]-1)<1e-14
tdown=abs(fnor[1,:]+1)<1e-14
ftop=np.compress(ttop,flst,axis=1)
fdown=np.compress(tdown,flst,axis=1)
m.set_region(top,ftop)
m.set_region(down,fdown)
#
md=gf.Model('real')
md.add_fem_variable('u',mf,2)
#
k=5.
md.add_initialized_data('k',k)
e=md.add_generic_elliptic_brick(mim,'u','k')
#
bc_top_q=np.ones((mf.nbdof(),2),dtype=float)*10
md.add_fem_data('NeumannData1',mf,2,2)
md.set_variable('NeumannData1',bc_top_q,0)
md.set_variable('NeumannData1',bc_top_q,1)
qn=md.add_normal_source_term_brick(mim,'u','NeumannData1',top)
#
bc_down=mf.eval('25')
md.add_initialized_fem_data('DirichletData2',mf,bc_down)
md.add_Dirichlet_condition_with_multipliers(mim,'u',mf,down,'DirichletData2')
#
dt=5.
md.add_initialized_data('dt',dt)
md.add_basic_d_on_dt_brick(mim,'u','dt')
#
theta=.5
md.add_initialized_data('theta',theta)
md.add_theta_method_dispatcher(e,'theta')
md.add_theta_method_dispatcher(qn,'theta')
#
u0=np.ones(mf.nbdof(),dtype=float)*20.
md.set_variable('u',u0,1)
md.first_iter()
md.solve()
#
#.... time loop

Best Regards,
Arkadiusz Wlodarczyk