Abstract The singular integral equation, obtained from applying Green's identity to the stationary part of the linear diffusion equation operator (the Laplace operator), is implemented in an element-by-element fashion for the solution of nonlinear heat conduction in heterogeneous media so that both the temperature and its normal directional distribution (the flux) are calculated directly at every nodal point. This formulation, referred to as the flux Green element method (flux GEM), is used to solve linear and nonlinear transient heat transfer problems in two-dimensional heterogeneous domains. The flux-based Green element formulation takes advantage of the ability of the boundary element theory to correctly calculate the fluxes so that the solution in each element is complete in the sense that any solution (temperature or its directional derivative) required at any other point other than the grid points is achieved by implementing the integrations only in the element in which the point is located. One apparent drawback of the formulation is increased number of degrees of freedom (unknowns) that are generated, but this is compensated by the fact that high accuracy is achieved with coarse discretization. Furthermore, this paper presents a novel way of addressing the closure problem associated with having more unknowns than discretized equations at internal nodes. Here, an additional compatibility equation is presented for the internal normal directional derivatives of the temperature. This equation has a universal appeal in the sense that it applies not only to this heat equation but to any situation where such directional derivatives exist at a point. For the nonlinear transient heat conduction, the Picard and Newton–Raphson algorithms are applied in linearizing the discrete equations. In all examples examined in this paper, the Newton–Raphson algorithm shows superior convergence characteristics. The flux GEM is applied to four examples, and in all cases, the flux GEM achieves comparable accuracy with between 20% and 33% of the grid of the earlier Green element formulation which differences the internal fluxes.