Consider the set $U = (-1,1) \times \{ 0\} \subset \mathbb R^2$ and a continuous function $f : U \rightarrow \mathbb R$.
Then for any domain $\Omega \subset \mathbb R^2$ such that $U \subset \partial\Omega$, there exist many functions $u$ harmonic in $\Omega$ such that $u|_U = f$ (as the Dirichlet problem is not completely determined).
My question is what happens if we want to prescribe the value of $u$ on $U$ and also the value of $\partial_\nu u$ on $U$? Does there exist a (maximal) domain $\Omega$ and function $u$ harmonic in $\Omega$ with the desired behavior in $U$?
I think, in general, there cannot exist a function $u$ and domain $\Omega$ satisfying this. For example, we could ask for $\partial_\nu u|_U=0$ and $u|_U(x)=\max(0,\vert x \vert -\frac 1 2)$. Thus, by unique continuation at the boundary for harmonic functions, since $u$ and its gradient are both $0$ in an open set of the boundary then $u$ must be identically zero. But the prescribed value at $u$ is not identically zero, so there cannot exist any $u$.
I was wondering whether we can find some compatibility conditions on the prescribed values of $u$ and its normal derivative such that we can ensure the existence of some solution. Also, if there exists a solution, what can we say about the maximal domain where it is defined? Does anybody know about any reference on this kind of problem?
In short, if you're interested in this problem, I would encourage you to do some research on the no-sign obstacle problem. I'll expand a bit below.
You're correct that the problem of harmonic continuation leads to a type of free boundary problem. Let $\Omega \subset \mathbb{R}^d$ be bounded and open. Let $f$ be a bounded measurable function on $\Omega$. Consider the Newtonian potential $$w(x) = (N*(f\chi_\Omega))(x) = \int_\Omega N(x-y)f(y) \ dy,$$ where $N$ is the fundamental solution to Laplace's equation in $\mathbb{R}^d$. Then, $w \in W_{\mathrm{loc}}^{2,p}(\mathbb{R}^d)$ for any $1<p<\infty$. Further, $w$ satisfies $$\Delta w = f\chi_\Omega,$$ in the sense of distributions and, in particular, $w$ is harmonic in $\mathbb{R}^d - \overline{\Omega}$.
Now, let $P \in \partial\Omega$ and suppose that, for some small $r>0$, a harmonic continuation of $w$ into $\Omega$ at $P$ exists, call it $v$. Then, $v$ is harmonic in the ball $B_r(P)$ and $v=w$ on $B_R(P)-\Omega$. Then, the difference $u = w - v$ satisfies $$\begin{cases} \Delta u = f\chi_\Omega &\text{ in } B_r(P)\\ u = |\nabla u| = 0 &\text{ in } B_r(P) - \Omega \end{cases}. \tag{*}$$
By Cauchy-Kovalevskaya, the harmonic continuation exists if both $\partial\Omega$ and $f$ are (real) analytic in some neighborhood of $P$. There are, however, a number of interesting questions one might consider regarding such problems.
The above problem $(*)$ is basically the no-sign obstacle problem and so the second question is really just a question about the regularity of the free boundary. I'd have to do some reading to recall the free boundary regularity results for the no-sign obstacle problem. I believe that the free boundary will be $C^{1,\alpha}$ except possible for a singular set $\Sigma$ where $\mathcal{H}^{d-1}(\Sigma)=0$. There are almost certainly sharper results out there. Further, in the case $d=2$, the singular set may be empty. I'll update if I come across any of these results or anyone is actually interested in talking about them.
Edit: I don't know of any simple solutions to the no-sign obstacle problem. Free boundary problems tend to be very difficult to solve. We can sometimes analytically solve very simple examples, but doing so, even in this case, tends to be very nontrivial. I am sure that you can find some examples on the internet.
I do know of some examples to a slightly different, but very related problem. We want to consider a general version of the mean-value property for harmonic functions. Recall that, if $h$ is harmonic in $B_1$, then $$\frac{1}{\omega_d}\int_{B_1} h(y) \ dy = h(0),$$ where $\omega_d$ is the volume of the unit ball in $\mathbb{R}^d$.
Consider a distribution $\mu$ with compact support. Let $\Omega$ be a bounded domain in $\mathbb{R}^d$ containing $\mathrm{supp}(\mu)$ such that, for any integrable harmonic function $h$ on $\Omega$, we have $$\int_\Omega h(y) \ dy = \int h \ d\mu.$$ Then, $\Omega$ is what is called a quadrature domain.
Let $x \in \complement\Omega$ and let $h(y) = |x-y|^{2-d}$ if $d \geq 3$ and $h(y) = -\log|x-y|$ if $d=2$. Then, we want $$\int_\Omega h(y) \ dy = \int h(y) \ d\mu(y).$$ Set $$u(x) := c_d\int_\Omega h(y) \ dy - c_d\int h(y) \ d\mu(y),$$ for a suitable constant $c_d$. Then, on $\mathbb{R}^d$, $u$ satisfies $$\Delta u = \chi_\Omega - \mu, \ u=|\nabla u| = 0 \text{ in } \complement\Omega.$$ So, for a point $P \in \partial\Omega$ and a small $r>0$, we will have $$\Delta u = \chi_\Omega \text{ in } B_r(P), \ u = |\nabla u| = 0 \text{ in } B_r(P) - \Omega.$$
There are an abundance of examples of quadrature domains. One example is the cardioid. We can represent the cardioid as the image of the unit circle under an appropriate conformal mapping ($z = w + \frac{1}{2}w^2$). In this case, the measure is $\mu = \frac{\pi}{2}(3\delta + \partial_x\delta)$.