An even more general form of the result can be found in section 7 of this paper: the result applies to Sobolev extension domains (which include Lipschitz domains) and the measure on the boundary is not necessarily the surface measure. This short paper is not self-contained: some of the heavy lifting is done in other sources such as Function spaces and potential theory by Adams and Hedberg (an expensive book, unfortunately). But it does a good job of putting the results together, and outlining two ways of obtaining the desired trace theorem: with and without Besov spaces.
Warning: if you plan to investigate this subject in any depth, you will not avoid (a) Besov spaces and (b) proofs that are not very accessible. They come with the territory.
Terminological remark: I never saw compact trace theorems being named after Rellich and Kondrachov. Their names are normally attached to compact embedding theorems between function spaces on the same domain.
To every second order elliptic PDE $$Lu = f$$on $\Omega \subset \mathbb{R}^n$, where $L$ is an elliptic operator and $f$ a measurable function, is associated a positive-definite, bounded, symmetric coefficient matrix $$Q(x) = [q_{ij}(x)], i,j = 1,\ldots n.$$ That is, the quadratic form $$\mathcal{Q}(x,\xi) = \xi^\top Q(x) \xi$$ satisfies, for some $c,C > 0$ $$c|\xi|^2 \leq \mathcal{Q}(x,\xi) \leq C|\xi|^2$$ for almost every $x \in \Omega$ and every $\xi \in \mathbb{R}^n$. It has been shown by various authors (but most fundamentally De Giorgi, Nash and Moser) that in the presence of local Sobolev and Poincare inequalities, and with the existence of an accumulating sequence of Lipschitz cutoff functions, that weak solutions to $Lu = f$ exist and are Holder continuous. As is well known, weak solutions live in Sobolev spaces, hence their importance to the subject.
Now, suppose we relax the condition on the quadratic form. Call $Lu = f$ degenerate elliptic if we have only that $$0 \leq \mathcal{Q}(x,\xi) \leq C|\xi|^2.$$ In particular, we say that $Q$ degenerates at $x \in \Omega$ if there exists $\xi \neq 0$ such that $\mathcal{Q}(x,\xi) = 0$. As it turns out, allowing that the quadratic form vanish causes major difficulties in adapting the theory of weak solutions. Weak solutions still exist, but the do NOT live, any longer, in Sobolev spaces. Instead, as has been proved by Sawyer and Wheeden (2009) and Rodney (2012), they reside in degenerate Sobolev spaces, which are quite a bit more difficult to deal with.
These spaces are defined with reference to the particular matrix $Q$ with which we are working. So given such a matrix $Q$, define the (possibly infinite) norm $$||w||_{QH^{1,p}(\Omega)} = \left( ||w||_p^p + \int_\Omega |\nabla w^\top Q(x) \nabla w|^\frac{p}{2} dx\right)^\frac{1}{p}$$ on $Lip_{loc}(\Omega)$, the space of locally Lipschitz functions. Note that the gradients of such functions exist almost everywhere by the Rademacher-Stepanov theorem. We define the degenerate Sobolev space $QH^{1,p}(\Omega)$ as the completion of the $$\{w \in Lip_{loc}(\Omega) : ||w||_{QH^{1,2}(\Omega)} < \infty\}$$ analogously to the (secondary) definition of classical Sobolev spaces, but the the gradient norm weighted to the matrix $Q$.
As a remark, when you are reading about this subject in Sawyer and Wheeden's paper, for example, they call the spaces as defined this way $W^{1,p}_Q (\Omega)$, even though this goes against the names given to the classical spaces.
Best Answer
Actually I would say the first thing to remark, is that if a function is in $L^p$ then it is only defined almost everywhere. Therefore, you just cannot in general define its trace since it would mean to get the values of the function on a set of measure $0$ (since of dimension smaller). However, if the function is continuous, you see that you can easily define the trace of your function and it will be continuous.
From this preliminary analysis, you deduce that in general, you need some regularity assumptions to define the trace of a function.
Now look at a function with a local singularity such as $$ f(x) = \frac{1}{|x|^a} $$ You can see that this function is locally in $L^p(\mathbb{R}^d)$ if $p<d/a$, but if you take the trace on a set of smaller dimension and containing $0$, you see that the trace will only locally be in $L^q$ with $q<d/a - 1/a$, so you loose a part of the integrability when you take the trace. This is from my point of view a way to understand intuitively why starting from a function with a certain regularity, you loose a part of the regularity when taking the trace.
The fractional Sobolev spaces created by real interpolation were investigated a lot by Jacques-Louis Lions, and actually were sometimes called trace spaces. A good reference is the book by Luc Tartar, An Introduction to Sobolev Spaces and Interpolation Spaces. Chapter 16 treat the case of the $L^2$ based $H^s$ Sobolev spaces and Chapter 40 of the more general case of $L^p$ based Sobolev spaces $W^{s,p}$.
An interesting part is also Chapter 33 about the space $H^\frac{1}{2}_{00}$, which in some sense the critical case where one can still define a trace on the border (since $H^s_0(\Omega) = H^s(\Omega)$ when $s\leq 1/2$).