Derivative function has no attribute 'subs'

variational formulation
1

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!

2

Note that F is a variational formulation (symbolical) written in the UFL form language.
To replace a coefficient in a ufl-form, you need to use ufl.replace. i.e.

from ufl import replace

F_ = replace(F, {y:y_new})

Note that this returns a ufl-form, a symbolical representation of F. To get a numerical value, you need to assemble it:

F_assembled = assemble(F_)

This returns a dolfin.Vector which you can take the norm of:

F_norm = F_assembled.norm("l2")

For your specific problem, this would end up as:

from ufl import replace
while  assemble(replace(F, {y:y_new})).norm("l2") >1e-6:
     # Add code

The same hold for this assignment, you need to use replace and assembly.

Please also note that you are constantly reassigning values in your loop.
You should use the assign function :

y_new = Function(V)
y_new.assign(project(y_BC, V))
while ...:
      if i!=0:
          y_initial.assign(y_new)
     y_new.vector()[:] = y_initial.vector()[:] - assemble(.....)[:]

Note that J is not defined outside of W(y).

2 Likes
3

Maybe I am skipping some details, or the formulation is meant to be a simplified version of another more involved problem, but the energy should be mu/2. * (Ic - 2) for 2D right? You may also want to take a look at Elasto-plastic analysis of a 2D von Mises material — Numerical tours of continuum mechanics using FEniCS master documentation or a post from the old forum Newton method programmed manually - FEniCS Q&A

Both of them write the Newton’s method manually, but that should be easily modifiable for gradient descent.

In any case, a simplified version of how I would go about implementing gradient descent would be something like (I have used displacement u instead of the deformation mapping y as in your case)

from dolfin import *
from copy import deepcopy

mesh = UnitCubeMesh(3,3,3)
V = VectorFunctionSpace(mesh, "CG", 1)
u = Function(V)
v = TestFunction(V)
du = TrialFunction(V)
B = Constant((0., -0.5, 0.))

left = CompiledSubDomain("near(x[0], 0.)")
bcs = [
    DirichletBC(V, Constant((0., 0., 0.)), left)
]

E, nu = Constant(10.), Constant(0.3)
mu = E/2./(1+nu)
lmbda = 2 * mu * nu / (1. - 2* nu)

F = Identity(len(u)) + grad(u)

Pi = mu/2. * (inner(F, F) - 3.) - mu * ln(det(F)) + lmbda/2. * (det(F) - 1.)**2
Pi *= dx
Pi -= inner(B, u) * dx

residual = derivative(Pi, u, v)
jacobian = derivative(residual, u, du)

A, b = assemble_system(jacobian, -residual, bcs)

res_initial = b.norm("l2")
res = deepcopy(res_initial)

alph = 0.05
max_iters = 5000
iters = 1

while res/res_initial > 1.e-6 and iters < max_iters:
    u.vector()[:] += alph * b[:]
    _, b = assemble_system(jacobian, -residual, bcs)
    res = b.norm("l2")
    print(f"Norm: {res}, iter: {iters}")
    iters += 1

It would be terribly slow, but then again the problem itself could warrant such a method in some cases.

2 Likes
4

Thank you so much for your very helpful responses!

5

Thanks for catching that. Yes, it should be mu/2. * (Ic - 2) for my 2D problem.

6

Hello,
Regarding the code that you have provided, I was wondering what is the initial u vector that assemble_system(jacobian, -residual, bcs) uses to calculate the b and res_initial before entering the while loop?
Also, did you change the y vector to the u vector on purpose?

Thanks a lot!

7

It assembles the residual with the initial u (which in this case is simply zeros) after applying the dirichlet boundary conditions bcs and, in this case, the body force term (B). For this specific case, it is simply the contribution from the last term

I changed the formulation to the displacement field u since that is more common practice, but should work for the deformation field too. The example that I linked nicely illustrates the Newton-iteration loop embedded within a global step-loop.

1 Like