Hello all,

I am new to Fenics and I want to implement a gradient descent algorithm to minimize the energy in the hyperelasticity demo. I couldn’t find any readily available solver for gradient descent in Fenics, and I tried to implement the algorithm by myself. However, when I used “subs” to substitute y_initial in the derivative of energy (F in this code), I got this error:

`

AttributeError: ‘Form’ object has no attribute ‘subs’

`

Can anyone help me with solving this problem?

Following is my minimal code:

```
import numpy as np
import matplotlib.pyplot as plt
from dolfin import *
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["representation"] = "uflacs"
# Create mesh and define function space
mesh = RectangleMesh(Point(-1,-1),Point(1,1), 50, 50)
V = VectorFunctionSpace(mesh, "Lagrange", 1)
################################################################################
# Define boundary condictions
################################################################################
F11=1
F12=0.0 + 0.1
F21=0.0 - 0.1
F22=1.0
y_BC = Expression(("F11*x[0]+F12*x[1]", "F21*x[0]+F22*x[1]"), F11=F11,F12=F12,F21=F21,F22=F22,degree=2)
tol = 1E-3
def x0_boundary(x, on_boundary):
return on_boundary and near(x[0], -1, tol) or near(x[0], 1, tol) or \
near(x[1], -1, tol) or near(x[1], 1, tol)
bcs = [DirichletBC(V, y_BC, x0_boundary)]
################################################################################
# Define functions
################################################################################
dy = TrialFunction(V) # Incremental displacement
v = TestFunction(V) # Test function
y = Function(V)
B = Constant((0, 0)) # Body force per unit area
T = Constant((0, 0)) # Traction force on the boundary
y_reference = Expression(("x[0]","x[1]"),degree=3)
################################################################################
# Elasticity model: Stored strain energy density (compressible neo-Hookean model)
################################################################################
E, nu = 10.0, 0.3
mu, blk = Constant(E/(2*(1 + nu))), Constant(E/(3*(1 - 2*nu)))
def W(y): #Stored strain energy density
F = grad(y) # Deformation gradient
C = F.T*F # Right Cauchy-Green tensor
Ic = tr(C)
J = det(F)
return (mu/2)*(Ic - 3) - mu*ln(J) + (blk/2 - mu/3)*(J-1)**2
################################################################################
# Total potential energy and solve
################################################################################
Pi = W(y)*dx -dot(B, (y-y_reference))*dx + dot(T, (y-y_reference))*ds
F = derivative(Pi, y, v) # first variation (derivative at u in the direction of v)
dF = derivative(F, y, dy) # Jacobian of F
####gradient descent solver:
alpha_ = 0.01
i=0
y_initial = project(y_BC, V)
y_new = Function(V)
y_new = project(y_BC, V)
while (np.abs(F.subs(y, y_new) ) >1e-6):
if i!=0:
y_initial = y_new
y_new = y_initial - alpha_ * J.subs(y, y_initial)
i = i+1
y = y_new
```

Thank you so much for your help!