Change in neumann boundary conditions through coordinate transformation of elliptic PDE, weak formulation

Question

The standard weak formulation of the Neumann problem for the Poisson equation is to find $u \in H^1 ( \Omega)$ such that for every $v \in H^1 ( \Omega)$:

$$ \int_{\Omega} \nabla u \nabla v d x = \int_{\Omega} fv d x + \int_{\partial \Omega} gv d s $$

for given $f \in L^2 ( \Omega)$ and $g \in H^{- 1 / 2} ( \Omega)$.

My goal is to introduce an arbitrary coordinate change and arrive at an analogous weak formulation.

I start with the classical formulation $- \Delta u ( x) = f ( x)$ and $\partial_n u ( x) = g (x)$ with everything fine and smooth. Let now $\phi : \Omega \rightarrow \tilde{\Omega} : x \mapsto \phi ( x) = : \tilde{x}$ be $C^{\infty}$ and bijective and define $\tilde{u} \circ \phi = u$. Applying the chain rule to $ \Delta_x ( \tilde{u} \circ \phi)$ we find the differential operator $L$:

$$ L \tilde{u} ( \tilde{x}) = a_{i j} ( x) \partial_{i j} \tilde{u} ( \tilde{x}) + b_i ( x) \tilde{u} ( \tilde{x}), $$

where $a_{i j} ( x) = \partial_k \phi_i ( x) \partial_k \phi_j ( x)$, $b_i ( x) = \Delta \phi_i ( x)$ and we used the usual summation convention over repeated indices. Yes the coordinates are all mixed up: I left out a composition with $\phi^{- 1}$ because it's messy.

It now holds that $\Delta_x u ( x) = L \tilde{u} (\phi ( x))$ for all $x \in \Omega$.

Question 1: My operator is not in divergence form, so deriving the weak formulation is going to be a real mess. What would be a better approach?

Question 2: How do I transform the boundary conditions? I'd like something like

$$ \partial_{\nu} \tilde{u} ( \tilde{x}) = \tilde{g}( \tilde{x}) . $$

The chain rule applied to $\partial_n ( \tilde{u} \circ \phi)$ results in $\partial_{\nu} \tilde{u} ( x) = \nabla \tilde{u} ( x) \nu ( x) = \nabla \tilde{u} ( x) D \phi ( x) n ( x)$, which is partly natural since the normal vector field is transformed with the differential, but is far from what I'd need during the partial integration deriving the weak formulation, which would rather be (assuming $L$ were in divergence form)

$$ \partial_{\nu} \tilde{u} = a_{i j} \partial_i \tilde{u} \nu_j .$$

This must be pretty basic, but I'm a bit confused here...

Answer 1

Actually the operator is in divergence form. You do not see it in those coordinates because you are missing the composition with $\phi^{-1}$ and because the computation is indeed a mess.

I think it is best to write the variational formulation and transform the functional. Then the divergence form will be apparent. Let me take $g=0$ to simplify the computation (you can work out the general case). The function $u$ minimizes the functional $$J(u) = \int \frac{|\nabla u|^2} 2 + f u \ dx $$ from all functions with average zero.

Let $\varphi = \phi^{-1}$, and let us write the functional in terms of $\tilde u$. $$ \tilde J(\tilde u) = \int \left( \frac{|\nabla \tilde u(y) D\varphi(y)^{-1}|^2} 2 + f(\varphi(y)) \tilde u \right) \det D \varphi \ dy $$

So, $\tilde u$ satisfies a divergence from equation like $\partial_i a_{ij} \partial_j \tilde u = \tilde f$ where $$ \{a_{ij}\} = \det D\varphi \ \left(D\varphi D\varphi^T \right)^{-1} \qquad \text{ and } \qquad \tilde f(y) = f(\varphi(y)) \det(D\varphi(y)).$$

with the natural Neumann condition $\nu_i a_{ij} \partial_j = 0$ on $\partial \phi(\Omega)$.

The case of $g \neq 0$ is messier because you have to compute its Jacobian restricted to the surface $\partial \Omega$.