Update (for reference): I have solved my problems using method 4 (above), with the small python function below. So I'm happy. Problems 1 and 3 may still be worth fixing, though.
==
def delta(V, pt):
"""Unit area delta function in discrete space V"""
q = Function(V)
PointSource(V, pt).apply(q.vector())
q /= assemble(q*dx)
return q
Update (for reference): I have solved my problems using method 4 (above), with the small python function below. So I'm happy. Problems 1 and 3 may still be worth fixing, though.
==
def delta(V, pt): q.vector( ))
"""Unit area delta function in discrete space V"""
q = Function(V)
PointSource(V, pt).apply(
q /= assemble(q*dx)
return q
...
q_s = delta(S, Point(0.0))