A more basic approach to solve a problem like this is to check for local maxima for our function in three general areas:
- Interior to the domain.
- Along the smooth boundary.
- At the corners of the boundary.
We then compare values of $f$ at corners and at critical points on boundary/interior. The largest of these will be global maximum.
A numerical approach like simplex method might be applicable to larger variety of problems, but for your example, because your function is differentiable and the boundary of the domain is smooth and easy to parameterize, we can solve it completely by basic calculus techniques.
Along boundary: $y = 1-x$. Substituting,
$$f(x,1-x) = -x^2-(1-x)^2 = -2x^2+2x-1.$$
This function has critical point at $x = 1/2$. Now,
$$f(1/2,1-1/2) = f(1/2,1/2) = -1/2.$$
There is no other boundary to the domain, so now we look at interior of the domain, asking when gradient of $f$ is zero.
$$f_x = -2x = 0,\quad f_y = -2y = 0,$$
which gives us a critical point of $(0,0)$ (you should verify that (0,0) is in the interior of the domain) and $f(0,0) = 0 > -1/2.$
So, the solution is $\max[f] = f(0,0) = 0.$