Converting Python code for stream function calculation for lid-driven cavity flow problem to cpp

I am trying to reconstruct the graphs in the FEniCS manual, mainly that of chapter 21. I have done so already, however I still would like to compare some of the different solver's accuracy with respect to another. The problem I am doing the comparison about is the lid-driven cavity flow problem in two dimensions -- also chapter 21.

My issue is that in that chapter the small snippet for the python code (page 407 -- section 21.4.2) is confusing to convert it to the cpp syntax -- at least for someone beginner as me momentarily. Below I have provided basically all what I am trying to do, but not working and in the main.cpp file, actually getting an error(s). Can someone please guide me through this ? I do not anticipate it being complicated, just a syntax thing.

I have tried writing a separate ufl file that includes the following:

# define function space
element  = FiniteElement("Lagrange", triangle, 2)
Velement = VectorElement("Lagrange", triangle, 2)

# define trial function
psi = TrialFunction(element)

# define test function
q = TestFunction(element)

# define coefficients
u = Coefficient(Velement)

# define bilinear form
a = inner(grad(psi), grad(q))*dx

# define linear form
L = inner( u[1].dx(0) - u[0].dx(1), q )*dx

from my expectation that should take the two-dimensional solution vector u and return the streamfunction variable psi -- after solving the system above, mind you.

Next, I went on to the main.cpp file which I wrote in cpp that calculates the lid-driven cavity flow in two-dimensions and produces the output correctly. I then inserted the following mess (sigh) after computing the solution at the final time. It is as follows:

  StreamFunction::BilinearForm a4(V, V);
  StreamFunction::LinearForm   L4(V);

  L4.u = u1;

  Matrix A4;
  Vector b4;

  auto psi = std::make_shared<Function>(V); // stream function

  assemble(A4, a4);
  assemble(b4, L4);

  DirichletBC bc_psi(V, zero, DomainBoundary());
  std::vector<DirichletBC*> bcPSI = { &bc_psi};

  for (std::size_t i=0; i<bcPSI.size(); i++){
    bcPSI[i]->apply(A4, b4);
  }
  solve(A4, *psi->vector(), b4);

  plot(psi, "streamfunction");

here, in the above, the parameter u1 is the velocity component consisting of both x and y coordinates after the final time. V is the function space I already used in the original code to determine the final velocity, its definition is:

  auto V = std::make_shared<VelocityUpdate::FunctionSpace>(mesh);

Thank you for all your help...