How to construct functions of calculated solutions?

2

In this case, f is relatively simple, and the easiest approach might be to just use a Calculator filter in Paraview. If f is more complicated, or needed elsewhere in the FEniCS script, you might consider a few options, as demonstrated below:

from dolfin import *
mesh = UnitIntervalMesh(10)
el = FiniteElement("CG",mesh.ufl_cell(),1)
el_x_el = MixedElement(2*[el,])
V = FunctionSpace(mesh,el_x_el)
u = project(Constant((1,2)),V)
m,c = u.split()
scaling = 1.0

# Output using Expression.  This directly interpolates $f$ into a CG
# FunctionSpace of degree 1.
f_expr = Expression("scaling*c*pow(m+1.0,2)/(2*m)+c-1.0",
                    scaling=scaling,m=m,c=c,degree=1)
V_out = FunctionSpace(mesh,"CG",1)
f_out = Function(V_out)
f_out.interpolate(f_expr)
File("f_expr.pvd") << f_out

# Output using $L^2$ projection.  This is more flexible (if $f$ involves,
# e.g., derivatives, conditionals, etc.), but requires solving a linear system.
f = scaling*c*pow(m+1.0,2)/(2*m)+c-1.0
File("f_proj.pvd") << project(f, V_out)

# Lumped-mass projection.  This remains flexible and avoids solving a linear
# system, but is only an approximation of the $L^2$ projection, and requires
# an extra user-defined function.  This is also more robust than full $L^2$
# projection if $f$ is discontinuous (e.g., from a UFL conditional).
def lumped_project(f):
    v = TestFunction(V_out)
    lhs = assemble(inner(Constant(1.0),v)*dx)
    rhs = assemble(inner(f,v)*dx)
    u = Function(V_out)
    as_backend_type(u.vector())\
        .vec().pointwiseDivide(as_backend_type(rhs).vec(),
                               as_backend_type(lhs).vec())
    return u
File("f_lump.pvd") << lumped_project(f)