Remembering my time as a postgraduate student I think you have to handle this kind of problem using "direct methods of variational calculus". Here you try to minimize the functional itself, i.e. the integrals, by minimizer-sequences contained in special Banach spaces (mostly equipped with weak*-topology). These minimizer-sequences converge to the desired solution usually in the weak-sense. Hence the solution is then called "weak solution" in contrast to the classical "strong solutions" which are smooth/differentiable and satisfy the Euler-Lagrange equations. I have attached some nicely readable lecture notes giving some inside on these direct methods.

Remembering my time as a postgraduate student I think you have to handle this kind of problem using "direct methods of variational calculus". Here you try to minimize the functional itself, i.e. the integrals, by minimizer-sequences contained in special Banach spaces (mostly equipped with weak*-topology). These minimizer-sequences converge to the desired solution usually in the weak-sense. Hence the solution is then called "weak solution" in contrast to the classical "strong solutions" which are smooth/differentiable and satisfy the Euler-Lagrange equations. I have attached some nicely readable lecture notes giving some inside on these direct methods.

The notes by Yan which Kreer has provided are quite nice and are relevant if appropriately used, noting that they refer only to the steady-state situation so do not immediately give equations of motion as desired.  It is easier to see the issues involved when looking at a scalar problem, rather than a continuum so suppose one has L=T-V with T=(m/2)y'^2 (usual kinetic energy with mass m) and the potential V=y^2.  One then gets the equation of motion (my')'=-2y (with force -dV/dy) -- a standard harmonic oscillator. [Using (my')' rather than my'' follows Newton and allows for a possibility of time-varying mass but we now ignore this and take m=1.]  If we would have V=|y|, then dV/dy is discontinuous and we need the considerations of Filippov to see what happens, especially in higher-dimensional cases with discontinuity of dV/dt across a more general switching surface.  For V=|y| a solution starting with y(0)=-1 would have to be y(t)=t-1 for t

Depends on which of the variables it is discostinuous. If L=L(x,u,p) with p the argument for Du(x), then just measurable in x is fine. If L is convex wrt p, which is necessary, it is OK in the scalar case, you just replace derivatives by sub-differentials. Non-smoothness wrt the u argument might be an issue, unless you have some asumption for this one too. Check standard sources on non-smooth problem: Clarke, Rockafellar, Temam, etc...

Experiments with a simple pendulum problem with additional jump in the potential energy show that the question actually makes sense; see the attached example in Mathematica. The computed solution (and especially the phase portrait) lead to a conclusion that a proper jump condition needs to be formulated for the trajectories, which cross the discontinuities in the potential energy. Integrating such equations over time would require detecting jumps and continuing the solution into the subsequent region with appropriate initial conditions. This can, however, be avoided by regularizing the problem using a small parameter in a fashion, which is similar to the provided example (provided that a suitable time integration scheme is used).