This is a read only copy of the old FEniCS QA forum. Please visit the new QA forum to ask questions

Error during assembly while using determinant of the deformation tensor F

+1 vote

I am trying a solve a simple hyperelastic model where the energy is given by J*\bs{E}:\bs{E}, where E is the small strain tensor and J is determinant of F = \bs{E} + \bs{I}.

from dolfin import *
mesh = UnitSquareMesh(2, 2)
S = FunctionSpace(mesh, 'CG', 1)
V_u = VectorFunctionSpace(mesh, 'CG', 1)
T = TensorFunctionSpace(mesh, 'CG', 1)

#----------------------------------------------------------------------------
# Define boundary sets for boundary conditions
#----------------------------------------------------------------------------
class Left(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1E-5 # tolerance for coordinate comparisons
        return on_boundary and abs(x[0] - 0) < tol

class Right(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1E-5 # tolerance for coordinate comparisons
        return on_boundary and abs(x[0] - 1) < tol
# Initialize sub-domain instances
left = Left() 
right = Right()
# define meshfunction to identify boundaries by numbers
boundaries = FacetFunction("size_t", mesh)
boundaries.set_all(0)
left.mark(boundaries, 1) # mark left as 1
right.mark(boundaries, 2) # mark right as 2
# Define new measure including boundary naming 
ds = Measure("ds")[boundaries] # left: ds(1), right: ds(2)   
def eps(v):
    return sym(grad(v))
# Define function spaces and test, trial functions
u, du, v = Function(V_u), TrialFunction(V_u), TestFunction(V_u)

# Define boundary condiition                                  
zero = Constant(0.0)
zero2 = Constant((0.0,0.0))
# Dirichlet boundary condition for a traction test boundary
u_L = zero2
u_R = Expression(("0.0","0.02"),degree=1) 
# Dispalcement boundary conditions
bc1 = DirichletBC(V_u, u_L, boundaries, 1) 
bc2 = DirichletBC(V_u, u_R, boundaries, 2) 
bc_u = [bc1, bc2]
# Formulation
E = sym(grad(u))
F = Identity(2) + E
J = det(F)
total_energy = J*(inner(E,E))*dx 
E_u = derivative(total_energy,u,v)
E_u_u = derivative(E_u,u,du)
E_du = replace(E_u,{u:du})
bc_u = []    
# Variational problem for the displacement
problem_u = LinearVariationalProblem(lhs(E_du), rhs(E_du), u, bc_u)

solver_u = LinearVariationalSolver(problem_u)
file_u = File("u.pvd")
solver_u.solve()
file_u << u

If i define the
total_energy = inner(E,E)dx there is no error
while if i define
total_energy = J
inner(E,E)*dx then i get error in the assembly

A part of the error is as follows

Calling FFC just-in-time (JIT) compiler, this may take some time.
Traceback (most recent call last):
  File "Post.py", line 54, in <module>
    problem_u = LinearVariationalProblem(lhs(E_du), rhs(E_du), u, bc_u)
  File "/usr/local/lib/python2.7/dist-packages/dolfin/fem/solving.py", line 82, in __init__
    a = Form(a, form_compiler_parameters=form_compiler_parameters) 

I feel that using determinant of a tensor in the formulation is causing the error. I would appreciate any help on this.

Thanks
Kaushik

asked Feb 3, 2017 by kaushikv123 FEniCS Novice (390 points)

1 Answer

+2 votes

Since the energy functional is not a quadratic form you should use the classes NonlinearVariationalProblem and NonlinearVariationalSolver.

answered Feb 3, 2017 by umberto FEniCS User (6,440 points)

Thanks umberto. It works when i use nonlinearsolver.

...