On two of the three boundaries, the function is just zero. On the last one, you know $x + y = 1$, so you can rewrite your function as a one-variable function
$$f(x, y) = f(x, 1 - x) = x^a e^{-x} (1 - x)^b e^{-(1 - x)} = \frac 1 e x^a (1 - x)^b$$
which isn't too hard to maximize.
Then for the interior, try to identify the critical points where $f_x = f_y = 0$.
We first find any critical point in the interior of the domain, which is $x^2 + y^2 < 1.$ So we have $f'(x, y, z) = (y + z, x - z, x - y)$ and $f'(x, y, z) = 0$ can only happen when $x = y = z$ and $y = -z,$ so $x = y = z = 0,$ but this is not in the interior of the domain. Hence, there are no critical points.
Consider now one part of the frontrier, the one defined by $x^2 + y^2 = z = 1.$ Hence we have to maximise $u(x, y) = f(x, y, 1) = x + xy - y$ with $x^2 + y^2 = 1.$ Set $g(x, y) = x^2 + y^2 - 1.$ Since $g'(x, y) = (2x, 2y),$ the derivative of $g$ is never zero and the Lagrange multipliers method apply. We have to solve $u'(x, y) = \dfrac{\lambda}{2} g'(x, y),$ which signifies $(1 + y, x - 1) = \lambda (x, y).$ Observe that $x = 1$ implies $y = 0$ and so we can consider the point $(1, 0)$ as one of the critical points in the restriction. Putting it aside, we have (by division), $\dfrac{1 + y}{x - 1} = \dfrac{x}{y}$ or else $y + y^2 = x^2 - x = 1 -y^2 - x$ or $x = 1 - 2y^2 - y.$ From the relation $x^2 + y^2 = 1$ we get $1 + 4y^4 + y^2 - 4y^2-2y+y^2 = 1,$ that is, $4y^4-2y^2-2y=0$ which leads to $y = 0$ (and hence $x = 1$) or else $2y^3-y-1=0.$ By inspection, we can factor $(y - 1)(2y^2+2y+1) = 0$ and the only solution is $y = 1$ with $x = -3,$ which is not in the restriction. Hence, we got only one point so far $(1, 0, 1).$
We finally tackle the other part of the restriction, which is $x^2 + y^2 - z = 0.$ Here we solve $f'(x, y, z) = \lambda (2x, 2y, -1)$ and this leads to (upon summing) $y + x = 2\lambda(x + y)$ and $y - x = \lambda.$ If $x = -y,$ then $z = 2x^2 \leq 1$ and we are optimising $h(x) = f(x, -x, 2x^2) = 4x^3 - x^2$ for $|x| \leq \dfrac{1}{2}$ and the optimisers for $h$ are $x = -\dfrac{1}{\sqrt{2}}$ and $x = \sqrt{2} - \dfrac{1}{2},$ this allow deducing the critical points for $f$ to be $\left(-\dfrac{1}{\sqrt{2}}, \dfrac{1}{\sqrt{2}}, 1\right)$ and $\left(\dfrac{1}{\sqrt{2}}, -\dfrac{1}{\sqrt{2}}, 1\right).$ If $x \neq -y,$ then $\lambda = \dfrac{1}{2}$ and $y = x + \dfrac{1}{2}.$ Again, substituting $x^2 + y^2 = z \leq 1$ leads to the point (details ommited) $\left(-\dfrac{1+\sqrt{7}}{4}, \dfrac{1 - \sqrt{7}}{4}, 1\right).$
To terminate the problem, recall the region is compact and evaluate in each of the four points.
Take my calculations with a grain of salt, I did them and checked them twice, but with calculus of several variables, I always commit arithmetic errors.
Best Answer
To find the local extrema of a function subject to inequality constraints, we can't just carelessly split the problem into the boundary and the interior as stated by @Alexandra. To demonstrate why, let's consider one of the points found by @Alexanda, namely $P = \left(\frac{1}{\sqrt{2}}, \frac{1}{\sqrt{2}}\right)$. If this was a local maxima subject to the constraints, then there should be some neighbourhood around this point that's in the feasible region decribed by the constraints where $f$ is always less than it's value at $P$. Notice that at point $P$ the gradient of $f$ is: $$\nabla f_P = - \left[\begin{matrix}e^{- \frac{3}{2}}\\e^{- \frac{3}{2}}\end{matrix}\right]$$ which points towards the interior of the circle (since it points radially inwards). Therefore, we can move in that direction to increase $f$ and so $f$ does not attain a local maximum at $P$.
You can always turn inequality constraints to equality constraints and use regular Lagrange multipliers. In your example, you have the objective function $f(x,y) = ye^{-x^2-2y^2}$ and constraints $g_1(x,y):= x^2 + y^2 -1 \leq 0$ and $g_2(x,y):= -y \leq 0$. If we bump up dimensions from say 2 dimensions to 4 dimensions, we can look at the constraints as $g_1(x,y,a,b):= x^2 + y^2 -1 + a^2 = 0$ and $g_2(x,y,a,b):= -y + b^2 = 0$. Now we can use the regular Lagrange multiplier method to get our solutions in this $4$-D space:
$$\Lambda(x,y,a,b,\lambda_1,\lambda_2) = f - \lambda_1g_1 - \lambda_2g_2$$ $\nabla \Lambda = 0$ gives us the solutions: $(x,y,a,b,\lambda_1,\lambda_2) = (-1, 0, 0, 0, 0, -e^{-1}) \left( 0, \frac{1}{2}, \pm \frac{\sqrt{3}}{2}, \mp \frac{\sqrt{2}}{2}, 0, 0\right), \left(0, 1, 0, \pm1, - \frac{3}{2 e^{2}}, 0\right), \left( 1,0, 0, 0, 0, - \frac{1}{e}\right)$
You can use the bordered hessian to check minima or maxima. Turns out first $(-1,0)$ is a saddle, $\left(0,\frac{1}{2}\right)$ is a max, $(0,1)$ is a min and $(1,0)$ is a saddle.