Periodic changes in a thermal soft X-ray flux of a rotating neutron star indicate a non-uniform distribution of the surface temperature. A possible cause of this phenomenon is a suppression of the heat flux across the magnetic field lines in a crust and an envelope of magnetized neutron stars. In this paper we study three-dimensional effects, associated with non-axisymmetric magnetic fields in neutron stars. We calculate the surface temperature distribution by solving numerically a three dimensional heat transfer equation in a magnetized neutron star crust. We adopt an anisotropic (tensorial) electron thermal conductivity coefficient, which is derived as an analytical solution of the Boltzmann equation with a Chapman-Enskog method. To calculate the surface temperature distribution, we construct a local one-dimensional plane-parallel model ("Ts-Tb"-relationship) of a magnetized neutron star envelope. We then use it as an outer boundary condition for the three-dimensional problem in the crust to find the self-consistent solution. To study possible observational manifestations from anisotropic temperature distributions we calculate light curves with a composite black-body model. Our calculations show, that a non-axisymmetric magnetic field distribution can lead to the irregular non-sinusoidal shape of a pulse profile as well as in some cases a significant amplification of pulsations of the thermal flux in comparison to the pure-dipolar magnetic field configurations.