Dear Lawrence,

Thanks so much! I have done as you suggested and now I have a different error:

Traceback (most recent call last):
  File "Elasticity2D_Firedrake.py", line 163, in <module>
    solve(lhs_varf == rhs_varf, sol, bcs=[bc1,bc2])
  File "/home/simone/firedrake/src/firedrake/firedrake/solving.py", line 125, in solve
    _solve_varproblem(*args, **kwargs)
  File "/home/simone/firedrake/src/firedrake/firedrake/solving.py", line 153, in _solve_varproblem
    solver.solve()
  File "/home/simone/firedrake/src/firedrake/firedrake/variational_solver.py", line 249, in solve
    dbc.apply(self._problem.u)
  File "</home/simone/firedrake/lib/python3.6/site-packages/decorator.py:decorator-gen-48>", line 2, in apply
  File "/home/simone/firedrake/src/PyOP2/pyop2/profiling.py", line 60, in wrapper
    return f(*args, **kwargs)
  File "/home/simone/firedrake/src/firedrake/firedrake/bcs.py", line 432, in apply
    r.assign(self.function_arg, subset=self.node_set)
  File "</home/simone/firedrake/lib/python3.6/site-packages/decorator.py:decorator-gen-40>", line 2, in assign
  File "/home/simone/firedrake/src/firedrake/firedrake/utils.py", line 61, in wrapper
    return f(*args, **kwargs)
  File "/home/simone/firedrake/src/firedrake/firedrake/function.py", line 377, in assign
    assemble_expressions.Assign(self, expr), subset)
  File "</home/simone/firedrake/lib/python3.6/site-packages/decorator.py:decorator-gen-45>", line 2, in evaluate_expression
  File "/home/simone/firedrake/src/firedrake/firedrake/utils.py", line 58, in wrapper
    return f(*args, **kwargs)
  File "/home/simone/firedrake/src/firedrake/firedrake/assemble_expressions.py", line 507, in evaluate_expression
    e, args, _ = ExpressionWalker().walk(tree)
  File "/home/simone/firedrake/src/firedrake/firedrake/assemble_expressions.py", line 346, in walk
    return (self.visit(expr),
  File "/home/simone/firedrake/src/ufl/ufl/algorithms/transformer.py", line 114, in visit
    r = h(o)
  File "/home/simone/firedrake/src/firedrake/firedrake/assemble_expressions.py", line 413, in operator
    operands = [operands[0], self.visit(operands[1])]
  File "/home/simone/firedrake/src/ufl/ufl/algorithms/transformer.py", line 114, in visit
    r = h(o)
  File "/home/simone/firedrake/src/firedrake/firedrake/assemble_expressions.py", line 417, in operator
    operands = list(map(self.visit, o.ufl_operands))
  File "/home/simone/firedrake/src/ufl/ufl/algorithms/transformer.py", line 114, in visit
    r = h(o)
  File "/home/simone/firedrake/src/firedrake/firedrake/assemble_expressions.py", line 383, in coefficient
    raise ValueError("Constant has mismatched shape for expression function space")
ValueError: Constant has mismatched shape for expression function space

I kept the function spaces as the original:

BDM = VectorFunctionSpace(mesh, "BDM", 1)
DGv = VectorFunctionSpace(mesh, "DG", 0)
DGs = FunctionSpace(mesh, "DG", 0)

Vh = BDM * DGv * DGs

Is this error due to the boundary conditions that I prescribed?

u0 = Constant((0.0, 0.0))                                                                           # zero displacement left, right and bottom

zero_stress = as_vector([Constant((0., 0.)), Constant((0., 0.))])                    # zero stress at top

Thank you so much for your help, I really appreciate it.

Best regards,

Simone