You are right indeed when saying (or writing) that, for digital signals, the function representing them is undefined in between the samples. When sampling an analogue signal you end up with a discrete version of it.
In the purely mathematical sense of the process, one would really like to use a brick-wall anti-aliasing filter during the analogue-to-digital conversion to get rid of all the frequencies above the Nyquist frequency. In the time-domain, this filter has an Infinite Impulse Response (IIR) ranging from $-\infty$ to $+\infty$, which means that it is a non-realisable filter.
Nevertheless, to follow the theoretical approach we will consider the analytical approach and skip the fact that this is impractical for real-life application. One more thing to keep in mind is that we do not consider the effects that quantisation has on the signal. Quantisation is a non-linear non-invertible process and to simplify things here we will not consider its results (this is most often the case when dealing with signals in an analytical way).
Now, consider a perfect sampling process (only sampling, not quantisation). This is to first band-limit the analogue signal with the perfect brick-wall anti-aliasing filter and then passing it through an ideal sampler, which gathers samples of the signal in infinitesimal time duration every $T_{s} = \frac{1}{f_{s}}$ seconds ($T_{s}$ is the sampling interval and $f_{s}$ the sampling frequency). After this process, you end up with the ideal sampled signal. If you consider the whole process in the frequency domain, all you have done so far is to filter the signal, resulting in a limited frequency range and then sampling it, which is like convolving the signal with an impulse train (with frequency $f_{s}$). This will replicate the frequency content of the band-limited signal at every $f_{s}$ infinitely many times towards $+\infty$ and $-\infty$. If you think about it, infinitely sharp (of infinitely short time) impulses will have an infinite bandwidth (all the possible frequencies, in order to create such a sharp function, for more info look at Fourier Transform). This coincides with the fact that after the sampling process you have frequencies in the range $\left[ -\infty, +\infty \right]$.
At the reconstruction stage of the signal, you should make sure that the extra frequencies introduced by the sampling process are removed since they are not part of the original signal. This can be achieved with a filter, most often termed reconstruction filter. Again, this should ideally be a brick-wall filter removing everything outside the $\left[0, f_{s}\right]$ frequency interval. Once more, in the time-domain, this would act as a smoothing filter to the impulses representing the samples of the digitised signal. In the simplest case where the digital signal has undergone no extra processing, the resulting analogue waveform/signal would have the exact same frequency content as the original band-limited one. If you now go back to the time-domain, according to the Fourier theorem, due to the linearity of the operator (partly also reflected in Parseval's theorem), two functions with the exact same frequency representation (watch out here, the phase must also be taken into account!), will have also the same time-domain representation. Thus, the reconstructed signal is identical to the original one.
As mentioned before, the original signal must be band-limited, but it can be shown that it does not have to be a base-band signal (this is to have frequencies in the $\left[ 0, \frac{f_{s}}{2} \right]$ range).
Conceptually, the reconstruction filter can be thought of as a smoothing filter, as you very well mention in your question. If one looks at the time-domain representation of the ideal brick-wall filters (both anti-aliasing and reconstruction) will realise that they are represented by a $sinc$ function (which is given by $sinc \left(t \right) = \frac{\sin \left(\pi t \right)}{\pi t}$) shown in the figure below.
Filtering in the time-domain is represented by the convolution operation which is pretty much like placing a scaled (with the scaling coefficient being the value of the sample) $sinc$ function at each sample. The zero crossings of the function coincide with the positions where the samples occur (at each sampling interval) resulting in the function not acting on the sample values but the sum of the $sinc$ functions centred on the samples will give the waveform of the reconstructed analogue signal (the $sinc$ functions have non-zero values at every point on the time axis other than where the samples occur, which sum up to create the waveform).
Since we have not considered any specific analogue signal for the analysis, this is valid for any signal. This also applies to a band-limited rectangular function, which by nature exhibits oscillations and overshoot close to the discontinuities of the non-band-limited function (which is a piecewise function).
Now, of course, the brick-wall filters are not realisable (non-causal filters) which means that alternatives have to be found. The best alternative so far (and maybe just because it is the easiest to implement in both hardware and software) is the zero-order hold reconstructor, which results in a staircase with "stair-values" equal to the sample amplitudes. This is then followed by a reconstruction filter. Since the zero-hold reconstructor does not possess a flat frequency response a digital correction filter is usually implemented before the digital signal goes through the reconstructor. For more information on that, you can refer to textbooks such as "Introduction to Signal Processing" by Orfanidis (which is freely distributed) or "Digital Signal Processing" by Proakis & Manolakis.
