The assessment of the temperature distribution on windows and façade components is crucial for safety and durability of building skins, as thermal stress from solar heating often causes breakage. Cast shadows further enhance temperature irregularities. Approximated methods of analysis, often proposed in standards, are used in the practice to determine the temperature distribution, but these do not always allow for precise calculation of the resulting thermal stress. On the other hand, advanced thermal analysis software is not always accessible to structural engineers. A “flux-based” approach was recently proposed to specifically evaluate the time-dependent 3D temperature distribution in monolithic and layered glazing. This theory expands upon Biot’s variational method from the 1950s, which is grounded in the definition of the heat displacement field, whose time derivative represents heat flux. The solar heat absorbed within the material, which can create a complex spatial distribution due to the refractive properties of the materials, is accounted for as an additional component of the heat displacement. The flux-based approach has been implemented in an in-house finite element code, taking advantage of the variational form of the equations. For layered glazing, this framework is in general more efficient than the traditional temperature-based approaches since the variable is the heat displacement, which enjoys a higher regularity than the temperature field. Therefore, it is particularly suitable for problems entailing non-smooth temperatures fields, for example at the interfaces between the different layers of the laminate and between shadowed and exposed regions of the glass surface. As a paradigmatic example, the developed code is used to determine the effect of size and shape of shadows in monolithic and laminated glass panes, obtaining results of general value that can be used to support the use of simplified models for the thermal design of glass façades.