Hello
I am solving a problem in hyperelasticity where I have a new term in strain energy function which is defined as:
In the above equation:
Where Ei represents the integral exponential function and C are the components of the right Cauchy stress tensor.
As the Ei function does not exist in the UFL, so I used function approximation based on this formula: 
. Gamma is a constant and Ein represents entire function that is defined as: 
. (The whole approximation could be found here. I implemented (and tested) this function approximation into my code. I also used trapezoidal rule for solving the first integral (G(\theta)). I am using 10 terms in order to approximate the Ei function and also 10 intervals for the integral calculation but this code is unbelievably slow although there is only 1 element!
Here is the full implementation:
from dolfin import *
import numpy as np
mesh = UnitCubeMesh.create(1,1,1,CellType.Type.hexahedron)
# 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] - 1.0) < DOLFIN_EPS and on_boundary
Element1 = FiniteElement("CG", mesh.ufl_cell(), 1)
Element2 = VectorElement("CG", mesh.ufl_cell(), 2)
W_elem = MixedElement([Element1, Element2])
W = FunctionSpace(mesh, W_elem)
left = Left()
right = Right()
boundaries = MeshFunction('size_t', mesh, mesh.topology().dim()-1)
subdomains = MeshFunction('size_t', mesh, mesh.topology().dim())
boundaries.set_all(0)
left.mark(boundaries, 1)
right.mark(boundaries, 2)
dx = Measure('dx', domain=mesh, subdomain_data=subdomains, metadata={'quadrature_degree': 10})
ds = Measure('ds', domain=mesh, subdomain_data=boundaries, metadata={'quadrature_degree': 10})
w = Function(W)
p,u = split(w) # p corresponds to hydrostatic pressure
d = u.geometric_dimension()
I = Identity(d)
F = I + grad(u)
C = F.T*F
Ic = tr(C)
J = det(F)
def trapezoidal(f, a, b, n):
h = float(b - a) / n
s = 0.0
s += f(a)/2.0
for i in range(1, n):
s += f(a + i*h)
s += f(b)/2.0
return s * h
c_1 = 10.
gamma = 0.5772156649
def Ei(X, k, b):
LN = ln(X)
for s in range(1, k + 1):
fact = 1
for i in range(1, s + 1):
fact = fact * i
b += (pow(-1., s + 1) * pow(-X, s) / (s * fact))
return gamma + LN - b
def G(theta):
LMDA = ((pow(sin(theta), 2) * C[0, 0] + pow(cos(theta), 2) * C[1, 1] )) ** (0.5)
return Ei(LMDA, 10, 0)
R = trapezoidal(G, -np.pi / 2., np.pi / 2., 10)
# Total potential energy
psi = (c_1 * (Ic - 3.) + R + p * (J - 1)) * dx
F1 = derivative(psi, w, TestFunction(W))
Jac = derivative(F1, w, TrialFunction(W))
bc_left = DirichletBC(W.sub(1), Constant((0., 0.0, 0.0)), boundaries, 1)
bc_right = DirichletBC(W.sub(1), Constant((0.5, 0.0, 0.0)), boundaries, 2)
bc_all = [bc_left,bc_right]
problem = NonlinearVariationalProblem(F1, w, bc_all, Jac)
solver = NonlinearVariationalSolver(problem)
solver.solve()
(p, u) = w.split(True)
Is there any way to speed up this code?
Thanks in advance for your help.
