I am working with a viscoelastic model which requires updating the internal variable (C^v) at each time step inside the domain. Moreover the viscosity is a complex nonlinear function of the applied deformation C = F^T F where F is the deformation gradient tensor. Keeping it short, I have something like this

```
# Defining the weak form a_uv and the problem
...
...
Jac = derivative(a_uv, u, utrial)
F = NonlinearVariationalProblem(a_uv, u, bcs, J=Jac)
solver = NonlinearVariationalSolver(F)
# Loop over time steps:
for i, t in enumerate(timeArray):
# Assign the loading
solver.solve()
# Update the internal variable
Cv = Cv + some_nonlinear_function(u, Cv)
...
# Using a local projection after the update
Cv = local_project(Cv, V) # this process is extremely slow
```

I am using @bleyerj 's idea of `LocalSolver`

as mentioned here and @kamensky 's suggestions here . I have tried using `Quadrature`

as well as `DG`

elements for my internal variable but there is no significant performance gain (even after using `LocalSolver`

). Even a simple `UnitCubeMesh(2,2,2)`

takes > 30 minutes in a serial computation.

Remark:

- If I use constant viscosity (a simpler update for the internal variable) then the code runs pretty fast (as expected).
- I cannot use a monolithic scheme as suggested here as I need a different update rule than implicit euler.

Any pointers as to how the projection could be made faster/more efficient?

Thanks in advance.