Hi,
For some reason, I am not sure why FeNicS is taking to run for a mixed field hyperelasticity problem. There are only 6 cells and the number of dofs shouldn't be a lot. Could be due to the way I wrote the code (see below) ? Thanks!
from dolfin import *
import numpy as np
from matplotlib import pylab as plt
from scipy import optimize
import os
parameters["form_compiler"]["quadrature_degree"] = 4
parameters["form_compiler"]["cpp_optimize"] = True
ffc_options = {"optimize": True, \
"eliminate_zeros": True, \
"precompute_basis_const": True, \
"precompute_ip_const": True}
nx = 1
ny = 1
nz = 1
mesh = BoxMesh(0, 0, 0, 10, 10, 1, nx, ny, nz)
def analytical_soln_disp(u):
C1 = 0.4#*0.0000980665 #MPa
x = [u, u]
x3 = 1.0/(x[0]*x[1])**2.0
sigma1 = 2*C1*(x[0]**2.0 - x3**2.0)
sigma2 = 2*C1*(x[1]**2.0 - x3**2.0)
return sigma1
def Fmat(u):
d = u.geometric_dimension()
I = Identity(d) # Identity tensor
F = I + grad(u) # Deformation gradient
return F
def Emat(u):
F = Fmat(u)
d = u.geometric_dimension()
I = Identity(d) # Identity tensor
return 0.5(F.TF - I)
############ Mooney Rivilin
def psiE_neohookean(E, p):
C1 = 0.4
d = u.geometric_dimension()
I = Identity(d) # Identity tensor
C = 2*E + I
J = pow(det(C), 0.5)
Cdev = J**(-2.0/3.0)*C
Ic = tr(C)
Ic_dev = pow(J, -float(2)/3)*Ic
SE = C1*(Ic_dev - 3) + p*(J - 1)
return SE
#
def Smat(u, p):
C1 = 0.4
E = Emat(u)
E = variable(E)
d = u.geometric_dimension()
I = Identity(d) # Identity tensor
C = 2*E + I
J = pow(det(C), 0.5)
Ic = tr(C)
Ic_dev = pow(J, -float(2)/3)*Ic
Psi = C1*(Ic_dev - 3) + p*(J - 1)
return diff(Psi,E)
def Sigma_mat(u, Kappa):
F = Fmat(u)
S = Smat(u, Kappa)
J = det(F)
return 1/J*F*S*F.T
Define boundary
class Left(SubDomain):
def inside(self, x, on_boundary):
return abs(x[0] - 0.0) < DOLFIN_EPS and on_boundary
class Right(SubDomain):
def inside(self, x, on_boundary):
return abs(x[0] - 10.0) < DOLFIN_EPS and on_boundary
class Bottom(SubDomain):
def inside(self, x, on_boundary):
return abs(x[2] - 0.0) < DOLFIN_EPS and on_boundary
class Top(SubDomain):
def inside(self, x, on_boundary):
return abs(x[2] - 1.0) < DOLFIN_EPS and on_boundary
class Lower(SubDomain):
def inside(self, x, on_boundary):
return abs(x[1] - 0.0) < DOLFIN_EPS and on_boundary
class Upper(SubDomain):
def inside(self, x, on_boundary):
return abs(x[1] - 10.0) < DOLFIN_EPS and on_boundary
V = VectorFunctionSpace(mesh, 'CG', 2)
Q = FunctionSpace(mesh, 'DG', 0)
W = MixedFunctionSpace([V,Q]) # mixed space
Initialize sub-domain instances
left = Left()
right = Right()
bottom = Bottom()
top = Top()
lower = Lower()
upper = Upper()
Initialize mesh function for boundary domains
boundaries = FacetFunction("size_t", mesh)
boundaries.set_all(0)
left.mark(boundaries, 1)
right.mark(boundaries, 2)
bottom.mark(boundaries, 3)
top.mark(boundaries, 4)
lower.mark(boundaries, 5)
upper.mark(boundaries, 6)
ds = Measure("ds")[boundaries]
n = FacetNormal(mesh)
Set up boundary condition at left end
c = Constant((0.0))
bcleft0 = DirichletBC(W.sub(0).sub(0), c, boundaries, 1)
bclower1 = DirichletBC(W.sub(0).sub(1), c, boundaries, 5)
bcbottom2 = DirichletBC(W.sub(0).sub(2), c, boundaries, 3)
uspecified = Expression("u", u = 1.0)
bc_right0 = DirichletBC(W.sub(0).sub(0), uspecified, boundaries, 2)
bc_upper1 = DirichletBC(W.sub(0).sub(1), uspecified, boundaries, 6)
Set up boundary conditions
bcs = [bcleft0, bclower1, bc_right0, bc_upper1, bcbottom2]
Define variational problem
du = TrialFunction(V)
u = Function(V)
v = TestFunction(V)
dup = TrialFunction(W)
up = Function(W)
vp = TestFunction(W)
up = Function(W)
u, p = split(up)
Psi = psiE_neohookean(Emat(u), p)
Total potential energy
Pi = Psi*dx
Pi = Psidx #- dot(pn, u)ds(2) - dot(pn, u)*ds(6)
Compute first variation of Pi (directional derivative about u in the direction of v)
F = derivative(Pi, up, vp)
Compute Jacobian of F
a = derivative(F, up, dup)
problem = NonlinearVariationalProblem(F, up, bcs, a)
solver = NonlinearVariationalSolver(problem)
prm = solver.parameters
prm['newton_solver']['absolute_tolerance'] = 1E-8
prm['newton_solver']['relative_tolerance'] = 1E-7
prm['newton_solver']['maximum_iterations'] = 25
prm['newton_solver']['relaxation_parameter'] = 1.0
usol, psol = up.split()
vtk_file = File("nonlinearelasticity.pvd")
vtk_file << usol
lambda_array = np.empty((0))
sigma_array = np.empty((0))
sigma_analytical = np.empty((0))
u_array = np.linspace(0,0.8,10)
p_array = np.linspace(0,100,20)
for udisp in u_array:
uspecified.u = udisp
solver.solve()
usol, psol = up.split()
## Compute solution
F = project(Fmat(usol),TensorFunctionSpace(mesh, "DG", 0))
S = project(Sigma_mat(usol, psol),TensorFunctionSpace(mesh, "DG", 0))
print S(10,5,0.5)[0], F(10,5,0.5)[0], F(10,5,0.5)[4], F(10,5,0.5)[8]
lbda = F(10,5,0.5)[0]
analytical_s11 = analytical_soln_disp(lbda)
sigma_analytical = np.append(sigma_analytical, analytical_s11)
lambda_array = np.append(lambda_array, F(10,5,0.5)[0])
sigma_array = np.append(sigma_array, S(10,5,0.5)[0])
vtk_file << usol
plt.plot(lambda_array, sigma_array, '--*')
plt.plot(lambda_array, sigma_analytical, '-')
plt.show()
