A computational framework is presented for the continuum modeling of cellular biomolecular diffusion influenced by electrostatic driving forces. This framework is developed from a combination of numerical methods, geometric meshing and computer visualization tools. In particular, a hybrid of (adaptive) finite element and boundary element methods is adopted to solve the Smoluchowski equation (SE), the Poisson equation (PE), and the Poisson-Nernst-Planck equation (PNPE) in order to describe electrodiffusion processes. The finite element method is used because of its flexibility in modeling irregular geometries and complex boundary conditions. The boundary element method is used due to the convenience of treating the singularities in the source charge distribution and its accurate solution to electrostatic problems on molecular boundaries. Nonsteady-state diffusion can be studied using this framework, with the electric field computed using the densities of charged small molecules and mobile ions in the solvent. A solution for mesh generation for biomolecular systems is supplied, which is an essential component for the finite element and boundary element computations. The partially coupled Smoluchowski equation and Poisson-Boltzmann equation (PBE) are considered as special cases of the PNPE in numerical algorithm, and therefore can be solved in this framework as well. The possible extensions of the physical model in this frame are also discussed. Some example computations are reported for: reaction-diffusion rate coefficient, ion density distribution, time-dependent diffusion process of the neurotransmitter consumption.