Hi, everyone.
I have a convergence problem when tried to solve a simple electro-elastic equations. The energy functional is
\Psi = (\mu/2.0)(Ic - d) - \muln(J) + 0.5 \lambda \ln(J))^2+c1 I: (E \otimes E)+c2 C:(E \otimes E). The functional needs to be minimized with respect to the two unknowns - the displacment u and the electric potential \phi It seems qutie a simple problem; yet I have been stuck with it for days. I want to get the displacement u and the electric potential \phi.
Thank you.
====Error message =====================================
Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 2.291e+01 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 1: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
Newton iteration 2: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
Newton iteration 3: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
Newton iteration 4: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
Newton iteration 5: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
Newton iteration 6: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
Newton iteration 7: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
Newton iteration 8: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
Newton iteration 9: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
Newton iteration 10: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
2dtest4.geo ==========================================================
Point(1) = {0, 0, 0, 1};
Point(2) = {10, 0, 0, 1};
Point(3) = {20, 0, 0, 1};
Point(4) = {0, 10, 0, 1};
Point(5) = {0, 20, 0, 1};
Line(1) = {2, 3};
Line(2) = {4, 5};
Circle(3) = {4, 1, 2};
Circle(4) = {5, 1, 3};
Line Loop(7) = {2, 4, -1, -3};
Plane Surface(7) = {7};
Physical Line(1) = {2};
Physical Line(2) = {1};
Physical Line(3) = {3};
Physical Line(4) = {4};
Physical Surface(1) = {7};
=====================================================================
from dolfin import *
import numpy as np
Optimization options for the form compiler
parameters["form_compiler"]["quadrature_degree"] =2
parameters["form_compiler"]["cpp_optimize"] =True
parameters["form_compiler"]["optimize"] =True
ffc_options = {"optimize": True, \
"eliminate_zeros": True, \
"precompute_basis_const": True, \
"precompute_ip_const": True, \
"quadrature_degree": 2}
Create mesh and define function space
mesh = Mesh('2dtest4.xml')
boundaries = MeshFunction('size_t', mesh, '2dtest4_facet_region.xml')
V1 = VectorFunctionSpace(mesh, "Lagrange",1)
V2 = FunctionSpace(mesh,"Lagrange",1)
V= MixedFunctionSpace([V1,V2])
Define functions
V_ux=V1.sub(0)
V_uy=V1.sub(1)
vq=TestFunction(V)
dup=TrialFunction(V)
up=Function(V)
(u,phi)=split(up)
Define Dirichlet boundary
zero=Constant(0.0)
volta=Constant(5.0)
voltb=Constant(-5.0)
T=Constant((0.1,0.0))
bcul = DirichletBC(V_ux,zero,boundaries,1)
bcur = DirichletBC(V_uy,zero,boundaries,2)
bcpl = DirichletBC(V2,volta,boundaries,3)
bcpr = DirichletBC(V2,voltb,boundaries,4)
bcs = [bcul,bcur,bcpl,bcpr]
Kinematics
d = u.geometric_dimension()
I = Identity(d) # Identity tensor
F = I + grad(u) # Deformation gradient
C = F.T*F # Right Cauchy-Green tensor
Invariants of deformation tensors
Ic = tr(C)
Jf = det(F)
Elasticity parameters
mu, k = 5.0, 10.0
lmbda = Constant(k-2.0*mu/3.0)
c1, c2 = 10.0,6.0
EE=-grad(phi)
Stored strain energy density (incompressible neo-Hookean model)
psi = (mu/2.0)(Ic - d) - muln(Jf) + (lmbda/2.0)*(ln(Jf))**2\
+c1*inner(I,outer(EE,EE))\
+c2*inner(C,outer(EE,EE))
ds = Measure("ds")[boundaries]
Total potential energy
Pi=psidx-dot(T, u)ds(4)
Compute first variation of Pi (directional derivative about u in the direction of v)
Func = derivative(Pi, up, vq)
Compute Jacobian of F
J = derivative(Func, up, dup)
Solve variational problem
problem=NonlinearVariationalProblem(Func,up,bcs,J=J,\
form_compiler_parameters=ffc_options)
solver=NonlinearVariationalSolver(problem)
solver.solve()
Save solution in VTK format
(u,phi)=up.split()
Plot and hold solution
plot(u,title="displacement")
plot(phi,title="potential")
interactive()
