In the below MWE, I have a unit square. I am applying an exponential growth only in the X-direction with the help of the growth tensor `Fg = as_tensor([[exp(r*T),0.0],[0.0,1.0]])`

, where r and T are the rate of growth and time, respectively. I am implementing a hyperelastic model to solve \eta\frac{\partial \mathbf{u}}{\partial t}=S \nabla.\mathbf{T}, where \eta is the friction coefficient arising from interaction of the growing square with a substrate it is immersed in, \mathbf{u} is a displacement vector, \mathbf{T} is the first Piola-Kirchoff stress tensor, and S is a constant that has dimension of area.

I solve the variational form subject to the condition that the Y-displacement is 0 for all time steps (otherwise, it generates a Y-displacement and I want the square to grow only in the X-direction).

When I calculate the Cauchy and the PK1 stress, I find that the stress in the Y-direction is nonzero. I think it should be zero.

Here is the MWE:

```
from ufl import *
from ufl.classes import *
from dolfin import *
import numpy as np
mesh = UnitSquareMesh(8,8)
ds_custom = Measure("ds", domain=mesh) # Boundary measure
dx_custom = Measure("dx", domain=mesh) # Surface measure
FS = FunctionSpace(mesh, 'P', 1)
VS = VectorFunctionSpace(mesh, 'P', 1)
TS = TensorFunctionSpace(mesh, 'P', 1)
u = Function(VS)
u_n = Function(VS)
v = TestFunction(VS)
NumStepg = 5
Time = 2 # Units: mins
dt = Time/NumStepg
j = 0
eta = 5.0 # Friction coefficient
r = 0.02 # Rate of growth
mu1 = 500.0/(1.0*(0.1**3)/12.0) # Flexural rigidity
nu1 = 0.499 # poisson ratio
csa3 = 0.01
vtk1 = File("displacement.pvd")
vtk2 = File("PK1stress.pvd")
vtk3 = File("CauchyStress.pvd")
vtk4 = File("Fe.pvd")
vtk5 = File("dW_dFe.pvd")
for T in np.linspace(0,Time,NumStepg):
d = mesh.geometry().dim()
I = Identity(d) # Identity tensor
F = I + grad(u) # Deformation gradient
Fg = as_tensor([[exp(r*T),0.0],[0.0,1.0]]) # Growth tensor only in the X-direction
Fg_inv = as_tensor([[exp(-r*T),0.0],[0.0,1.0]]) # Inverse of the growth tensor
Fe = variable(F*Fg_inv)
Ce = variable(Fe.T*Fe)
Je = det(Fe)
# Invariants of deformation tensors
I1 = tr(Ce)
W = ((mu1/2)*(I1 - 3) + ((nu1)*(mu1)/(1-(2*nu1)))*(Je-1)**2 - mu1*(ln(Je)))
TT = ((det(Fg)*diff(W,Fe)*(Fg_inv.T))) # 1st PK Stress for compressible morphoelasticity
PE = (eta)*dot(u,v)*dx_custom + (csa3)*(dt)*inner(TT,grad(v))*dx_custom - (eta)*dot(u_n, v)*dx_custom
all_domains = MeshFunction("size_t", mesh, mesh.topology().dim() - 1, 1)
bcs2 = DirichletBC(VS.sub(1), Constant(0.0), all_domains, 1) # making Y- component of displacement zero
bcs = [bcs2]
# Solve variational problem
solve(PE == 0, u , bcs,
solver_parameters={"newton_solver": {"relative_tolerance": 5e-9,
"absolute_tolerance": 5e-10,"maximum_iterations": 50}})
print(f"Number of iterations: {j}", flush=True)
print(f"Time = {T} mins", flush=True)
j += 1
u_n.assign(u)
cauchystress = (2/Je)*((mu1/2)*(Fe*Fe.T - I) + (nu1*mu1*(Je-1)*Je*I)/(1-2*nu1)) # Cartesian Cauchy Stress for compressible morphoelasticity
PK1 = project(TT, TS)
disp = project(u, VS)
CS = project(cauchystress, TS)
dwdfe = project(diff(W,Fe), TS)
fe = project(Fe, TS)
PK1.rename('PK1','PK1')
disp.rename('disp','disp')
CS.rename('CS','CS')
dwdfe.rename('dwdfe','dwdfe')
fe.rename('fe','fe')
vtk1<<(disp, T)
vtk2<<(PK1, T)
vtk3<<(CS, T)
vtk4<<(fe, T)
vtk5<<(dwdfe, T)
```

Furthermore, when I calculate \frac{\partial W}{\partial \mathbf{F}_e} manually, it gives me zero values in the (1,1) position when \mathbf{F}_e has the (1,1) component 1 for all time points. Since W is a scalar function and \mathbf{F}_e has the value 1 in the (1,1) position for all time points, \frac{\partial W}{\partial \mathbf{F}_e} must have 0 at the (1,1) position since W can be rewritten as a function of the first (0,0) component of \mathbf{F}_e only.

Please let me know if there is anything wrong with the code or my understanding.

Thanks in advance!