# Setting Up Boundary Conditions for Magnetostatic PDE

Bug introduced in 11.2.0 and fixed in 11.3.0

The system is a hollow cylinder (thin solenoid) with a current density $\text{J}$ and I'm looking to solve the magnetic potential ($\text{A}$) inside the radius of the cylinder using FEA. I'm struggling with the boundary conditions for the open ends (and in general). Suggestions are appreciated.

u0 = 4 Pi 1*^-7;

j = 100; (*current density parallel to z-axis*)

solution =
NDSolveValue[{
-Div[(1/u0) Grad[u[r, z], {r, th, z}, "Cylindrical"],
{r, th, z}, "Cylindrical"] == NeumannValue[j, r == 1],
DirichletCondition[u[r, z] == 0, z^2 == 1]},
u, {r, 0, 1}, {z, -1, 1},
Method -> {"PDEDiscretization" -> {"FiniteElement",
"MeshOptions" -> {"MaxBoundaryCellMeasure" -> 0.01}}}]

ContourPlot[solution[r, z], {r, z} ∈ solution["ElementMesh"],
PlotRange -> All, Contours -> 19, ImageSize -> Medium]


To compute the magnetic density field ($\text{B}$):

bField = Curl[solution[r, z], {r, th, z}, "Cylindrical"]


Update:

I've done additional research and the following is an approach by FEMM to solve such problems. The cross section is still in cylindrical coordinates but the open boundary is spherical.

bound = Table[{boundRadius Sin[theta],
boundRadius Cos[theta]}, {theta, 0, Pi, Pi/100}];
boundIndex = Partition[Last[FindShortestTour[bound]], 2, 1];

coil = {{1, -1}, {1, 1}};
coilIndex = {{Length[bound] + 1, Length[bound] + 2}};

reg = ImplicitRegion[r^2 + z^2 < boundRadius^2 && r >= 0, {r, z}];

bmesh = ToBoundaryMesh["Coordinates" -> Join[bound, coil],
"BoundaryElements" -> {LineElement[Join[boundIndex, coilIndex]]}];
mesh = ToElementMesh[bmesh];

Show[
RegionPlot[reg, Epilog -> Point[bound], AspectRatio -> 2,
ImageSize -> Medium],
bmesh["Wireframe"]
] I need to apply the following mixed boundary condition to the spherical surface:

$\frac{1}{\mu_0}\frac{\partial A}{\partial r}+\frac{A}{\mu_0 R} = 0$

I have the following failed attempt (NDSolveValue spits out bcnop warning and the internal b.c. is ignored):

u0 = 4 Pi 1*^-7;
j = 100;

solution =
NDSolveValue[{-Div[(1/u0) Grad[u[r, z], {r, th, z}, "Cylindrical"],
{r, th, z}, "Cylindrical"] - j ==
NeumannValue[0, r == 0] -
DirichletCondition[u[r, z] == u0 j, r == 1 && z^2 <= 1]},
u, {r, z} ∈ mesh,
Method -> {"PDEDiscretization" -> {"FiniteElement"}}]

ContourPlot[solution[r, z], {r, z} ∈ solution["ElementMesh"],
PlotRange -> All, Contours -> 19, ImageSize -> Medium,
AspectRatio -> 2]

• Just a precison for beginners in magnetostatic : this problem can be solved anatically. See here – andre314 Sep 22 '17 at 18:27
• That's what I was thinking. My message was for "visitors" – andre314 Sep 22 '17 at 18:30
• The problem of open ends (=boundaries at infinity) arrises very often in magnetostatics and electrostatics. Example : In electrostatics in 2 dimensions, on may use a conformal mapping to map infinity to finite boundaries, and then use FEM. Otherwise, FEM solvers have special complex methods for dealing with some "infinite mesh elements" (?). I'm afraid that's not a simple affair. Would be happy to be wrong. – andre314 Sep 22 '17 at 20:03
• OK, @xzczd there is a regression in V11.2 that triggered the message and that will need to get fixed. Never the less I think the wire should have some extension. – user21 Oct 1 '17 at 12:07
• @xzczd, I just merged a fix for this problem, the next release will work as expected. Thanks for digging this out. – user21 Oct 17 '17 at 0:57

This is fixed in version 11.3:

Needs["NDSolveFEM"]
boundRadius Cos[theta]}, {theta, 0, Pi, Pi/100}];
boundIndex = Partition[Last[FindShortestTour[bound]], 2, 1];

coil = {{1, -1}, {1, 1}};
coilIndex = {{Length[bound] + 1, Length[bound] + 2}};

reg = ImplicitRegion[r^2 + z^2 < boundRadius^2 && r >= 0, {r, z}];

bmesh = ToBoundaryMesh["Coordinates" -> Join[bound, coil],
"BoundaryElements" -> {LineElement[Join[boundIndex, coilIndex]]}];
mesh = ToElementMesh[bmesh];

Show[RegionPlot[reg, Epilog -> Point[bound], AspectRatio -> 2,
ImageSize -> Medium], bmesh["Wireframe"]] Show[Graphics[
GraphicsComplex[mesh["Coordinates"],
Point[Flatten[ElementIncidents[mesh["PointElements"]]]]]]] Note, how the interior points are now (again) part of the boundary mesh. And the solution:

u0 = 4 Pi 1*^-7;
j = 100;

solution =
NDSolveValue[{-Div[(1/u0) Grad[u[r, z], {r, th, z},
"Cylindrical"], {r, th, z}, "Cylindrical"] - j ==
NeumannValue[0, r == 0] -
DirichletCondition[u[r, z] == u0 j, r == 1 && z^2 <= 1]},
u, {r, z} \[Element] mesh];

ContourPlot[solution[r, z], {r, z} \[Element] solution["ElementMesh"],
PlotRange -> All, Contours -> 19, ImageSize -> Medium,
AspectRatio -> 2, PlotPoints -> 50] Well, I'm not quite familiar with electromagnetism and it's not immediately clear to me how to compare the numeric solution with the analytic solution, but I've circumvented the bcnop warning anyway, so let me post an answer.

The idea is creating a very narrow slit in the domain to approximate the wire i.e. the LineElement in your failed attempt. Notice I've also enlarge boundRadius because after some trial I noticed boundRadius = 2 is too small:

boundRadius = 4;
u0 = 4 Pi 1*^-7;
j = 100;

eps = 10^-4;

reg = ImplicitRegion[
r^2 + z^2 < boundRadius^2 && r >= 0 && ! (1 - eps < r < 1 + eps && z^2 < 1), {r, z}];

solution = NDSolveValue[{-Div[(1/u0) Grad[u[r, z], {r, th, z}, "Cylindrical"], {r, th,
z}, "Cylindrical"] -
DirichletCondition[u[r, z] == u0 j, 1 - eps < r < 1 + eps && z^2 < 1]},
u, {r, z} ∈ reg,
Method -> {FiniteElement, MeshOptions -> MaxCellMeasure -> 0.001}]

interestedreg =
ImplicitRegion[
r^2 + z^2 < interestedRadius^2 && ! (1 - eps < r < 1 + eps && z^2 < 1), {r, z}];

ContourPlot[
Evaluate@PiecewiseExpand@If[r > 0, solution[r, z], solution[-r, z]], {r, z} ∈
interestedreg, PlotRange -> All, PlotPoints -> 100, Contours -> 19] • I think this is the way to go - also from a physics point of view. You may want to create the mesh manually like in the question for very 'thin' wires. – user21 Oct 1 '17 at 8:41
• I eventually obtained a result close to your graph but when you take the Curl of solution to obtain the magnetic field it doesn't match the analytical solution. – Young Oct 1 '17 at 17:32
• I think there is something wrong with how I formed the differential equation. – Young Oct 1 '17 at 17:39
• @Young Er…sorry, but the $A$ i.e. u in your equation is a scalar, how can you calculate its Curl? – xzczd Oct 2 '17 at 2:50
• @xzczd It's a little bit subtle : The value u[x,y] calculated by Young is only a component of the potential vector (the component Aphi in cylindrical coordinate, that is to the component perpendicular to the paper). The other components are 0 in this case (because of the symetries of the problem). – andre314 Oct 2 '17 at 5:51

Extended comment

The result of the OP's code are not the same on Mathematica 11.0 and Mathematica 11.2

Here is the code I have used :

<<NDSolveFEM

boundRadius Cos[theta]}, {theta, 0, Pi, Pi/100}];
boundIndex = Partition[Last[FindShortestTour[bound]], 2, 1];

coil = {{1, -1}, {1, 1}};
coilIndex = {{Length[bound] + 1, Length[bound] + 2}};

reg = ImplicitRegion[r^2 + z^2 < boundRadius^2 && r >= 0, {r, z}];

bmesh = ToBoundaryMesh["Coordinates" -> Join[bound, coil],
"BoundaryElements" -> {LineElement[Join[boundIndex, coilIndex]]}];
mesh = ToElementMesh[bmesh];

(*Show[
RegionPlot[reg, Epilog -> Point[bound], AspectRatio -> 2,
ImageSize -> Medium],
bmesh["Wireframe"]
]*)

u0 = 4 Pi 1*^-7;
j = 100;

solution =
NDSolveValue[{-Div[(1/u0) Grad[u[r, z], {r, th, z}, "Cylindrical"],
{r, th, z}, "Cylindrical"] - j ==
NeumannValue[0, r == 0] -
DirichletCondition[u[r, z] == u0 j, r == 1 && z^2 <= 1]},
u, {r, z} \[Element] mesh,
Method -> {"PDEDiscretization" -> {"FiniteElement"}}]

ContourPlot[solution[r, z], {r, z} \[Element] solution["ElementMesh"],
PlotRange -> All, Contours -> 19, ImageSize -> Medium,
AspectRatio -> 2]


Mathematica 11.2 Mathematica 11.0 - first evaluation (on a fresh kernel)

It doesn't work : Mathematica 11.0 - second evaluation (= reevaluation of the same code) Edit

same as before, but with PlotPoints->50added to ContourPlot: • Will the last result be improved if you add e.g. PlotPoints->50 to ContourPlot? – xzczd Oct 1 '17 at 10:55
• @xzczd see my edit. Nice ! – andre314 Oct 1 '17 at 11:07
• It helps if you define the boundRadius = 2; before the bound. – user21 Oct 1 '17 at 11:54
• OK, thanks Andre (@xzczd) this is a regression in 11.2 and I'll see to get it fixed for the next release. – user21 Oct 1 '17 at 12:05
• A precision for visitors : the comment of @user21 above (boundRadius =2 before the bound) is destinated to fix the problem of the failing first evaluation on Mathematica version 11.0 – andre314 Oct 1 '17 at 12:08

Extended comment:

Something like this

epsi = 10^-1;
coil = {{1 - epsi, -1}, {1 - epsi, 1}, {1 + epsi,
1}, {1 + epsi, -1}};
coilIndex = {{1, 2}, {2, 3}, {3, 4}, {4, 1}} + Length[bound];