Solving the 3D Poisson equation on the surface of a cube using a relaxation method

The code below is from another question asked by someone else about solving the 3D Laplace equation for truncated octahedron in a cube matrix, How to solve Laplace equation in 3D?. I wanna know, how would you solve the 3D Poisson equation (which is basically the Laplace equation with a source function), on the surface of a cube, meaning with no boundary conditions, using a relaxation method such as Successive Over Relaxation or Gauss-Seidel or any other relaxation method?

 TruncatedOctahedron = {x + y + z <= 10 && x + y - z <= 10 &&
x - y + z <= 10 && -x + y + z <= 10 && x + y + z >= -10 &&
x + y - z >= -10 && x - y + z >= -10 && -x + y + z >= -10 &&
-6 <= x <= 6 && -6 <= y <= 6 && -6 <= z <= 6}；
Cube = Cuboid[{-100, -100, -100}, {100, 100, 100}];
Subscript[Γ, D] = {DirichletCondition[u[x, y, z] == 200.,
{x, y, z} ∈ TruncatedOctahedron],
DirichletCondition[u[x, y, z] == 15., {x, y, z} ∈ Cube]};
Ω = RegionDifference[Cube, TruncatedOctahedron]
sol = NDSolveValue[{Inactive[Laplacian][u[x, y, z], {x, y, z}] == 0,
Subscript[Γ, D]}, u, {x, y, z} ∈ Ω];

• Please add link to previous question, mentioned in the question. – bbgodfrey Dec 24 '18 at 20:56
• I don't get what your major concern is. Is it really important to you to obtain the solution by a relaxation method? For this simple PDE on a surface, direct factorizations are usually faster and more robust. – Henrik Schumacher Dec 25 '18 at 10:18
• Moreover, without boundary conditions, the Laplacian as a one-dimensional null space (the constant functions) and since it is selfadjoint, the Laplacian is not surjective (the image is the $L^2$-orthogonal complement of the constant functions). So in general, you can only expect least-squares solutions. – Henrik Schumacher Dec 25 '18 at 10:23
• The equation I really wanna solve is: Laplacian[ u[x,y,z], {x,y,z} ] = Exp[u[x,y,z]]. This is more complicated and requires a a relaxation method. – Thando Dec 26 '18 at 22:26
• The exponential on the right-hand side makes it more complicated. Do you have any ideas about how I can solve it on the surface of a cube or a sphere without boundary conditions? Thanks Henrik. – Thando Dec 26 '18 at 22:30

Here is an approach that uses the direct solver Pardiso along with a vanishing-average constraint in order to obtain the least-squares solution.

Since surface finite elements are not built into Mathematica (yet), I use the functions getSurfaceLaplaceBeltrami, getSurfaceMass, getSurfaceLaplacianCombinatorics, SurfaceLaplaceBeltrami, and SurfaceMassMatrix from this post. In order to be self-contained, here is their code:

Quiet[Block[{xx, x, PP, P, UU, U, VV, V, f, Df, u, Du, v, Dv, g,
xx = Table[x[[i]], {i, 1, 2}];
PP = Table[P[[i, j]], {i, 1, 3}, {j, 1, 3}];
UU = Table[U[[i]], {i, 1, 3}];
VV = Table[V[[i]], {i, 1, 3}];
(*local affine parameterization of the surface with respect to the \
"standard triangle"*)
f = x \[Function]
PP[[1]] + x[[1]] (PP[[2]] - PP[[1]]) + x[[2]] (PP[[3]] - PP[[1]]);
Df = x \[Function] Evaluate[D[f[xx], {xx}]];
(*the Riemannian pullback metric with respect to f*)
g = x \[Function] Evaluate[Df[xx]\[Transpose].Df[xx]];
(*two affine functions u and v and their derivatives*)
u = x \[Function]
UU[[1]] + x[[1]] (UU[[2]] - UU[[1]]) + x[[2]] (UU[[3]] - UU[[1]]);
Du = x \[Function] Evaluate[D[u[xx], {xx}]];
v = x \[Function]
VV[[1]] + x[[1]] (VV[[2]] - VV[[1]]) + x[[2]] (VV[[3]] - VV[[1]]);
Dv = x \[Function] Evaluate[D[v[xx], {xx}]];
integrant =
x \[Function]
Evaluate[D[D[v[xx] u[xx] Sqrt[Abs[Det[g[xx]]]], {UU}, {VV}]]];
(*since the integrant is quadratic over each triangle,
we use a three-
point Gauss quadrature rule (for the standard triangle)*)
quadraturepoints = {{0, 1/2}, {1/2, 0}, {1/2, 1/2}};
quadratureweights = {1/6, 1/6, 1/6};
getSurfaceMass =
With[{code =
Part -> CompileGetElement},
Compile[{{P, _Real, 2}}, code, CompilationTarget -> "C",
RuntimeAttributes -> {Listable}, Parallelization -> True]];
integrant =
x \[Function]
Evaluate[
D[D[Dv[xx].Inverse[g[xx]].Du[xx] Sqrt[
Abs[Det[g[xx]]]], {UU}, {VV}]]];
(*since the integrant is constant over each trianle,we use a one-
point Gauss quadrature rule (for the standard triangle)*)
quadraturepoints = {{1/3, 1/3}};
getSurfaceLaplaceBeltrami =
With[{code =
Part -> CompileGetElement},
Compile[{{P, _Real, 2}}, code, CompilationTarget -> "C",
RuntimeAttributes -> {Listable}, Parallelization -> True]]]];

getSurfaceLaplacianCombinatorics =
Quiet[Module[{ff},
With[{code =
Flatten[Table[Table[{ff[[i]], ff[[j]]}, {i, 1, 3}], {j, 1, 3}],
1] /. Part -> CompileGetElement},
Compile[{{ff, _Integer, 1}}, code, CompilationTarget -> "C",
RuntimeAttributes -> {Listable}, Parallelization -> True]]]];

SurfaceLaplaceBeltrami[pts_, flist_, pat_] :=
With[{spopt = SystemOptions["SparseArrayOptions"],
vals = Flatten[
getSurfaceLaplaceBeltrami[Partition[pts[[flist]], 3]]]},
InternalWithLocalSettings[
SetSystemOptions[
"SparseArrayOptions" -> {"TreatRepeatedEntries" -> Total}],
SparseArray[Rule[pat, vals], {Length[pts], Length[pts]}, 0.],
SetSystemOptions[spopt]]];

SurfaceMassMatrix[pts_, flist_, pat_] :=
With[{spopt = SystemOptions["SparseArrayOptions"],
vals = Flatten[getSurfaceMass[Partition[pts[[flist]], 3]]]},
InternalWithLocalSettings[
SetSystemOptions[
"SparseArrayOptions" -> {"TreatRepeatedEntries" -> Total}],
SparseArray[Rule[pat, vals], {Length[pts], Length[pts]}, 0.],
SetSystemOptions[spopt]]];


Creating a triangle mesh on the cube:

R = BoundaryDiscretizeRegion[
Scale[PolyhedronData["Cube", "MeshRegion"], 2],
MaxCellMeasure -> {1 -> 0.05}];


Assembling mass matrix M and stiffness matrix A:

pts = MeshCoordinates[R];
triangles = MeshCells[R, 2, "Multicells" -> True][[1, 1]];
pat = Flatten[getSurfaceLaplacianCombinatorics[triangles], 1];
A = SurfaceLaplaceBeltrami[pts, Flatten@triangles, pat]
M = SurfaceMassMatrix[pts, Flatten@triangles, pat]


Assembling the saddle point matrix L and computing its factorization S (the factorization takes a few seconds):

a = SparseArray[{Total[M]}];
L = ArrayFlatten[{{A, a\[Transpose]}, {a, 0.}}];
S = LinearSolve[L, Method -> "Pardiso"];


Choosing a source function f and sampling it in order to obtain a right-hand side b for the equation L.u == b.

f = X \[Function] Cos[15/3 Indexed[X, 1]] + 1/2 Cos[7/3 Indexed[X, 2]] Cos[13/3 Indexed[X, 3]];
b = Join[M.(f /@ pts), {0.}];


Solving the linear equation:

u = Most@S[b];


Plotting the result:

Graphics3D[{
EdgeForm[],
GraphicsComplex[
pts,
Polygon[triangles],
VertexColors -> ColorData["SunsetColors"] /@ Rescale[u]
]
},
Lighting -> "Neutral"
]


Alternatively, you may use an iterative solver such as "ConjugateGradient":

vals = (# - a[[1]].#/Total[a[[1]]]) &[(f /@ pts)];
c = M.vals;
v = LinearSolve[A, c,
Method -> {"Krylov",
"Preconditioner" -> "ILU0", "Tolerance" -> .0001,
"MaxIterations" -> 300}
];

Max[Abs[A.v - c]]/Max[Abs[c]]
Sqrt[Dot[u - v, M, u - v]]

Max[Abs[A.v - c]]/Max[Abs[c]] (* relative residual *)
Sqrt[Dot[u - v, M, u - v]] (* L^2-distance between the solutions *)


0.000294991

0.000087751

Liouville–Bratu–Gelfand equation

As mentioned in this post, OP wants to solve the following Liouville–Bratu–Gelfand-type equation

$$\Delta u = \exp(u + h ) -1 = \exp(u) \, \exp(h ) -1,$$

with a prescribed function $$h$$. We may write down the system as a function F for which we search a root. Then we employ FindRoot and supply it with the precomputed Jacobian of the systen (implemented as DF) and a not-to-bad initial guess. FindRoot applies Newton's method with line search (the line search seems to be essential here; Newton's method with step size 1 tends to diverge). One can use an arbitrary function h, but OP wants

$$h(x) = \sum_{i=1}^n \log \| x - p_i\|^2$$

with few source points $$p_1,\dotsc,p_n$$ on the surface of the cube. So, let's set it up:

n = 20;
sources = RandomSample[MeshCoordinates[R], n];
Exphvals = Times @@ DistanceMatrix[sources, MeshCoordinates[R]]^2;

ClearAll[F, DF];
F[u_?VectorQ] := A.u + M.(Exp[u] Exphvals - 1.);
DF[u_?VectorQ] := A + M.DiagonalMatrix[SparseArray[Exp[u] Exphvals]];
u0 = ConstantArray[1., Length[Exphvals]];

u = u /. FindRoot[F[u], {u, u0}, Jacobian :> DF[u]];
Max[Abs[F[u]]]


4.58583*10^-14

Graphics3D[{
Sphere[sources, 0.05],
EdgeForm[],
GraphicsComplex[pts, Polygon[triangles],
VertexColors -> ColorData["SunsetColors"] /@ Rescale[u]]
},
Lighting -> "Neutral"
]


• Thanks Henrik, this looks great! However I don't understand the various parts of the code. Would you please explain some parts of the code to me. Is the code solving the equation: Laplacian[ u[x,y,z], {x,y,z} ] = Exp[u[x,y,z]]? Is the size of the cube 3x3x3? Thanks for this code, it's very helpful. I would like to make some changes to it to suit the specific problem I'm working on, which is: 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}]. – Thando Dec 27 '18 at 8:34
• The point {x_i, y_i, z_i} represents the position of the i-th vortex on the surface of the cube , and {i, 1, N} means that the number of vortices on the surface of the cube range from 1 to N. I have to adjust the code to include all of these things, so that I can vary N and see what happens when I do. Thanks. – Thando Dec 27 '18 at 8:39
• Well, the first code block is about assembling the mass and stiffness matrix; I don't like to go too much into detail, but I can say that a standard Ritz-Galerkin scheme with piecewise-linear testfunctions is employed. I used the standard cube as surface R, but you can use any MeshRegion or BoundaryMeshRegion of dimension 2, immersed into $\mathbb{R}^3$ with triangular faces as R. You may replace f also by any other source function. With the nonlinear lower-order term, I would use Newton's method. Give me a few minutes I'll see what I can do about it. – Henrik Schumacher Dec 27 '18 at 9:12
• Just to be sure: You want to solve a PDE on the surface of the cube, right? In that case, the code example in your question is very misleading. – Henrik Schumacher Dec 27 '18 at 9:12
• Sure I want to solve the PDE on the surface of the cube. Sorry about the code example I gave. Thanks. – Thando Dec 27 '18 at 10:41

In version 11.3, no tricks are needed. All scales and amplitudes can be reduced to 1. Code for the Laplace equation with data on a cube and truncated octahedron

Needs["NDSolveFEM"]
L = 3/2; L1 = 1; L2 = 3/5; truncatedOctahedron =
ImplicitRegion[
x + y + z <= L1 && x + y - z <= L1 &&
x - y + z <= L1 && -x + y + z <= L1 && x + y + z >= -L1 &&
x + y - z >= -L1 &&
x - y + z >= -L1 && -x + y + z >= -L1 && -L2 <= x <= L2 && -L2 <=
y <= L2 && -L2 <= z <= L2, {x, y, z}];
cube = ImplicitRegion[-L <= x <= L && -L <= y <= L && -L <= z <=
L, {x, y, z}];

\[CapitalOmega] = RegionDifference[cube, truncatedOctahedron];

sol = NDSolveValue[{Laplacian[u[x, y, z], {x, y, z}] == 0,
DirichletCondition[
u[x, y, z] ==
0, -L2 <= x <= L2 && -L2 <= y <= L2 && -L2 <= z <= L2],
DirichletCondition[u[x, y, z] == 1, True]},
u, {x, y, z} \[Element] \[CapitalOmega]];
{ContourPlot[sol[x, y, 0], {x, -L, L}, {y, -L, L}, Contours -> 20,
ColorFunction -> Hue, FrameLabel -> Automatic,
PlotLegends -> Automatic],
ContourPlot[sol[0, y, z], {y, -L, L}, {z, -L, L}, Contours -> 20,
ColorFunction -> Hue, FrameLabel -> Automatic,
PlotLegends -> Automatic],
ContourPlot[sol[x, 0, z], {x, -L, L}, {z, -L, L}, Contours -> 20,
ColorFunction -> Hue, FrameLabel -> Automatic,
PlotLegends -> Automatic]}


Code for the Poisson equation with data on a cube and truncated octahedron

sol1 = NDSolveValue[{Laplacian[u[x, y, z], {x, y, z}] ==
Sin[2*x]*Sin[2*y] + Cos[3*z],
DirichletCondition[
u[x, y, z] ==
0, -L2 <= x <= L2 && -L2 <= y <= L2 && -L2 <= z <= L2],
DirichletCondition[u[x, y, z] == 1, True]},
u, {x, y, z} \[Element] \[CapitalOmega]];

{ContourPlot[sol1[x, y, 0], {x, -L, L}, {y, -L, L}, Contours -> 20,
ColorFunction -> Hue, FrameLabel -> Automatic,
PlotLegends -> Automatic],
ContourPlot[sol1[0, y, z], {y, -L, L}, {z, -L, L}, Contours -> 20,
ColorFunction -> Hue, FrameLabel -> Automatic,
PlotLegends -> Automatic],
ContourPlot[sol1[x, 0, z], {x, -L, L}, {z, -L, L}, Contours -> 20,
ColorFunction -> Hue, FrameLabel -> Automatic,
PlotLegends -> Automatic]}


• Thanks Alex. Would the code work for the problem: 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}] ? – Thando Dec 27 '18 at 9:08
• You can solve the problem, but you need to formulate the problem first and not at the end. So far everything is very abstract. Do you want to solve a problem in a cuboid, in a sphere or in $\Omega$? – Alex Trounev Dec 27 '18 at 12:40
• @Thando Equation Laplacian[ u[x,y,z], {x,y,z} ] = (Exp[ u[x,y,z] + h[x,y,z] ]) -1` is called Liouville–Bratu–Gelfand equation. – Alex Trounev Dec 27 '18 at 12:56
• Cool, I thought it was a Poisson equation. Yes sure I want to solve it on the surface of a cube or the surface of a sphere, which ever is simpler to code. – Thando Dec 27 '18 at 13:51
• @Thando Then what is truncated octahedron and 3D for? – Alex Trounev Dec 27 '18 at 14:03