# 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. – Alex Trounev Nov 19 '19 at 17:05
• This link might be interesting for you. – yarchik Nov 19 '19 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 – user2757771 Nov 19 '19 at 19:16
• @AlexTrounev Asymptotically, the potential should be $1/r$ as follows from the series expansion of the elliptic functions, not $\log(r)$ – yarchik Nov 19 '19 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). – andre314 Nov 19 '19 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 ;) – CA Trevillian Nov 20 '19 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. – Alex Trounev Nov 20 '19 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. – user2757771 Nov 20 '19 at 14:07
• This link is unavailable. – Alex Trounev Nov 20 '19 at 20:37
• A line break character was inserted use "static.vecteezy.com/system/resources/previews/000/552/501/…" – user2757771 Nov 28 '19 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.