Unknown in space R defined by an integral

Hello,

I solve a two-phase 1D Poiseuille flow written in cylindrical coordinates in a tube with radius R_t. The regions are [0,\lambda] (fluid1) & [\lambda,R_t] (fluid2). \lambda is equal to 0.8 and R_t is equal to 1. The unknowns of the problem are the velocity (v) and the shear stress (trz). The two regions have different viscosity only. The sample of my code is the following:



# subdomains
dx     = Measure("dx")(subdomain_data=subdomains)
dx_f1  = dx(fluid1,subdomain_data=subdomains)
dx_f2  = dx(fluid2,subdomain_data=subdomains)

phi = []

FE1    = FiniteElement('CG',mesh.ufl_cell(),1)
FER    = FiniteElement('R' ,mesh.ufl_cell(),0)


FEM    = FunctionSpace(mesh,MixedElement([FE1,FE1]))  **Line 1**
phi    = TestFunctions(FEM)

# Solution Vector
uknown     = Function(FEM)
uknown0    = Function(FEM)
v,trz      = split(uknown) **Line 2**


#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
# 							Initial Conditions
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-


v_0   = Expression('0.0'  , degree=0)
trz_0 = Expression('0.1'  , degree=0)

v0     = interpolate(v_0   , FEM.sub(0).collapse())
trz0   = interpolate(trz_0 , FEM.sub(1).collapse())


#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
# 								BCs
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-


def on_left ( x, on_boundary ):
	tol = 1E-14
	return on_boundary and near(x[0], xl, tol)

def on_right ( x, on_boundary ):
	tol = 1E-14
	return on_boundary and near(x[0], xr, tol)


ns0      = Constant(0.0)

vBC      = DirichletBC(FEM.sub(0), ns0, on_right  )
trzBC    = DirichletBC(FEM.sub(1), ns0, on_left   )


bcs    = [vBC,trzBC]
bcs_   = [vBC,trzBC]


#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
#             Weak forms
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+


# Fluid 1
Res1   = ( Re*(1./dt)*(v-v0)*phi[0] - phi[0] - (1./x[0])*trz*phi[0] + trz*phi[0].dx(0) ) * dx_f1
Res2   = (trz - 2.*v.dx(0)) * phi[1] * dx_f1

# Fluid 2
Res3   = ( Re*(1./dt)*(v-v0)*phi[0] - phi[0] - (1./x[0])*trz*phi[0] + trz*phi[0].dx(0) ) * dx_f2
Res4   = (trz - v.dx(0)) * phi[1] * dx_f2

**Line 3** 

Everything works fine and I get the right solution. Now I want to solve for the extra unknown h in the R space (probably with a Lagrangian multiplier) the equation of which is defined as:

h \int_{0}^{R_t} u r dr = \int_{0}^{\lambda} u r dr

So I tried to add this in the previous code by modifying specific lines

FEM      = FunctionSpace(mesh,MixedElement([FE1,FE1,FER]))  **Line 1**
v,trz,h  = split(uknown) **Line 2**
Res5     = h*v*x[0]*phi[2]*dx - v*x[0]*phi[2]*dx_f1 **Line 3**

But I get NaN in the solution. The problem is at the first term of Res5 and probably the multiplication by the velocity v. If I write Res5 like this h*x[0]*phi[2]*dx - v*x[0]*phi[2]*dx_f1 (without v in the first term), everything is good.

Thanks