An efficient way to calculate u = (rho_u)/rho

Question

Given VectorElements u and rho_u, and a FiniteElement rho, whats an efficient way to calculate u = rho_u / rho when running with a large number of processes?

I've tried interpolating an expression, which is convenient to code, but it takes just as much time as my pressure solve to perform the interpolation.

I've tried setting up a linear system and solving for it, which also works but is even slower. I didn't expect this one to be fastest, though I was hopeful that having the system pre-assembled would make solving it relatively fast.

I've also tried the vector operation directly with the local portion of the global vector,

rho->vector()->get_local(rho_val);
rhou->vector()->get_local(rhou_val);
u->vector()->get_local(u_val);

std::pair<std::size_t, std::size_t> V1range = V1->dofmap()->ownership_range()


for(size_t i=0; i < rhou_val.size(); ++i)
{
    u_val[i] = rhou_val[i] / rho_val[i/dimensions];
}

u->vector()->set_local(u_val);

which works if I make assumptions on how the vector space's degrees of freedom are mapped relative to the scalar's. It also requires them to have collocated degrees of freedom which I generally don't have.

I've also noticed dividing by rho in my variational forms results in nearly doubling the quadrature degree required. Pre-computing rho_inv = 1 / rho on the vector directly has been helpful there.