Combining continuous and discontinuous elements and receiveing an error "All terms in form must have same rank."

Question

Hello,

I'm working on a variant of the an elastostatic problem where I would like to include the effect of residual stress in the formulation. I've based my UFL on examples in the book and elsewhere on the web. I have introduced a new material property for the residual strain, and incorporated it into the equations. Unfortunately, I'm getting an error, which is:

All terms in form must have same rank.

Looking at the equation for epsilon, the new term is completely on the right. When added, the error is that the rank does not match. However, I have verified that the rank does indeed match. I suspect that the error in fact has to do with the automatic indexing used with the trial and test functions, but it is still not clear how the error should be corrected.

The following compiles without error:

cell = tetrahedron
element = VectorElement("Lagrange", cell, 1 )
property = FiniteElement("DG", cell, 0 )

v = TestFunction(element)
u = TrialFunction(element) 
f = Coefficient(element)

E = Coefficient(property)
nu = Coefficient(property)
epsr = Coefficient(property)

mu = E / (2*(1+nu))
lmbda = E*nu / ((1+nu)*(1-2*nu))

def epsilon(v):
    return 0.5*(grad(v) + transpose(grad(v)))

def sigma(v):
    return 2*mu*epsilon(v) + lmbda*(tr(epsilon(v)) * Identity(cell.d))

a = inner( sigma(u), grad(v) )*dx
L = dot( f, v ) * dx

Changing the definition of epsilon to the following causes the error:

def epsilon(v):
    return 0.5*(grad(v) + transpose(grad(v))) + epsr*Identity(cell.d)

Answer 1

Hi, the problem is that as a result of the new definition of epsilon there will be terms in a which have no trial function but only test function. You could isolate them and add to L but it is simpler to let UFL figure this out for you. The following compiles okay

cell = tetrahedron
element = VectorElement("Lagrange", cell, 1 )
property = FiniteElement("DG", cell, 0 )

v = TestFunction(element)
u = TrialFunction(element)
f = Coefficient(element)

E = Coefficient(property)
nu = Coefficient(property)
epsr = Coefficient(property)

mu = E / (2*(1+nu))
lmbda = E*nu / ((1+nu)*(1-2*nu)) 

def epsilon(v):
    return 0.5*(grad(v) + transpose(grad(v))) + epsr*Identity(cell.d)

def sigma(v):
    return 2*mu*epsilon(v) + lmbda*(tr(epsilon(v)) * Identity(cell.d))

F = inner( sigma(u), grad(v) )*dx + dot( f, v ) * dx
a = lhs(F)
L = rhs(F)

Comments

Thank-you. That was very useful.