Electrostatics: Finite Elements

I'm messing around with FEM in mathematica and am having trouble solving a very simple problem of the electric field around a unifromly charged sphere.

Here is my workflow.

Needs["NDSolveFEM"]
Needs["FEMAddOns"]

(* Define Boundary *)
domain = ToBoundaryMesh[Rectangle[{-3, -3}, {3, 3}]];
circle = ToBoundaryMesh[Disk[]];
bmesh = BoundaryElementMeshJoin[domain, circle]

(* Define Elements *)
air = {0, 2};
sphere = {0, 0};
markerSpecification = {{air, 1}, {dielectric, 2}};
mesh = ToElementMesh[bmesh, "RegionMarker" -> markerSpecification,
"MeshOrder" -> 1, "NodeReordering" -> True];
mesh["Wireframe"["MeshElementStyle" -> FaceForm /@{White, LightRed}]]


(* Solve Laplaces equation with charge density = 1 inside the unit disk *)
usol = NDSolveValue[{-Laplacian[v[x, y], {x, y}] ==
Piecewise[{{1, x^2 + y^2 <= 1}}, 0]
(* Set potential to 1 on the boundary of the disk *)
, DirichletCondition[v[x, y] == 0, x^2 + y^2 == 1]}
, v, {x, y} \[Element] mesh]
ContourPlot[usol[x, y], {x, y} \[Element] mesh]


Thats really not the result I expected since the solution should be Piecewise[{{r^2/4, r < 1}}, Log[r]/2 + 1/4]

I also get even worse results if I make the potential = 0 on the boundary.

• You solve the 2D problem for the disk, but you are trying to compare it with the solution for the sphere. Piecewise[{{r, r <= 1}}, 1/r^2] it a solution for a charged sphere for radial component of the electric field, not for potential. Nov 19, 2019 at 17:05
• This link might be interesting for you. Nov 19, 2019 at 17:17
• @AlexTrounev Thats a good point it's the cross section of a cylinder so should be log(r) for the potential and 1/r for the field Nov 19, 2019 at 19:16
• @AlexTrounev Asymptotically, the potential should be $1/r$ as follows from the series expansion of the elliptic functions, not $\log(r)$ Nov 19, 2019 at 19:46
• As soon as the total charge is 0, the sum of the contributions is 1/r in the expansion, but a charge alone contributes as Log[r] (in 2 D). Nov 19, 2019 at 22:24

We solve the problem for the disk and compare the numerical solution and the analytical

Needs["NDSolveFEM"];
reg = Rectangle[{-3., -3.}, {3., 3.}];
mesh = ToElementMesh[reg,
MeshRefinementFunction ->
Function[{vertices, area},
area > 0.0001 (0.1 + 10 Norm[Mean[vertices]])]]
mesh["Wireframe"]

eq = Laplacian[u[x, y], {x, y}] == Piecewise[{{1, x^2 + y^2 <= 1}}, 0];
bc = DirichletCondition[u[x, y] == Log[x^2 + y^2]/4 + 1/4, True];
U = NDSolveValue[{eq, bc}, u, {x, y} \[Element] mesh]


Here, on the left, the two solutions almost coincide; on the right, the difference of solutions

{Plot[{U[x, 0],
Piecewise[{{x^2/4, -1 <= x <= 1}, {Log[x^2]/4 + 1/4,
True}}]}, {x, -3, 3}, PlotLegends -> {"FEM", "Analytical"}],
Plot[{U[x, 0] -
Piecewise[{{x^2/4, -1 <= x <= 1}, {Log[x^2]/4 + 1/4,
True}}]}, {x, -3, 3}]}


Potential and electric field

f = Evaluate[-Grad[U[x, y], {x, y}]];

ContourPlot[U[x, y], {x, y} \[Element] mesh, Contours -> 20,
ColorFunction -> "Rainbow"]

StreamDensityPlot[f, {x, y} \[Element] mesh,
ColorFunction -> "Rainbow", PlotLegends -> Automatic]


• Alex, please provide some clarity for me with this answer, can you use these same techniques to do magnetic fields? Your answers seem so simple when I see them, and it would be great to see an example with a uniformly magnetized sphere ;) Nov 20, 2019 at 4:19
• @CATrevillian Another issue is discussed here. But if you open the topic and show the code, then perhaps we can solve your problem. For the magnetic field there are other algorithms. Nov 20, 2019 at 12:10
• Thanks, any insight why my approach failed? What if we had a more complex object then a sphere and wanted to specify the charge density based on a mesh? For example I would like to see worflow for solving the efeild around this Import["https://static.vecteezy.com/system/resources/previews/000/552/\ 501/large_2x/vector-heart-romantic-love-graphic.jpg"] // ColorNegate // ImageMesh // RegionResize[#, {{-1, 1}, {-1, 1}}] & // ToBoundaryMesh assuming it is uniformely charged. Nov 20, 2019 at 14:07
• This link is unavailable. Nov 20, 2019 at 20:37
• A line break character was inserted use "static.vecteezy.com/system/resources/previews/000/552/501/…" Nov 28, 2019 at 21:03

Here is a way to get much better accuracy for u in much less time.

Because this has a discontinuity it is better to make a mesh that respects that.

Needs["NDSolveFEM"];
reg = RegionDifference[Rectangle[{-3., -3.}, {3., 3.}], Disk[]];
mesh = ToElementMesh[reg, "RegionHoles" -> None,
"RegionMarker" -> {{{0, 0}, 1}, {{2, 0}, 2}}]
mesh["Wireframe"]


I also added markers which makes the equation set up a bit easier I find:

eq = Laplacian[u[x, y], {x, y}] == If[ElementMarker == 1, 1, 0];
bc = DirichletCondition[u[x, y] == Log[x^2 + y^2]/4 + 1/4, True];
sol = NDSolveValue[{eq, bc}, u, {x, y} \[Element] mesh];


This computes the error to the analytical solution.

Plot[{
sol[x, 0] -
Piecewise[{{x^2/4, -1 <= x <= 1}, {Log[x^2]/4 + 1/4, True}}],
U[x, 0] -
Piecewise[{{x^2/4, -1 <= x <= 1}, {Log[x^2]/4 + 1/4,
True}}]}, {x, -3, 3}, PlotRange -> All]
`

The orange line is the result from Alex the blue line is this computation.