Newton iteration 47: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton iteration 48: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton iteration 49: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton iteration 50: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Traceback (most recent call last):
File "demo_leaf.py", line 83, in
solve(F == 0, u_total)
File "/usr/lib/python2.7/dist-packages/dolfin/fem/solving.py", line 300, in solve
_solve_varproblem(*args, **kwargs)
File "/usr/lib/python2.7/dist-packages/dolfin/fem/solving.py", line 349, in _solve_varproblem
solver.solve()
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: 2016.1.0
*** Git changeset: unknown
*** -------------------------------------------------------------------------
from dolfin import *
import numpy as np
import scipy as sp
Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
ffc_options = {"optimize": True, \
"eliminate_zeros": True, \
"precompute_basis_const": True, \
"precompute_ip_const": True}
Sub domain for Periodic boundary condition
class PeriodicBoundary(SubDomain):
# Left boundary is "target domain" G
def inside(self, x, on_boundary):
return bool(x[0] < DOLFIN_EPS and x[0] > -DOLFIN_EPS and on_boundary)
# Map right boundary (H) to left boundary (G)
def map(self, x, y):
y[0] = x[0] - 20.0
y[1] = x[1]
computation domain
L = 5.0
W = 1.0
h = 0.02
growth epsilon parameters
n = 10
beta = 0.065
epsilon_g = Expression("0.065*pow(x[1]/1, 10)")
create mesh
mesh = RectangleMesh(Point(0.0, 0.0), Point(L, W), 20, 4, "right/left")
print("Plotting a RectangleMesh")
plot(mesh, title="Rectangle")
define function space
P2 = VectorElement("Lagrange", mesh.ufl_cell(), 2)
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
TH = P2 * P1
W = FunctionSpace(mesh, TH)
Mark boundary subdomians
left = CompiledSubDomain("near(x[0], side) && on_boundary", side = 0.0)
right = CompiledSubDomain("near(x[0], side) && on_boundary", side = 20.0)
define functions
du_total = TrialFunction(W)
u_total = Function(W)
v = TestFunction(W)
u, zeta = split(u_total)
Kinematics
epsilon = sym(grad(u)) + grad(grad(zeta))/2
epsilon_gg = as_matrix([[epsilon_g, 0], [0, 0]])
epsilon_e = epsilon - epsilon_gg
Elasticity parameters
E = 1 # the Young 's modulus
mu = 0.3 # the Poisson 's ratio
the bending energy density
F_b = Epow(h, 3)/12/(1 - pow(mu, 2))(pow(tr(grad(grad(zeta))), 2) + 2(1 - mu)det(tr(grad(grad(zeta)))))*dx
the stretching energy density
F_s = Eh/(1 - pow(mu, 2))(pow(tr(epsilon_e), 2) + 2(1 - mu)det(epsilon_e))*dx
total potential energy
Psi = F_s + F_b
compute first variation of Pi
F = derivative(Psi, u_total, v)
computer jacobian of F
J = derivative(F, u_total, du_total)
solver variational problem
solve(F == 0, u_total)
save solution in VTK format
file = File("output/leaf.pvd");
file << zeta
Plot and hold solution
plot(zeta, mode = "displacement")
interactive()
