Quadrature elements are currently a bit buggy, but can be used by using some deprecated functionality as follows:

```
from dolfin import *
from ufl import Max
# The following shuts off a deprecation warning for quadrature representation:
import warnings
from ffc.quadrature.deprecation \
import QuadratureRepresentationDeprecationWarning
warnings.simplefilter("ignore", QuadratureRepresentationDeprecationWarning)
# A consistent quadrature degree must be used throughout, to avoid mismatch
# between the quadrature space and quadature of the variational form:
QUAD_DEG = 4
dx = dx(metadata={"quadrature_degree":QUAD_DEG})
mesh = UnitSquareMesh(10,10)
# Function space with DoFs at quadrature points:
Qelement = FiniteElement("Quadrature", mesh.ufl_cell(),
degree=QUAD_DEG,
quad_scheme="default")
Q = FunctionSpace(mesh, Qelement)
# Function to store previous max strain energy:
psi_max_old = Function(Q)
# Placeholder, to be replaced with whatever the expression for strain
# energy is in the problem:
psi = SpatialCoordinate(mesh)[0]
# Use this in the variational problem where psi_max is needed:
psi_max = Max(psi,psi_max_old)
# After solving in each load step, update psi_max_old. When the unknowns of
# a variational problem are in a Quadrature-type space, you need to use the
# deprecated quadrature representation.
psi_max_old.assign(project(psi_max,Q,
form_compiler_parameters
={"quadrature_degree":QUAD_DEG,
"representation":"quadrature"}))
# Plot (projection onto linears) to verify:
from matplotlib import pyplot as plt
plot(project(psi_max_old,FunctionSpace(mesh,"Lagrange",1),
form_compiler_parameters={"quadrature_degree":QUAD_DEG}))
plt.show()
```

Also check out Jeremy Bleyerâ€™s plasticity tutorial, which uses quadrature spaces:

https://comet-fenics.readthedocs.io/en/latest/demo/2D_plasticity/vonMises_plasticity.py.html