finite element method - Non-linear Poisson equation over non-rectangular domain

I need to solve non-linear Poisson equation

Laplacian[u[x, y], {x, y}] == u[x, y]^2

Over a non-rectangular domain

The problem in short: non-linear Poisson equation over rectangular domain runs OK, and linear Poisson equation over non-rectangular domain runs OK, but not the non-linear over non-rectangular.

The domain is

boundaries = {-y, .25^2 - (x)^2 - y^2, -x, y - 1, x - 1};

\[CapitalOmega]in = 
ImplicitRegion[And @@ (# <= 0 & /@ boundaries), {x, y}];


Show[RegionPlot[\[CapitalOmega]in], 
ContourPlot[
Evaluate[Thread[boundaries == 0]], {x, 0., 1}, {y, 0, 1.}, 
ContourStyle -> {Purple, Green, Red, Blue, Purple}], 
PlotRange -> {{0.0, 1}, {0., 1.}}, AspectRatio -> Automatic] 

with simple boundary conditions

  Conditions = {DirichletCondition[u[t, x, y] == 1, 
  boundaries[[1]] == 0.],
  DirichletCondition[u[t, x, y] == 1, boundaries[[2]] == 0],

  DirichletCondition[u[t, x, y] == 1, boundaries[[3]] == 0.],
  DirichletCondition[u[t, x, y] == 1, boundaries[[4]] == 0.],
  DirichletCondition[u[t, x, y] == 1, boundaries[[5]] == 0.],
  u[0, x, y] == 1};

I try to run a relaxation scheme

 Eq = Laplacian[u[t, x, y], {x, y}] - u[t, x, y]^2

sol = NDSolveValue[{Eq == Derivative[1, 0, 0][u][t, x, y], 
Conditions}, u, {t, 0, 1}, {x, y} \[Element] \[CapitalOmega]in, 

Method -> {"MethodOfLines", Method -> "Automatic", 
"DifferentiateBoundaryConditions" -> {True, "ScaleFactor" -> 1}}]

The problem is "Nonlinear coefficients are not supported in this version of NDSolve".

The linear Poisson equation runs OK

Eq = Laplacian[u[t, x, y], {x, y}] - u[t, x, y]

sol = NDSolveValue[{Eq == Derivative[1, 0, 0][u][t, x, y], 
 Conditions}, u, {t, 0, 1}, {x, y} \[Element] \[CapitalOmega]in, 
Method -> {"MethodOfLines", Method -> "Automatic", 

"DifferentiateBoundaryConditions" -> {True, "ScaleFactor" -> 1}}]

ContourPlot[sol[1, x, y], {x, y} \[Element] \[CapitalOmega]in, 
 ColorFunction -> "TemperatureMap", Contours -> 50, 
 AspectRatio -> Automatic]

Also, non-linear over a rectangular domain runs OK:

boundaries = {-y, -x, y - 1, x - 1};
\[CapitalOmega]in = 
ImplicitRegion[And @@ (# <= 0 & /@ boundaries), {x, y}];


 Show[RegionPlot[\[CapitalOmega]in], 
 ContourPlot[
 Evaluate[Thread[boundaries == 0]], {x, 0., 1}, {y, 0, 1.}, 
 ContourStyle -> {Purple, Green, Red, Blue, Purple}], 
 PlotRange -> {{0.0, 1}, {0., 1.}}, AspectRatio -> Automatic] 

Conditions = {DirichletCondition[u[t, x, y] == 1, 
 boundaries[[1]] == 0.],
 DirichletCondition[u[t, x, y] == 1, boundaries[[2]] == 0],

 DirichletCondition[u[t, x, y] == 1, boundaries[[3]] == 0.],
 DirichletCondition[u[t, x, y] == 1, boundaries[[4]] == 0.],
 u[0, x, y] == 1};

 sol = NDSolveValue[{Eq == Derivative[1, 0, 0][u][t, x, y], 
  Conditions}, u, {t, 0, 1}, {x, 0, 1}, {y, 0, 1}, 
 Method -> {"MethodOfLines", Method -> "Automatic", 
 "DifferentiateBoundaryConditions" -> {True, "ScaleFactor" -> 1}}]

 ContourPlot[sol[1, x, y], {x, y} \[Element] \[CapitalOmega]in, 

 ColorFunction -> "TemperatureMap", Contours -> 50, 
 AspectRatio -> Automatic]