Need suggestions on how to implement this time stepping for wave equation

Question

I have the following system of equations obtained by implementing Sympletic Euler time scheme to wave equation. I want to model this in Fenics. Here 'u' is the displacement and 'p' is corresponding velocity.

$$ \dot{u} = p $$
$$ u(x, 0) = u_i$$
$$ p(x, 0) = 0$$
$$ u(x, t) = 0 \quad \text{on} \quad \partial \Omega$$

I read in a documentation that the velocity 'p' should be obtained first and thereafter 'u' should be computed. Also, that it must be ensured that we never have to explicitly compute and store the inverse of the mass matrix M.

Here, $M = uvdx$ and $K = \nabla u \nabla v dx$ .

Is this the correct way to implement the above system of equations? Can anyone kindly guide on how this can be corrected / refined.

u = interpolate(ui,V)
u_p = function(V)
u_new = function(V)
p = interpolate(Constant(0.0),V)

while t =< T
    bc.apply(M,u)        # apply boundary conditions    
    solve(M, u_p.vector(), u.vector)    
    p_new = p.vector() - dt*K*u_p.vector() 
    t += dt              # solve for next timestep
    u_new.vector()[:] = u.vector() + dt*p_new.vector()
    u.assign(u_new)      # update old vector
    p.vector()[:] = p_new

Comments

Just tried implementing it.....error: bc.apply(M,u.vector) # apply boundary conditions
TypeError: in method 'DirichletBC_apply', argument 3 of type 'dolfin::GenericVector &'