The first results of a new three-dimensional, finite temperature Skyrme-Hartree-Fock+BCS study of the properties of inhomogeneous nuclear matter at densities and temperatures leading to the transition to uniform nuclear matter are presented. Calculations are carried out in a cubic box representing a unit cell of the locally periodic structure of the matter. A constraint is placed on the two independent components of the quadrupole moment of the neutron density in order to investigate the dependence of the total energy-density of matter on the geometry of the nuclear structure in the unit cell.This approach allows self-consistent modeling of effects such as (i) neutron drip, resulting in a neutron gas external to the nuclear structure, (ii) shell effects of bound and unbound nucleons, (iii) the variety of exotic nuclear shapes that emerge, collectively termed 'nuclear pasta' and (iv) the dissolution of these structures into uniform nuclear matter as density and/or temperature increase.In part I of this work the calculation of the properties of inhomogeneous nuclear matter in the core collapse of massive stars is reported. Emphasis is on exploring the effects of the numerical method on the results obtained; notably, the influence of the finite cell size on the nuclear shapes and energy-density obtained. Results for nuclear matter in beta-equilibrium in cold neutrons stars are subject of part II. The calculation of the band structure of unbound neutrons in neutron star matter, yielding thermal conductivity, specific heat and entrainment parameters, will be outlined in part III.Calculations are performed at baryon number densities of n b = 0.04 -0.12 fm −3 , a proton fraction of y p = 0.3 and temperatures in the range 0 -7.5 MeV. A wide variety of nuclear shapes are shown to emerge. It is suggested that thermodynamical properties change smoothly in the pasta regime up to the transition to uniform matter; at that transition, thermodynamic properties of the matter vary discontinuously, indicating a phase transition of first or second order. The calculations are carried out using the SkM * Skyrme parameterization; a comparison with calculations using Sly4 at n b = 0.08 fm −3 , T = 0 MeV is made.