Abstract Current work attempts to study the heatline patterns during natural convection for different types of Dirichlet heatfunction boundary conditions. The enclosures with various shapes (square, curved, trapezoidal, tilted square and parallelogrammic) are considered with various thermal boundary conditions such as (a) case 1: hot left wall, cold right wall and adiabatic horizontal walls, (b) case 2: hot bottom wall, cold left and right walls and adiabatic top wall and (c) case 3: hot bottom wall with other cold walls. Traditionally, the reference of heatfunction (Π=0) is assumed at the adiabatic wall and the implementation of reference (Π=0) may be non-trivial for the case with zero or multiple adiabatic wall(s). Various heatfunction boundary conditions have been formulated based on locations of Π=0 for systems with more than one adiabatic walls (case 1) or no adiabatic wall (case 3). As test problems, Π=0 is considered at the junctions of isothermal walls (cases 2 and 3) or on the isothermal wall (case 3). The governing equations are solved via the Galerkin finite element method at various Rayleigh numbers (103 and 105) and Prandtl numbers (Pr=0.015 and 7.2). The magnitudes of the heatfunctions change drastically with the location of the datum of Π (Π=0) whereas, the heat flow patterns remain same irrespective of the heatfunction boundary conditions. The gradients of heatfunctions or the heat flux along the active walls (hot/cold) are invariant of the choice of the reference (Π=0). The local and average Nusselt numbers are also independent of the choice of Π=0 and the Nusselt numbers are found to be identical with heatfunction gradients obtained with various locations of Π=0. Current work may be useful for heat flow visualization in various thermal systems involving complex thermal boundary conditions. © 2015 Elsevier Ltd.