Hi everyone,

in my last issue there was great help with one equation of the system which I am trying to solve (Problem with nonlinear variational problem). The full system which i try to solve is:

\begin{align}
\int u\psi -\int u^n\psi +\Delta t\left(\int\nabla p\psi +\int \sigma\cdot\nabla\psi + \mu\int\nabla u\cdot \nabla\psi\right)&=0\\
\sigma &=C-1\\
\int \left(C-C^n+\frac{1}{\lambda}(C-1)\right)v -\Delta t\begin{pmatrix}
\frac{1}{\lambda}\int\partial_1u v_1\\
\frac{1}{\lambda}\int\partial_2u v_2\\
2\int (\partial_1 uC_1+\partial_2 u C_2)v_3
\end{pmatrix}&= 0
\end{align}

Where \sigma and C are 3 dimensional vectors and u is scalar. I tried to solve this task by using this piece of code (only the things relevant for the variational form):

```
domain=mshr.Cylinder(Point(0,0,0),Point(2,0,0),1.,1.)
mesh = mshr.generate_mesh(domain, nCells)
StressTensorElement1 = FiniteElement('P', tetrahedron, 1)
StressTensorElement2 = FiniteElement('P', tetrahedron, 1)
StressTensorElement3 = FiniteElement('P', tetrahedron, 1)
VelocityElements = FiniteElement('P', tetrahedron, 1)
element = MixedElement(StressTensorElement1, StressTensorElement2, StressTensorElement3, VelocityElements)
V = FunctionSpace(mesh, element)
# stuff with boundary conditions was here but not relevant saved in "bc"
dt = Constant(0.1)
mu = Constant(0.1)
gradP = Constant(5)
Lambda = Constant(2.)
C1, C2, C3, u = TrialFunctions(V)
CV1, CV2, CV3, v = TestFunctions(V)
u_n = interpolate(Constant((0,0,1,0)), V) # previous timestep and initial condition at start
un1, un2, un3, un4 = split(u_n)
def getSigma(self, C1, C2, C3, Lambda):
"""get sigma from C"""
return 0.5/Lambda*as_vector([C1-1, C2-1, C3-1])
F = u * v * dx - un4 * v * dx + dt*(gradP * v * dx\
+ dot(getSigma(uk1, uk2, uk3, Lambda), grad(v)) * dx)\
+ Constant(mu) * dot(grad(u), grad(v)) * dx \
+ (C1 - un1 + 1 / Lambda * (C1 - 1)) * CV1 * dx\
+ (C2 - un2 + 1 / Lambda * (C2 - 1)) * CV2 * dx \
+ (C3 - un3 + 1 / Lambda * (C3 - 1)) * CV3 * dx \
- dt *(1/Lambda*u.dx(0) * CV1 * dx \
+ 1 / Lambda * u.dx(1) * CV2 * dx + 2 * (u.dx(0) * C1 + u.dx(1) * C2) * CV3 * dx)
U = Function(V)
solve(F == 0, U, bc)
```

I omitted the timeloop here. The error I get is

`adding expression with non-matching form arguments`

But this problem is definitly nonlinear and if i call it the linear way it will of course still not work.

Thank you for your time to read this!