A model fractal forcing for direct numerical simulations (DNS) of stationary, homogeneous and isotropic turbulence is proposed. The power spectrum of this forcing is a power-law function of wavenumber with a positive exponent that is an increasing function of $D_{f}$, the fractal dimension of the fractal stirrer. The following results are obtained for DNS turbulence subjected to fully self-similar fractal forcing of fractal dimension $D_f$. (i) The Taylor and Kolmogorov microscales are proportional to each other and to the smallest length scale of the fractal forcing. (ii) The integral length scale is much smaller than the size $L_b$ of the DNS box and a decreasing function of the extent of the fractal forcing range because fractal forcing generates very irregular velocity fields. (iii) In qualitative (but not quantitative) agreement with renormalization group (RG) theories of turbulence, higher values of $D_f$ lead to increased energy at the highest wavenumbers. (iv) The energy and energy input rate spectra and the inter-scale energy transfer $T(k)$ all scale with the turbulence r.m.s. velocity $u'$ and the Taylor microscale $\lambda$. (v) More than 80% of the total dissipation occurs in the fractal forcing range of scales which extends from $L_b$ to about one to two times the Taylor microscale. In that range, $T(k)$ is negligible. (vi) $\lambda \,{\sim}\, \nu/u'$ (where $\nu$ is the kinematic viscosity), the kinetic energy dissipation rate per unit mass $\epsilon \,{\sim}\, u'^{3}/\lambda \,{\sim}\,u'^{4}/\nu$ and the velocity derivative skewness $S$ is independent of Reynolds number in the limit where the Reynolds number and the fractal forcing range are increased together. (vii) The intermediate eigenvalue of the strain rate tensor is on average positive, and the negligible values of $T(k)$ in the fractal forcing range are accompanied by lower values of $S$ and a significant reduction in local compression by comparison to turbulence forced only at the large scales. (viii) The geometrical alignments between vorticity, strain rate tensor eigenvectors and vortex stretching vector are qualitatively as in turbulence forced only at the large scales but significantly weakened. (ix) The p.d.f.s of velocity increments are approximately, though not exactly, Gaussian at all scales between the Kolmogorov and integral length scales. A few preliminary DNS results are also given for the case of turbulence generated by a fractal forcing that is discretely, as opposed to fully, self-similar.