Hi,
I am solving a nearly-incompressible hyperelastic problem. I am using NeoHookean material and the Cauchy stress for this material is given as:
This is my code. I am comparing Fenics result with analytical solution. When I pick the material parameter value c = 1e6, Fenics and analytical solution match with each other.
This is the comparison of Fenics and analytical with c=1e6 (Near the end of the code I have a code line to specify this value, i.e. solve_elasticity(geometry_3d(), 1e6, 0.1, 1.0, ‘Results’))
However, when I choose c=1e1, Fenics fails to solve. This is the error message.
---------------------------------------------------------------------------
RuntimeError Traceback (most recent call last)
<ipython-input-57-26ef847f64a6> in <module>
146 if __name__ == '__main__':
147 parameters['std_out_all_processes'] = False
--> 148 solve_elasticity(geometry_3d(), 1e1, 0.1, 1.0, 'Results')
149 print('Done !!!')
<ipython-input-57-26ef847f64a6> in solve_elasticity(facet_function, c1, dt, T_end, output_dir)
119 # Prepare to solve and solve
120 w0.assign(w)
--> 121 solver.solve()
122
123 # Save DefGrad to file in VTK format
RuntimeError:
*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
*** fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error: Unable to solve nonlinear system with NewtonSolver.
*** Reason: Newton solver did not converge because maximum number of iterations reached.
*** Where: This error was encountered inside NewtonSolver.cpp.
*** Process: 0
***
*** DOLFIN version: 2019.1.0
*** Git changeset: 42aa1842a10aa10b95547cade0497565f05abd49
*** ------------------------------------------------------------------------
Could anyone please help me out.
This is my code:
from dolfin import *
import matplotlib.pyplot as plt
import os
import dolfin as dlf
import numpy as np
import xml.etree.ElementTree as ET
import math
def geometry_3d():
"""Prepares 3D geometry. Returns facet function with 1, 2 on parts of
the boundary."""
mesh = UnitCubeMesh(2, 2, 2)
gdim = mesh.geometry().dim()
boundary_parts = MeshFunction('size_t', mesh, mesh.topology().dim()-1)
x0 = AutoSubDomain(lambda x: near(x[0], 0))
x1 = AutoSubDomain(lambda x: near(x[0], 1))
y0 = AutoSubDomain(lambda x: near(x[1], 0))
z0 = AutoSubDomain(lambda x: near(x[2], 0))
x0 .mark(boundary_parts, 1)
y0 .mark(boundary_parts, 2)
z0 .mark(boundary_parts, 3)
x1.mark(boundary_parts, 4)
return boundary_parts
def solve_elasticity(facet_function,c1, dt, T_end, output_dir):
# Get mesh and prepare boundary measure
mesh = facet_function.mesh()
gdim = mesh.geometry().dim()
dx = Measure("dx")
ds = Measure("ds", subdomain_data=facet_function, subdomain_id=4)
# Limit quadrature degree
dx = dx(degree=4)
ds = ds(degree=4)
# Build function space
element_v = VectorElement("P", mesh.ufl_cell(), 1)
element_s = FiniteElement("P", mesh.ufl_cell(), 1)
mixed_element = MixedElement([element_v, element_v, element_s])
W = FunctionSpace(mesh, mixed_element)
# def grad and stress output
W_DFnStress = TensorFunctionSpace(mesh, "DG", degree=0)
defGrad = Function(W_DFnStress, name='F')
PK1_stress = Function(W_DFnStress, name='PK1')
# Define constitutive law
def stress(u, p):
"""Returns 1st Piola-Kirchhoff stress and (local) mass balance
for given u, p."""
F = I + grad(u)
J = det(F)
B = F * F.T
devB = B - (1/3)*dlf.tr(B)*I
T = -p*I + (2*c1/J**(5.0/3.0))*devB # Cauchy stress
S = J*T*inv(F).T # 1st Piola-Kirchhoff stress
pp = -dlf.ln(J)/J
return S, pp
# Timestepping theta-method parameters
q = Constant(0.5)
dt = Constant(dt)
# Unknowns, values at previous step and test functions
w = Function(W)
u, v, p = split(w)
w0 = Function(W)
u0, v0, p0 = split(w0)
_u, _v, _p = TestFunctions(W)
I = Identity(W.mesh().geometry().dim())
# Balance of momentum
S, pp = stress(u, p)
S0, pp0 = stress(u0, p0)
F1 = (1.0/dt)*inner(u-u0, _u)*dx \
- ( q*inner(v, _u)*dx + (1.0-q)*inner(v0, _u)*dx )
F2a = inner(S, grad(_v))*dx + pp*_p*dx
F2b = inner(S0, grad(_v))*dx + pp0*_p*dx
F2 = (1.0/dt)*inner(v-v0, _v)*dx + q*F2a + (1.0-q)*F2b
# Whole system and its Jacobian
F = F1 + F2
J = derivative(F, w)
# Extract solution components
u, v, p = w.split()
# vector to store stretch and sterss
stress_vec=[]
stretch_vec=[]
# Time-stepping loop
t = 0
print('Start solving...')
while t <= T_end:
print("Time: {}".format(t))
# Prepare BCs
bc0 = DirichletBC(W.sub(0).sub(0), Constant(0.), facet_function, 1)
bc1 = DirichletBC(W.sub(0).sub(1), Constant(0.), facet_function, 2)
bc2 = DirichletBC(W.sub(0).sub(2), Constant(0.), facet_function, 3)
tDirBC = Expression(('0.1*time_'),time_=t , degree=0)
bc3 = DirichletBC(W.sub(0).sub(0), tDirBC, facet_function, 4)
bcs = [bc0,bc1,bc2,bc3]
# Initialize solver
problem = NonlinearVariationalProblem(F, w, bcs=bcs, J=J)
solver = NonlinearVariationalSolver(problem)
solver.parameters['newton_solver']['relative_tolerance'] = 1e-6
solver.parameters['newton_solver']['linear_solver'] = 'mumps'
# Prepare to solve and solve
w0.assign(w)
solver.solve()
# Save DefGrad to file in VTK format
DF = I + grad(u)
defGrad.assign(project(DF, W_DFnStress))
stretch_vec.append(defGrad(0.5,0,0)[0])
# Save Stress to file in VTK format
PK1_s,pp_t = stress(u, p)
PK1_stress.assign(project(PK1_s, W_DFnStress))
stress_vec.append(PK1_stress(0.5,0,0)[0])
# increment time
t += float(dt)
# compare FE and analytical solution
stress_vec=np.array(stress_vec,dtype=np.float64)
stretch_vec=np.array(stretch_vec,dtype=np.float64)
stress_ana = 2*c1*(stretch_vec-1.0/stretch_vec**2) # 1st PK stress
f = plt.figure()
plt.plot(stretch_vec,stress_vec,'r-')
plt.plot(stretch_vec,stress_ana,'b*')
if __name__ == '__main__':
parameters['std_out_all_processes'] = False
solve_elasticity(geometry_3d(), 1e2, 0.1, 1.0, 'Results')
print('Done !!!')