  Initial data corresponding to spacetimes containing black holes are
considered in the time symmetric case. The solutions are obtained by matching
across the apparent horizon different, conformally flat, spatial metrics. The
exterior metric is the vacuum solution obtained by the well known conformal
imaging method. The interior metric for every black hole is regular everywhere
and corresponds to a positive energy density. The resulting matched solutions
cover then the whole initial (Cauchy) hypersurface, without any singularity,
and can be useful for numerical applications. The simpler cases of one black
hole (Schwarzschild data) or two identical black holes (Misner data) are
explicitly solved. A procedure for extending this construction to the multiple
black hole case is also given, and it is shown to work for all time symmetric
vacuum solutions obtained by the conformal imaging method. The numerical
evolution of one such 'stuffed' black hole is compared with that of a pure
vacuum or 'plain' black hole in the spherically symmetric case.
