Note that the weak form \int_\Omega k\nabla u\cdot\nabla v\,d\mathbf{x} corresponds to the strong form -\nabla\cdot(k\nabla u), which is not the same as -k\nabla^2 u when k varies in space. Likewise, when c_p^{-1} varies in space, it cannot be pulled inside the divergence operator, so
where the left side is what you get by dividing the original conservation law through by c, but the right side is what the weak form in the code corresponds to when normal == False.
While k and c_p^{-1} are piecewise-constant, there are Dirac deltas in their spatial derivatives at the discontinuity. (You could pull c_p^{-1} inside the divergence operator on either side of the internal interface where c_p and k are discontinuous, but then when you integrate by parts on either side to get a weak form, there would be non-cancelling boundary terms at the interface.)