0
$\begingroup$

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} ∈ Ω];
$\endgroup$
  • 1
    $\begingroup$ Please add link to previous question, mentioned in the question. $\endgroup$ – bbgodfrey Dec 24 '18 at 20:56
  • $\begingroup$ 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. $\endgroup$ – Henrik Schumacher Dec 25 '18 at 10:18
  • $\begingroup$ 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. $\endgroup$ – Henrik Schumacher Dec 25 '18 at 10:23
  • $\begingroup$ 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. $\endgroup$ – Thando Dec 26 '18 at 22:26
  • $\begingroup$ 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. $\endgroup$ – Thando Dec 26 '18 at 22:30
8
$\begingroup$

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, 
    integrant, quadraturepoints, quadratureweights}, 
   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 = 
       N[quadratureweights.Map[integrant, quadraturepoints]] /. 
        Part -> Compile`GetElement}, 
     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}};
   quadratureweights = {1/2};
   getSurfaceLaplaceBeltrami = 
    With[{code = 
       N[quadratureweights.Map[integrant, quadraturepoints]] /. 
        Part -> Compile`GetElement}, 
     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 -> Compile`GetElement}, 
     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]]]}, 
   Internal`WithLocalSettings[
    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]]]}, 
   Internal`WithLocalSettings[
    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"
 ]

enter image description here

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",
     Method -> "ConjugateGradient",
     "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"
 ]

enter image description here

$\endgroup$
  • $\begingroup$ 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}]. $\endgroup$ – Thando Dec 27 '18 at 8:34
  • $\begingroup$ 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. $\endgroup$ – Thando Dec 27 '18 at 8:39
  • $\begingroup$ 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. $\endgroup$ – Henrik Schumacher Dec 27 '18 at 9:12
  • $\begingroup$ 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. $\endgroup$ – Henrik Schumacher Dec 27 '18 at 9:12
  • $\begingroup$ Sure I want to solve the PDE on the surface of the cube. Sorry about the code example I gave. Thanks. $\endgroup$ – Thando Dec 27 '18 at 10:41
2
$\begingroup$

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["NDSolve`FEM`"]
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]}

fig1 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]}

fig2

$\endgroup$
  • $\begingroup$ 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}] ? $\endgroup$ – Thando Dec 27 '18 at 9:08
  • $\begingroup$ 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 $? $\endgroup$ – Alex Trounev Dec 27 '18 at 12:40
  • $\begingroup$ @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. $\endgroup$ – Alex Trounev Dec 27 '18 at 12:56
  • $\begingroup$ 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. $\endgroup$ – Thando Dec 27 '18 at 13:51
  • $\begingroup$ @Thando Then what is truncated octahedron and 3D for? $\endgroup$ – Alex Trounev Dec 27 '18 at 14:03

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.