# Penalty
k = 1
d = 1
# Other parameters 
n = FacetNormal(mesh) 
h = CellSize(mesh) 
vn = 0.5*(dot(vel,n) + abs(dot(vel,n))) 
# Weak form
a = v*u/dt*dx - dot(grad(v),vel*u)*dx + dot(jump(v),vn('+')*u('+')-vn('-')*u('-'))*dS + dot(v, vn*u)*ds \
  + dot(grad(v),D*grad(u))*dx - dot(jump(v,n),avg(D*grad(u)))*dS - dot(avg(D*grad(v)),jump(u,n))*dS \
  + gamma*dot(jump(v,n),jump(u,n))*dS  
L = v*f*dx + v*u0/dt*dx
...
I have a couple questions:
1) If say I am solving for the transport of 4 chemical species, Normally I would create a mixed function space like V_mixed = V*V*V*V, split u into u1, u2, u3, and u4, write bilinear and linear forms for each of these individual components and add them
 up. That is, a = a1 + a2 + a3 + a4, L = L1 + L2 + L3 + L4. This can become quite a lot of code if I end up having to solve ~20 chemical species eventually. The catch is that the velocity and diffusivity will be the same for all species but the forcing function
 and boundary conditions are going to be different. Is there a more efficient of writing the bilinear and linear form instead of having to copy and paste the a = ... and L = ... multiple times and changing the trial/test functions?