Can I safely evaluate an integral (ufl assemble_scalar) of a quantity with 2nd order spatial derivatives?

As mentioned in this post, I am solving for “temp” and “volt” which are defined with CG polynomials of degree 4. My weak form does not involve 2nd order spatial derivatives, as it would be a variational crime, which I remedied by integrating by parts.

I am, however, in needs to evaluate an integral over the mesh and it involves 2nd order spatial derivatives. In code, it would be assemble_scalar(form(- temp * (S_xx * J_vector[0].dx(0) + S_yy * J_vector[1].dx(1)) * dx) where J_vector = -sigma * grad(volt) - sigma * Seebeck_tensor * grad(temp) involves the gradient of both temp and volt, therefore the .dx(0) and .dx(1) are actually 2nd order spatial derivatives. I know terms like this should not enter the weak form, but can they enter the computation of an integral via ufl.assemble_scalar?

As long as the space is sufficiently high order, second order derivatives can enter assemble_scalar.
If you have a fourth order Lagrange space the second order derivative is a second order discontinuous Lagrange function, which can be integrated over for instance cells, with no issues.