 # Strain energy function with additional integral term

Hi
I am solving a Hyperelasticity problem with a strain energy function in this form:

The \mu and \lambda are Lame constants. The last integral is an additional term. The term p(\theta) is defined as: where \alpha is a constant. The parameter \beta is defined as \beta=a_0 \cdot C \cdot a_0

C is the right Cauchy deformation tensor and a_0 is function of \theta: The parameter \beta could be rewritten:

The \beta in integrand (last term in the first equation) includes the components of C tensor . If I ignore the \beta and keep only p(\theta) in the integrand and use quad from scipy.integrate, it works fine but I need to have \beta included! In this case scipy.integrate does NOT work anymore.

Is there any idea how this integration should be implemented based on the above explanation?

Here is my implementation:

from dolfin import *
import numpy as np
import matplotlib.pyplot as plt

n = 1
mesh = BoxMesh(Point(0, 0, 0), Point(1, 1, 1), n, n, n)

#define boundaries
class left(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x, 0.0)

class right(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-14
return (on_boundary and (abs(x - 1.0) < tol))

Left = left()
Right = right()

boundaries = MeshFunction('size_t', mesh, mesh.topology().dim()-1)
boundaries.set_all(0)

#Marking boundaries
Left.mark(boundaries, 1)
Right.mark(boundaries, 2)

domains = MeshFunction('size_t', mesh, mesh.topology().dim())
ds = Measure("ds", subdomain_data=boundaries)
dx = Measure('dx',subdomain_data=domains)

n = FacetNormal(mesh)
V = VectorFunctionSpace(mesh, 'CG', degree=1)

bc = DirichletBC(V, Constant((0., 0.,0.)),boundaries, 1)

du  = TrialFunction(V)           # Trial function
v  = TestFunction(V)             # Test function
u  = Function(V)                 # Displacement field

d = u.geometric_dimension()
I = Identity(d)             # Identity tensor
C = F.T*F                   # Right Cauchy-Green tensor
b = F * F.T                 # Left Cauchy-Green tensor
# Invariants of deformation tensors
Ic = tr(C)
J  = det(F)

E, nu = 1E6, 0.46
mu = 0.4e6
lmbda= Constant(E*nu/((1 + nu)*(1 - 2*nu)))

force  = Constant((1.0,  0.0, 0.0))
alpha = 1.26

# DefinE BETA
a_0 = as_vector([[1, 0,0]])
BETA = inner(dot(a_0,C),a_0)

def integrand(Theta):

BETA = pow(np.cos(Theta),2)*C[0,0] + pow(np.sin(Theta),2)*C[1,1] + np.cos(Theta)*np.sin(Theta)*C[0,1] + np.cos(Theta)*np.sin(Theta)*C[1,0]

return ((1. / (np.pi * alpha)) * np.exp(np.cos(2. * (Theta))))*BETA  # It works if BETA is removed

P_int = quad(integrand, - np.pi/2., np.pi/2.)
VAL = list(P_int)
#print (VAL)

psi = ((lmbda / 2.)*pow(ln(J),2)+0.5*mu*(Ic-3) + VAL)*dx - dot(force, u)*ds(2)
F1 = derivative(psi, u, v)

# Compute Jacobian
Jac = derivative(F1, u, du)
problem = NonlinearVariationalProblem(F1, u, bc, Jac)

solver = NonlinearVariationalSolver(problem)
solver.solve()

You might try SciPy’s fixed_quad function instead of quad, and find a suitable quadrature degree empirically. The following quick test looks promising:

from dolfin import *

# Non-numerical function returning UFL, for testing:
def func(theta):
return theta*Constant(1.0)

# For unclear reasons, fixed_quad returns a tuple, where the second component
# is always(?) None, but the first argument is what one expects, and does
# not try to cast the output of func or its own return value to type float.

# The result remains a UFL type, not a numerical value:
print(int_func)

1 Like

Thank you. It was really promising and solved this issue.
I was wondering if there is a way to define the parameter \beta in terms of an exponential integral function. To be more specific I am trying to include Ei (\beta) in the strain energy function. The exponential integral function is available in scipy by calling: from scipy.special import expi
Here is an example:

from dolfin import *
import numpy as np
from scipy.special import expi

n = 1
mesh = BoxMesh(Point(0, 0, 0), Point(1, 1, 1), n, n, n)

#define boundaries
class left(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x, 0.0)

class right(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-14
return (on_boundary and (abs(x - 1.0) < tol))

Left = left()
Right = right()

boundaries = MeshFunction('size_t', mesh, mesh.topology().dim()-1)
boundaries.set_all(0)

#Marking boundaries
Left.mark(boundaries, 1)
Right.mark(boundaries, 2)

domains = MeshFunction('size_t', mesh, mesh.topology().dim())
ds = Measure("ds", subdomain_data=boundaries)
dx = Measure('dx',subdomain_data=domains)

n = FacetNormal(mesh)
V = VectorFunctionSpace(mesh, 'CG', degree=1)

bc = DirichletBC(V, Constant((0., 0.,0.)),boundaries, 1)

du  = TrialFunction(V)           # Trial function
v  = TestFunction(V)             # Test function
u  = Function(V)                 # Displacement field

d = u.geometric_dimension()
I = Identity(d)             # Identity tensor
C = F.T*F                   # Right Cauchy-Green tensor
b = F * F.T                 # Left Cauchy-Green tensor
# Invariants of deformation tensors
Ic = tr(C)
J  = det(F)

E, nu = 1E6, 0.46
mu = 0.4e6
lmbda= Constant(E*nu/((1 + nu)*(1 - 2*nu)))

force  = Constant((1.0,  0.0, 0.0))
alpha = 1.26

def func(Theta):
BETA = expi(np.cos(Theta)*C[0,0])  # I have simplified BETA just as an example
return ((1. / (np.pi * alpha)) * np.exp(np.cos(2. * (Theta))))*BETA


It returns: ufunc 'expi' not supported for the input types