# Slow code (function approximation in hyperelasticity)

**URL:** <https://fenicsproject.discourse.group/t/slow-code-function-approximation-in-hyperelasticity/2419>\
**Category:** variational formulation\
**Created:** [February 12, 2020, 7:28am UTC](https://fenicsproject.discourse.group/t/slow-code-function-approximation-in-hyperelasticity/2419 "2020-02-12T07:28:06Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![Leo](https://avatars.discourse-cdn.com/v4/letter/l/c2a13f/32.png) [@Leo](https://fenicsproject.discourse.group/u/Leo)\
**Post date:** [February 12, 2020, 7:28am UTC](https://fenicsproject.discourse.group/t/slow-code-function-approximation-in-hyperelasticity/2419/1 "2020-02-12T07:28:06Z")

</div>

Hello  
I am solving a problem in hyperelasticity where I have a new term in strain energy function which is defined as:

![image](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/1X/e9306a9db6a6166ae7493f76570c05271cf1d676.png)  
In the above equation:

![image](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/1X/755780968ba1d43631e3c72069ebf2c7e6c3ecda.png)

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: ![image](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/1X/5c64aa6a3b08d6a4b3c7e0870a3032bbe84f10c1.png). Gamma is a constant and **Ein** represents _ **entire function** _ that is defined as: ![image](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/1X/a140ad98831c6ede23683f6ce89713878348ef14.png) . (The whole approximation could be found [here](https://en.wikipedia.org/wiki/Exponential_integral). 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.

---

<div class="post-metadata">

**Author:** ![kamensky](https://avatars.discourse-cdn.com/v4/letter/k/e95f7d/32.png) [@kamensky](https://fenicsproject.discourse.group/u/kamensky)\
**Post date:** [February 12, 2020, 3:41pm UTC](https://fenicsproject.discourse.group/t/slow-code-function-approximation-in-hyperelasticity/2419/2 "2020-02-12T15:41:41Z")

</div>

First, I doubt there is much benefit from using 10-th degree quadrature in space (i.e., the degree of `dx` and `ds`, not the integrals over \theta at each stress point) with P2–P1 Taylor–Hood elements, even if the constitutive model is quite complicated. You can get some savings by reducing to degree 4 quadrature in space.

However, the much bigger cost saving is from switching to the experimental [TSFC](https://epubs.siam.org/doi/10.1137/17M1130642) form representation. TSFC is not installed by default, but you can install Jan Blechta’s FEniCS port of it like

```bash
pip3 install git+https://github.com/blechta/tsfc.git@2018.1.0
pip3 install git+https://github.com/blechta/COFFEE.git@2018.1.0
pip3 install git+https://github.com/blechta/FInAT.git@2018.1.0
pip3 install singledispatch networkx pulp

```

in my experience, this is usually much faster and more robust for quad/hex elements, and also gets around some bugs in the default UFLACS representation related to the use of quadrature function spaces. Once it’s installed you can invoke TSFC in your script with

```python
parameters["form_compiler"]["representation"] = "tsfc"

```

right after importing DOLFIN. I also find that, for some complicated forms, I need to manually increase Python’s recursion limit when using TSFC. You can do that with

```python
import sys
sys.setrecursionlimit(10000)

```

---

<div class="post-metadata">

**Author:** ![Leo](https://avatars.discourse-cdn.com/v4/letter/l/c2a13f/32.png) [@Leo](https://fenicsproject.discourse.group/u/Leo)\
**Post date:** [February 13, 2020, 12:44am UTC](https://fenicsproject.discourse.group/t/slow-code-function-approximation-in-hyperelasticity/2419/3 "2020-02-13T00:44:50Z")

</div>

Thank you for the solution. After installing TSFC the speed of code increased significantly. However changing the quadrature degree from 10 to 4 affects the results a little bit.

---

<div class="post-metadata">

**Author:** ![Leo](https://avatars.discourse-cdn.com/v4/letter/l/c2a13f/32.png) [@Leo](https://fenicsproject.discourse.group/u/Leo)\
**Post date:** [February 13, 2020, 6:52am UTC](https://fenicsproject.discourse.group/t/slow-code-function-approximation-in-hyperelasticity/2419/4 "2020-02-13T06:52:56Z")

</div>

I have 1 more question if you do not mind me asking. I changed the above code a little bit by adding 1 more term to the strain energy:

```
K = 1000000
psi = (c_1 * (Ic - 3.) + R + p * (J - 1)) * dx+ K * pow(ln(J), 2) * dx

```

and also I defined the boundary conditions in a pointwise manner. I want to fix the nodes on the left side in the x-direction and apply displacement in x-direction on the nodes on the right side. So instead of:

```
bc_left = DirichletBC(W.sub(1), Constant((0., 0.0, 0.0)), boundaries, 1)
bc_right = DirichletBC(W.sub(1), Constant((0.6, 0.0, 0.0)), boundaries, 2)
bc_all = [bc_left, bc_right]

```

I want to do this:

```
def point_left_1(x, on_boundary):
    tol = DOLFIN_EPS
    return (abs(x[0] - 0.0) < tol) and (abs(x[1] - 0) < tol) and (abs(x[2] - 0) < tol)

BC_POINT_left_1 = DirichletBC(W.sub(1), Constant((0., 0., 0.)), point_left_1, method="pointwise")

def point_left_2(x, on_boundary):
    tol = DOLFIN_EPS
    return (abs(x[0] - 0.0) < tol) and (abs(x[1] - 0) < tol) and (abs(x[2] - 1) < tol)

BC_POINT_left_2 = DirichletBC(W.sub(1), Constant((0., 0., 0.)), point_left_2, method="pointwise")

def point_left_3(x, on_boundary):
    tol = DOLFIN_EPS
    return (abs(x[0] - 0.0) < tol) and (abs(x[1] - 1) < tol) and (abs(x[2] - 0) < tol)

BC_POINT_left_3 = DirichletBC(W.sub(1), Constant((0., 0., 0.)), point_left_3, method="pointwise")

def point_left_4(x, on_boundary):
    tol = DOLFIN_EPS
    return (abs(x[0] - 0.0) < tol) and (abs(x[1] - 1) < tol) and (abs(x[2] - 1) < tol)

BC_POINT_left_4 = DirichletBC(W.sub(1), Constant((0., 0., 0.)), point_left_4, method="pointwise")

def point_right_1(x, on_boundary):
    tol = DOLFIN_EPS
    return (abs(x[0] - 1.0) < tol) and (abs(x[1] - 0) < tol) and (abs(x[2] - 0) < tol)

BC_POINT_right_1 = DirichletBC(W.sub(1), Constant((1.0, 0., 0.)), point_right_1, method="pointwise")

def point_right_2(x, on_boundary):
    tol = DOLFIN_EPS
    return (abs(x[0] - 1.0) < tol) and (abs(x[1] - 1.0) < tol) and (abs(x[2] - 0) < tol)

BC_POINT_right_2 = DirichletBC(W.sub(1), Constant((1., 0, 0.)), point_right_2, method="pointwise")

def point_right_3(x, on_boundary):
    tol = DOLFIN_EPS
    return (abs(x[0] - 1.0) < tol) and (abs(x[1] - 0) < tol) and (abs(x[2] - 1.) < tol)

BC_POINT_right_3 = DirichletBC(W.sub(1), Constant((1.0, 0., 0.)), point_right_3, method="pointwise")

def point_right_4(x, on_boundary):
    tol = DOLFIN_EPS
    return (abs(x[0] - 1.0) < tol) and (abs(x[1] - 1.) < tol) and (abs(x[2] - 1.) < tol)

BC_POINT_right_4 = DirichletBC(W.sub(1).sub(2), Constant((0.)), point_right_4, method="pointwise")

bc_all = [BC_POINT_left_1, BC_POINT_left_2, BC_POINT_left_3, BC_POINT_left_4 ,BC_POINT_right_1, BC_POINT_right_2, BC_POINT_right_3, BC_POINT_right_4]

```

But it diverges.I played with different parts of the code to figure it out. I found that If I use `K = 100` or setting `'quadrature_degree': 2` it converges (However this is not what I want as the results are wrong). Is there any way to fix this convergence issue by any chance?  
Thanks again!

---

<div class="post-metadata">

**Author:** ![kamensky](https://avatars.discourse-cdn.com/v4/letter/k/e95f7d/32.png) [@kamensky](https://fenicsproject.discourse.group/u/kamensky)\
**Post date:** [February 13, 2020, 3:32pm UTC](https://fenicsproject.discourse.group/t/slow-code-function-approximation-in-hyperelasticity/2419/5 "2020-02-13T15:32:39Z")

</div>

I believe the K term in this formulation is misguided. Note that there is already a Lagrange multiplier, p, enforcing the constraint J=1. If that constraint is satisfied, \ln J = 0, so the K term has no influence on the exact solution to the PDE system. However, for large K, this term is trying to enforce J=1 _pointwise_ in the discrete solution, whereas the Lagrange multiplier p is enforcing the constraint weighted by every test function in the \text{CG}\_1 pressure space. Trying to enforce pointwise (near-)incompressibility in the discrete problem is typically unstable. In the linearized problem, this instability can be formalized in terms of the inf-sup/Babuška–Brezzi condition on the variational problem. For example, the Taylor–Hood element you are using is well-known to be inf-sup stable for the Stokes flow problem, which is formally identical to linearized incompressible isotropic elasticity. Typically, the pressure space with respect to which incompressibility is enforced needs to be lower-order than the displacement/velocity space. For more mathematical detail on inf-sup conditions, some thorough lecture notes can be found [here](https://www.oden.utexas.edu/media/reports/2006/0608.pdf).

A classic observation that you’ve re-discovered is that reduced quadrature seems to improve stability when penalizing compression pointwise. For some cases, it turns out that using reduced quadrature for just the dilational terms in a formulation (i.e., “selective reduced quadrature”) can be re-framed as discretizing pressure in a lower-order space ([https://doi.org/10.1016/0045-7825(78)90005-1](https://doi.org/10.1016/0045-7825(78)90005-1)), which is the mathematical explanation for why it helps.

---

<div class="post-metadata">

**Author:** ![Leo](https://avatars.discourse-cdn.com/v4/letter/l/c2a13f/32.png) [@Leo](https://fenicsproject.discourse.group/u/Leo)\
**Post date:** [February 13, 2020, 4:06pm UTC](https://fenicsproject.discourse.group/t/slow-code-function-approximation-in-hyperelasticity/2419/6 "2020-02-13T16:06:06Z")

</div>

Thank you for this clarification. Strain energy I want to use comes from [here (section 4.1.2.17)](https://help.febio.org/FEBio/FEBio_um_2_8/) and it is in this form: ![ss](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/1X/df690f77c6e65ff8fae971434bf67137221b2657.png)  
I tried to simplify it in my question. It is absolutely right that in case of incompressibility, the **ln(J)** should be equal to zero. So this is also my question why they have it in the strain energy.  
In another paper, they say : **"To enforce incompressibiliy we took K = 100MPa"**  
Does it mean by keeping the last term in the strain energy function (Large K), we can remove the **p \* (J - 1)**?

---

<div class="post-metadata">

**Author:** ![bhaveshshrimali](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/bhaveshshrimali/32/124_2.png) [@bhaveshshrimali](https://fenicsproject.discourse.group/u/bhaveshshrimali)\
**Post date:** [February 13, 2020, 4:16pm UTC](https://fenicsproject.discourse.group/t/slow-code-function-approximation-in-hyperelasticity/2419/7 "2020-02-13T16:16:09Z")

</div>

I believe that when you use a mixed formulation the term K/2 (\ln J)^2 or equivalently K/2(J-1)^2 shouldn’t be present. The constraint J=1 is taken care of by the term with the lagrange multiplier p(J-1). For more clarity, you can start with a “pure-displacement” strain energy (with the K term) and take a partial legendre transform of that to arrive at a mixed formulation.

Edit: You can keep the K/2(\ln J)^2 term and instead of using “Taylor-Hood” elements use higher order lagrange elements (p=2,3) to prevent volumetric locking and (may) use K=1000. Among other things, this might be a good check for your displacements coming from a mixed formulation