A new and effective computational approach is presented for analyzing transient heat conduction problems. The approach consists of a meshless Fragile Points Method (FPM) being utilized for spatial discretization, and a Local Variational Iteration (LVI) scheme for time discretization. Anisotropy and nonhomogeneity do not give rise to any difficulties in the present implementation. The meshless FPM is based on a Galerkin weak-form formulation and thus leads to symmetric matrices. Local, very simple, polynomial and discontinuous trial and test functions are employed. In the meshless FPM, Interior Penalty Numerical Fluxes are introduced to ensure the consistency of the method. The LVIM in the time domain is a combination of the Variational Iteration Method (VIM) and a collocation method in each finitely large time interval. The present methodology represents a considerable improvement to the current state of science in computational transient heat conduction in anisotropic nonhomogeneous media. In this first part of the two-paper series, we focus on the theoretical formulation and implementation of the proposed methodology. Numerical results and validation are then presented in Part II.