Here is a brief sketch of a simple way to do this. For convenience, I will consider an open subset $U$ of the plane. We cover the latter in the natural way with squares of side length $2^{-n}$ and let $U_n$ be the union of all of those inside $U$. Then $U$ is the union of the $U_n$. This reduces to the case of such $U_n$ and this can easily be done explicitly.

Here is a brief sketch of a simple way to do this. For convenience, I will consider an open subset $U$ of the plane. We cover the latter in the natural way with squares of side length $2^{-n}$ and let $U_n$ be the union of all of those inside $U$. Then $U$ is the union of the $U_n$. This reduces to the case of such $U_n$ and this can easily be done explicitly.

Edit: Since I can't comment, here is a reply to the criticism below, not in complete detail since I am working with a pad. Since my $U_n$ is a finite Union of squares (well hypercubes in higher dimensions) it suffices to approximate these externally by smooth sets. In one dimension, we have a compact interval inside an open set. Then we can take a suitable bump function which is $1$ on the interval and zero outside the open set (the convolution of the characteristic function of the interval with a very thin bell function. For a hypercube we simply take tensor products of such bump functions. These functions yield smooth approximations to the hyper cubes in a standard way.