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
from scipy.integrate import quad
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.0)
class right(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-14
return (on_boundary and (abs(x[0] - 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
F = I + grad(u) # Deformation gradient
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)[0]
#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()
```