The first is really easy--- you make a lattice approximation, put a lattic-spacing dependent coefficient in front of $\delta^2(x)$ to make it $c(\epsilon) (\delta_\epsilon(x))^2$ and take the continuum limit $\epsilon\rightarrow 0$, with $c\propto\epsilon$, and from this you learn that $c \delta(x)^2$ becomes $ c'\delta(x)$ in the limit, where c' is the renormalized coefficient. In other words, squaring the delta function gives a delta function with an infinite coefficient.
This is regularization dependent. If you point-split, so that the product is for delta-functions at separate points, you can get 0 as the answer. You have to define the lattice thing, and show how to take the limit. That's what people do in quantum field theory anyway.
For more complicated cases, like the product $f(x) \delta''(x)\delta'(x)$, no point splitting, you can get a $\delta$ function (if the two are at the same point), or you can get a linear combination of $\delta$ and $\delta'$ if the two approximate delta-functions are nearby but not right on top of each other.
The main issue is that distributions are in the continuum limit from the start, and field theory is defined so that you take the continuum limit at the end, because the continuum limit is complicated.
This is why distribution language is disfavored for modern field theory formulations, in favor of an explicit regulator and explicit renormalization scheme. This was the Feynman-Pauli-Villars regulator approach, which is universally used today.
Your second example is more complicated, because to define the delta-metric, you have to imagine a series of manifolds that approximate the delta-metric, and there is no guarantee that the limiting conception is sensible.
Consider a 2d metric which is
$$ g_{\mu\nu} = (1+ \delta_\epsilon(x)\delta_\epsilon(y)) (-l_{00} + l_{11})$$
This is conformally flat (which is not saying much in 2d). The result will have opposite sign stresses right near the origin, which come from differentiating the approximate delta functions, and it is not going to be physical. The requirements that the energy is positive, or that the vacuum equations are satisfied, each impose an ellipticity constraint that eliminates such localized curvatures.
It is not possible to smoosh too much curvature in a small region, because the hoop bound tells you that mass M requires an extension of about M to be visible without a horizon.
So one can ask: what kind of localized singularities are allowed in vacuum or energy-condition restricted GR? The only case I know of is the 2d cone-metric, for a 2d space, where there is a delta-function in the curvature at the cone point. This generalizes to a codimension 2 cone point in higher dimensions, and to orbifold singularities of various dimensions in string theory, where you can T-dualize the orbifold to lower dimensions.
Aside from orbifolds (where the metric is not singular, only the curvature), I don't believe there are other examples where Einstein's equations allow a singular metric (the centers of black holes are not part of the manifold).