Hi,

I have been trying to implement a hyperelastic material model defined here using an axisymmetric formulation for a cylinder (a 2 dimensional slice of). I have gone through what I can to account for missing terms in the 2D representation, changing

```
Pi = psi*df.dx
```

to

```
Pi = psi*2*np.pi*X[1]*df.dx
```

for example (revolving around X[0]).

FEniCS returns a reasonable looking solution for the displacement field, but I’m struggling to work out if I need to change the lines

```
residual = df.derivative(Pi, u, v)
Jacobian = df.derivative(residual, u, du)
```

to account for missing terms in the 2D representation. I’m not sure how I could change them very easily if I do need to. I have left my reduced code below and I’m using FEniCS version 2019.2.0.dev0

Thanks,

Kit

```
import dolfin as df, numpy as np
# Define mesh and displacement of RHS
L, R, density = 10, 3, 10
mesh = df.RectangleMesh(df.Point(0, 0), df.Point(L, R), L*density, R*density)
displacement = 2
# Define function space
V = df.VectorFunctionSpace(mesh, 'CG', 1)
du = df.TrialFunction(V)
v = df.TestFunction(V)
u = df.Function(V)
# Define material parameters
k0 = 1e3
mu0 = 50
k1 = 1e3
k2 = 2
a04 = df.Constant((1.0, 0))
# Kinematics
X = df.SpatialCoordinate(mesh) # X[0] = Z, X[1] = R
dimension = mesh.geometry().dim()
I = df.Identity(dimension) # Identity tensor
F = I + df.grad(u) # Deformation gradient
C = F.T * F # Right Cauchy-Green tensor
B = F * F.T # Left Cauchy-Green tensor
J_naive = df.det(F) # Naive determinant of F
J = (1 + u[1]/X[1]) * J_naive # Corrected determinant of F
I1_naive = df.tr(C) # Naive first invariant of C
I1 = I1_naive + (1 + u[1]/X[1]) ** 2 # Corrected first invariant of C
I1_ico = J ** (-2/3) * I1 # Isochoric first invariant of C
I4 = df.dot(a04, C * a04) # Anisotropic invariant
# Strain energy density
psi_vol = 1/2 * k0 * (J - 1) ** 2
psi_ico = 1/2 * mu0 * (I1_ico - 3)
psi_ani = (k1 / (2 * k2)) * (df.exp(k2 * (I4 - 1) ** 2) - 1)
psi = psi_vol + psi_ico + psi_ani
# Total potential energy
Pi = psi*2*np.pi*X[1]*df.dx
# Define boundaries: x[0] = Z, x[1] = R
x = np.unique(mesh.coordinates()[:,0])
y = np.unique(mesh.coordinates()[:,1])
class Left(df.SubDomain):
def inside(self, x, on_boundary):
return df.near(x[0], min(x)) and on_boundary
class Right(df.SubDomain):
def inside(self, x, on_boundary):
return df.near(x[0], max(x)) and on_boundary
class Axis(df.SubDomain):
def inside(self, x, on_boundary):
return df.near(x[1], min(y)) and on_boundary
class CurvedSurface(df.SubDomain):
def inside(self, x, on_boundary):
return df.near(x[1], max(y)) and on_boundary
# Mark facets
facets = df.MeshFunction("size_t", mesh, 1)
facets.set_all(0)
Left().mark(facets, 1)
Right().mark(facets, 2)
Axis().mark(facets, 3)
CurvedSurface().mark(facets, 4)
# Clamped at left, sliding at right, sliding on axis
bcs = [df.DirichletBC(V, df.Constant((0., 0.)), facets, 1),
df.DirichletBC(V.sub(0), df.Constant(displacement), facets, 2),
df.DirichletBC(V.sub(1), df.Constant(0.), facets, 3)]
# Displacement field solution
residual = df.derivative(Pi, u, v)
Jacobian = df.derivative(residual, u, du)
df.solve(residual == 0, u, bcs, J = Jacobian,
solver_parameters={"newton_solver":{"linear_solver":"lu"}})
df.plot(u, mode='displacement')
```