if it is urgent then you can compute the RHS by a matrix product
use a couple Finite element
fespace Wh(Th,[P1dc,P1dc,P1dc]; );
Vh [ux, uy,uz],[ uu1, uu2,uu3], [dux,duy,duz], [vx, vy, vz];
varf vM([ uu1, uu2,uu3],[vx, vy, vz]) = -intallfaces(Th)(
((1./dt)*k *(mean(dx(uu1))*N.x*jump(vx)+mean(dy(uu1))*N.y*jump(vx)+mean(dz(uu1))*N.z*jump(vx)+mean(dx(uu2))*N.x*jump(vy)+mean(dy(uu2))*N.y*jump(vy)+mean(dz(uu2))*N.z*jump(vy)+mean(dx(uu3))*N.x*jump(vz)+mean(dy(uu3))*N.y*jump(vz)+mean(dz(uu3))*N.z*jump(vz)
+mean(dx(vx))*N.x*jump(uu1)+mean(dy(vx))*N.y*jump(uu1)+mean(dz(vx))*N.z*jump(uu1)+mean(dx(vy))*N.x*jump(uu2)+mean(dy(vy))*N.y*jump(uu2)+mean(dz(vy))*N.z*jump(uu2)+mean(dx(vz))*N.x*jump(uu3)+mean(dy(vz))*N.y*jump(uu3)+mean(dz(vz))*N.z*jump(uu3)
+ (pena*NN)*(jump(uu1)*jump(vx)+jump(uu2)*jump(vy)+jump(uu3)*jump(vz))
))/ nElementonB
)
matrix M(Wh,Wh);
real[int] RHS = M*uu1[];
use RHS as you rhs of you equation !
like in Time dependent schema optimization for heat equations