Solving the Liouville-Bratu-Gelfand equation on the surface of a cube

Question

How do I solve the equation:

Laplacian[ u[x,y,z], {x,y,z} ] = (Exp[ u[x,y,z] + h[x,y,z] ]) -1, where h[x,y,z]= Sum [ { Log |{x,y,z} - {x_i, y_i, z_i}|^2 }, { i, 1, N}], on the surface of a cube, which is a compact surface without a boundary?

Answer 1

I will give an explanation to the solution of the problem. First, $u[x,y]=-h[x,y]$ is the exact solution of the problem on a plane for h[x_, y_] := Sum[Log[(x - xn[[i]])^2 + (y - yn[[i]])^2], {i, 1, n}]. This raises the question of how the influence of the source given on one face extends to all other faces? If it is through 3D, then this is one problem, and if through 2D, then this is another problem and it has a different solution. Secondly, it is possible to cut and unfold the surface of the cube on a plane. Therefore, the problem is reduced to 2D. But it is necessary to sew the solution on all cuts. Third, we can solve the problem on the surface of a cube in 3D, as did Henrik Schumacher in this post. I wrote code for 2D sources that are a solution to the Laplace equation in 2D. I did not continue the sources through the cuts. If the author needs a continuation, I will add a couple of lines. It shows what the solution looks like on a plane.

c1 = ImplicitRegion[-1 <= x <= 0 && 0 <= y <= 1, {x, y}];
c2 = ImplicitRegion[0 <= x <= 1 && 0 <= y <= 1, {x, y}];
c3 = ImplicitRegion[1 <= x <= 2 && 0 <= y <= 1, {x, y}];
c4 = ImplicitRegion[2 <= x <= 3 && 0 <= y <= 1, {x, y}];
c5 = ImplicitRegion[0 <= x <= 1 && 1 <= y <= 2, {x, y}];
c6 = ImplicitRegion[0 <= x <= 1 && -1 <= y <= 0, {x, y}];
cubS = RegionUnion[c1, c2, c3, c4, c5, c6];

n = 20; n1 = RandomInteger[{0, n}]; n2 = 
 RandomInteger[{0, n - n1}]; n3 = 
 RandomInteger[{0, n - n2 - n1}]; n4 = 
 RandomInteger[{0, n - n3 - n2 - n1}]; n5 = 
 RandomInteger[{0, n - n4 - n3 - n2 - n1}];
n6 = n - n5 - n4 - n3 - n2 - n1;
t0 = 1/50; k = 10; x1 = Table[RandomReal[{-1, 0}], {n1}]; y1 = 
 Table[RandomReal[{0, 1}], {n1}]; x2 = 
 Table[RandomReal[{0, 1}], {n2}]; y2 = Table[RandomReal[{0, 1}], {n2}];
x3 = Table[RandomReal[{1, 2}], {n3}]; y3 = 
 Table[RandomReal[{0, 1}], {n3}];
x4 = Table[RandomReal[{2, 3}], {n4}]; y4 = 
 Table[RandomReal[{0, 1}], {n4}];
x5 = Table[RandomReal[{0, 1}], {n5}]; y5 = 
 Table[RandomReal[{1, 2}], {n5}];
x6 = Table[RandomReal[{0, 1}], {n6}]; y6 = 
 Table[RandomReal[{-1, 0}], {n6}];
xn = Join[x1, x2, x3, x4, x5, x6]; yn = 
 Join[y1, y2, y3, y4, y5, y6]; pn = 
 Table[{xn[[i]], yn[[i]]}, {i, 1, Length[xn]}];
Show[RegionPlot[cubS, AspectRatio -> Automatic], 
 Graphics[{Red, PointSize[.01], Point[pn]}]]

h[x_, y_] := Sum[Log[(x - xn[[i]])^2 + (y - yn[[i]])^2], {i, 1, n}]

U[0][x_, y_] := -h[x, y]
Do[U[i] = 
    NDSolveValue[{(u[x, y] - U[i - 1][x, y])/t0 == 
       Laplacian[u[x, y], {x, y}] - Exp[U[i - 1][x, y] + h[x, y]] + 1,
       DirichletCondition[u[x, y] == U[i - 1][-1, y], 
       x == 3 && 0 <= y <= 1], 
      DirichletCondition[u[x, y] == U[i - 1][0, x], 
       y == 0 && -1 <= x <= 0], 
      DirichletCondition[u[x, y] == U[i - 1][1, -x + 1], 
       y == 0 && 1 <= x <= 2], 
      DirichletCondition[u[x, y] == U[i - 1][0, -x + 1], 
       y == 1 && -1 <= x <= 0], 
      DirichletCondition[u[x, y] == U[i - 1][1, x], 
       y == 1 && 1 <= x <= 2], 
      DirichletCondition[u[x, y] == U[i - 1][x - 2, -1], 
       y == 0 && 2 <= x <= 3], 
      DirichletCondition[u[x, y] == U[i - 1][x - 2, 2], 
       y == 1 && 2 <= x <= 3]}, u, {x, y} \[Element] cubS];, {i, 1, 
   k}];
Table[ContourPlot[U[i][x, y], {x, y} \[Element] cubS, Contours -> 20, 
  ColorFunction -> Hue, PlotLegends -> Automatic, PlotLabel -> i, 
  PlotRange -> All, AspectRatio -> Automatic], {i, 1, k}]

The distribution of sources on the surface and solution evolution starting from i= 0 and up toi = 10. Note that the solution $u(x,y)=-h(x,y)$ at i=0 contains logarithmic singularities.