Hi

I am working on a hyperelasticity problem.The material model is same as the hyperelasticity demo.

Now I want to add a new term to the strain energy function which is : cos(\theta) .tr(C)

Where C is the Right Cauchy-Green tensor and \theta is defined in this figure:

In my implementation the strain energy function is defined as:

`psi = ((lmbda / 2.)*pow(ln(J),2)-mu*ln(J)+0.5*mu*(Ic-3))*x[0]*dx - inner(-p*n, u)*x[0]*ds(3) + cos(theta) * Ic * dx`

The **last term** is what I want to add to the \psi. As you can see, \theta is calculated with respect to the **center**.

I want to add the term cos(\theta) .tr(C) to the strain energy function **in a way that it is calculated with the orientation of the elements**. For example I want \theta to be zero at the bottom edge, 45 degree at the middle of the geometry and 90 degree at the left edge, etc. I mean it is supposed to change with orientation of the elements (With respect to the center).

What is the best way to do it in FEniCS?

Here is my minimal work:

**Note:** The geometry and problem definition has been taken from Numerical tours of continuum mechanics

```
from dolfin import *
from mshr import *
Re = 11.
Ri = 9.
rect = Rectangle(Point(0., 0.), Point(Re, Re))
domain = Circle(Point(0., 0.), Re, 100) - Circle(Point(0., 0.), Ri, 100)
domain = domain - Rectangle(Point(0., -Re), Point(-Re, Re)) \
- Rectangle(Point(0., 0.), Point(Re, -Re))
mesh = generate_mesh(domain, 10)
class Bottom(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], 0) and on_boundary
class Left(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 0) and on_boundary
class Outer(SubDomain):
def inside(self, x, on_boundary):
return near(sqrt(x[0]**2+x[1]**2), Re, 1e-1) and on_boundary
facets = MeshFunction("size_t", mesh, 1)
facets.set_all(0)
Bottom().mark(facets, 1)
Left().mark(facets, 2)
Outer().mark(facets, 3)
ds = Measure("ds", subdomain_data=facets)
x = SpatialCoordinate(mesh)
def GRAD(v):
return as_tensor([[v[0].dx(0), 0, v[0].dx(1)],
[0, v[0]/x[0], 0],
[v[1].dx(0), 0, v[1].dx(1)]])
E = Constant(1e6)
nu = Constant(0.3)
mu = E/2/(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)
n = FacetNormal(mesh)
p = Constant(10.)
V = VectorFunctionSpace(mesh, 'CG', degree=1)
du = TrialFunction(V) # Trial function
v = TestFunction(V) # Test function
u = Function(V) # Displacement field
F = Identity(3) + GRAD(u)
C = F.T*F # Right Cauchy-Green tensor
Ic = tr(C)
J = det(F)
#How should the last term (theta) be defined?
psi = ((lmbda / 2.)*pow(ln(J),2)-mu*ln(J)+0.5*mu*(Ic-3))*x[0]*dx - inner(-p*n, u)*x[0]*ds(3) + cos(theta) * Ic * dx
bcs = [DirichletBC(V.sub(1), Constant(0), facets, 1),
DirichletBC(V.sub(0), Constant(0), facets, 2)]
F1 = derivative(psi, u, v)
Jac = derivative(F1, u, du)
problem = NonlinearVariationalProblem(F1, u, bcs, Jac)
solver = NonlinearVariationalSolver(problem)
solver.solve()
```