Dacorogna, Bernard (2004). Introduction to the Calculus of Variations. London: Imperial College Press.
Or
Gilbarg, David; Neil S. Trudinger (1988). Elliptic Partial Differential Equations of Second Order. Springer
It seems to me that a definition of "regularity" really depends on where it is being applied, for instance the regularity theory of weak solutions refers to the extra differentiability a function may have, due to it solving a PDE. The regularity a solution can inherit depends on the properties of the problem, i.e., the smoothness of a domain boundary, the smoothness of boundary conditions and initial data. Notice here that the smoothness of the initial data falls under the category of "regularity of initial data".
Let's consider the example you gave, let $u$ be a harmonic function, which I assume you mean in the classical sense, i.e., $u\in C^2(\Omega)\cap C(\overline{\Omega})$, and $$-\Delta u=0\quad\mbox{in}\quad\Omega,$$
even in this classical setting we can gain "regularity" by plugging a derivative of $u$ into the equation and noticing that $-\Delta$ commutes with the derivative. We then find that $u\in C^\infty(\Omega)-$it is in fact more regular than a general function in $u\in C^2(\Omega)\cap C(\overline{\Omega})$. If we wanted to extend this smoothness to the boundary, we would need $\partial\Omega$ to satisfy some smoothness conditions.
You mentioned that if we multiply the equation by a test function $\varphi\in C^\infty_c(\Omega)$ and integrate by parts, we obtain
$$\int_\Omega u\Delta\varphi=0\quad\forall\varphi\in C^\infty_c(\Omega),$$
which indeed tells us that $u\in L^1_{loc}(\Omega)$. But we already knew this, since $u\in C^2(\Omega)\cap C(\overline{\Omega})$.
Regularity theory works the other way: (to make my job easier) take $u\in H^1(\Omega)$, and assume
$$\int_\Omega u\Delta\varphi=0\quad\forall\varphi\in C^\infty_c(\Omega),$$
this now tells us more about the nature of $u$. We can integrate by parts once to find
$$\int_\Omega \nabla u\cdot\nabla\varphi=0\quad\forall\varphi\in C^\infty_c(\Omega),$$
now by the choice of a suitable test function (using difference quotients), we can actually deduce that $u\in H^2_{loc}(\Omega)$, in which case we can justify integrating by parts once more to find
$$\int_\Omega \varphi\Delta u=0\quad\forall\varphi\in C^\infty_c(\Omega).$$
This holds a.e. and so $-\Delta u=0$ a.e.
We can bootstrap the difference quotients - test function method to deduce that $u\in H^k_{loc}(\Omega)$ for all $k\in\Bbb N$, i.e., $u\in C^\infty(\Omega)$. Again if the boundary is more regular, we can extend this notion of improved differentiability to the boundary.
Overall we have used properties of the Laplace equation to deduce that weak solutions behave in a nicer way than they "need to" to be considered a solution, i.e., we go from $H^1(\Omega)$ to $C^\infty(\overline{\Omega})$.
We see that the solution is more "regular", since really the word "regular" could be coined with the notion of "well behaved", and indeed in general a $C^\infty(\overline{\Omega})$ function is more well behaved than a $H^1(\Omega)$ function.
Best Answer
As suggested in the comments, you want to integrate by parts. From what I can tell, your worry is that $\Omega$ is only assumed to be open and bounded, and integration by parts requires some kind of regularity assumption on the boundary. You're right, but this is not an issue here because $\varphi $ has compact support in $\Omega_\varepsilon$. Indeed, for each $\varphi \in C^\infty_0(\Omega_\varepsilon)$, choose a set $U \subset \Omega$ such that $\operatorname{supp} \varphi \subset \subset U \subset \subset \Omega_\varepsilon$ and $\partial U$ is smooth. Then, by Green's identity\begin{align*} \int_{\Omega_\varepsilon} \varphi \Delta u_\varepsilon \, dx &=\int_U \varphi \Delta u_\varepsilon \, dx \\ &= \int_U u_\varepsilon \Delta \varphi \, dx + \int_{\partial U}\bigg(u_\varepsilon \frac{\partial \varphi }{\partial \nu} - \varphi \frac{\partial u_\varepsilon}{\partial \nu} \bigg)\, d\mathcal H^{n-1}_x\\ &=0. \end{align*} Note that all the derivative on $u_\varepsilon$ make sense because $u_\varepsilon$ is smooth. Since $\varphi$ was arbitrary, it follows that $\Delta u_\varepsilon =0$ in $\Omega_\varepsilon$ as required.