Nonrectangular region for NDSolve

Question

I have a PDE with mixed boundaries (Neumann and Dirichlet on some sides) in the region

$(t,x,y) \in \left( 0, T\right) \times\left\{ -L \leq x \leq L, 0 \leq y \leq h(x) \right\}$

where $h(x)$ is something like $\exp\left\{ -(x-x^*)^2\right\}$, doesn't matter. And I have Neumann boundary condition on the curve $\left(x, h(x)\right)$.

How can I provide such region to NDSolve?

Answer 1

This problem can be easily solved using V10's new FEM functionality. For concreteness, let's suppose we want to solve the heat equation $$u_t - \Delta u = 0$$ over the region $$\left\{(x,y): -1 \leq x \leq 1, \; 0 \leq y \leq e^{-x^2}\right\}.$$ We'll take the initial temperature distribution to be identically 1, i.e. $u(x,y,0)=1$; we'll suppose the bottom edge is held at 1 while the left and right edges are held at 0, i.e. $u(x,0,t)=1$ and $u(-1,y,t)=u(1,y,t)=0$; and we'll suppose the curved top is insulated, i.e. the normal derivate of $u$ is zero along the curve $(x,e^{-x^2})$.

Needs["NDSolve`FEM`"];
Clear[u];
omega = ImplicitRegion[-1 <= x <= 1 && 0 <= y <= Exp[-x^2], {x, y}];
mesh = ToElementMesh[omega];
gamma1 = DirichletCondition[u[t, x, y] == 0, x == 1 || x == -1];
gamma2 = DirichletCondition[u[t, x, y] == 1, y == 0];
u = NDSolveValue[{D[u[t, x, y], t] - Laplacian[u[t, x, y], {x, y}] == 
  NeumannValue[0, y == Exp[-x^2]], gamma1, gamma2,
  u[0, x, y] == 1}, u, Element[{x, y}, mesh], {t, 0, 3}];

Note that the NeumannValue is specified as part of the differential equation itself.

We can now plot the solution:

pics = Table[Plot3D[u[t, x, y], Element[{x, y}, omega],
  BoundaryStyle -> Thick, ViewPoint -> {2.4, 2.25, 0.9},
  ColorFunction -> "TemperatureMap", ColorFunctionScaling -> False, 
  PlotRange -> {0, 1.01}], {t, 0, 3, 0.05}];
ListAnimate[pics]

Answer 2

NDSolve requires a rectangular domain, so you have to make a change of coordinates.

$(t,x,y) \in \left( 0, T\right) \times\left\{ a \leq x \leq b, g(x) \leq y \leq h(x) \right\}$

Since you have explicit expressions for your boundaries ($x=a$, $x=b$, $y=g(x)$, and $y=h(x)$), we can use a linear interpolation

coords={u,v};
x[u_, v_] = (1 - u)*a + u*b
y[u_, v_] = With[{x=x[u,v]},(1 - v)*g[x] + v*h[x]]
$Assumptions = {Element[{g[_], h[_]}, Reals], h[a_] > g[a_]}

In terms of these coordinates your domain becomes

$(t,u,v) \in \left( 0, T\right) \times\left\{ 0 \leq u \leq 1, 0 \leq v \leq 1 \right\}$

In my other answer I described how to transform the PDE. The general steps are

1. write your metric
2. find its inverse and determinant
3. find the basis vectors and basis covectors
4. transform components of vectors into the new coordinate bases
5. find expressions for operations (directional derivatives, laplacians, etc.)
6. transform PDE and IV/BVs

NDSolve should then be able to give you a solution in terms of $u$ and $v$.

sol=F/.First@NDSolve[eqns,F,{t,0,T},{u,0,1},{v,0,1}]

You then have to map your result onto your original domain.

Clear[x,y];
u[x_,y_]:= (x-a)/(b-a)
v[x_,y_]:= (y-g[x])/(h[x]-g[x])
F[t_,x_,y_]:= sol[t,u[x,y],v[x,y]]