Least effort to handle a point source inside the domain of PDE(s)

The good news is that yes, there is an easy way to put your problem into NDSolve by using the new finite element functionality in v10. The bad news is that it seems the specific problem you're trying to solve is ill-posed.

---

NDSolve can now handle internal boundaries; see e.g. the first figure under "Details" for DirichletCondition. Generating a mesh with such internal boundaries is described in the "Element Mesh Generation" tutorial. I don't know if a single constrained point technically counts as a "boundary", but it seems to work.

Create a spatial mesh with a node at the point source:

Needs["NDSolve`FEM`"];
bmesh = ToBoundaryMesh[
   "Coordinates" -> {{-1, -1}, {-1, 1}, {1, 1}, {1, -1}, {0, 0}}, 
   "BoundaryElements" -> {LineElement[{{1, 2}, {2, 3}, {3, 4}, {4, 1}}]}];
mesh = ToElementMesh[bmesh];
Show[mesh["Wireframe"], Graphics[{Red, PointSize[Large], Point[{0, 0}]}]]

(One could also use the not-really-documented "IncludePoints" option, as in this other answer.)

Direcly specifying u[t, 0, 0] doesn't work, as you already know, but DirichletCondition does:

sol = NDSolve[{
   D[u[t, x, y], t, t] == D[u[t, x, y], x, x] + D[u[t, x, y], y, y], 
   DirichletCondition[u[t, x, y] == Sin[10 t], x == 0 && y == 0], 
   u[0, x, y] == 0,
   Derivative[1, 0, 0][u][0, x, y] == 0, 
   u[t, -1, y] == 0,
   u[t, 1, y] == 0,
   u[t, x, -1] == 0, 
   u[t, x, 1] == 0},
  u, {t, 0, 3}, {x, y} ∈ mesh];

It complains that "NDSolve has computed initial values that give a zero residual for the differential-algebraic system, but some components are different from those specified", which is to be expected. But it gives a solution anyway.

frames = Table[
   Plot3D[u[t, x, y] /. sol, {x, -1, 1}, {y, -1, 1}, 
    PlotRange -> {-1, 1}, Mesh -> None, PlotStyle -> White], {t, 0, 3, 0.05}];
Export["a.gif", frames];

---

We start to see a problem if we change the resolution of the mesh. To avoid the massive memory requirements of a uniformly refined mesh, it's better to refine only where the solution changes rapidly, i.e. in the neighbourhood of the point source. One can use

mesh = ToElementMesh[bmesh, 
   "MeshRefinementFunction" -> Function[{vertices, area}, 
     area > Max[a, 1*^-2 Min[Norm /@ vertices]]]];

which smoothly refines the mesh to have elements of area $a$ near the point source at the origin. Here are some meshes with $a=10^{-2}$, $10^{-4}$, and $10^{-6}$, followed by zooms to $[-0.1,0.1]\times[-0.1,0.1]$:

And here are the corresponding solutions at $t=1$:

The solutions seem to be getting weaker the finer we make the mesh. What's going on? I don't know for absolutely certain, but I'm guessing that the problem is ill-posed and the solution we've computed is essentially an artifact of the numerical discretization. As an analogy, consider the Laplace problem on a punctured domain with Dirichlet boundary conditions: $$\begin{align} \nabla^2f(x,y)&=0&\text{for }&x\in\Omega\setminus\{(0,0)\},\\ f(x,y)&=0&\text{for }&x\in\partial\Omega,\\ f(0,0)&=1. \end{align}$$ You can solve this numerically and obtain a reasonable-looking numerical solution, but it is an illusion because a one-point set has zero capacity for the Laplacian, and if you refine the mesh the solution goes to zero. I believe the same thing is happening here. Numerically, the energy that the source imparts to the system is mesh-dependent, being related to the area of its neighbouring elements. Theoretically, I guess there is no solution.

---

So yeah. Can you use NDSolve for this problem? You can. But... maybe you shouldn't.

Disclaimer: I am not a functional analyst and this is not mathematical advice. Consult your friendly neighbourhood applied mathematician.

The periodic driving at one point doesn't seem to be compatible with the boundary conditions expected by NDSolve, so I modified the problem in two ways: first, broaden the point source into a Gaussian, and then incorporate this driving as a source term in the actual differential equation.

So we're actually solving the inhomogeneous wave equation here. For the initial condition, I chose something that looks like the driving term spatially. The result is still not very satisfactory because NDSolve is much slower than your hand-coded finite-difference scheme.

σ = .1;

sol = u /. First@NDSolve[{
     D[u[t, x, y], t, t] == 
      D[u[t, x, y], x, x] + D[u[t, x, y], y, y] + 
       Sin[10 t] Exp[-(x^2 + y^2)/(2 σ^2)],
     Derivative[1, 0, 0][u][0, x, y] == 0,
     u[0, x, y] == .004 Exp[-(x^2 + y^2)/(2 σ^2)],
     u[t, -1, y] == 0,
     u[t, 1, y] == 0,
     u[t, x, -1] == 0,
     u[t, x, 1] == 0},
    u,
    {t, 0, 3},
    {x, -1, 1},
    {y, -1, 1}, 
    Method -> {"MethodOfLines", 
      "SpatialDiscretization" -> {"TensorProductGrid", 
        "MaxPoints" -> 100}}];

frames = 
  Table[Plot3D[sol[t, x, y], {x, -1, 1}, {y, -1, 1}, 
    PlotRange -> .02 {-1, 1}, PlotPoints -> 20, Mesh -> False], {t, 0,
     3, .05}];

ListAnimate[frames]

The parameter $\sigma$ dictates the spatial width of the Gaussian in the driving term, and I added the "Method" options explicitly so that we can control the value of "MaxPoints" if you decide to make $\sigma$ smaller than what I chose. The smaller $\sigma$, the larger "MaxPoints" needs to be. This makes NDSolve painfully slow if you really want to approach a well-localized point source. Of course, the plotting of the resulting InterpolatingFunction eats up some time, too. So in terms of speed, your discretized calculation is faster both in the actual differential-equation solving and in the plotting steps.

Nevertheless, the above at least shows how we can make NDSolve produce the desired result in principle.