TSFC compiler (AssertionError when using "ds")

variational formulation
1

Hello
I am solving a hyperelasticity problem. For the sake of speed, I want to use the TSFC compiler. I am getting an error after including ds when I want to apply a load on the right boundary. Here is the code:

from dolfin import *

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

# 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

mesh = UnitCubeMesh.create(1,1,1,CellType.Type.hexahedron)

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

left = Left()
right = Right()

du  = TrialFunction(V)
v  = TestFunction(V)
u  = Function(V)

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': 2})
ds = Measure('ds', domain=mesh, subdomain_data=boundaries, metadata={'quadrature_degree': 2})

d = u.geometric_dimension()
I = Identity(d)
F = I + grad(u)
C = F.T*F
Ic = tr(C)
J  = det(F)

mu = 20
K = 1E6

force  = Constant((0.1,  0.0, 0.0))
psi = (mu* (Ic - 3.) ) * dx + 0.5 *  K * pow(ln(J),2) * dx - dot(force, u)*ds(2)

F1 = derivative(psi, u, v)
# Compute Jacobian of F
Jac = derivative(F1, u, du)

bc_left = DirichletBC(V, Constant((0,0,0)), boundaries , 1)
problem = NonlinearVariationalProblem(F1, u, bc_left, Jac)

solver = NonlinearVariationalSolver(problem)
solver.solve()

The above code fails with “AssertionError”. If I comment out the second line (parameters["form_compiler"]["representation"] = "tsfc") it works. It also works if I keep the “tsfc” and replace ds(2) with for example dx but I want to use the “tsfc” and keeping "ds(2) simultaneously.
Does anybody know where the problem is?

2

Although maybe not a “solution”, the code works with the following modification:

#mesh = UnitCubeMesh.create(1,1,1,CellType.Type.hexahedron)
mesh = UnitCubeMesh.create(1,1,1,CellType.Type.tetrahedron)

The error you’re getting may just be a manifestation of the current incomplete support for hex elements in FEniCS.

3

The strange thing is that by removing the second line the code works (even by keeping the hex element). However the code becomes very slow.
Anyways, thanks for your help.