I am using fenics version 2019.2.0.dev0 in python to solve a drift-diffusion PDE for current flow in semiconductors and have a solution for the PDE. I am trying to extract the current flow across a few internal and external boundaries by integrating the current density function over the boundary surface. I followed a few examples from the discourse to define such surface integrals however when trying to integrate the current density function the code fails with a

`` assert expr.*ufl_is_terminal_modifier*

AssertionError"

when performing an assemble. My code looks as follows

```
from dolfin import *
import numpy as np
# Define the mesh and function space
mesh = RectangleMesh(Point(0,0),Point(40,40),40,10)
element = VectorElement('P',triangle,3,dim=3)
V = FunctionSpace(mesh, element)
# Define subdomains of interest
SubDomain1 = CompiledSubDomain( '(x[0]>= -1e-14)&&(x[0]<=20+1e-14)' )
SubDomain2 = CompiledSubDomain( '(x[0]>= 20-1e-14)&&(x[0]<=40+1e-14)' )
# Define the boundaries of interest
Boundary1 = CompiledSubDomain( 'near(x[0],0)' )
Boundary2 = CompiledSubDomain( 'near(x[0],20)' )
Boundary3 = CompiledSubDomain( 'near(x[0],40)' )
# Assign Id to boundaries
BoundaryId = MeshFunction('size_t',mesh,mesh.topology().dim()-1)
Boundary1.mark(BoundaryId,1)
Boundary2.mark(BoundaryId,2)
Boundary3.mark(BoundaryId,3)
# create measures to integrate over the marked boundaries (internal and external boundaries use ds and dS respectively)
dMeas_int = ds(domain = mesh,subdomain_data=BoundaryId)
dMeas_ext = dS(domain = mesh,subdomain_data=BoundaryId)
# Define the PDE to solve
############# Initial condition #####################
ic = Expression(('p0','n0','2.0-0.05*x[0]'),degree=3,p0=1e-18,n0=1e-4)
u0 = project(ic,V)
pi_p0, pi_n0, phi_0 = split(u0)
############# Dirichlet boundary conditions #########
DC1 = DirichletBC( V.sub(2),Constant(2.0),Boundary1 )
DC2 = DirichletBC( V.sub(2),Constant(0.0),Boundary3 )
############## Trial and solution functions ############
u = Function(V)
pi_p, pi_n, phi_i = split(u)
vi_p,vi_n, vi_phi = TestFunctions(V)
Rip = 0.
Rin=0.
gamma = 0.02585
dt = 0.01
############### Potential PDE ###################
ai_phi = (0.6454775 * dot(grad(phi_i),grad(vi_phi)) - (pi_p-pi_n)*vi_phi)*dx
Li_phi = 1e-4 * vi_phi * dx
######## Drift Diffusion Current ###############
Jip = gamma * grad(pi_p)+pi_p * grad(phi_i)
Jin = gamma * grad(pi_n)-pi_n * grad(phi_i)
######## Define the weak pde form #############
ai_p = (pi_p * vi_p + dot(Jip, grad(vi_p) ) * dt ) * dx
ai_n = (pi_n * vi_n + dot(Jin, grad(vi_n) ) * dt ) * dx
Lip = (pi_p0 * vi_p + Rip * vi_p * dt) * dx
Lin = (pi_n0 * vi_n + Rin * vi_n*dt) * dx
F = ai_p+ai_n+ai_phi - Lip - Lin- Li_phi
# Solve the PDE
solve(F==0,u,[DC1,DC2])
```

The code upto this point works fine and a solution to the problem is obtained in u.

```
# Integrate the current density to get total current flowing (the problematic part that failed)
nx = Constant(np.array((1.,0.)))
J = Jin-Jip
Jx = dot(J('+'),nx)
########## Surface integral forms to get the total current #############
I1 = Jx * dMeas_int(1)+Jx * dMeas_ext(1) # integration over Boundary1
I2 = Jx * dMeas_int(2)+Jx * dMeas_ext(2) # integration over Boundary2
I3 = Jx * dMeas_int(3)+Jx * dMeas_ext(3) # integration over Boundary3
Curr1 = assemble(I1)
Curr2 = assemble(I2)
Curr3 = assemble(I3)
```

Here on the line of ``assemble(I1)" I get a very long run-time error saying

```
Calling FFC just-in-time (JIT) compiler, this may take some time.
Traceback (most recent call last):
...
...
...
in balance_modified_terminal
assert expr._ufl_is_terminal_modifier_
AssertionError
```

The error gives no clue as to what is wrong with the code and I am stuck trying to evaluate the total current integral

Any idea what might be going wrong?