Time-dependent Point Heat Source on Heat Equation

variational formulation
1

I try to model a situation that there are a number of heaters on the 2D rectangle room. At each time t, each heater produces the heat and the heat diffuses through the 2D space. Since the heater should be modeled as a point heat source, this is my current code:

from fenics import *
from mshr import *
import numpy as np

dt = 0.1
num_steps = 10

domain = Rectangle(Point(-1, -1), Point(1, 1))
mesh = generate_mesh(domain, 30)
V = FunctionSpace(mesh, 'P', 1)

def boundary(x, on_boundary):
    return on_boundary
bc = DirichletBC(V, Constant(0), boundary)

u_0 = Constant(0)
u_n = interpolate(u_0, V)

epsilon = 0.03
v = TestFunction(V)
t = 0

for i in range(num_steps):
    u = TrialFunction(V)
    f = Expression('0.0+1000*(i+1)*(-0.6 - epsilon<=x[0] & x[0]<= -0.6 + epsilon & -0.7 - epsilon<=x[1] & x[1]<=-0.7 + epsilon)', degree=2, i=i, epsilon=epsilon)
    F = u*v*dx + dt*dot(grad(u), grad(v))*dx - (u_n + dt*f)*v*dx
    a, L = lhs(F), rhs(F)
    u = Function(V)
    t += dt
    solve(a == L, u, bc)
    plot(u)
    u_n.assign(u)

Which means that there is one heater at (-0.6, -0.7) on [-1, 1]^2 2D box and the heater produce 1000*t heat at each time t. I heard that the class PointSource can handle such kind of situation, but I don’t clearly understand how to use that class. Is there any better idea to model this situation? Thanks!