I have some general questions about FE. On phone so can’t easily use latex.
1.I notice that in general, when we solve the heat equation for instance, we obtain the temperature distribution, possibly as a CG of order 1 FE. From it, we can compute the heat flux which involves grad(T). No matter how we post process T, it seems like we get a DG heat flux of order n1. Is this true?

Can fenicsx directly solve for the heat flux rather than the temperature? In this case the weak form looks simpler to me. What was Neumann b.c. (constant fluxes) becomes Dirichlet b.c. on the boundaries. And we can define the heat flux in a CG space, thereby completely defeating all approximations in 1., can this really be that easy? I.e. when we are interested in computing as accurately as possible both the temperature and the thermal flux, we could solve both PDEs separately? I have never seen this being done in any document online, I guess I am missing something.

What would be the drawbacks to solving for the flux instead of primary variable? What advantages?

If I also need to compute terms wtih spatial derivative of the flux, am I better to solve for.the flux defined in a CG space and then obtain the quantity in a DG space, or would two consecutive spatial derivatives of the temperature be acceptable, provided the temperature is defined in a FE space of high enough degree?