# Boundary conditions issue with Laplace equation

I'm trying to solve the Laplace equation for two spheres of radius rs raised to a potential of 1 a distance l away from each other. The potential at the space boundary is 0. This is what I have now:

l = 40;
rs = 15;
r = 50;
\[CapitalOmega] = DiscretizeRegion[RegionDifference[Cuboid[{-r, -r, -r}, {r, r, r}],
Ball[{{l/2, 0, 0}, {-l/2, 0, 0}}, rs]]];
NDSolveValue[{Laplacian[V[x, y, z], {x, y, z}] == 0,
DirichletCondition[V[x, y, z] == 1,
{x, y, z} \[Element] Sphere[{{-l/2, 0, 0}, {l/2, 0, 0}}, rs]],
V[r, y, z] == V[-r, y, z] == V[x, r, z] == V[x, -r, z] ==
V[x, y, r] == V[x, y, -r] == 0}, V, {x, y,
z} \[Element] \[CapitalOmega]]

and I get this error:

NDSolveValue::bcnop: No places were found on the boundary where {x,y,z}\[Element]
Sphere[{{-20,0,0},{20,0,0}},15] was True, so DirichletCondition[V==1,{x,y,z}\[Element]Sphere[{{-20,0,0},{20,0,0}},15]]
will effectively be ignored.

Any help would be appreciated.

• These conditions are not clear V[r, y, z] == V[-r, y, z] == V[x, r, z] == V[x, -r, z] == V[x, y, r] == V[x, y, -r] == 0. Are they periodic boundary conditions? Commented Jul 13 at 2:41
• @AlexTrounev I want the potential to be 0 at infinity in the actual problem, so I make it 0 at my boundary. Commented Jul 13 at 8:40
• The main problem is that the multi Ball definition is not recognized. For the FEM code I just committed a fix, that will be available in the version after 14.1 (i.e. 14.2 or 15.0). Unit then you can use the workarounds presented here. Commented Jul 15 at 6:20

We can use OpenCascadeLink to generate mesh for this problem as follows (note, we rescale l, r, rs because the solution does not depend on the choice of scale)

<< NDSolveFEM; Needs["OpenCascadeLink"]
L = 4; r = 5; rs = 3/2;

ball1 = Ball[{-L/2, 0, 0}, rs]; ball2 = Ball[{L/2, 0, 0}, rs]; cube1 =
Cuboid[{-r, -r, -r}, {r, r, r}]; shape1 = OpenCascadeShape[ball1];

mesh = ToElementMesh[bmesh];

mesh["Wireframe"]

Solution

sol = NDSolveValue[{-Laplacian[u[x, y, z], {x, y, z}] == 0,
DirichletCondition[
u[x, y, z] ==
1, (x - L/2)^2 + y^2 + z^2 == rs^2 || (x + L/2)^2 + y^2 + z^2 ==
rs^2], DirichletCondition[u[x, y, z] == 0,
x^2 == r^2 || y^2 == r^2 || z^2 == r^2]}, u,
Element[{x, y, z}, mesh]];

Visualization

DensityPlot[sol[x, y, 0], {x, -r, r}, {y, -r, r},
ColorFunction -> "Rainbow", PlotPoints -> 50,
PlotLegends -> Automatic]

Update 1. Unfortunately, the quality of this solution is not very high. For example, if we plot sol[x,0,0], we will see that the boundary conditions are not met. We can try to improve solution by adding elements and holes as in user21 answer

mesh1 = ToElementMesh[bmesh, MaxCellMeasure -> .025,
"RegionHoles" -> {{L/2, 0, 0}, {-L/2, 0, 0}}]
mesh1["Wireframe"]

Solution

sol1 = NDSolveValue[{-Laplacian[u[x, y, z], {x, y, z}] == 0,
DirichletCondition[
u[x, y, z] ==
1, (x - L/2)^2 + y^2 + z^2 == rs^2 || (x + L/2)^2 + y^2 + z^2 ==
rs^2], DirichletCondition[u[x, y, z] == 0,
x^2 == r^2 || y^2 == r^2 || z^2 == r^2]}, u,
Element[{x, y, z}, mesh1]]

Visualization

DensityPlot[sol1[x, y, 0], {x, -r, r}, {y, -r, r},
ColorFunction -> "Rainbow", PlotPoints -> 50,
PlotLegends -> Automatic]

So, we see that the holes are in place. But if we now build a solution on the line y=z=0, then we will see that the boundary conditions are not met,

Plot[sol1[x, 0, 0], {x, -r, r}, PlotPoints -> 200, PlotRange -> All]

If we compare 2 solutions, we will see that solution sol looks more preferable,

Show[Plot[sol[x, 0, 0], {x, -r, r}, PlotPoints -> 200,
PlotRange -> All, PlotStyle -> {Red, Dashed}],
Plot[sol1[x, 0, 0], {x, -r, r}, PlotPoints -> 200, PlotRange -> All]]

Update 2 We can improve solutions with holes and without holes, if we implement boundary conditions on spheres using ElementMarker as sugested by user21. With our first code we have

sol = NDSolveValue[{-Laplacian[u[x, y, z], {x, y, z}] == 0,
DirichletCondition[u[x, y, z] == 1,
ElementMarker == 7 || ElementMarker == 8],
DirichletCondition[u[x, y, z] == 0,
x^2 == r^2 || y^2 == r^2 || z^2 == r^2]}, u,
Element[{x, y, z}, mesh]]

And with our second code we have

sol1 = NDSolveValue[{-Laplacian[u[x, y, z], {x, y, z}] == 0,
DirichletCondition[u[x, y, z] == 1,
ElementMarker == 7 || ElementMarker == 8],
DirichletCondition[u[x, y, z] == 0,
x^2 == r^2 || y^2 == r^2 || z^2 == r^2]}, u,
Element[{x, y, z}, mesh1]]

Now we plot sol, sol1 on the line y=z=0and see very good quality

Show[Plot[sol1[x, 0, 0], {x, -r, r}, PlotPoints -> 200,
PlotRange -> All, PlotStyle -> Blue],
Plot[sol[x, 0, 0], {x, -r, r}, PlotPoints -> 200, PlotRange -> All,
PlotStyle -> {Red, Dashed}]]

• You should be able to directly use OpenCascadeShape on the complete RegionDifference expression. Infact you should be able to call ToElementMesh on the complete Region expression. Commented Jul 14 at 12:36
• @user21 Can you show in code what actually I should be able? :) Commented Jul 14 at 12:39
• @AlexTrounev Thank you for your solution. Is there any way I can improve the accuracy? If you plot Plot[sol[x, 0, 0], {x, -r, r}] you can see that the potential at the spheres isn't really 1 and not quite constant inside. I tried MaxCellMeasure and MeshQualityGoal when defining the mesh and some options in NDSolve but they seem not to have any effect. Commented Jul 14 at 16:17
• Actually your result is not quite correct, I you'd need to specify region holes. Commented Jul 15 at 5:19
• Added an answer. Have a look. Commented Jul 15 at 5:38

In Version 14.0 this works:

Needs["NDSolveFEM"]

l = 40;
rs = 15;
r = 50;
reg=RegionDifference[Cuboid[{-r,-r,-r},{r,r,r}],Ball[{{l/2,0,0},{-l/\
2,0,0}},rs]];

You'd need to specify the region holes:

\[CapitalOmega] =
ToElementMesh[reg, "RegionHoles" -> {{l/2, 0, 0}, {-l/2, 0, 0}}]

A slightly more accurate way to do it:

reg = RegionDifference[Cuboid[{-r, -r, -r}, {r, r, r}],
RegionUnion[Ball[#, rs] & /@ {{l/2, 0, 0}, {-l/2, 0, 0}}]];
\[CapitalOmega] =
ToElementMesh[reg, "RegionHoles" -> {{l/2, 0, 0}, {-l/2, 0, 0}}]

In Version 14.1 (you no longer need the "RegionHoles" specification in this case)

You can then use

sol = NDSolveValue[{Laplacian[V[x, y, z], {x, y, z}] == 0,
DirichletCondition[
V[x, y, z] == 1, {x, y, z} \[Element]
Sphere[{{-l/2, 0, 0}, {l/2, 0, 0}}, rs]],
V[r, y, z] == V[-r, y, z] == V[x, r, z] == V[x, -r, z] ==
V[x, y, r] == V[x, y, -r] == 0},
V, {x, y, z} \[Element] \[CapitalOmega]]

For the DirichletCondition predicate a more robust version is to use:

sol = NDSolveValue[{Laplacian[V[x, y, z], {x, y, z}] == 0,
DirichletCondition[
V[x, y, z] ==
1, (Abs[(x - l/2)^2 + y^2 + z^2 - rs^2] <=
10^-5) || (Abs[(x + l/2)^2 + y^2 + z^2 - rs^2] <= 10^-5)],
V[r, y, z] == V[-r, y, z] == V[x, r, z] == V[x, -r, z] ==
V[x, y, r] == V[x, y, -r] == 0},
V, {x, y, z} \[Element] \[CapitalOmega]]

DensityPlot[sol[x, y, 0], {x, -r, r}, {y, -r, r},
ColorFunction -> "Rainbow", PlotPoints -> 50,
PlotLegends -> Automatic]

Note, that in Alexei's answer the region has elements inside the Balls.

You may have to play a bit with the tolerance of the boundary conditions. If you do not want to do that, then using markers removes this uncertainty entirely. Let me show you how:

Needs["NDSolveFEM`"]

l = 40;
rs = 15;
r = 50;
reg = RegionDifference[Cuboid[{-r, -r, -r}, {r, r, r}],
RegionUnion[Ball[#, rs] & /@ {{l/2, 0, 0}, {-l/2, 0, 0}}]];
\[CapitalOmega] =
ToElementMesh[reg, "RegionHoles" -> {{l/2, 0, 0}, {-l/2, 0, 0}}];
\[CapitalOmega]["BoundaryElementMarkerUnion"]

(* {1, 2, 3, 4, 5, 6, 7, 8} *)

Each continuous surface is attributed with a marker (See documentation, for example the ElementMesh generation tutorial), Here we use the face markers, other markers exist but are not relevant here.

The following code from the ElementMesh Visualization tutorial helps us to scroll through the boundaries.

bmesh = \[CapitalOmega];
bIDs = bmesh["BoundaryElementMarkerUnion"];
edgeframe = bmesh["Edgeframe"];
outline = bmesh["Wireframe"["MeshElementStyle" ->
Directive[Opacity[0.2], FaceForm[LightBlue], EdgeForm[]]]];
Manipulate[Show[
outline,
edgeframe,
bmesh["Wireframe"[ElementMarker == bIDs[[id]],
"MeshElementStyle" ->
Directive[FaceForm[LightGray], EdgeForm[]]]]
], {{id, 1, "ElementMarker ID"}, 1, Length[bIDs], 1,
Appearance -> "Open"}, SaveDefinitions -> True]

From this we find that the boundary element markers have IDs 7 and 8 and we use them in the predicate:

sol = NDSolveValue[{Laplacian[V[x, y, z], {x, y, z}] == 0,
DirichletCondition[V[x, y, z] == 1,
ElementMarker == 7 || ElementMarker == 8],
V[r, y, z] == V[-r, y, z] == V[x, r, z] == V[x, -r, z] ==
V[x, y, r] == V[x, y, -r] == 0},
V, {x, y, z} \[Element] \[CapitalOmega]];

DensityPlot[sol[x, y, 0], {x, -r, r}, {y, -r, r},
ColorFunction -> "Rainbow", PlotPoints -> 50,
PlotLegends -> Automatic]

• Thank you very much (+1). Unfortunately, this solution has less quality than that computed without holes. Commented Jul 15 at 7:38
• @AlexTrounev can you quantify that? Commented Jul 15 at 8:02
• Please, see update to my answer. Commented Jul 15 at 8:31
• Yes, last solution is very good. Commented Jul 15 at 10:23

It works if you remove one of the two balls. For example : :

l = 40;
rs = 15;
r = 50;
domain =
DiscretizeRegion[
RegionDifference[Cuboid[{-r, -r, -r}, {r, r, r}],
Ball[(*{{l/2,0,0},{-l/2,0,0}}*)   {l/2, 0, 0}   , rs]]
, PlotRange -> {{-50, 50}, {0, 50}, {0, 50}}]

sol = NDSolveValue[{Laplacian[V[x, y, z], {x, y, z}] == 0,
DirichletCondition[
V[x, y, z] == 1, {x, y, z} \[Element]
Sphere[{{-l/2, 0, 0}, {l/2, 0, 0}}, rs]],
V[r, y, z] == V[-r, y, z] == V[x, r, z] == V[x, -r, z] ==
V[x, y, r] == V[x, y, -r] == 0},
V, {x, y, z} \[Element] domain]
ContourPlot3D[sol[x, y, z], {x, y, z} \[Element] domain,
PlotRange -> {{-50, 50}, {0, 50}, {0, 50}}, PlotPoints -> 10,
BoxRatios -> {2, 1, 1}] // EchoTiming

Timing : 347.187 seconds (!)

It seems, at least on my version 12.2 that the problem comes from de definition of the domain. Your code with two balls gives :

l = 40;
rs = 15;
r = 50;
DiscretizeRegion[
RegionDifference[Cuboid[{-r, -r, -r}, {r, r, r}],
Ball[{{l/2, 0, 0}, {-l/2, 0, 0}}, rs]]
,PlotRange -> {{0, 50}, {0, 50}, {0, 50}}]

The holes are missing.
I have tried some other syntaxes, for example RegionDifference[RegionDifference[...,Ball[...] ],Ball[...]] and RegionDifference[...,RegionUnion[Ball[...],Ball[...]] ]. It doesn't work neither.
When the family Region/MeshRegion/RegionDifference... was introduced several years ago, there were bugs. Since then the situation get better, so maybe your code works in more recent versions of Mathematica.

• I'm on 13.3 and it doesn't draw them out either. Strangely, if they overlap (l/2 < R) everything works fine. Also, if I don't use DiscretizeRegion when defining the domain, it draws them out, but then NDSolve gives a mesh error (mesh cannot be generated). Commented Jul 12 at 19:02
• In a recent version you should be able to call ToElementMesh on the region and it should work Commented Jul 14 at 12:37